C
C Program MATINV to compute symmetric matrix M2, which is inverse to
C input symmetric matrix M1.
C
C Version: 6.20
C Date: 2008, June 12
C
C Coded by Petr Bulant
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C
C.......................................................................
C                                                    
C Description of data files:
C
C Input data read from the standard input device (*):
C     The data are read by the list directed input (free format) and
C     consist of a single string 'SEP':
C     'SEP'...String in apostrophes containing the name of the input
C             SEP parameter or history file with the input data.
C     No default, 'SEP' must be specified and cannot be blank.
C
C                                                     
C Input data file 'SEP':
C     File 'SEP' has the form of the SEP
C     parameter file.  The parameters, which do not differ from their
C     defaults, need not be specified in file 'SEP'.
C Filenames of the header files of the matrices:
C     MATIN1='string' .. Name of the header file of the input matrix M1.
C             Matrix M1 must be symmetric.
C             No default, MATIN1 must be specified and cannot be blank.
C     MATOUT='string' . Name of the header file of the output matrix M2.
C             No default, MATOUT must be specified and cannot be blank.
C     For general description of the files with matrices refer to file
C     forms.htm.
C Form of the output file with the matrix M2:
C     FORMM='string' ... Form of the output file with the matrix.
C             Possible values of FORMM are:
C               FORMM='FORMATTED'   ... formatted file
C               FORMM='UNFORMATTED' ... unformatted file
C             Default: FORMM='FORMATTED'
C     SPARSE=integer ... Identifies whether the matrix should be sparse.
C             Possible values of SPARSE are:
C               SPARSE=1  ... sparse matrix
C               SPARSE=0  ... automatic selection of the sparseness:
C                 matrix with half or more zero elements will be sparse
C                 in the index format, otherwise non-sparse matrix
C               SPARSE=-1 ... normal (not sparse) matrix
C             Default: SPARSE=0 (automatic selection of the sparseness)
C     For detailed description of the forms of the files with matrices
C     refer to file forms.htm.
C Optional parameter specifying the form of the output formatted
C matrix data files:
C     NUMLINM=positive integer ... Number of reals or integers in one
C             line of the output file.  See the description in file
C             mat.for.
C
C For detailed description of storage of matrices in the memory
C refer to file mat.for.
C
C-----------------------------------------------------------------------
C Subroutines and external functions required:
      EXTERNAL SMBLO,IND,ERROR,RSEP1,RSEP3I,RSEP3T,UPPER,
     * RMATH,RMATD,WMATH,WMATD,GSMAT,GSMATR,SGMATR,NELMAT,ISYM,SINV
      INTEGER IND,NELMAT,ISYM
C     SMBLO,IND ... This file.
C     ERROR ... File error.for.
C     RSEP1,RSEP3I,RSEP3T ...
C     File sep.for.
C     UPPER ... File
C     length.for.
C     RMATH,RMATD,WMATH,WMATD,GSMAT,GSMATR,SGMATR,NELMAT,ISYM ...
C     File mat.for.
C     SINV ... File sinv.for.
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C
      CHARACTER*80 FILSEP,FILE1,FILE3,FILED1,FILED3,SYM1,FORM1,
     *             FORM3,SPARS1
      INTEGER M1,M2,NEL1,IA,NA,NSPAR1,ISPAR3,LU1,I1,I2,I3,I4,
     * N,NN,IANEW,NB,IBMI,IBMA,IER,ISYM1
      REAL RSPAR1,EPS
      CHARACTER*72 TXTERR
      PARAMETER (LU1=1)
C
C-----------------------------------------------------------------------
C
      EPS=0.000001
C
C     Reading a name of the file with the input data:
      WRITE(*,'(A)') '+MATINV: Enter input filename: '
      FILSEP=' '
      READ(*,*) FILSEP
C
C     Reading all the data from the SEP file into the memory:
      IF (FILSEP.NE.' ') THEN
        CALL RSEP1(LU1,FILSEP)
      ELSE
C       MATINV-01
        CALL ERROR('MATINV-01: SEP file not given.')
      ENDIF
C
C     Reading the names of matrices header files:
      CALL RSEP3T('MATIN1',FILE1,' ')
      IF (FILE1.EQ.' ') THEN
C       MATINV-02
        CALL ERROR('MATINV-02: File MATIN1 not given.')
      ENDIF
      CALL RSEP3T('MATOUT',FILE3,' ')
      FILED3=' '
      IF (FILE3.EQ.' ') THEN
C       MATINV-03
        CALL ERROR('MATINV-03: File MATOUT not given.')
      ENDIF
C     Reading the header file of the input matrix:
      CALL RMATH(LU1,FILE1,FILED1,M1,M2,SPARS1,NEL1,SYM1,FORM1)
      ISYM1=ISYM(SYM1)
      IF (SPARS1.EQ.' ') THEN
        NSPAR1=-1
      ELSEIF (SPARS1.EQ.'CSC') THEN
        NSPAR1=NELMAT(M1,M2,ISYM1)-NEL1
      ELSE
C       MATINV-06
        CALL ERROR('MATINV-06: Wrong sparseness of the input matrix.')
C       In this version only dense or sparse CSC matrix is allowed.
      ENDIF
      IF (SYM1.NE.'sym') THEN
C       MATINV-07
        CALL ERROR('MATINV-07: Input matrix is not symmetric.')
      ENDIF
C     Reading the properties of the output file:
      CALL RSEP3T('FORMM',FORM3,'FORMATTED')
      CALL UPPER(FORM3)
      CALL RSEP3I('SPARSE',ISPAR3,0)
C
      N=M1
      NN=M1*(M1+1)/2
      IF (3*NN+N+1.GT.MRAM) THEN
C       MATINV-04
        WRITE(TXTERR,'(A,I9,A)')
     *  'MATINV-04: Array RAM too small,',3*NN+N+1-MRAM,' units missing'
        CALL ERROR(TXTERR)
      END IF
C
C     Reading input matrix:
      IF (NSPAR1.GE.0) THEN
        IA=1
        NA=2*NEL1
      ELSE
        IA=MRAM-NN+1
        NA=NEL1
      ENDIF
      CALL RMATD(LU1,FILED1,M2,SPARS1,NEL1,FORM1,IA)
      IF (NSPAR1.GE.0) THEN
C       Changing input matrix to non-sparse:
        CALL SGMATR(M1,M2,ISYM1,NSPAR1,NEL1,1,MRAM,IA,NA)
        IANEW=MRAM-NN+1
        CALL MSHIFT(M1,M2,ISYM1,NSPAR1,NEL1,1,MRAM,IA,NA,IANEW)
      ENDIF
C
      WRITE(*,'(A)') '+MATINV: Working...            '
C
C     Identification of blocks:
      CALL SMBLO(RAM(MRAM-NN+1),N,NN,NB,IRAM(MRAM-NN-N))
C
C     Preparing array for results:
      DO 5, I1=MRAM-2*NN-N,MRAM-NN-N-1
        RAM(I1)=0.
   5  CONTINUE
C
C     Computing the inverse matrix:
      DO 20, I1=1,NB
C       Current block:
        IBMI=IRAM(MRAM-NN-N+I1-1)+1
        IBMA=IRAM(MRAM-NN-N+I1)
C       Moving the block to RAM(1):
        I4=0
        DO 10, I2=IBMI,IBMA
          DO 8, I3=IBMI,I2
            I4=I4+1
            RAM(I4)=RAM(MRAM-NN+IND(I3,I2))
   8      CONTINUE
  10    CONTINUE
        I2=IBMA-IBMI+1
        I3=I2*(I2+1)/2
C       Inverting matrix:
        CALL SINV(RAM,I2,EPS,IER)
        IF(IER.NE.0) THEN
C         MATINV-05
          WRITE(TXTERR,'(A,I5,A)')
     *      'MATINV-05: Error',IER,' in subroutine SINV'
          CALL ERROR(TXTERR)
        ENDIF
C       Moving the block to RAM(MRAM-NN-N-1-NN+1):
        I4=0
        DO 14, I2=IBMI,IBMA
          DO 12, I3=IBMI,I2
            I4=I4+1
            RAM(MRAM-2*NN-N-1+IND(I3,I2))=RAM(I4)
  12      CONTINUE
  14    CONTINUE
  20  CONTINUE
C
C
      IA=MRAM-2*NN-N
      IF (ISPAR3.EQ.1) THEN
C       Changing non-sparse matrix to sparse matrix:
        CALL GSMATR(M1,M2,ISYM1,NSPAR1,NEL1,1,MRAM,IA,NA)
      ELSEIF (ISPAR3.EQ.0) THEN
C       Automatic selection of the sparseness:
C       Matrix is non-sparse, it will be changed to sparse if its
C       sparseness is 0.5 or more:
        RSPAR1=0.5
        CALL GSMAT(M1,M2,ISYM1,NSPAR1,NEL1,1,MRAM,IA,NA,RSPAR1)
      ENDIF
C
C     Writing output matrix:
      SPARS1=' '
      IF (NSPAR1.GE.0) SPARS1='CSC'
      CALL WMATH(LU1,FILE3,FILED3,M1,M2,SPARS1,NEL1,SYM1,FORM3)
      CALL WMATD(LU1,FILED3,M1,M2,SPARS1,NEL1,FORM3,IA)
C
      WRITE(*,'(A)') '+MATINV: Done.                 '
C
      STOP
      END
C
C=======================================================================
C
      SUBROUTINE SMBLO(SM,N,NN,NB,IB)
C
C=======================================================================
C     Subroutine to find blocks of a symmetric matrix
      INTEGER N,NN,NB,IB(0:N)
      REAL SM(NN)
C     Input: SM ... symmetric matrix
C            N  ... number of rows (and columns) of the matrix
C            NN ... N*(N+1)/2
C     Output: NB ... number of blocks of the matrix
C             IB ... array with numbers of lines (and rows), at which
C                    the blocks finish
      INTEGER I1,I2
      EXTERNAL IND
      INTEGER IND
C-----------------------------------------------------------------------
C     Initiating first block:
      NB=1
      IB(0)=0
      IB(1)=1
      I1=0
  10  CONTINUE
C       Loop for lines of the matrix:
        I1=I1+1
        DO 20, I2=N,IB(NB)+1,-1
C         Loop for rows of the line I1:
          IF (SM(IND(I1,I2)).NE.0.) THEN
C           The block is larger:
            IB(NB)=I2
            GOTO 21
          ENDIF
  20    CONTINUE
  21    CONTINUE
        IF (IB(NB).EQ.N) THEN
C         End of the search for blocks.
          RETURN
        ENDIF
        IF (IB(NB).EQ.I1) THEN
C         Block finished, continuing with the next block.
          NB=NB+1
          IB(NB)=I1+1
        ENDIF
      GOTO 10
      END
C=======================================================================
      INTEGER FUNCTION IND(I,J)
      INTEGER I,J
C     I ... Index of a line.
C     J ... Index of a row.
      IND=J*(J-1)/2+I
      RETURN
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'mat.for'
C     mat.for
      INCLUDE 'sinv.for'
C     sinv.for
      INCLUDE 'mfsd.for'
C     mfsd.for
C
C=======================================================================
C