C
C Subroutine file 'mat.for' with subroutines dealing with matrices.
C
C Version: 6.70
C Date: 2012, November 7
C
C.......................................................................
C
C This file consists of the following external procedures:
C     WMATH... Subroutine designed to write a matrix header file.
C             WMATH
C     RMATH... Subroutine designed to read a matrix header file.
C             RMATH
C     WMATD... Subroutine designed to write a matrix data file.
C             WMATD
C     RMATD... Subroutine designed to read a matrix data file.
C             RMATD
C     WMATR... Auxiliary subroutine to WMATD, designed to write
C             the array of matrix elements into a file.
C             WMATR
C     TMATR... Subroutine to transpose a matrix.
C             TMATR
C     SMATRE... Subroutine to extend sparse symmetric or antisymmetric
C             matrix into sparse general matrix.
C             SMATRE
C     GSMATR... Subroutine to change normal (not sparse) matrix into
C             the sparse matrix.
C             GSMATR
C     SGMATR... Subroutine to change sparse matrix into the normal
C             (not sparse) matrix.
C             SGMATR
C     GSMAT... Subroutine to change normal (not sparse) matrix into
C             the sparse matrix if the matrix is sparser than given
C             limit.
C             GSMAT
C     VELEM... Real function to find the value of an element
C             of a matrix.
C             VELEM
C     NELMAT... Integer function to calculate total number of elements
C             of a non-sparse matrix.
C             NELMAT
C     NSPMAT... Integer function to calculate number of zero elements
C             of a non-sparse matrix.
C             NSPMAT
C     ISYM... Integer function to assign the integer corresponding
C             to the symmetry.
C             ISYM
C     MSHIFT... Subroutine designed to shift the matrix in memory.
C             MSHIFT
C Subroutines to enable functionality of some programs which
C do not use above subroutines directly:
C     OMAT... Subroutine designed to set the form of the file for
C             reading or writing a matrix, and to open the file.
C             OMAT
C     WMAT... Subroutine designed to write a given matrix into the given
C             file.
C             WMAT
C     RMAT... Subroutine designed to read a matrix from the given file.
C             RMAT
C
C=======================================================================
C                                                     
C Storage of the matrices in the memory:
C     The matrices are stored in real array RAM of common block RAMC
C     defined in file ram.inc.
C     To store integers in this array, integer array IRAM coinciding
C     with real array RAM is also defined.
C     Normal (not sparse) matrices:
C             Normal (dense) matrices are identified by the value of
C             SPARSE=' ' in the matrix header file.  For normal matrices
C             all values of the matrix are stored columnwise, starting
C             from the first value of the first column to the last
C             stored value of the first column, then the second column,
C             and so on to the value of the last stored element of the
C             last column.
C             If the matrix is other than general, only N elements
C             corresponding to the given symmetry are stored, i.e. for
C             a symmetric matrix, just elements from the first row to
C             the diagonal are stored for each column (N=M1*(M1+1)/2);
C             for a skew matrix, just elements above the diagonal
C             are stored for each column (N=M1*(M1-1)/2); for a diagonal
C             matrix, just diagonal elements are stored (N=M1).
C     Sparse matrices:
C             Sparse matrices are identified by the value of
C             SPARSE='CSC'. NELEM is the number of nonzero elements
C             of the matrix, calculating only the elements corresponding
C             to the given symmetry, see above.
C             The first M2 storage locations contain the addresses
C             of the beginnings of all columns of the matrix.
C             The nonzero elements of column J of the matrix are thus
C             stored in the positions from (I)RAM(IRAM(IMAT+J-1))
C             to (I)RAM(IRAM(IMAT+J)-1), where IMAT is the address
C             of the first storage location of the matrix.
C             The (M2+1)th storage location contains the address in the
C             array just after the last stored element.  This address
C             may be understood as the address of the beginning of a
C             fictitious column number M2+1.
C             Next 2*NELEM storage locations contain the nonzero
C             elements of the matrix stored columnwise. Each element
C             occupies two storage locations, the first one containing
C             index of the row of the element and written in the integer
C             array, the second one containing the value of the element
C             and written in the real array.  In this form the matrix
C             occupies 1+M2+2*NELEM storage locations.
C
C For the description of the forms of the disk files with matrices
C and of the parameters describing the matrices in the matrix header
C files refer to file forms.htm.
C
C=======================================================================
C
C     
C
      SUBROUTINE WMATH(LU,FILEH,FILED,M1,M2,SPARSE,NELEM,SYMM,FORM)
      INTEGER LU,M1,M2,NELEM
      CHARACTER*(*) FILEH,FILED,SPARSE,SYMM,FORM
C
C Subroutine designed to write the matrix header file.
C
C Input:
C     LU...   Logical unit number of the matrix header output file.
C             The output matrix header file is opened, written, and
C             closed.
C     FILEH... String containing the name of the matrix header file.
C             The matrix header file has the form of the SEP parameter
C             file.
C             If FILEH=' ', no action is done.
C     FILED... String containing the name of the matrix data file.
C             If FILED=' ', default filename is determined and returned.
C             The default matrix data filename is constructed from the
C             matrix header filename by either replacing its last
C             character by '@' if the filename has a three-character
C             extension, or by adding '@' in other cases.
C     M1...   Number of rows of the matrix.
C     M2...   Number of columns of the matrix.
C     SPARSE... Sparseness of the matrix data file:
C             SPARSE=' ':  Matrix is dense (not sparse).
C             SPARSE='CSC' (case insensitive): Matrix is a sparse matrix
C               stored in the compressed column format.
C     NELEM... Number of nonzero elements of the sparse matrix,
C             calculating only the elements corresponding to the
C             symmetry of the matrix.
C             For SPARSE=' ', NELEM has no meaning and is not used.
C     SYMM... Symmetry of the matrix data file:
C             SYMM=' ': General matrix.
C             SYMM='sym': Symmetric matrix.
C             SYMM='skew': Skew matrix.
C             SYMM='diag': Diagonal matrix.
C     FORM... Form of the matrix data file:
C             FORM='FORMATTED': Formatted file.
C             FORM='UNFORMATTED': Unformatted file.
C             If FORM=' ', default FORM is determined according to the
C             SEP parameter FORMM and returned.  Default for FORMM is
C             FORMM='FORMATTED'.
C
C Output:
C     FILED... If FILED=' ' on input, default filename is determined
C             and returned.
C             Otherwise, the input value of FILED is returned.
C     FORM... If FORM=' ' on input, default FORM for writing
C             is determined, converted to uppercase, and returned.
C
C Input SEP parameter:
C     FORMM='string'... Form of the matrix data file:
C             FORMM='FORMATTED': Formatted file.
C             FORMM='UNFORMATTED': Unformatted file.
C             Default: FORMM='FORMATTED'
C
C Subroutines and external functions required:
      EXTERNAL LENGTH,UPPER,WSEP3I,WSEP3T,RSEP3T
      INTEGER  LENGTH
C     LENGTH,UPPER... File length.for.
C     WSEP3I,WSEP3T,RSEP3T... File sep.for.
C
C Coded by Ludek Klimes and Petr Bulant
C
C-----------------------------------------------------------------------
C
      INTEGER L
      CHARACTER*72 TXTERR
C
      IF(FILEH.EQ.' ') THEN
        RETURN
      END IF
C
C     Opening the matrix header file:
      OPEN(LU,FILE=FILEH,ERR=100)
C
C     Default data file name:
      IF(FILED.EQ.' ') THEN
        L=LENGTH(FILEH)
        IF (L.GE.4) THEN
          IF (FILEH(L-3:L-3).EQ.'.') THEN
            L=L-1
          ENDIF
        ENDIF
        L=L+1
        FILED=FILEH
        FILED(L:L)='@'
      END IF
C
C     Default form of the data file:
      IF(FORM.EQ.' ') THEN
        CALL RSEP3T('FORMM',FORM,'FORMATTED')
        CALL UPPER(FORM)
      END IF
C
C     Writing the values of the parameters describing the output matrix:
      CALL WSEP3T(LU,'IN',FILED)
      CALL WSEP3I(LU,'M1',M1)
      CALL WSEP3I(LU,'M2',M2)
      CALL WSEP3T(LU,'SPARSE',SPARSE)
      IF(SPARSE.NE.' ') THEN
        CALL WSEP3I(LU,'NELEM',NELEM)
      END IF
      CALL WSEP3T(LU,'SYMMETRY',SYMM)
      CALL WSEP3T(LU,'FORM',FORM)
C
      CLOSE(LU)
      RETURN
  100 CONTINUE
C     MAT-01
      WRITE(TXTERR,'(A,A,A)') 'MAT-01: Error when opening file ''',
     *FILEH(1:MIN0(LENGTH(FILEH),37)),'''.'
      CALL ERROR(TXTERR)
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE RMATH(LU,FILEH,FILED,M1,M2,SPARSE,NELEM,SYMM,FORM)
      INTEGER LU,M1,M2,NELEM
      CHARACTER*(*) FILEH,FILED,SPARSE,SYMM,FORM
C
C Subroutine designed to read the matrix header file.
C
C Input:
C     LU...   Logical unit number of the input matrix header file.
C             The input matrix header file is opened, read, and
C             closed.
C     FILEH... String containing the name of the matrix header file.
C             The matrix header file has the form of the SEP parameter
C             file.
C             If FILEH=' ', no file is read and defaults are returned.
C The input parameters are not altered.
C
C Output:
C     FILED... String containing the name of the matrix data file.
C             The default matrix data filename is constructed from the
C             matrix header filename by either replacing its last
C             character by '@' if the filename has a three-character
C             extension, or by adding '@' in other cases.
C     M1...   Number of rows of the matrix.  Default: M1=1
C     M2...   Number of columns of the matrix.  Default: M2=1
C     SPARSE... Sparseness of the matrix data file:
C             SPARSE=' ':  Matrix is dense (not sparse).
C             SPARSE='CSC': Matrix is a sparse matrix stored in the
C               compressed column format.
C             Default: SPARSE=' '
C     NELEM... Number of elements of the matrix stored in the matrix
C             data file. I.e. either the number of nonzero elements
C             of the sparse matrix, or the number of elements
C             of dense matrix.
C     SYMM... Symmetry of the matrix data file:
C             SYMM=' ': General matrix.
C             SYMM='sym': Symmetric matrix.
C             SYMM='skew': Skew matrix.
C             SYMM='diag': Diagonal matrix.
C             Default: SYMM=' '
C     FORM... Form of the matrix data file:
C             FORM='FORMATTED': Formatted file.
C             FORM='UNFORMATTED': Unformatted file.
C             Default: FORM='FORMATTED'
C     Note: the values of SPARSE and FORM are converted to uppercase
C     on output, and the value of SYMM to lowercase, in order to
C     simplify comparing of the strings in the calling program.
C
C Subroutines and external functions required:
      EXTERNAL LENGTH,LOWER,UPPER,SSEP,RSEP1,RSEP3I,RSEP3T,ISYM,NELMAT
      INTEGER  LENGTH,ISYM,NELMAT
C     LENGTH,LOWER,UPPER ... File 'length.for'.
C     SSEP,RSEP1,RSEP3I,RSEP3T... File 'sep.for'.
C
C Coded by Ludek Klimes and Petr Bulant
C
C-----------------------------------------------------------------------
C
      CHARACTER*80 FILEDD
      INTEGER L,ISYMM
C     Indices of the sets of SEP parameters
      INTEGER ISEP,IOLD,ITMP
      SAVE ISEP
      DATA ISEP/0/
C
C     Switching to a new parameter set and reading the matrix header:
      CALL SSEP(ISEP,IOLD)
      CALL RSEP1(LU,FILEH)
C
C     Reading the values of the parameters describing the input matrix:
      L=LENGTH(FILEH)
      IF (L.GE.4) THEN
        IF (FILEH(L-3:L-3).EQ.'.') THEN
          L=L-1
        ENDIF
      ENDIF
      L=L+1
      FILEDD=FILEH
      FILEDD(L:L)='@'
      CALL RSEP3T('IN',FILED,FILEDD)
      CALL RSEP3I('M1',M1,1)
      CALL RSEP3I('M2',M2,1)
      CALL RSEP3T('SPARSE',SPARSE,' ')
      CALL RSEP3T('SYMMETRY',SYMM,' ')
      CALL RSEP3T('FORM',FORM,'FORMATTED')
C
C     Converting input strings to lowercase:
      CALL LOWER(SYMM)
      CALL UPPER(FORM)
      CALL UPPER(SPARSE)
C
C     Number NELEM of stored matrix elements:
      ISYMM=ISYM(SYMM)
      NELEM=NELMAT(M1,M2,ISYMM)
      L=NELEM
      IF(SPARSE.NE.' ') THEN
        CALL RSEP3I('NELEM',NELEM,L)
      ENDIF
C
      ISEP=-ISEP
      CALL SSEP(ISEP,ITMP)
      CALL SSEP(IOLD,ISEP)
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE WMATD(LU,FILED,M1,M2,SPARSE,NELEM,FORM,IMAT)
      CHARACTER*(*) FILED,FORM,SPARSE
      INTEGER LU,M1,M2,NELEM,IMAT
C
C Subroutine designed to write a given matrix into the matrix data file.
C The matrix to be written is assumed to be stored in common block RAMC
C of file 'ram.inc'.
C
C Input:
C     LU...   Logical unit number to be used for the output.
C     FILED... Destination filename.  If not blank, the file will be
C             opened and closed.  If blank, the file is assumed to be
C             already open, and will not be closed in this subroutine.
C     M1...   Number of rows of the matrix.
C             Used only for SPARSE='CSC'.
C     M2...   Number of columns of the matrix.
C             Used only for SPARSE='CSC'.
C     SPARSE... Sparseness of the matrix data file:
C             SPARSE=' ':  Matrix is dense (not sparse).
C             SPARSE='CSC' (case insensitive): Matrix is a sparse matrix
C               stored in the compressed column format.
C     NELEM... Number of elements of the matrix to be written.
C     FORM... Form of the matrix data file:
C             FORM='FORMATTED': Formatted file.
C             FORM='UNFORMATTED': Unformatted file.
C     IMAT... Address of the first storage location of the matrix
C             in array RAM (or IRAM) of common block RAMC.
C No output.
C
C     The description of storage of matrices in the memory may be found
C     above.
C     For description of storage of matrices in disk files refer to file
C     forms.htm.
C
C Subroutines and external functions required:
      EXTERNAL LENGTH,UPPER,WARRAI,WMATR
      INTEGER LENGTH
C     LENGTH,UPPER... File 'length.for'.
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
      INTEGER IRAM(MRAM)
      EQUIVALENCE (RAM,IRAM)
C
C Coded by Ludek Klimes and Petr Bulant
C
C-----------------------------------------------------------------------
C
C     Local storage locations:
      INTEGER O,I1,I2
      CHARACTER*72 TXTERR
      CHARACTER*11 FORMU,SPARSU
C     O...    Position of the first matrix element in (I)RAM.
C
C.......................................................................
C
C     Form of the file with the matrix, opening the file:
      FORMU=FORM
      CALL UPPER(FORMU)
      IF(FILED.NE.' ') THEN
        WRITE(*,'(2A)') '+Writing: ',FILED(1:MIN0(LENGTH(FILED),70))
        OPEN(LU,FILE=FILED,FORM=FORMU,ERR=100)
      END IF
C
C     Writing the matrix:
      O=IMAT
      SPARSU=SPARSE
      CALL UPPER(SPARSU)
      IF(SPARSU.EQ.'CSC') THEN
C       Rewriting the array of pointers into the disk form:
        I2=(IMAT-1)+(M2+1)+1
        DO 10, I1=IMAT,IMAT+M2
          IRAM(I1)=(IRAM(I1)-I2)/2
  10    CONTINUE
C       Writing the pointers:
        CALL WARRAI(LU,' ',FORMU,.FALSE.,0,.FALSE.,0,M2,IRAM(IMAT+1))
C       Rewriting the array of pointers into the memory form:
        DO 20, I1=IMAT,IMAT+M2
          IRAM(I1)=IRAM(I1)*2+I2
  20    CONTINUE
        O=IMAT+M2+1
      END IF
      CALL WMATR(LU,M1,M2,SPARSE,NELEM,FORMU,O)
C
C     Closing output file:
      IF(FILED.NE.' ') THEN
        CLOSE(LU)
        WRITE(*,'(1A)')
     *  '+                                                             '
      END IF
      RETURN
  100 CONTINUE
C     MAT-02
      WRITE(TXTERR,'(A,A,A)') 'MAT-02: Error when opening file ''',
     *FILED(1:MIN0(LENGTH(FILED),37)),'''.'
      CALL ERROR(TXTERR)
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE RMATD(LU,FILED,M2,SPARSE,NELEM,FORM,IMAT)
      CHARACTER*(*) FILED,FORM,SPARSE
      INTEGER LU,M2,NELEM,IMAT
C
C Subroutine designed to read a matrix from the matrix data file
C and to store it in common block RAMC of file 'ram.inc'.
C
C Input:
C     LU...   Logical unit number to be used for the input.
C     FILED... Name of the matrix data file. If not blank, the file will
C             be opened and closed.  If blank, the file is assumed to be
C             already open, and will not be closed in this subroutine.
C     M2...   Number of columns of the matrix.
C             Used only for SPARSE='CSC'.
C     SPARSE... Sparseness of the matrix data file:
C             SPARSE=' ':  Matrix is dense (not sparse).
C             SPARSE='CSC' (case insensitive): Matrix is a sparse matrix
C               stored in the compressed column format.
C     NELEM... Number of elements of the matrix to be read.
C     FORM... Form of the matrix data file:
C             FORM='FORMATTED': Formatted file.
C             FORM='UNFORMATTED': Unformatted file.
C     IMAT... Address of the first storage location in array RAM
C             (or IRAM) of common block RAMC where to store the matrix.
C
C No output.
C
C Subroutines and external functions required:
      EXTERNAL LENGTH,UPPER,RARRAI,RARRAY
      INTEGER LENGTH
C     LENGTH,UPPER... File 'length.for'.
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
      INTEGER IRAM(MRAM)
      EQUIVALENCE (RAM,IRAM)
C
C Coded by Ludek Klimes and Petr Bulant
C
C-----------------------------------------------------------------------
C
C     Local storage locations:
      INTEGER O,I,I1,I2
      CHARACTER*72 TXTERR
      CHARACTER*11 FORMU,SPARSU
C     O...    Position of the first matrix element in (I)RAM.
C
C.......................................................................
C
C     Opening the file with matrix data:
      FORMU=FORM
      CALL UPPER(FORMU)
      IF(FILED.NE.' ') THEN
        WRITE(*,'(2A)') '+Reading: ',FILED(1:MIN0(LENGTH(FILED),70))
        OPEN(LU,FILE=FILED,FORM=FORMU,STATUS='OLD',ERR=100)
      END IF
C
C     Reading the matrix:
      O=IMAT
      SPARSU=SPARSE
      CALL UPPER(SPARSU)
      IF(SPARSU.EQ.'CSC') THEN
C       Sparse matrix:
C       Reading the pointers:
        IRAM(O)=0
        CALL RARRAI(LU,' ',FORMU,.FALSE.,0,M2,IRAM(O+1))
C       Rewriting the array of pointers into the memory form:
        I2=(IMAT-1)+(M2+1)+1
        DO 20, I1=IMAT,IMAT+M2
          IRAM(I1)=IRAM(I1)*2+I2
  20    CONTINUE
C       Reading the indices and values:
        O=IMAT+M2+1
        IF(FORMU.EQ.'FORMATTED') THEN
          READ(LU,*) (IRAM(I),RAM(I+1),I=O,O+2*NELEM-1,2)
        ELSE
          READ(LU)   (IRAM(I),RAM(I+1),I=O,O+2*NELEM-1,2)
        END IF
      ELSE
        CALL RARRAY(LU,' ',FORMU,.FALSE.,0.,NELEM,RAM(O))
      END IF
C
C     Closing input file:
      IF(FILED.NE.' ') THEN
        CLOSE(LU)
        WRITE(*,'(1A)')
     *  '+                                                             '
      END IF
      RETURN
  100 CONTINUE
C     MAT-03
      WRITE(TXTERR,'(A,A,A)') 'MAT-03: Error when opening file ''',
     *FILED(1:MIN0(LENGTH(FILED),37)),'''.'
      CALL ERROR(TXTERR)
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE WMATR(LU,M1,M2,SPARSE,NELEM,FORM,IMAT)
      CHARACTER*(*) SPARSE,FORM
      INTEGER LU,M1,M2,NELEM,IMAT
C
C Subroutine designed to write a given array of matrix elements into the
C matrix data file.
C
C Input:
C     LU...   Logical unit number to be used for the output.
C             The file should be already opened, and is not closed after
C             writing.
C     M1...   Number of rows of the matrix.
C             Used only for SPARSE='CSC'.
C     M2...   Number of columns of the matrix.
C             Used only for SPARSE='CSC'.
C     SPARSE... Sparseness of the matrix data file:
C             SPARSE=' ':  Matrix is dense (not sparse).
C             SPARSE='CSC' (case insensitive): Matrix is a sparse matrix
C               stored in the compressed column format.
C     NELEM... Number of matrix elements to be written.
C     FORM... Form of the matrix data file in lowercase:
C             FORM='FORMATTED': Formatted file.
C             FORM='UNFORMATTED': Unformatted file.
C     IMAT... Address of the first storage location in array RAM
C             (or IRAM) of common block RAMC where the matrix data
C             are stored (i.e. first location after pointers for
C             sparse matrix).
C
C No output.
C
C                                                 
C Input SEP parameter:
C     NUMLINM=positive integer... Number of the numbers to be written
C             to each line of the output file. For sparse matrices
C             number of pairs of row indices and values of matrix
C             elements to be written to each line.
C             NUMLINM must be less than 100 (99 at most).
C             Default: NUMLINM=5
C
C Subroutines and external functions required:
      EXTERNAL UPPER
C     UPPER... File 'length.for'.
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
      INTEGER IRAM(MRAM)
      EQUIVALENCE (RAM,IRAM)
C
C Coded by Ludek Klimes and Petr Bulant
C
C-----------------------------------------------------------------------
C
C     Local storage locations:
      INTEGER I
      CHARACTER*20 FORMAT
      CHARACTER*11 FORMU,SPARSU
C     FORMAT... String containing the output format.
C
      INTEGER NUMLIN
      SAVE NUMLIN
      DATA NUMLIN/-1/
C
C.......................................................................
C
C     Setting output format:
      IF(NUMLIN.EQ.-1) THEN
        CALL RSEP3I('NUMLINM',NUMLIN,5)
      ENDIF
      FORMAT='(00(E13.7,1X))'
      FORMAT(3:3)=CHAR(ICHAR('0')+MOD(NUMLIN,10))
      FORMAT(2:2)=CHAR(ICHAR('0')+    NUMLIN/10 )
C
C     Writing the array:
      FORMU=FORM
      CALL UPPER(FORMU)
      SPARSU=SPARSE
      CALL UPPER(SPARSU)
      IF(SPARSU.EQ.'CSC') THEN
C       Sparse matrix:
        IF(FORMU.EQ.'FORMATTED') THEN
C         FORMAT='(00(I0,1X,E13.7,1X))'
          FORMAT(5:20)='I0,1X,E13.7,1X))'
          IF     (M1.GT.999999999) THEN
C           MAT-04
            CALL ERROR
     *      ('MAT-04: Insufficient space in format specification')
          ELSE IF (M1.GT.99999999) THEN
            FORMAT(6:6)='9'
          ELSE IF (M1.GT.9999999) THEN
            FORMAT(6:6)='8'
          ELSE IF (M1.GT.999999) THEN
            FORMAT(6:6)='7'
          ELSE IF (M1.GT.99999) THEN
            FORMAT(6:6)='6'
          ELSE IF (M1.GT.9999) THEN
            FORMAT(6:6)='5'
          ELSE IF (M1.GT.999) THEN
            FORMAT(6:6)='4'
          ELSE IF (M1.GT.99) THEN
            FORMAT(6:6)='3'
          ELSE IF (M1.GT.9) THEN
            FORMAT(6:6)='2'
          ELSE
            FORMAT(6:6)='1'
          END IF
          WRITE(LU,FORMAT)
     *         (IRAM(I),RAM(I+1),I=IMAT,IMAT+2*NELEM-1,2)
        ELSE
          WRITE(LU)
     *         (IRAM(I),RAM(I+1),I=IMAT,IMAT+2*NELEM-1,2)
        END IF
      ELSE
C       Dense (not sparse) matrix:
        IF(FORMU.EQ.'FORMATTED') THEN
          WRITE(LU,FORMAT) (RAM(I),I=IMAT,IMAT+NELEM-1)
        ELSE
          WRITE(LU)        (RAM(I),I=IMAT,IMAT+NELEM-1)
        END IF
      END IF
C
      RETURN
      END
C=======================================================================
C
C     
C
      SUBROUTINE TMATR(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA)
      INTEGER M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA
C
C Subroutine designed to transpose a matrix.
C
C Input:
C     M1 ...  Number of rows of the matrix.
C     M2 ...  Number of columns of the matrix.
C     ISYM... Index of symmetry of the matrix.
C             ISYM=1: SYM='diag'
C             ISYM=2: SYM='sym'
C             ISYM=3: SYM='skew'
C             ISYM=4: SYM=' '
C     NSPAR... Sparseness of the matrix.
C             NSPAR.LT.0: Matrix is stored element by element.
C             NSPAR.GE.0: Matrix is stored as a sparse matrix,
C               NSPAR is the number of zero matrix elements.
C     NELEM... Number of elements of the matrix stored in array RAM.
C     MTMP... Dimension of auxiliary arrays IATMP and ATMP.
C     MIA,MAA... Minimum and maximum address of arrays (I)RAM
C             available for the subroutine.  Entire arrays (I)RAM
C             from MIA to MAA may be used for temporary storage.
C             For NSPAR.LT.0 and general matrix (symmetry=' '),
C               MAA-MIA+1 must be at least 2*NELEM.  For NSPAR.LT.0 and
C               other symmetries, MAA-MIA+1 may equal NELEM.
C             For NSPAR.GE.0, MAA-MIA+1 must be at least NA+NAout.
C     IA...   Address of the first storage location in array (I)RAM
C             used for the matrix.
C     NA...   Number of storage locations for the input matrix.
C
C Output:
C     M1...   Number of rows of the transposed matrix.
C     M2...   Number of columns of the transposed matrix.
C     IA...   If possible, IAout equals IAin. For sparse matrix stored
C             at the end of (I)RAM and NAout.GT.NAin, IAout is set to
C             MAA-NAout+1.
C     NA...   Number of storage locations for the output matrix.
C
C Coded by Petr Bulant
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL ERROR
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C     IRAM,RAM... Indices of elements and elements of the matrix, see
C             the detailed description of storage of matrices in
C             the memory given above.
C
C     Local storage location:
      INTEGER I1,I2,I3,I4,IAOLD,IANEW,IB,NB
      CHARACTER*72 TXTERR
C
C.......................................................................
C
      IF (ISYM.EQ.1.OR.ISYM.EQ.2) THEN
C       Symmetric or diagonal matrix
        RETURN
      ENDIF
      IF (ISYM.EQ.3) THEN
C       Antisymmetric matrix
        IF (NSPAR.GE.0) THEN
C         Sparse matrix
          DO 1, I1=IA+M2+1+1,IA-1+NA,2
            RAM(I1)=-RAM(I1)
  1       CONTINUE
        ELSE
C         Normal (not sparse) matrix
          DO 10, I1=IA,IA-1+NA
            RAM(I1)=-RAM(I1)
  10      CONTINUE
        ENDIF
        RETURN
      ENDIF
C     General matrix:
      IF (NSPAR.GE.0) THEN
C       Transposing sparse matrix:
        NB=NA-M2+M1
        IF (NA+NB.GT.MAA-MIA+1) THEN
C         MAT-05
          WRITE(TXTERR,'(A,I9,A)')
     *    'MAT-05: Array RAM too small,',NA+NB-(MAA-MIA+1),
     *    ' units missing.'
          CALL ERROR(TXTERR)
        ENDIF
        IAOLD=IA
        IF (MAA-(IA+NA-1).LT.NB) THEN
          IANEW=MAA-(NA+NB)+1
          CALL MSHIFT(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,IANEW)
        ENDIF
        IB=IA+NA
C       Computing lengths of columns of transp(A):
        DO 20, I1=0,M1
           IRAM(IB+I1)=0
  20    CONTINUE
        DO 24, I1=1,M2
           DO 22, I2=IRAM(IA-1+I1),IRAM(IA+I1)-2,2
              I3=IB+IRAM(I2)
              IRAM(I3)=IRAM(I3)+1
  22       CONTINUE
  24    CONTINUE
C       Computing pointers of transp(A) from lengths:
        IRAM(IB)=IB+M1+1
        DO 26, I1=1,M1
          IRAM(IB+I1)=IRAM(IB+I1-1)+2*IRAM(IB+I1)
  26    CONTINUE
C       Creating transp(A):
        DO 30, I1=1,M2
          DO 28, I2=IRAM(IA-1+I1),IRAM(IA+I1)-2,2
            I3=IRAM(I2)
            I4=IRAM(IB-1+I3)
            IRAM(I4)=I1
            RAM(I4+1)=RAM(I2+1)
            IRAM(IB-1+I3)=I4+2
  28      CONTINUE
  30    CONTINUE
        DO 32, I1=M1,1,-1
           IRAM(IB+I1)=IRAM(IB-1+I1)
  32    CONTINUE
        IRAM(IB)=IB+M1+1
C       Shifting transposed matrix to the original position:
        IA=IB
        NA=NB
        I1=M1
        M1=M2
        M2=I1
        IF (IAOLD-1+NA.LE.MAA) THEN
          CALL MSHIFT(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,IAOLD)
        ENDIF
      ELSE
C       Transposing normal (not sparse) matrix:
        IF (2*NELEM.GT.MAA-MIA+1) THEN
C         MAT-06
          WRITE(TXTERR,'(A,I9,A)')
     *    'MAT-06: Array RAM too small,',2*NELEM-(MAA-MIA+1),
     *    ' units missing.'
          CALL ERROR(TXTERR)
        ENDIF
        IAOLD=IA
        IF (MAA-(IA+NELEM-1).LT.NELEM) THEN
          IANEW=MAA-2*NELEM+1
          CALL MSHIFT(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,IANEW)
        ENDIF
        I3=IA-1
        I4=IA-1+NELEM
        DO 36, I2=1,M2
          DO 34, I1=1,M1
            RAM(I4+(I1-1)*M2+I2)=RAM(I3+(I2-1)*M1+I1)
  34      CONTINUE
  36    CONTINUE
        IA=I4+1
        I1=M1
        M1=M2
        M2=I1
        CALL MSHIFT(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,IAOLD)
      ENDIF
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE SMATRE(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA)
      INTEGER M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA
C
C Subroutine designed to extend sparse symmetric or antisymmetric
C matrix to the sparse general matrix.
C Sparse matrix to be extended must be stored in form "as on a disk",
C i.e. matrix indices in IARRAY(1) to IARRAY(NELEM) and values of matrix
C elements in ARRAY(NELEM+1) to ARRAY(NELEM+NELEM).
C
C Input:
C     M1...   Number of rows of the matrix.
C     M2...   Number of columns of the matrix.
C     ISYM... Index of symmetry of the matrix.
C             ISYM=1: SYM='diag'
C             ISYM=2: SYM='sym'
C             ISYM=3: SYM='skew'
C             ISYM=4: SYM=' '
C     NSPAR... Sparseness of the matrix. Must be NSPAR.GE.0 on input.
C     NELEM... Number of elements of the matrix stored in array RAM.
C     MIA,MAA... Minimum and maximum address of arrays (I)RAM
C             available for the subroutine.  Entire arrays (I)RAM
C             from MIA to MAA may be used for temporary storage.
C             At least NA+2*NNE storage locations must be available,
C             where NNE is the number of nonzero elements of the
C             extended matrix.
C     IA...   Address of the first storage location in array IRAM
C             used for the matrix.
C     NA...   Number of storage locations for the input matrix.
C
C Output:
C     ISYM... ISYM=4, matrix is extended to general matrix.
C     NSPAR... Sparseness of the extended matrix (NSPAR.GE.0).
C     NELEM... Number of elements of the extended matrix.
C     IA...   Address of the first storage location in array IRAM
C             used for the extended matrix.
C             If IA-1+M2+1+2*NNE.LE.MAA, IA is not changed.
C             Otherwise IA is set to MAA-(M2+1+2*NNE)+1.
C     NA...   Number of storage locations for the extended matrix.
C
C Coded by Petr Bulant
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL ERROR,NELMAT,MSHIFT,VELEM
      INTEGER  NELMAT
      REAL VELEM
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C     IRAM,RAM... Indices of elements and elements of the matrix, see
C             the detailed description of storage of matrices in
C             the memory given above.
C
C     Local storage location:
      INTEGER I1,I2,I3,I4,NNE,IANEW,I2END
      REAL SYMMUL,VEL
      CHARACTER*72 TXTERR
C
C.......................................................................
C
      IF ((NSPAR.LT.0).OR.((ISYM.NE.2).AND.(ISYM.NE.3))) THEN
C       MAT-07
        CALL ERROR('MAT-07: Wrong invocation of SMATRE.')
      ENDIF
C     Calculating the new number of elements in the extended matrix:
      NNE=NELEM
      DO 3, I2=1,M2
        DO 2, I3=IRAM(IA-1+I2),IRAM(IA+I2)-2,2
          I1=IRAM(I3)
          IF (I1.NE.I2) NNE=NNE+1
   2    CONTINUE
   3  CONTINUE
      IF (M2+1+2*NNE.GT.MAA-MIA+1) THEN
C       MAT-08
        WRITE(TXTERR,'(A,I9,A)') 'MAT-08: Array RAM too small,',
     *  M2+1+2*NNE-(MAA-MIA+1),' units missing.'
        CALL ERROR(TXTERR)
      ENDIF
C     Multiplication switch for symmetric or antisymmetric matrix:
      IF (ISYM.EQ.2) THEN
        SYMMUL=1.
      ELSE
        SYMMUL=-1.
      ENDIF
C     Checking the available space for the extended matrix,
C     changing IA and shifting the matrix:
      IF (IA-1+M2+1+2*NNE.GT.MAA) THEN
        IANEW=MAA-(M2+1+2*NNE)+1
        CALL MSHIFT(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,IANEW)
      ENDIF
C     Constructing the extended matrix:
      I4=IA+M2+1+2*NNE
      DO 14, I2=M2,1
C       Address of the old end of the column I2:
        I2END=IRAM(IA-1+I2+1)-2
C       Updating the address of the beginning of column I2+1:
        IRAM(IA-1+I2+1)=I4
C       Recording the values under the diagonal:
        DO 12, I1=M1,I2+1,-1
          VEL=SYMMUL*VELEM(M1,M2,ISYM,NSPAR,IA,I2,I1)
          IF (VEL.NE.0.) THEN
            I4=I4-2
            IRAM(I4)=I1
            RAM(I4+1)=VEL
          ENDIF
  12    CONTINUE
C       Copying the values on and above the diagonal:
        DO 13, I3=I2END,IRAM(IA-1+I2),-2
          I4=I4-2
          IRAM(I4)=IRAM(I3)
          RAM(I4+1)=RAM(I3+1)
  13    CONTINUE
  14  CONTINUE
C
      ISYM=4
      NSPAR=NELMAT(M1,M2,ISYM)-NNE
      NELEM=NNE
      NA=M2+1+2*NNE
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE GSMATR(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA)
      INTEGER M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA
C
C Subroutine designed to change normal (not sparse) matrix to the
C sparse matrix. The sparse matrix is written to the end
C of arrays I(RAM).
C
C Input:
C     M1...   Number of rows of the matrix.
C     M2...   Number of columns of the matrix.
C     ISYM... Index of symmetry of the matrix.
C             ISYM=1: SYM='diag'
C             ISYM=2: SYM='sym'
C             ISYM=3: SYM='skew'
C             ISYM=4: SYM=' '
C     NSPAR... Sparseness of the matrix. Must be NSPAR.LT.0 on input.
C     NELEM... Number of elements of the matrix stored in array RAM.
C     MIA,MAA... Minimum and maximum address of arrays (I)RAM
C             available for the subroutine.  Entire arrays (I)RAM
C             from MIA to MAA may be used for temporary storage.
C             MAA-MIA should be at least M2+1+NAout, where NAout
C             stays for the value of NA on output, but higher value
C             of MAA-MIA is recommended.  MAA-MIA equal to
C             NAin+M2+1+NAout is always sufficient.
C     IA...   Address of the first storage location in array RAM
C             used for the matrix.
C     NA...   Number of storage locations for the input matrix.
C
C Output:
C     NSPAR... Sparseness of the matrix.  NSPAR.GE.0 on output.
C     NELEM... Number of nonzero elements of the matrix
C             in arrays (I)RAM.
C     IA...   New address of the first storage location in array RAM
C             used for the matrix. IA equals MAA-NA+1 on output.
C     NA...   New number of storage locations for the sparse matrix.
C
C Coded by Petr Bulant
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL ERROR,NSPMAT
      INTEGER NSPMAT
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C     IRAM,RAM... Indices of elements and elements of the matrix, see
C             the detailed description of storage of matrices in
C             the memory given above.
C
C     Local storage location:
      INTEGER I1,I2,I11,I12,J1,J2,J3,NN,NSPARN,IANEW
      CHARACTER*72 TXTERR
C
C.......................................................................
C
      IF ((NSPAR.GE.0).OR.(NELEM.NE.NA)) THEN
C       MAT-09
        CALL ERROR('MAT-09: Wrong invocation of GSMATR.')
      ENDIF
C
C     Number NSPARN of zero and NN of nonzero elements:
      NSPARN=NSPMAT(NA,RAM(IA))
      NN=NA-NSPARN
C
C     Moving the matrix to the optimum position:
      IF (IA.NE.MIA+M2+1) THEN
        IANEW=MAX0(MIA+M2+1,(MAA-(2*NN+NA)+1))
        CALL MSHIFT(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,IANEW)
      ENDIF
C
      IF ((IA+2*NN).GT.(MAA)) THEN
C       MAT-10
        I1=IA+2*NN-MAA
        WRITE(TXTERR,'(A,I9,A)')
     *  'MAT-10: Array RAM too small,',I1,' units missing.'
        CALL ERROR(TXTERR)
      ENDIF
C     Moving the nonzero elements of the matrix to the end of RAM,
C     storing the corresponding row indices of the elements,
C     storing the column pointers:
      J1=MAA
      J2=IA-1+NA
      J3=IA-1-(M2+1)
C     J1 points where to write the next element of the sparse matrix
C     J2  I3 points to the next element of the dense matrix
C     J3 points before the array of pointers
C     Pointer for column M2+1:
      IRAM(J3+M2+1)=MAA+1
C     Loop over indices of columns:
      DO 20, I2=M2,1,-1
        IF (ISYM.EQ.1) THEN
          I11=I2
          I12=I2
        ELSEIF (ISYM.EQ.2) THEN
          I11=I2
          I12=1
        ELSEIF (ISYM.EQ.3) THEN
          I11=I2-1
          I12=1
        ELSE
          I11=M1
          I12=1
        ENDIF
C       Loop over indices of rows of the column I2:
        DO 18, I1=I11,I12,-1
          IF (RAM(J2).NE.0.) THEN
            IF (J1-1.LT.J2) THEN
C             MAT-11
              CALL ERROR('MAT-11: Array RAM too small')
C             An attemp to change normal matrix to sparse matrix
C             "in place" failed, more memory is required.
            ENDIF
C           Moving the value, recording the row index:
            RAM(J1)=RAM(J2)
            IRAM(J1-1)=I1
            J1=J1-2
          ENDIF
          J2=J2-1
  18    CONTINUE
C       Writing the pointer for column I2:
        IRAM(J3+I2)=J1+1
  20  CONTINUE
      IF (J2.NE.IA-1) THEN
C       MAT-12
        CALL ERROR('MAT-12: Disorder in GSMATR')
C       This error should not appear.
      ENDIF
C     Updating IA:
      IA=J1+1-(M2+1)
C     Moving the array of pointers to the proper position:
      IF (IA.NE.J3+1) THEN
        DO 30, I1=M2+1,1,-1
          IRAM(IA-1+I1)=IRAM(J3+I1)
  30    CONTINUE
      ENDIF
C     Recording the numbers corresponding to the sparse matrix:
      NA=2*NN+M2+1
      NSPAR=NSPARN
      NELEM=NN
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE SGMATR(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA)
      INTEGER M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA
C
C Subroutine designed to change sparse matrix to the normal (not sparse)
C matrix.
C
C Input:
C     M1...   Number of rows of the matrix.
C     M2...   Number of columns of the matrix.
C     ISYM... Index of symmetry of the matrix.
C             ISYM=1: SYM='diag'
C             ISYM=2: SYM='sym'
C             ISYM=3: SYM='skew'
C             ISYM=4: SYM=' '
C     NSPAR... Sparseness of the matrix.
C     NELEM... Number of elements of the matrix stored in array RAM.
C     MIA,MAA... Minimum and maximum address of arrays (I)RAM
C             available for the subroutine.  Entire arrays (I)RAM
C             from MIA to MAA may be used for temporary storage.
C             There must be always enough memory for columns 1 to K
C             of sparse matrix plus K to M2 of dense matrix.
C             (MAA-MIA+1).GE.(NAin+NAout) is always sufficient.
C     IA...   Address of the first storage location in array IRAM
C             used for the matrix.
C     NA...   Number of storage locations for the input matrix.
C
C Output:
C     NSPAR... Sparseness of the matrix. NSPAR=-1 on output.
C     NELEM... Number of elements of the matrix stored in array RAM.
C     IA...   New address of the first storage location in array RAM.
C             If IAin-1+NAout.LE.MAA, IAout equals IAin. Otherwise
C             IAout equals maximum of (MIA,MAA-NAout-NAin+1).
C     NA...   New number of storage locations for the matrix.
C
C Coded by Petr Bulant
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL ERROR,NELMAT,MSHIFT
      INTEGER NELMAT
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C     IRAM,RAM... Indices of elements and elements of the matrix, see
C             the detailed description of storage of matrices in
C             the memory given above.
C
C     Local storage location:
      INTEGER I1,I2,I3,J1,J2,NNE,IANEW,IAOLD
      CHARACTER*72 TXTERR
C
C.......................................................................
C
      IF (NSPAR.LT.0) RETURN
C
C     Rewriting the matrix from the sparse form
C     to the normal (not sparse) form:
C     Number of elements of the dense matrix:
      NNE=NELMAT(M1,M2,ISYM)
C     Checking the memory, shifting the matrix if necessary:
      IAOLD=IA
      IF (IA-1+NA+NNE.GT.MAA) THEN
        IF (IA.GT.MIA) THEN
          IANEW=MAX0(MIA,MAA-NNE-NA+1)
          CALL MSHIFT(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,IANEW)
        ENDIF
      ENDIF
C     J2 points to the current position in the dense matrix:
      J2=MIN0(MAA,IA-1+NA+NNE)
      DO 20, I2=M2,1,-1
C       J1 points before the beginning of column I2 in the dense matrix:
        IF     (ISYM.EQ.1) THEN
          J1=J2-1
        ELSEIF (ISYM.EQ.2) THEN
          J1=J2-I2
        ELSEIF (ISYM.EQ.3) THEN
          J1=J2-I2+1
        ELSE
          J1=J2-M1
        ENDIF
        IF (J1.LT.IRAM(IA+I2)-1) THEN
C         MAT-13
          WRITE(TXTERR,'(A,I9,A)')
     *    'MAT-13: Array RAM too small,',IRAM(IA+I2)-1-J1,
     *    ' units missing.'
          CALL ERROR(TXTERR)
        ENDIF
        DO 8, I3=J1+1,J2
          RAM(I3)=0.
   8    CONTINUE
C       Rewriting column I2:
        DO 10, I3=IRAM(IA+I2)-1,IRAM(IA+I2-1)+1,-2
          I1=IRAM(I3-1)
          IF (ISYM.EQ.1) THEN
            RAM(J1+1)=RAM(I3)
          ELSE
            RAM(J1+I1)=RAM(I3)
          ENDIF
  10    CONTINUE
        J2=J1
  20  CONTINUE
      IA=J2+1
      NA=NNE
      NELEM=NA
      NSPAR=-1
C     If possible, shifting the matrix to the original IA:
      IF (IA.NE.IAOLD) THEN
        IF (IAOLD-1+NA.LE.MAA) THEN
          CALL MSHIFT(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,IAOLD)
        ENDIF
      ENDIF
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE GSMAT(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,RSPAR)
      INTEGER M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA
      REAL RSPAR
C
C Subroutine designed to calculate the number NSPAR of zero elements
C of a normal (not sparse) matrix, and to change the matrix into the
C sparse matrix if NSPAR/NELEM.GE.RSPAR and the number of storage
C locations required for the sparse matrix is lower than NA.
C If the matrix is changed to sparse, it is written to the end
C of arrays I(RAM).
C
C Input:
C     M1...   Number of rows of the matrix.
C     M2...   Number of columns of the matrix.
C     ISYM... Index of symmetry of the matrix.
C             SYM='diag' ... ISYM=1
C             SYM='sym'  ... ISYM=2
C             SYM='skew' ... ISYM=3
C             SYM=' '    ... ISYM=4
C     NSPAR... Sparseness of the matrix. Must be NSPAR.LT.0 on input.
C     NELEM... Number of elements of the matrix stored in array RAM.
C     MIA,MAA... Minimum and maximum address of arrays (I)RAM
C             available for the subroutine.  Entire arrays (I)RAM
C             from MIA to MAA may be used for temporary storage.
C     IA...   Address of the first storage location in array RAM
C             used for the matrix.
C     NA...   Number of storage locations for the input matrix.
C     RSPAR... Minimum rate of sparseness to change the matrix
C             into the sparse matrix.
C
C Output:
C   For NSPAR.LT.0 on output, only the value of RSPAR is
C   calculated, and there are no other changes on output:
C     RSPAR... Rate of sparseness of the matrix. (Number of zero elements
C             of the matrix divided by the number of all elements.)
C   Otherwise:
C     NSPAR... Number of zero elements of the matrix. NSPAR.GE.0.
C     NELEM... Number of nonzero elements of the matrix.
C     IA...   New address of the first storage location in array IRAM.
C     NA...   New number of storage locations for the matrix.
C     RSPAR... Rate of sparseness of the matrix. (Number of zero elements
C             of the matrix divided by the number of all elements.)
C
C Coded by Petr Bulant
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL ERROR,NSPMAT,GSMATR,MSHIFT
      INTEGER NSPMAT
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C     IRAM,RAM... Indices of elements and elements of the matrix, see
C             the detailed description of storage of matrices in
C             the memory given above.
C
C     Local storage location:
      INTEGER NN
      REAL RSP0
C
C.......................................................................
C
      IF ((NSPAR.GE.0).OR.(NELEM.NE.NA)) THEN
C       MAT-14
        CALL ERROR('MAT-14: Wrong invocation of GSMAT.')
      ENDIF
C     Number of zero elements:
      NN=NSPMAT(NA,RAM(IA))
      RSP0=RSPAR
      RSPAR=FLOAT(NN)/FLOAT(NA)
      IF ((RSPAR.LT.RSP0).OR.(M2+1+2*(NA-NN).GT.NA).OR.
     *    (M2+1+NA.GT.MAA-MIA)) THEN
C       Matrix will stay non-sparse.
        RETURN
      ENDIF
C     Changing the matrix into the sparse matrix:
      CALL GSMATR(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA)
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      REAL FUNCTION VELEM(M1,M2,ISYM,NSPARS,IA,IROW,ICOL)
      INTEGER M1,M2,ISYM,NSPARS,IA,IROW,ICOL
C
C Subroutine designed to calculate Value of ELEMent of a matrix.
C
C Input:
C     M1...   Number of rows of the matrix.
C     M2...   Number of columns of the matrix.
C     ISYM... Index of the symmetry of the matrix.
C     NSPARS... Sparseness of the matrix.
C     IA...   Address of the first storage location of the matrix.
C     IROW... Number of the row of the matrix element.
C     ICOL... Number of the column of the matrix element.
C     For the detailed description of storage of matrices in the memory
C     refer to above.
C
C Output:
C     VELEM... Value of the matrix element.
C
C Coded by Petr Bulant
C
C-----------------------------------------------------------------------
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C     IRAM,RAM... Indices of elements and elements of the matrix, see
C             the detailed description of storage of matrices in
C             the memory given above.
C
C     Local storage location:
      INTEGER IRO,ICO,II1,II2,II3,II21,IR1,IR2,IR3
      REAL RMUL
C.......................................................................
C
      IF (NSPARS.GE.0) THEN
C       Sparse matrix, searching for VELEM by halving intervals:
        VELEM=0.
C       For diagonal of antisymmetric matrix:
        IF ((ISYM.EQ.3).AND.(IROW.EQ.ICOL)) RETURN
        RMUL=1.
        IRO=IROW
        ICO=ICOL
        IF (((ISYM.EQ.2).OR.(ISYM.EQ.3)).AND.(IROW.GT.ICOL)) THEN
C         (Anti)symmetric matrix: Aij=RMUL*Aji
          IRO=ICOL
          ICO=IROW
          IF (ISYM.EQ.3) RMUL=-1.
        ENDIF
C       Searching between beginning and end of column ICO.
C       II1 and II2 are the adresses of the two matrix elements between
C       which we search, IR1 and IR2 are their indices of rows.
        II1=IRAM(IA+ICO-1)
        II2=IRAM(IA+ICO)-2
        IF (II2.GE.II1) THEN
          IR1=IRAM(II1)
          IR2=IRAM(II2)
          IF (IR1.EQ.IRO) THEN
            VELEM=RMUL*RAM(II1+1)
            RETURN
          ENDIF
          IF (IR2.EQ.IRO) THEN
            VELEM=RMUL*RAM(II2+1)
            RETURN
          ENDIF
          IF ((IR1.LT.IRO).AND.(IRO.LT.IR2)) THEN
C           IRO may be between IR1 and IR2, halving the interval:
  10        CONTINUE
            II21=II2-II1
            IF (II21.GT.2) THEN
C             II3 and IR3 is the half of the interval:
              II3=II1+(II21/4)*2
              IR3=IRAM(II3)
              IF (IR3.EQ.IRO) THEN
                VELEM=RMUL*RAM(II3+1)
                RETURN
              ENDIF
              IF (IR3.LT.IRO) THEN
                II1=II3
                IR1=IR3
                GOTO 10
              ELSEIF (IRO.LT.IR3) THEN
                II2=II3
                IR2=IR3
                GOTO 10
              ENDIF
            ENDIF
          ENDIF
        ENDIF
      ELSE
C       Non-sparse matrix:
        IF     (ISYM.EQ.2) THEN
C         'sym'
          IF (IROW.LE.ICOL) THEN
            VELEM=RAM(IA-1+(ICOL-1)*ICOL/2+IROW)
          ELSE
            VELEM=RAM(IA-1+(IROW-1)*IROW/2+ICOL)
          ENDIF
        ELSEIF (ISYM.EQ.3) THEN
C         'skew'
          IF (IROW.LT.ICOL) THEN
            VELEM=RAM(IA-1+(ICOL-1)*(ICOL-2)/2+IROW)
          ELSEIF (IROW.EQ.ICOL) THEN
            VELEM=0.
          ELSE
            VELEM=-RAM(IA-1+(IROW-1)*(IROW-2)/2+ICOL)
          ENDIF
        ELSEIF (ISYM.EQ.1) THEN
C         'diag'
          IF (IROW.EQ.ICOL) THEN
            VELEM=RAM(IA-1+ICOL)
          ELSE
            VELEM=0.
          ENDIF
        ELSE
C         ' '
          VELEM=RAM(IA-1+(ICOL-1)*M1+IROW)
        ENDIF
      ENDIF
      RETURN
      END
C
C=======================================================================
C
C     
C
      INTEGER FUNCTION NELMAT(M1,M2,ISYM)
      INTEGER M1,M2,ISYM
C
C Integer function to calculate number of elements of non-sparse matrix.
C
C Input:
C     M1...   Number of rows of the matrix.
C     M2...   Number of columns of the matrix.
C     ISYM... Index of the symmetry of the matrix.
C
C Output:
C     NELMAT... Number of elements of M1*M2 non-sparse matrix
C             of given symmetry. Note that for indices corresponding to
C             symmetries 'diag', 'sym' and 'skew' number M2 is used
C             to calculate NELMAT.
C
C Coded by Petr Bulant
C
C-----------------------------------------------------------------------
C
      IF     (ISYM.EQ.4) THEN
C       General matrix
        NELMAT=M1*M2
      ELSEIF (ISYM.EQ.2) THEN
C       Symmetric matrix
        NELMAT=M2*(M2+1)/2
      ELSEIF (ISYM.EQ.3) THEN
C       Skew matrix
        NELMAT=M2*(M2-1)/2
      ELSEIF (ISYM.EQ.1) THEN
C       Diagonal matrix
        NELMAT=M2
      ELSE
C       MAT-15
        CALL ERROR('MAT-15: Wrong index of matrix symmetry.')
C       Input argument ISYM should be equal to 1, 2, 3 or 4.
      ENDIF
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      INTEGER FUNCTION NSPMAT(NELEM,ARRAY)
      INTEGER NELEM
      REAL ARRAY(*)
C
C Integer function to calculate number of zero elements of non-sparse
c matrix.
C
C Input:
C     NELEM... Number of elements of the matrix.
C     ARRAY... Elements of the matrix.
C
C Output:
C     NSPMAT... Number of zero elements of the matrix.
C
C Coded by Petr Bulant
C
      INTEGER I
C-----------------------------------------------------------------------
C
      NSPMAT=0
      DO 10, I=1,NELEM
        IF (ARRAY(I).EQ.0.) NSPMAT=NSPMAT+1
  10  CONTINUE
C
      RETURN
      END
C
C
C=======================================================================
C
C     
C
      INTEGER FUNCTION ISYM(SYM)
      CHARACTER*(*) SYM
C
C Integer function to assign the integer corresponding to the symmetry.
C
C Input:
C     SYM...  Symmetry of a matrix.
C
C Output:
C     ISYM... Integer number corresponding to the symmetry SYM.
C             SYM='diag' ... ISYM=1
C             SYM='sym'  ... ISYM=2
C             SYM='skew' ... ISYM=3
C             SYM=' '    ... ISYM=4
C
C Coded by Petr Bulant
C
C-----------------------------------------------------------------------
C
      IF     (SYM.EQ.'diag') THEN
        ISYM=1
      ELSEIF (SYM.EQ.'sym' ) THEN
        ISYM=2
      ELSEIF (SYM.EQ.'skew') THEN
        ISYM=3
      ELSEIF (SYM.EQ.' '   ) THEN
        ISYM=4
      ELSE
C       MAT-16
        CALL ERROR('MAT-16: Wrong matrix symmetry.')
C       Input argument SYM should be equal
C       to one of strings ' ', 'sym', 'skew', 'diag'.
      ENDIF
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE MSHIFT(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,IANEW)
      INTEGER M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,IANEW
C
C Subroutine designed to shift the matrix in arrays (I)RAM
C to the new position IANEW.
C
C Input:
C     M1...   Number of rows of the matrix.
C     M2...   Number of columns of the matrix.
C     ISYM... Index of symmetry of the matrix.
C             ISYM=1: SYM='diag'
C             ISYM=2: SYM='sym'
C             ISYM=3: SYM='skew'
C             ISYM=4: SYM=' '
C     NSPAR... Sparseness of the matrix.
C     NELEM... Number of elements of the matrix stored in arrays (I)RAM.
C     MIA,MAA... Minimum and maximum address of arrays (I)RAM
C             available for the subroutine.
C             Only the storage locations between IA and IA+NA-1,
c             and between IANEW and IANEW+NA-1 are altered.
C     IA...   Address of the first storage location in arrays (I)RAM
C             used for the matrix.
C     NA...   Number of storage locations for the matrix.
C     IANEW... Address of the first storage location in arrays (I)RAM
C             where the matrix is to be shifted.
C
C Output:
C     IA...   IA is changed to the value of IANEW.
C
C Coded by Petr Bulant
C
C-----------------------------------------------------------------------
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C     IRAM,RAM... Indices of elements and elements of the matrix, see
C             the detailed description of storage of matrices in
C             the memory given above.
C
C     Local storage location:
      INTEGER II,ISHIFT
C
C.......................................................................
C
      IF (IANEW.EQ.IA) RETURN
C
      IF ((IANEW.LT.MIA).OR.(IANEW+NA-1.GT.MAA).OR.
     *    (IA.LT.MIA).OR.(IA+NA-1.GT.MAA)) THEN
C       MAT-17
        CALL ERROR('MAT-17: Wrong invocation of MSHIFT.')
C       Storage locations IA to IA+NA-1 and
C       IANEW to IANEW+NA-1 must fit into the available memory.
      ENDIF
C
      IF (NSPAR.LT.0) THEN
C       Non-sparse matrix:
        IF (IANEW.LT.IA) THEN
          DO 10, II=0,NA-1
            RAM(IANEW+II)=RAM(IA+II)
  10      CONTINUE
        ELSE
          DO 11, II=NA-1,0,-1
            RAM(IANEW+II)=RAM(IA+II)
  11      CONTINUE
        ENDIF
      ELSE
C       Sparse matrix, form CSC:
        ISHIFT=IANEW-IA
        IF (IANEW.LT.IA) THEN
          DO 30, II=0,M2
            IRAM(IANEW+II)=IRAM(IA+II)+ISHIFT
  30      CONTINUE
          DO 31, II=0,NELEM-1
            IRAM(IANEW+M2+1+II*2)  =IRAM(IA+M2+1+II*2)
             RAM(IANEW+M2+1+II*2+1)= RAM(IA+M2+1+II*2+1)
  31      CONTINUE
        ELSE
          DO 32, II=NELEM-1,0,-1
             RAM(IANEW+M2+1+II*2+1)= RAM(IA+M2+1+II*2+1)
            IRAM(IANEW+M2+1+II*2)  =IRAM(IA+M2+1+II*2)
  32      CONTINUE
          DO 33, II=M2,0,-1
            IRAM(IANEW+II)=IRAM(IA+II)+ISHIFT
  33      CONTINUE
        ENDIF
      ENDIF
C
      IA=IANEW
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE OMAT(LU,FILE,IRW,FORMM)
      CHARACTER*(*) FILE,FORMM
      INTEGER LU,IRW
C
C Subroutine for backward compatibility.
C
C Coded by Petr Bulant
C
C-----------------------------------------------------------------------
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE WMAT(LU,FILE,M1,M2,OUT)
      CHARACTER*(*) FILE
      INTEGER LU,M1,M2
      REAL OUT(*)
C
C Subroutine for backward compatibility.  WMATH and WMATD should be used
C instead of using WMAT !!!
C Subroutine designed to write a given matrix into the file.
C
C Input:
C     LU...   Logical unit number to be used for the output.
C     FILE... Destination header file name.  Must not be blank.
C     M1...   Number of rows of the given matrix.
C     M2...   M2=0 for a symmetric matrix,
C             M2=-1 for a diagonal matrix,
C             M2=number of columns for a general matrix.
C     OUT...  Components of the given matrix stored columnwise.
C             For a symmetric matrix, just components from the first row
C             to the diagonal are stored for each column, i.e., array
C             OUT has M1*(M1+1)/2 matrix components.
C             For a diagonal matrix, just M1 diagonal components are
C             stored.
C
C No output.
C
C Coded by Ludek Klimes and Petr Bulant
C
C-----------------------------------------------------------------------
C
C     Local storage locations:
C
      EXTERNAL ISYM,NELMAT,WMATD,WMATH,ERROR
      INTEGER  ISYM,NELMAT
      CHARACTER*13 FORMM
      INTEGER I1,I2,NELEM
C
C     FORMM ..Form of the files with matrices.
C
C.......................................................................
C
      CHARACTER*80 FILED
      CHARACTER*13 FORMAT
      CHARACTER*4 SYMM
C
      IF (FILE.EQ.' ') THEN
C       MAT-18
        CALL ERROR('MAT-18: Matrix header file not given')
C       In this version of WMAT, the matrix header file
C       must be specified and sequential writing of matrix data file
C       is not allowed.
      ENDIF
      FORMM=' '
      FILED=' '
      IF(M2.EQ.0) THEN
C       Symmetric matrix
        I2=M1
        SYMM='sym'
      ELSEIF(M2.EQ.-1) THEN
C       Diagonal matrix
        I2=M1
        SYMM='diag'
      ELSE
C       General matrix
        I2=M2
        SYMM=' '
      END IF
      NELEM=NELMAT(M1,I2,ISYM(SYMM))
      CALL WMATH(LU,FILE,FILED,M1,I2,' ',NELEM,SYMM,FORMM)
C
      FORMAT='(5(G14.7,1X))'
C     Writing the matrix:
      OPEN(LU,FILE=FILED,FORM=FORMM)
      IF(M2.EQ.0) THEN
C       Symmetric matrix
          IF (FORMM.EQ.'FORMATTED') THEN
            DO 11 I2=1,M1
              WRITE(LU,FORMAT) (OUT(I1),I1=I2*(I2-1)/2+1,I2*(I2+1)/2)
  11        CONTINUE
          ELSE
            WRITE(LU) (OUT(I1),I1=1,M1*(M1+1)/2)
          ENDIF
      ELSE
C       Diagonal or general matrix
        IF (FORMM.EQ.'FORMATTED') THEN
          DO 12 I2=M1,M1*IABS(M2),M1
            WRITE(LU,FORMAT) (OUT(I1),I1=I2-M1+1,I2)
  12      CONTINUE
        ELSE
          WRITE(LU) (OUT(I1),I1=1,M1*IABS(M2))
        ENDIF
      END IF
      CLOSE(LU)
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE RMAT(LU,FILE,M1,M2,ARRAY)
      CHARACTER*(*) FILE
      INTEGER LU,M1,M2
      REAL ARRAY(*)
C
C Subroutine for backward compatibility.  RMATH and RMATD should be used
C instead of using RMAT !!!
C Subroutine designed to read a matrix from the file.
C In this version of RMAT sequential reading of matrix data file
C is not allowed.
C
C Input:
C     LU...   Logical unit number to be used for the input.
C     FILE... Destination header file name.  Must not be blank.
C     M1...   Number of rows of the matrix.
C     M2...   M2=0 for a symmetric matrix,
C             M2=1 for a diagonal matrix,
C             M2=number of columns for a general matrix.
C
C Output:
C     ARRAY... Components of the given matrix stored columnwise.
C             For a symmetric matrix, just components from the first row
C             to the diagonal are stored for each column, i.e., array
C             ARRAY has M1*(M1+1)/2 matrix components.
C             For a diagonal matrix, just M1 diagonal components are
C             stored.
C
C Coded by Ludek Klimes and Petr Bulant
C
C-----------------------------------------------------------------------
C
      EXTERNAL RMATD,RMATH,ERROR
C     Local storage location:
      CHARACTER*13 FORMM,SPARSE
      CHARACTER*80 FILED
      CHARACTER*4 SYMM
      INTEGER I,M1R,M2R,NELEM
C
      IF (FILE.EQ.' ') THEN
C       MAT-19
        CALL ERROR('MAT-19: Matrix header file not given')
C       In this version of RMAT, the matrix header file
C       must be specified and sequential reading of matrix data file
C       is not allowed.
      ENDIF
      CALL RMATH(LU,FILE,FILED,M1R,M2R,SPARSE,NELEM,SYMM,FORMM)
      IF ((M1.NE.M1R).OR.(SPARSE.NE.' ').OR.(SYMM.EQ.'skew').OR.
     *    ((SYMM.EQ.'sym' ).AND.(M2.NE.0)).OR.
     *    ((SYMM.EQ.' '   ).AND.(M2.NE.M2R)).OR.
     *    ((SYMM.EQ.'diag').AND.(M2.NE.1))) THEN
C       MAT-20
        CALL ERROR('MAT-20: Unexpected values in matrix header file.')
C       Some of the values of the matrix header file differs from
C       the values estimated using input parameters M1 and M2.
C       For example, this version of RMAT cannot deal with sparse
C       matrices.
      ENDIF
C
      OPEN(LU,FILE=FILED,FORM=FORMM)
C     Reading the matrix:
      IF(M2.LE.0) THEN
C       Symmetric matrix
        IF (FORMM.EQ.'FORMATTED') THEN
          READ(LU,*) (ARRAY(I),I=1,M1*(M1+1)/2)
        ELSE
          READ(LU)   (ARRAY(I),I=1,M1*(M1+1)/2)
        ENDIF
      ELSE
C       Diagonal or general matrix
        IF (FORMM.EQ.'FORMATTED') THEN
          READ(LU,*) (ARRAY(I),I=1,M1*M2)
        ELSE
          READ(LU)   (ARRAY(I),I=1,M1*M2)
        ENDIF
      END IF
      CLOSE(LU)
C
      RETURN
      END
C
C=======================================================================
C