C 
C Subroutine file 'mat.for' with subroutines dealing with matrices.
C
C Version: 6.10
C Date: 2007, June 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 the array
C             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     SMATR...Subroutine to reorganize sparse matrix in the memory from
C             the form "as on a disk" to the form "as in the memory"
C             or back.
C             SMATR
C     RAMATR..Subroutine to reorganize indices and values of elements
C             of sparse matrix in the memory from the form "as on
C             a disk" to the form "as in the memory" or back.
C             Auxiliary subroutine to SMATR.
C             RAMATR
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     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     Normal (not sparse) matrices:
C             Normal matrices are usually identified by the value of
C             NSPAR=-1.   For normal matrices all the values
C             of the matrix are stored columnwise, starting with
C             from the first value of the first column to the last value
C             of the first column, then the second column, and so on
C             to the value of the last element of the last column.
C             If the matrix is other than general, only the elements
C             corresponding to the given symmetry are stored, e.g. for
C             diagonal matrix only the values of the diagonal elements
C             are stored.
C     Sparse matrices:
C             Sparse matrices are usually identified by the value of
C             NSPAR nonnegative. NSPAR is then the number of zero
C             elements of the matrix, NSPAR=0 means that the matrix
C             is stored as sparse, but it has all elements nonzero.
C             There are two forms of storage of the sparse matrices,
C             described usually by parameter ISTOR:
C       ISTOR=0 ... form "as on a disk": in this form the first NN
C             storage locations of the matrix contain matrix indices
C             of the nonzero elements of the matrix written in integer
C             array. Matrix index of an element of a matrix is
C             the sequence number, in which the element would be written
C             in the general matrix.  This means that the matrix index
C             is influenced only by the number of rows and of columns
C             of the matrix, not by the symmetry of the matrix.
C             The next NN storage locations contain the values of the
C             nonzero elements of the matrix written in real array.
C             In this form the matrix occupies 2*NN storage locations.
C       ISTOR=1 ... form "as in the memory": in this form the first
C             storage location contains the address in the array
C             just before the first nonzero element, this address
C             may be understood as the address of the end of column
C             number zero.  Next M2 storage locations contain
C             the addresses of the ends of all the columns
C             of the matrix.  The nonzero elements of column J
C             of the matrix are thus stored in
C             ARRAY(IARRAY(OARRAY+J-1)+1) to ARRAY(IARRAY(OARRAY+J)),
C             where OARRAY is the address of the first storage location
C             of the matrix.  Next 2*NN storage locations contain the
C             nonzero elements of the matrix stored columnwise. Each
C             element occupies two storage locations, the first one
C             containing index of the row of the element written
C             in the integer array, the second one containing the value
C             of the element written in the real array.  In this form
C             the matrix occupies 1+M2+2*NN storage locations.
C
C For detailed description of the forms of the files with matrices
C and of the parameters describing the matrices
C refer to file forms.htm.
C
C=======================================================================
C
C     
C
      SUBROUTINE WMATH(LU,FILEH,FILED,M1,M2,NSPARS,SYMM,FORM)
      INTEGER LU,M1,M2,NSPARS
      CHARACTER*(*) FILEH,FILED,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     M1...   Number of rows of the matrix.
C     M2...   Number of columns of the matrix.
C     NSPARS... Sparseness of the matrix data file:
C             NSPARS.LT.0: Matrix is stored element by element.
C             NSPARS.GE.0: Matrix is stored as a sparse matrix.
C               NSPARS is the number of zero matrix elements.
C     SYMM... Symmetry of the matrix data file:
C             SYMMETRY=' ': General matrix.
C             SYMMETRY='sym': Symmetric matrix.
C             SYMMETRY='skew': Skew matrix.
C             SYMMETRY='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 and returned.
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 and returned.
C             Otherwise, the input value of FORM is converted to
C             lowercase and returned.
C
C Subroutines and external functions required:
      EXTERNAL LENGTH,LOWER,WSEP3I,WSEP3T,RSEP3T
      INTEGER  LENGTH
C     LENGTH,LOWER... 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 LOWER(FORM)
      END IF
      IF (FORM.NE.'formatted'.AND.FORM.NE.'unformatted') THEN
C       MAT-01
        CALL ERROR('MAT-01: Wrong value of FORMM.')
C       Input parameter FORMM, if specified,
C       must equal either 'formatted' or 'unformatted'.
      ENDIF
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 WSEP3I(LU,'SPARSE',NSPARS)
      CALL WSEP3T(LU,'SYMMETRY',SYMM)
      CALL WSEP3T(LU,'FORM',FORM)
C
      CLOSE(LU)
      RETURN
  100 CONTINUE
C     MAT-27
      WRITE(TXTERR,'(A,A,A)') 'MAT-27: 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,NSPARS,SYMM,FORM,NELEM)
      INTEGER LU,M1,M2,NSPARS,NELEM
      CHARACTER*(*) FILEH,FILED,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     M1...   Number of rows of the matrix.
C     M2...   Number of columns of the matrix.
C     NSPARS... Sparseness of the matrix data file:
C             NSPARS.LT.0: Matrix is stored element by element.
C             NSPARS.GE.0: Matrix is stored as a sparse matrix.
C               NSPARS is the number of zero matrix elements.
C     SYMM... Symmetry of the matrix data file:
C             SYMMETRY=' ': General matrix.
C             SYMMETRY='sym': Symmetric matrix.
C             SYMMETRY='skew': Skew matrix.
C             SYMMETRY='diag': Diagonal matrix.
C     FORM... Form of the matrix data file:
C             FORM='formatted': Formatted file.
C             FORM='unformatted': Unformatted file.
C     NELEM.. Number of elements of the matrix.
C
C Subroutines and external functions required:
      EXTERNAL LENGTH,SSEP,RSEP1,RSEP3I,RSEP3T,ISYM,NELMAT
      INTEGER  LENGTH,ISYM,NELMAT
C     LENGTH... 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
C     Number of SEP parameter sets (NSEP=1 is the input history file):
      INTEGER NSEP
      SAVE NSEP
      DATA NSEP/1/
C
C     Switching to a new parameter set and reading the matrix header:
      NSEP=NSEP+1
      CALL SSEP(NSEP)
      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 RSEP3I('SPARSE',NSPARS,-1)
      CALL RSEP3T('SYMMETRY',SYMM,' ')
      CALL RSEP3T('FORM',FORM,'formatted')
C
C     Converting input strings to lowercase:
      CALL LOWER(SYMM)
      CALL LOWER(FORM)
C
C     Number NELEM of stored matrix elements:
      ISYMM=ISYM(SYMM)
      NELEM=NELMAT(M1,M2,ISYMM)
      IF (NSPARS.GE.0) NELEM=NELEM-NSPARS
C
      CALL SSEP(1)
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE WMATD(LU,FILED,NSPARS,FORM,NELEM,IARRAY,ARRAY)
      CHARACTER*(*) FILED,FORM
      INTEGER LU,NSPARS,NELEM,IARRAY(*)
      REAL ARRAY(*)
C
C Subroutine designed to write a given matrix into the matrix data file.
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     NSPARS... Sparseness of the matrix data file:
C             NSPARS.LT.0: Matrix is stored element by element.
C             NSPARS.GE.0: Matrix is stored as a sparse matrix.
C               NSPARS is the number of zero matrix elements.
C     FORM... Form of the matrix data file:
C             FORM='formatted': Formatted file.
C             FORM='unformatted': Unformatted file.
C     NELEM.. Number of elements of the matrix.
C     IARRAY..Indices of the non-zero matrix elements if the matrix is
C             to be stored as a sparse matrix.
C             The matrix elements are indexed from 1 to N.  For N refer
C             to the description of argument ARRAY below.
C             The number of non-zero matrix elements is NELEM=N-NSPARS.
C             Integer array IARRAY is not used if NSPARS.LT.0.
C     ARRAY.. Elements of the given matrix stored columnwise.
C             For NSPARS.LT.0, all N matrix elements will be stored:
C               For a general matrix, ARRAY has N=M1*M2 components.
C               For a symmetric matrix, just elements from the first row
C                 to the diagonal are stored for each column, i.e.,
C                 ARRAY has N=M1*(M1+1)/2 components.
C               For a skew matrix, just elements above the diagonal are
C                 stored for each column, i.e., ARRAY has N=M1*(M1-1)/2
C                 components.
C               For a diagonal matrix, just N=M1 diagonal elements are
C                 to be stored.
C             For NSPARS.GE.0, NELEM=N-NSPARS indices of non-zero matrix
C               elements from IARRAY, and the corresponding NELEM
C               matrix elements from ARRAY will be stored.
C               The matrix elements are assumed to be located in array
C               ARRAY in the positions from ARRAY(NELEM+1) to
C               ARRAY(NELEM+NELEM).
C
C No output.
C
C Subroutines and external functions required:
      EXTERNAL LENGTH,LOWER,WARRAI,WMATR
C     LENGTH,LOWER... File 'length.for'.
C
C Coded by Ludek Klimes and Petr Bulant
C
C-----------------------------------------------------------------------
C
C     Local storage locations:
      INTEGER N,O
      CHARACTER*72 TXTERR
C     N...    Number of matrix elements to be stored.
C     O...    Position of the first matrix element in ARRAY.
C
C.......................................................................
C
C     Form of the file with the matrix, opening the file:
      IF(FILED.NE.' ') THEN
        WRITE(*,'(2A)') '+Writing: ',FILED(1:MIN0(LENGTH(FILED),70))
        OPEN(LU,FILE=FILED,FORM=FORM,ERR=100)
      END IF
C
      N=NELEM
C
C     Writing the matrix:
      O=1
      IF(NSPARS.GE.0) THEN
        CALL WARRAI(LU,' ',FORM,.FALSE.,0,.FALSE.,0,N,IARRAY)
        O=N+1
      END IF
      CALL WMATR(LU,FORM,N,ARRAY(O))
C
C     Closing output file:
      IF(FILED.NE.' ') THEN
        CLOSE(LU)
        WRITE(*,'(1A)')
     *  '+                                                             '
      END IF
      RETURN
  100 CONTINUE
C     MAT-28
      WRITE(TXTERR,'(A,A,A)') 'MAT-28: Error when opening file ''',
     *FILED(1:MIN0(LENGTH(FILED),37)),'''.'
      CALL ERROR(TXTERR)
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE RMATD(LU,FILED,NSPARS,FORM,NELEM,IARRAY,ARRAY)
      CHARACTER*(*) FILED,FORM
      INTEGER LU,NSPARS,NELEM,IARRAY(*)
      REAL ARRAY(*)
C
C Subroutine designed to read a matrix from the matrix data file.
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     NSPARS... Sparseness of the matrix data file:
C             NSPARS.LT.0: Matrix is stored element by element.
C             NSPARS.GE.0: Matrix is stored as a sparse matrix.
C               NSPARS is the number of zero matrix elements.
C     FORM... Form of the matrix data file:
C             FORM='formatted': Formatted file.
C             FORM='unformatted': Unformatted file.
C     NELEM.. Number of elements of the matrix.
C
C Output:
C     IARRAY..Indices of the non-zero matrix elements if the matrix is
C             stored as a sparse matrix.
C             The matrix elements are indexed from 1 to N.  For N refer
C             to the description of argument ARRAY below.
C             The number of non-zero matrix elements is NELEM=N-NSPARS.
C             Integer array IARRAY is not used if NSPARS.LT.0.
C     ARRAY.. Elements of the given matrix stored columnwise.
C             For NSPARS.LT.0, all N matrix elements are stored:
C               For a general matrix, ARRAY has N=M1*M2 components.
C               For a symmetric matrix, just elements from the first row
C                 to the diagonal are stored for each column, i.e.,
C                 ARRAY has N=M1*(M1+1)/2 components.
C               For a skew matrix, just elements above the diagonal are
C                 stored for each column, i.e., ARRAY has N=M1*(M1-1)/2
C                 components.
C               For a diagonal matrix, just N=M1 diagonal elements are
C                 stored.
C             For NSPARS.GE.0, NELEM=N-NSPARS indices of non-zero matrix
C               elements are read to the array IARRAY, and the
C               corresponding NELEM matrix elements are read to
C               the array ARRAY. The matrix elements are read in array
C               ARRAY into the positions from ARRAY(NELEM+1) to
C               ARRAY(NELEM+NELEM).
C
C Subroutines and external functions required:
      EXTERNAL LENGTH,LOWER,RARRAI,RARRAY
C     LENGTH,LOWER... File 'length.for'.
C
C Coded by Ludek Klimes and Petr Bulant
C
C-----------------------------------------------------------------------
C
C     Local storage locations:
      INTEGER N,O
      CHARACTER*72 TXTERR
C     N...    Number of stored matrix elements.
C     O...    Position of the first matrix element in ARRAY.
C
C.......................................................................
C
C     Opening the file with matrix data:
      IF(FILED.NE.' ') THEN
        WRITE(*,'(2A)') '+Reading: ',FILED(1:MIN0(LENGTH(FILED),70))
        OPEN(LU,FILE=FILED,FORM=FORM,STATUS='OLD',ERR=100)
      END IF
C
C     Reading the matrix:
      N=NELEM
      O=1
      IF(NSPARS.GE.0) THEN
        CALL RARRAI(LU,' ',FORM,.FALSE.,0,N,IARRAY)
        O=N+1
      END IF
      CALL RARRAY(LU,' ',FORM,.FALSE.,0.,N,ARRAY(O))
C
C     Closing input file:
      IF(FILED.NE.' ') THEN
        CLOSE(LU)
        WRITE(*,'(1A)')
     *  '+                                                             '
      END IF
      RETURN
  100 CONTINUE
C     MAT-29
      WRITE(TXTERR,'(A,A,A)') 'MAT-29: Error when opening file ''',
     *FILED(1:MIN0(LENGTH(FILED),37)),'''.'
      CALL ERROR(TXTERR)
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE WMATR(LU,FORM,N,ARRAY)
      CHARACTER*(*) FORM
      INTEGER LU,N
      REAL ARRAY(N)
C
C Subroutine designed to write a given array of matrix elements into the
C 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     FORM... Form of the matrix data file in lowercase:
C             FORM='formatted': Formatted file.
C             FORM='unformatted': Unformatted file.
C     N...    Number of matrix elements.
C     ARRAY...Array of matrix elements.
C
C No output.
C
C                                                  
C Input SEP parameter:
C     NUMLIN=positive integer... Number of the numbers to be written
C             to each line of the output file.
C             NUMLIN must be less than 100 (99 at most).
C             Default: NUMLIN=5
C
C Coded by Ludek Klimes and Petr Bulant
C
C-----------------------------------------------------------------------
C
C     Local storage locations:
      CHARACTER*14 FORMAT
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('NUMLIN',NUMLIN,5)
      ENDIF
      FORMAT='(00(G13.7,1X))'
      FORMAT(3:3)=CHAR(ICHAR('0')+MOD(NUMLIN,10))
      FORMAT(2:2)=CHAR(ICHAR('0')+    NUMLIN/10 )
C
C     Writing the array:
      IF (FORM.EQ.'formatted') THEN
        WRITE(LU,FORMAT) ARRAY
      ELSE
        WRITE(LU)        ARRAY
      ENDIF
C
      RETURN
      END
C=======================================================================
C
C     
C
      SUBROUTINE TMATR(M1,M2,SYM,NSPAR,NELEM,IARRAY,ARRAY,
     *                 IATMP,ATMP,MTMP)
      CHARACTER*(*) SYM
      INTEGER M1,M2,NSPAR,NELEM,MTMP
      REAL     ARRAY(*),ATMP(*)
      INTEGER IARRAY(*),IATMP(*)
C
C Subroutine designed to transpose a matrix.
C Sparse matrix to be transposed must be stored in form "as on a disk".
C
C Input:
C     M1...   Number of rows of the matrix.
C     M2...   Number of columns of the matrix.
C     SYM ... Symmetry of the matrix.
C     NSPAR...Sparseness of the matrix.
C     NELEM...Number of elements of the matrix.
C     IARRAY..Components of the given matrix.
C     ARRAY...Components of the given matrix.
C     For detailed description of the parameters describing the matrices
C     refer to file forms.htm.
C     Detailed description of storage of matrices in the memory
C     is given above.
C     MTMP .. Dimension of auxiliary arrays IATMP and ATMP.
C
C Output:
C     M1...   Number of rows of the transposed matrix.
C     M2...   Number of columns of the transposed matrix.
C     IARRAY..Components of the transposed matrix.
C     ARRAY...Components of the transposed matrix.
C
C Coded by Petr Bulant
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL ERROR,INDEXI
C
C     Local storage location:
      INTEGER I1,I2,I3
      CHARACTER*72 TXTERR
C
C.......................................................................
C
      IF (SYM.EQ.'sym'.OR.SYM.EQ.'diag') THEN
C       Symmetric or diagonal matrix
        RETURN
      ENDIF
      IF (SYM.EQ.'skew') THEN
C       Antisymmetric matrix
        IF (NSPAR.GE.0) THEN
C         Sparse matrix
          DO 1, I1=NELEM+1,2*NELEM
            ARRAY(I1)=-ARRAY(I1)
  1       CONTINUE
        ELSE
C         Normal (not sparse) matrix
          DO 10, I1=1,NELEM
            ARRAY(I1)=-ARRAY(I1)
  10      CONTINUE
        ENDIF
        RETURN
      ENDIF
C     General matrix:
      IF (NSPAR.GE.0) THEN
C       Transposing sparse matrix:
        IF (2*NELEM.GT.MTMP) THEN
C         MAT-02
          WRITE(TXTERR,'(A,I9,A)')
     *    'MAT-02: Array RAM too small,',2*NELEM-MTMP,' units missing.'
          CALL ERROR(TXTERR)
        ENDIF
        DO 2, I3=1,NELEM
          I1=MOD((IARRAY(I3)-1),M1)+1
          I2=(IARRAY(I3)-1)/M1+1
          IARRAY(I3)=(I1-1)*M2+I2
          IATMP(I3)=IARRAY(I3)
          ATMP(NELEM+I3)=ARRAY(NELEM+I3)
   2    CONTINUE
        CALL INDEXI(NELEM,IATMP,IARRAY)
        DO 3, I3=1,NELEM
          ARRAY(NELEM+I3)=ATMP(NELEM+IARRAY(I3))
          IARRAY(I3)=IATMP(IARRAY(I3))
   3    CONTINUE
      ELSE
C       Transposing normal (not sparse) matrix:
        IF (NELEM.GT.MTMP) THEN
C         MAT-03
          WRITE(TXTERR,'(A,I9,A)')
     *    'MAT-03: Array RAM too small,',NELEM-MTMP,' units missing.'
          CALL ERROR(TXTERR)
        ENDIF
        DO 30, I2=1,M2
          DO 20, I1=1,M1
            ATMP((I1-1)*M2+I2)=ARRAY((I2-1)*M1+I1)
  20      CONTINUE
  30    CONTINUE
        DO 40, I1=1,NELEM
          ARRAY(I1)=ATMP(I1)
  40    CONTINUE
      ENDIF
      I1=M1
      M1=M2
      M2=I1
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE SMATRE(M1,M2,SYM,NSPAR,NELEM,
     *                  IARRAY,ARRAY,MARRAY,NARRAY,ISTOR)
      CHARACTER*(*) SYM
      INTEGER M1,M2,NSPAR,NELEM,NARRAY,MARRAY,ISTOR
      REAL     ARRAY(*)
      INTEGER IARRAY(*)
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
C Input:
C     M1...   Number of rows of the matrix.
C     M2...   Number of columns of the matrix.
C     SYM ... Symmetry of the matrix.
C     NSPAR...Sparseness of the matrix.
C     NELEM...Number of elements of the matrix.
C     NARRAY..Number of storage locations for the matrix.
C     IARRAY..Components of the given matrix.
C     ARRAY...Components of the given matrix.
C     ISTOR...Form of the given matrix.
C     For detailed description of the parameters describing the matrices
C     refer to file forms.htm.
C     Detailed description of storage of matrices in the memory
C     is given above.
C     MARRAY..Dimension of array ARRAY.
C
C Output:
C     SYM ... SYM=' '. New symmetry of the matrix (general matrix).
C     NSPAR...Sparseness of the extended matrix.
C     NELEM...Number of elements of the extended matrix.
C     NARRAY..Number of storage locations for the extended matrix.
C     IARRAY..Components of the extended matrix.
C     ARRAY...Components of the extended matrix.
C
C Coded by Petr Bulant
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL ERROR,INDEXI
C
C     Local storage location:
      INTEGER I1,I2,I3,I4,NNE
      REAL SYMMUL
      CHARACTER*72 TXTERR
C
C.......................................................................
C
      IF ((NSPAR.LT.0).OR.((SYM.NE.'sym').AND.(SYM.NE.'skew')).OR.
     *    (MARRAY.LT.NARRAY).OR.(ISTOR.NE.0)) THEN
C       MAT-04
        CALL ERROR('MAT-04: Wrong invocation of SMATRE.')
      ENDIF
C     Calculating the new number of elements in the extended matrix:
      NNE=NELEM
      DO 1, I3=1,NELEM
        I1=MOD((IARRAY(I3)-1),M1)+1
        I2=(IARRAY(I3)-1)/M1+1
        IF (I1.NE.I2) NNE=NNE+1
   1  CONTINUE
      IF (4*NNE.GT.MARRAY) THEN
C       MAT-05
        WRITE(TXTERR,'(A,I9,A)')
     *  'MAT-05: Array RAM too small,',4*NNE-MARRAY,' units missing.'
        CALL ERROR(TXTERR)
      ENDIF
C     Multiplication switch for symmetric or antisymmetric matrix:
      IF (SYM.EQ.'sym') THEN
        SYMMUL=1.
      ELSE
        SYMMUL=-1.
      ENDIF
C     Filling the auxiliary arrays:
      I4=0
      DO 2, I3=1,NELEM
        I4=I4+1
        I1=MOD((IARRAY(I3)-1),M1)+1
        I2=(IARRAY(I3)-1)/M1+1
        IARRAY(2*NNE+I4)=IARRAY(I3)
        ARRAY(3*NNE+I4)=ARRAY(NELEM+I3)
        IF (I1.NE.I2) THEN
          I4=I4+1
          IARRAY(2*NNE+I4)=(I1-1)*M2+I2
          ARRAY(3*NNE+I4)=SYMMUL*ARRAY(NELEM+I3)
        ENDIF
   2  CONTINUE
C     Sorting:
      CALL INDEXI(NNE,IARRAY(2*NNE+1),IARRAY)
C     Copying back the sorted matrix elements:
      DO 3, I3=1,NNE
        ARRAY(NNE+I3)=ARRAY(3*NNE+IARRAY(I3))
        IARRAY(I3)=IARRAY(2*NNE+IARRAY(I3))
   3  CONTINUE
C
      NELEM=NNE
      NARRAY=2*NNE
      SYM=' '
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE SMATR(M1,M2,NSPAR,NELEM,
     *                 IARRAY,ARRAY,MARRAY,OARRAY,NARRAY,ISTOR)
      INTEGER M1,M2,NSPAR,NELEM,OARRAY,NARRAY,MARRAY,ISTOR
      REAL     ARRAY(*)
      INTEGER IARRAY(*)
C
C Subroutine designed to change a sparse matrix from form "as on a disk"
C to form "as in the memory" or back.
C
C Input:
C     M1...   Number of rows of the matrix.
C     M2...   Number of columns of the matrix.
C     NSPAR...Sparseness of the matrix.
C     NELEM...Number of elements of the matrix.
C     OARRAY..Address of the first storage location in array ARRAY
C             to be used for the matrix.
C     NARRAY..Number of storage locations for the matrix.
C     IARRAY..Components of the given matrix.
C     ARRAY...Components of the given matrix.
C     ISTOR...Form of the given matrix.
C     For detailed description of the parameters describing the matrices
C     refer to file forms.htm.
C     Detailed description of storage of matrices in the memory
C     is given above.
C     MARRAY..Dimension of array ARRAY.
C
C Output:
C     OARRAY..New address of the first storage location in array ARRAY.
C     NARRAY..New number of storage locations for the matrix.
C     IARRAY..Components of the given matrix.
C     ARRAY...Components of the given matrix.
C     ISTOR...New form of the given matrix, opposite to the input.
C
C Coded by Petr Bulant
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL ERROR,RAMATR
C
C     Local storage location:
      INTEGER I1,I2,I3,I4,I20,NN,OAROLD,OARNEW
      CHARACTER*72 TXTERR
C
C.......................................................................
C
      IF (NSPAR.LT.0) THEN
C       MAT-06
        CALL ERROR('MAT-06: Wrong invocation of SMATR.')
      ENDIF
      OAROLD=OARRAY
      IF (ISTOR.EQ.0) THEN
C       Rearranging sparse matrix from the form 'as on a disk'
C       to the form 'as in the memory':
        NN=NELEM
        OARNEW=MARRAY-2*NN-M2
        IF (((OARNEW+M2).LT.(OAROLD+NN)).OR.(OARNEW.LT.1)) THEN
C         MAT-07
          WRITE(TXTERR,'(A,I9,A)')
     *    'MAT-07: Array RAM too small,',OAROLD+NN-OARNEW-M2,
     *    ' units missing.'
          CALL ERROR(TXTERR)
        ENDIF
C       Moving the indices and the values to the new positions:
        DO 2, I3=NN,1,-1
          ARRAY(OARNEW+M2+2*I3)=ARRAY(OAROLD+NN+I3-1)
          IARRAY(OARNEW+M2+2*I3-1)=IARRAY(OAROLD+I3-1)
   2    CONTINUE
C       Calculating end addresses of the columns,
C       changing matrix indices of elements to the row numbers:
        I20=M2+1
        DO 4, I3=NN,1,-1
          I1=MOD((IARRAY(OARNEW+M2+2*I3-1)-1),M1)+1
          I2=(IARRAY(OARNEW+M2+2*I3-1)-1)/M1+1
          IARRAY(OARNEW+M2+2*I3-1)=I1
          IF (I2.NE.I20) THEN
            DO 3, I4=I20-1,I2,-1
              IARRAY(OARNEW+I4)=OARNEW+M2+2*I3
   3        CONTINUE
            I20=I2
          ENDIF
   4    CONTINUE
        DO 5, I4=I20-1,0,-1
          IARRAY(OARNEW+I4)=OARNEW+M2
   5    CONTINUE
        NARRAY=1+M2+2*NN
        ISTOR=1
      ELSE
C       Rearranging sparse matrix from the form 'as in the memory'
C       to the form 'as on a disk':
        OARNEW=1
C       Changing indices of rows to matrix indices:
        DO 20, I2=1,M2
          DO 18, I3=IARRAY(OAROLD+I2-1)+1,IARRAY(OAROLD+I2)-1,2
            I1=IARRAY(I3)
            IARRAY(I3)=(I2-1)*M1+I1
  18      CONTINUE
  20    CONTINUE
        CALL RAMATR(ARRAY,IARRAY,OAROLD+M2+1,2*NELEM,OARNEW,ISTOR)
        NARRAY=2*NELEM
      ENDIF
      OARRAY=OARNEW
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE RAMATR(ARRAY,IARRAY,OA,NA,OANEW,ISTOR)
      INTEGER OA,NA,OANEW,ISTOR
      REAL     ARRAY(*)
      INTEGER IARRAY(*)
C
C Subroutine designed to Reorganize Array with a matrix from the form
C "as on a disk" to the form "as in the memory" or back.  This
C subroutine reorganizes indices and values of the matrix from the form
C "all indices, all values" to the form "index, value, index, value, .."
C or back, it does not change the indices.
C For ISTOR=1 on input the value OA must be less than OANEW.
C For ISTOR=0 on input the value OA must be greater than OANEW.
C
C Input:
C     OA ...  Address of the first storage location in array ARRAY
C             used for the index of the first element of the matrix.
C     OANEW...New address of the first storage location in array ARRAY
C             to be used for the index of the first element
C             of the reorganized matrix.
C     NA ...  Number of storage locations for the indices and values
C             of the matrix elements.
C     ISTOR...Form of the given matrix.
C     IARRAY..Components of the given matrix.
C     ARRAY...Components of the given matrix.
C     Detailed description of storage of matrices in the memory
C     is given above.
C
C Output:
C     IARRAY..Reorganized components of the given matrix.
C     ARRAY...Reorganized components of the given matrix.
C     ISTOR...New form of the given matrix, opposite to the input.
C
C Coded by Petr Bulant
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL ERROR
C
C     Local storage location:
      INTEGER I1,I2,I3,I4,J2,J3,J4,NN
      CHARACTER*72 TXTERR
C
C.......................................................................
C
      NN=NA/2
      IF (ISTOR.EQ.0) THEN
C       IF (I1.GT.J2) THEN
C         MAT-08
C         WRITE(TXTERR,'(A,I9,A)')
C    *    'MAT-08: Array RAM too small,',I1-J2,' units missing.'
C         CALL ERROR(TXTERR)
C       ENDIF
        CALL ERROR('MAT-08: Option ISTOR.EQ.0 is not yet coded.')
      ELSEIF (ISTOR.EQ.1) THEN
C       Rearranging the matrix to the form 'as on a disk':
        IF (OANEW.GE.OA) THEN
C         MAT-09
          CALL ERROR('MAT-09: Wrong invocation of RAMATR')
C         If ISTOR=1, OANEW must be less than OA.
        ENDIF
C       Moving the indices to the new positions, while the values stay
C       on place, and only if there is no place for the indices, then
C       the values are collected and shifted towards the end of ARRAY.
C       J2 points just before the collected values
C       J3 points to the end of the collected values
        J2=OA
        J3=OA+1
        DO 34, I3=1,NN
          I1=OANEW+I3-1
          I2=OA+2*I3-2
          IF (I1.GT.J2) THEN
            J4=I2-2
            DO 31, I4=I2-3,J3,-2
              ARRAY(J4)=ARRAY(I4)
              J4=J4-1
  31        CONTINUE
            DO 32, I4=J3-1,J2+1,-1
              ARRAY(J4)=ARRAY(I4)
              J4=J4-1
  32        CONTINUE
            J3=I2-1
            J2=J4
            IF (I1.GT.J2) THEN
C             MAT-10
              WRITE(TXTERR,'(A,I9,A)')
     *        'MAT-10: Array RAM too small,',I1-J2,' units missing.'
              CALL ERROR(TXTERR)
C             This error should not appear.
            ENDIF
          ENDIF
          IARRAY(I1)=IARRAY(I2)
  34    CONTINUE
C       Indices are moved, now collecting the rest of values
C       and moving them to the end of ARRAY:
        I2=OA+2*NN
        J4=I2-2
        DO 35, I4=I2-3,J3,-2
          ARRAY(J4)=ARRAY(I4)
          J4=J4-1
  35    CONTINUE
        DO 36, I4=J3-1,J2+1,-1
          ARRAY(J4)=ARRAY(I4)
          J4=J4-1
  36    CONTINUE
        J2=J4
C       Moving the values to the new positions:
        DO 38, I3=1,NN
          ARRAY(OANEW+NN+I3-1)=ARRAY(J2+I3)
  38    CONTINUE
        ISTOR=0
      ENDIF
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE GSMATR(M1,M2,SYM,NSPAR,NELEM,ARRAY,IARRAY,MARRAY,
     *                  OARRAY,NARRAY,ISTOR)
      CHARACTER*(*) SYM
      INTEGER M1,M2,NSPAR,NELEM,OARRAY,NARRAY,MARRAY,ISTOR
      REAL     ARRAY(*)
      INTEGER IARRAY(*)
C
C Subroutine designed to change normal (not sparse) matrix to the
C sparse matrix.
C
C Input:
C     M1...   Number of rows of the matrix.
C     M2...   Number of columns of the matrix.
C     SYM ... Symmetry of the matrix.
C     NSPAR...Sparseness of the matrix. Must be NSPAR.LT.0 on input.
C     NELEM...Number of elements of the matrix.
C     OARRAY..Address of the first storage location in array ARRAY
C             used for the matrix.
C     NARRAY..Number of storage locations for the matrix.
C     IARRAY..Components of the given matrix.
C     ARRAY...Components of the given matrix.
C     ISTOR...Form to be used for the given matrix.
C             ISTOR=0 ... form "as on a disk"
C             ISTOR=1 ... form "as in the memory"
C     For detailed description of the parameters describing the matrices
C     refer to file forms.htm.
C     Detailed description of storage of matrices in the memory
C     is given above.
C     MARRAY..Dimension of array ARRAY.
C
C Output:
C     NSPAR...Sparseness of the matrix.
C     NELEM...Number of elements of the matrix.
C     NARRAY..New number of storage locations for the matrix.
C     IARRAY..Components of the given matrix.
C     ARRAY...Components of the given matrix.
C
C Coded by Petr Bulant
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL ERROR,RAMATR,NSPMAT
      INTEGER NSPMAT
C
C     Local storage location:
      INTEGER I1,I2,I3,I4,I20,I11,I12,J1,NN
      CHARACTER*72 TXTERR
C
C.......................................................................
C
      IF ((NSPAR.GE.0).OR.(NELEM.NE.NARRAY)) THEN
C       MAT-11
        CALL ERROR('MAT-11: Wrong invocation of GSMATR.')
      ENDIF
C
C     Number NSPAR of zero and NELEM of nonzero elements:
      NSPAR=NSPMAT(NARRAY,ARRAY(OARRAY))
      NELEM=NARRAY-NSPAR
      NN=NELEM
      IF (((ISTOR.EQ.0).AND.((1+2*NN+1).GT.(MARRAY))).OR.
     *    ((ISTOR.EQ.1).AND.((1+M2+2*NN).GT.(MARRAY)))) THEN
C       MAT-12
        IF (ISTOR.EQ.0) THEN
          I1=1+2*NN+1-MARRAY
        ELSE
          I1=1+M2+2*NN-MARRAY
        ENDIF
        WRITE(TXTERR,'(A,I9,A)')
     *  'MAT-12: Array RAM too small,',I1,' units missing.'
        CALL ERROR(TXTERR)
      ENDIF
C     Moving the nonzero elements of the matrix to the end of ARRAY,
C     storing the corresponding matrix indices of the elements:
      J1=MARRAY
      I3=OARRAY-1+NARRAY
      DO 20, I2=M2,1,-1
        IF (SYM.EQ.'diag') THEN
          I11=I2
          I12=I2
        ELSEIF (SYM.EQ.'sym') THEN
          I11=I2
          I12=1
        ELSEIF (SYM.EQ.'skew') THEN
          I11=I2-1
          I12=1
        ELSE
          I11=M1
          I12=1
        ENDIF
        DO 18, I1=I11,I12,-1
          IF (ARRAY(I3).NE.0.) THEN
            IF (J1-1.LT.I3) THEN
C             MAT-13
              WRITE(TXTERR,'(A,I9,A)')
     *        'MAT-13: Array RAM too small,',I3-J1+1,' units missing.'
              CALL ERROR(TXTERR)
            ENDIF
            ARRAY(J1)=ARRAY(I3)
            IARRAY(J1-1)=(I2-1)*M1+I1
            J1=J1-2
          ENDIF
          I3=I3-1
  18    CONTINUE
  20  CONTINUE
      IF (I3.NE.OARRAY-1) THEN
C       MAT-14
        WRITE(TXTERR,'(A,I9,A)')
     *  'MAT-14: Array RAM too small,',OARRAY-1-I3,' units missing.'
        CALL ERROR(TXTERR)
      ENDIF
C     J1 now points just before the index of the first nonzero element
      IF (ISTOR.EQ.0) THEN
C       Rearranging the matrix to the form 'as on a disk':
        ISTOR=1
        CALL RAMATR(ARRAY,IARRAY,J1+1,2*NN,OARRAY,ISTOR)
        NARRAY=2*NN
      ELSE
C       Rearranging the matrix to the form 'as in the memory':
C       Moving the indices and the values to the new positions:
        DO 42, I3=1,NN
          ARRAY(OARRAY+M2+2*I3)=ARRAY(J1+2*I3)
          IARRAY(OARRAY+M2+2*I3-1)=IARRAY(J1+2*I3-1)
  42    CONTINUE
C       Calculating end addresses of the columns,
C       changing matrix indices of the elements to the row numbers:
        I20=M2+1
        DO 44, I3=NN,1,-1
          I1=MOD(IARRAY(OARRAY+M2+2*I3-1),M1)
          IF (I1.EQ.0) I1=M1
          I2=(IARRAY(OARRAY+M2+2*I3-1)-1)/M1+1
          IARRAY(OARRAY+M2+2*I3-1)=I1
          IF (I2.NE.I20) THEN
            DO 43, I4=I20-1,I2,-1
              IARRAY(OARRAY+I4)=OARRAY+M2+2*I3
  43        CONTINUE
            I20=I2
          ENDIF
  44    CONTINUE
        DO 45, I4=I20-1,0,-1
          IARRAY(OARRAY+I4)=OARRAY+M2
  45    CONTINUE
        NARRAY=1+M2+2*NN
      ENDIF
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE SGMATR(M1,M2,SYM,NSPAR,NELEM,ARRAY,IARRAY,MARRAY,
     *                  OARRAY,NARRAY,ISTOR)
      CHARACTER*(*) SYM
      INTEGER M1,M2,NSPAR,NELEM,OARRAY,NARRAY,MARRAY,ISTOR
      REAL     ARRAY(*)
      INTEGER IARRAY(*)
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     SYM ... Symmetry of the matrix.
C     NSPAR...Sparseness of the matrix. NSPAR.GE.0 is assumed on input.
C     OARRAY..Address of the first storage location in array ARRAY
C             used for the matrix.
C     NARRAY..Number of storage locations for the matrix.
C     IARRAY..Components of the given matrix.
C     ARRAY...Components of the given matrix.
C     ISTOR...Form of the given matrix.
C             ISTOR=0 ... form "as on a disk"
C             ISTOR=1 ... form "as in the memory"
C     For detailed description of the parameters describing the matrices
C     refer to file forms.htm.
C     Detailed description of storage of matrices in the memory
C     is given above.
C     MARRAY..Dimension of array ARRAY.
C
C Output:
C     NSPAR...Sparseness of the matrix. NSPAR=-1 on output.
C     NELEM...Number of elements of the matrix.
C     ISTOR...Form of the given matrix. ISTOR=0 on output.
C     OARRAY..New address of the first storage location in array ARRAY.
C     NARRAY..New number of storage locations for the matrix.
C     IARRAY..Components of the given matrix.
C     ARRAY...Components of the given matrix.
C
C Coded by Petr Bulant
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL ERROR,ISYM,NELMAT
      INTEGER ISYM,NELMAT
C
C     Local storage location:
      INTEGER ISY,I1,I2,I3,J1,J2,J3,I11,I12
      CHARACTER*72 TXTERR
C
C.......................................................................
C
      IF (NSPAR.LT.0) THEN
C       MAT-15
        CALL ERROR('MAT-15: Wrong invocation of SGMATR.')
      ENDIF
C
      ISY=ISYM(SYM)
C
      IF (ISTOR.EQ.0) THEN
C       Rewriting the matrix from the sparse form "as on the disk"
C       to the normal (not sparse) form:
        J2=NELMAT(M1,M2,ISY)
        J3=J2
        IF (OARRAY+NELEM+J2-1.GT.MARRAY) THEN
C         MAT-16
          WRITE(TXTERR,'(A,I9,A)')
     *    'MAT-16: Array RAM too small,',OARRAY+NELEM+J2-1-MARRAY,
     *    ' units missing.'
          CALL ERROR(TXTERR)
        ENDIF
        J1=NELEM
C       J1 is the index in sparse matrix, J2 in the non-sparse matrix
        DO 4, I2=M2,1,-1
C         I2 is the number of column
          I12=1
          IF     (ISY.EQ.1) THEN
            I11=I2
            I12=I2
          ELSEIF (ISY.EQ.2) THEN
            I11=I2
          ELSEIF (ISY.EQ.3) THEN
            I11=I2-1
          ELSE
            I11=M1
          ENDIF
          DO 2, I1=I11,I12,-1
C           I1 is the number of line
            I3=(I2-1)*M1+I1
C           I3 is the matrix index of the element [I1,I2]
            IF     (IARRAY(OARRAY-1+J1).LT.I3) THEN
              ARRAY(OARRAY-1+NELEM+J2)=0.
            ELSEIF (IARRAY(OARRAY-1+J1).EQ.I3) THEN
              ARRAY(OARRAY-1+NELEM+J2)=ARRAY(OARRAY-1+NELEM+J1)
              J1=J1-1
            ELSE
C             MAT-17
              CALL ERROR('MAT-17: Disorder in SGMATR')
C             This error should not appear.
            ENDIF
            J2=J2-1
   2      CONTINUE
   4    CONTINUE
        IF ((J1.NE.0).OR.(J2.NE.0)) THEN
C         MAT-18
          CALL ERROR('MAT-18: Disorder in SGMATR')
C         This error should not appear.
        ENDIF
        NSPAR=-1
        J2=NELEM
        NELEM=J3
        NARRAY=J3
        CALL MSHIFT(M1,M2,ISY,NSPAR,NELEM,IARRAY,ARRAY,MARRAY,
     *                  OARRAY+J2,NARRAY,ISTOR,OARRAY)
        RETURN
      ENDIF
C     Rewriting the matrix from the sparse form "as in the memory"
C     to the normal (not sparse) form:
      J2=MARRAY
      DO 20, I2=M2,1,-1
C       Initiating storage locations for column I2:
        IF (ISY.EQ.1) THEN
          J1=J2-1
        ELSEIF (ISY.EQ.2) THEN
          J1=J2-I2
        ELSEIF (ISY.EQ.3) THEN
          J1=J2-I2+1
        ELSE
          J1=J2-M1
        ENDIF
        IF (J1.LE.IARRAY(OARRAY+I2)) THEN
C         MAT-19
          WRITE(TXTERR,'(A,I9,A)')
     *    'MAT-19: Array RAM too small,',IARRAY(OARRAY+I2)-J1,
     *    ' units missing.'
          CALL ERROR(TXTERR)
        ENDIF
        DO 8, I3=J1+1,J2
          ARRAY(I3)=0.
   8    CONTINUE
C       Rewriting column I2:
        DO 10, I3=IARRAY(OARRAY+I2),IARRAY(OARRAY+I2-1)+2,-2
          I1=IARRAY(I3-1)
          IF ( (ISY.EQ.4).OR.
     *        ((ISY.EQ.1).AND.(I1.EQ.I2)).OR.
     *        ((ISY.EQ.2) .AND.(I1.LE.I2)).OR.
     *        ((ISY.EQ.3).AND.(I1.LT.I2))) THEN
            IF (ISY.EQ.1) THEN
              ARRAY(J1+1)=ARRAY(I3)
            ELSE
              ARRAY(J1+I1)=ARRAY(I3)
            ENDIF
          ENDIF
  10    CONTINUE
        J2=J1
  20  CONTINUE
      NARRAY=MARRAY-J2
C     Shifting the matrix to the beginning of ARRAY:
      J2=J2+1
      DO 30, I1=J2,MARRAY
        ARRAY(OARRAY+I1-J2)=ARRAY(I1)
  30  CONTINUE
      NELEM=NARRAY
      NSPAR=-1
      ISTOR=0
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE GSMAT(M1,M2,SYM,NSPAR,NELEM,ARRAY,IARRAY,MARRAY,
     *                  OARRAY,NARRAY,ISTOR,RSPAR)
      CHARACTER*(*) SYM
      INTEGER M1,M2,NSPAR,NELEM,OARRAY,NARRAY,MARRAY,ISTOR
      REAL     ARRAY(*),RSPAR
      INTEGER IARRAY(*)
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.
C If the matrix is changed to sparse, it is written to the beginning
C of arrays IARRAY and ARRAY.
C The subroutine uses whole arrays ARRAY and IARRAY.
C
C Input:
C     M1...   Number of rows of the matrix.
C     M2...   Number of columns of the matrix.
C     SYM ... Symmetry of the matrix.
C     NSPAR...Sparseness of the matrix.
C             NSPAR.LT.0 is assumed on input.
C     NELEM...Number of elements of the matrix.
C     OARRAY..Address of the first storage location in array ARRAY
C             used for the matrix.
C     NARRAY..Number of storage locations for the matrix.
C     ARRAY...Components of the given matrix.
C     For detailed description of the parameters describing the matrices
C     refer to file forms.htm.
C     Detailed description of storage of matrices in the memory
C     is given above.
C     MARRAY..Dimension of array ARRAY.
C     RSPAR...Minimum rate of sparseness to change the matrix
C             into the sparse matrix.
C
C Output:
C     NSPAR...Sparseness of the matrix.
C             For NSPAR.LT.0 on output, there are no changes on output.
C     NELEM...Number of elements of the matrix.
C     OARRAY..New address of the first storage location in array ARRAY.
C             For NSPAR.GE.0 on output OARRAY=1 on output.
C     NARRAY..New number of storage locations for the matrix.
C     IARRAY..Components of the given matrix.
C     ARRAY...Components of the given matrix.
C     ISTOR...Form of the matrix:
C             For NSPAR.GE.0 on output ISTOR=0 on output,
C             i.e. form "as on a disk".
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
      INTEGER NSPMAT
C
C     Local storage location:
      INTEGER I1,NN
      REAL RSP0
C
C.......................................................................
C
      IF ((NSPAR.GE.0).OR.(NELEM.NE.NARRAY)) THEN
C       MAT-20
        CALL ERROR('MAT-20: Wrong invocation of GSMAT.')
      ENDIF
C     Number of zero elements:
      NN=NSPMAT(NARRAY,ARRAY(OARRAY))
      RSP0=RSPAR
      RSPAR=FLOAT(NN)/FLOAT(NARRAY)
      IF (RSPAR.LT.RSP0) THEN
C       Matrix will stay non-sparse.
        RETURN
      ENDIF
C     Moving the matrix to the beginning of ARRAY:
      IF (OARRAY.NE.1) THEN
        DO 10, I1=1,NELEM
          ARRAY(I1)=ARRAY(OARRAY+I1-1)
  10    CONTINUE
        OARRAY=1
      ENDIF
C     Changing the matrix into the sparse matrix:
      ISTOR=0
      CALL GSMATR(M1,M2,SYM,NSPAR,NELEM,ARRAY,IARRAY,MARRAY,
     *                  OARRAY,NARRAY,ISTOR)
C
      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-21
        CALL ERROR('MAT-21: 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
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-22
        CALL ERROR('MAT-22: 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,IARRAY,ARRAY,MARRAY,
     *                  OARRAY,NARRAY,ISTOR,OARNEW)
      INTEGER M1,M2,ISYM,NSPAR,NELEM,OARRAY,NARRAY,MARRAY,ISTOR,OARNEW
      REAL     ARRAY(*)
      INTEGER IARRAY(*)
C
C Subroutine designed to shift the matrix in arrays IARRAY and ARRAY
C to the new position OARNEW.
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     NSPAR...Sparseness of the matrix.
C     NELEM...Number of elements of the matrix.
C     IARRAY..Components of the given matrix.
C     ARRAY...Components of the given matrix.
C     MARRAY..Dimension of arrays IARRAY and ARRAY.
C     OARRAY..Address of the first storage location in arrays ARRAY
C             and IARRAY used for the matrix.
C     NARRAY..Number of storage locations for the matrix.
C     ISTOR...Form used for the matrix if the matrix is sparse.
C             ISTOR=0 ... form "as on a disk"
C             ISTOR=1 ... form "as in the memory"
C     OARNEW..Address of the first storage location in arrays ARRAY
C             and IARRAY where the matrix is to be shifted.
C     For detailed description of the parameters describing the matrices
C     refer to file forms.htm.
C     Detailed description of storage of matrices in the memory
C     is given above.
C
C Output:
C     OARRAY..OARRAY is put to OARNEW.
C     IARRAY..Components of the shifted matrix.
C     ARRAY...Components of the shifted matrix.
C
C Coded by Petr Bulant
C
C-----------------------------------------------------------------------
C
C     Local storage location:
      INTEGER II,ISHIFT
C
C.......................................................................
C
      IF ((OARNEW.LT.1).OR.(OARNEW+NARRAY-1.GT.MARRAY).OR.
     *    (OARRAY.LT.1).OR.(OARRAY+NARRAY-1.GT.MARRAY)) THEN
C       MAT-23
        CALL ERROR('MAT-23: Wrong invocation of MSHIFT.')
C       Storage locations OARRAY to OARRAY+NARRAY and
C       OARNEW to OARNEW+NARRAY must fit into the memory.
C       This error should not appear.
      ENDIF
C
      IF (OARNEW.EQ.OARRAY) RETURN
C
      IF (NSPAR.LT.0) THEN
C       Non-sparse matrix:
        IF (OARNEW.LT.OARRAY) THEN
          DO 10, II=0,NARRAY-1
            ARRAY(OARNEW+II)=ARRAY(OARRAY+II)
  10      CONTINUE
        ELSE
          DO 11, II=NARRAY-1,0,-1
            ARRAY(OARNEW+II)=ARRAY(OARRAY+II)
  11      CONTINUE
        ENDIF
      ELSE
C       Sparse martix:
        IF (ISTOR.EQ.0) THEN
C         Form "as on disk":
          IF (OARNEW.LT.OARRAY) THEN
            DO 20, II=0,NELEM-1
              IARRAY(OARNEW+II)=IARRAY(OARRAY+II)
  20        CONTINUE
            DO 21, II=0,NELEM-1
              ARRAY (OARNEW+NELEM+II)=ARRAY(OARRAY+NELEM+II)
  21        CONTINUE
          ELSE
            DO 22, II=NELEM-1,0,-1
              ARRAY (OARNEW+NELEM+II)=ARRAY(OARRAY+NELEM+II)
  22        CONTINUE
            DO 23, II=NELEM-1,0,-1
              IARRAY(OARNEW+II)=IARRAY(OARRAY+II)
  23        CONTINUE
          ENDIF
        ELSE
C         Form "as in the memory":
          ISHIFT=OARNEW-OARRAY
          IF (OARNEW.LT.OARRAY) THEN
            DO 30, II=0,M2
              IARRAY(OARNEW+II)=IARRAY(OARRAY+II)+ISHIFT
  30        CONTINUE
            DO 31, II=0,NELEM-1
              IARRAY(OARNEW+M2+1+II*2)  =IARRAY(OARRAY+M2+1+II*2)
               ARRAY(OARNEW+M2+1+II*2+1)= ARRAY(OARRAY+M2+1+II*2+1)
  31        CONTINUE
          ELSE
            DO 32, II=NELEM-1,0,-1
               ARRAY(OARNEW+M2+1+II*2+1)= ARRAY(OARRAY+M2+1+II*2+1)
              IARRAY(OARNEW+M2+1+II*2)  =IARRAY(OARRAY+M2+1+II*2)
  32        CONTINUE
            DO 33, II=M2,0,-1
              IARRAY(OARNEW+II)=IARRAY(OARRAY+II)+ISHIFT
  33        CONTINUE
          ENDIF
        ENDIF
      ENDIF
C
      OARRAY=OARNEW
      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 fo 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 I2,NELEM
C
C     FORMM ..Form of the files with matrices.
C
C.......................................................................
C
      CHARACTER*80 FILED
      CHARACTER*4 SYMM
      INTEGER IOUT(1)
C
      IF (FILE.EQ.' ') THEN
C       MAT-24
        CALL ERROR('MAT-24: 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,-1,SYMM,FORMM)
      CALL WMATD(LU,FILED,-1,FORMM,NELEM,IOUT,OUT)
      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 fo using RMAT.
C Subroutine designed to read a matrix from the file.
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             out 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
      CHARACTER*80 FILED
      CHARACTER*4 SYMM
      INTEGER IOUT(1),M1R,M2R,NSPAR,NELEM
C
      IF (FILE.EQ.' ') THEN
C       MAT-25
        CALL ERROR('MAT-25: 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,NSPAR,SYMM,FORMM,NELEM)
      IF ((M1.NE.M1R).OR.(NSPAR.GE.0).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-26
        CALL ERROR('MAT-26: 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
      CALL RMATD(LU,FILED,NSPAR,FORMM,NELEM,IOUT,ARRAY)
      RETURN
      END
C
C=======================================================================
C