C
C Program MATLIN to compute matrix M3 as a linear combination
C M3=COEF1*M1+COEF2*M2 of matrices M1 and M2. The matrices M1 and M2
C must have the same number of rows and the same number of columns,
C and they must be of the same symmetry, resulting matrix M3 will have
C the same symmetry.
C
C Version: 6.10
C Date: 2007, June 8
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             No default, MATIN1 must be specified and cannot be blank.
C     MATIN2='string' .. Name of the header file of the input matrix M2.
C             Default: MATIN2=' ' ... no input matrix M2, matrix M1 will
C             be written on output. This option may be used to multiply
C             matrix M1 by a constant, or to change its form.
C     MATOUT='string' . Name of the header file of the output matrix M3.
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 M3:
C     MATFORM='string' ... Form of the output file with the matrix.
C             Possible values of MATFORM are:
C               MATFORM='formatted'   ... formatted file
C               MATFORM='unformatted' ... unformatted file
C             Default: MATFORM='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               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 Coefficients to multiply the input matrices:
C     COEF1=real ... coefficient to multiply matrix M1.
C             Default: COEF1=1.
C     COEF2=real ... Same as COEF1, but for matrix M2.
C             Default: COEF2=1.
C Optional parameter specifying the form of the output formatted
C matrix data files:
C     NUMLIN=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 ERROR,RSEP1,RSEP3I,RSEP3T,LOWER,RMATH,RMATD,WMATH,WMATD,
     *GSMAT,GSMATR,SGMATR,ISYM,NELMAT,SWITCH
      INTEGER ISYM,NELMAT
C     ERROR ... File error.for.
C     RSEP1,RSEP3I,RSEP3T ...
C             File sep.for.
C     RMATH,RMATD,WMATH,WMATD  ...
C     GSMAT,GSMATR,SGMATR,ISYM,NELMAT ...
C             File mat.for.
C     SWITCH ... This file.
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C
      CHARACTER*80 FILSEP,FILE1,FILE2,FILE3,FILED1,FILED2,FILED3
      CHARACTER*80 SYM1,SYM2,SYM3,SYMC,FORM1,FORM2,FORM3
      INTEGER M1,M2,M1B,M2B,NEL1,NEL2,NELC,OA,NA,
     *OB,NB,OC,NC,NSPAR1,NSPAR2,ISPAR3,NSPARC,ISYM1,ISYM2,ISYMC,
     *ISTORA,ISTORC,LU1,I1,I2,I3,I4,IR,IC
      REAL RSPAR3,COEF1,COEF2
      CHARACTER*72 TXTERR
      PARAMETER (LU1=1)
C
C-----------------------------------------------------------------------
C
C     Reading a name of the file with the input data:
      WRITE(*,'(A)') '+MATLIN: 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       MATLIN-01
        CALL ERROR('MATLIN-01: SEP file not given.')
      ENDIF
C
C     Reading the names of matrices header files:
      CALL RSEP3T('MATIN1',FILE1,' ')
      IF (FILE1.EQ.' ') THEN
C       MATLIN-02
        CALL ERROR('MATLIN-02: File MATIN1 not given.')
      ENDIF
      CALL RSEP3T('MATIN2',FILE2,' ')
      CALL RSEP3T('MATOUT',FILE3,' ')
      FILED3=' '
      IF (FILE3.EQ.' ') THEN
C       MATLIN-03
        CALL ERROR('MATLIN-03: File MATOUT not given.')
      ENDIF
C     Reading the multiplication coefficients:
      CALL RSEP3R('COEF1',COEF1,1.)
      CALL RSEP3R('COEF2',COEF2,1.)
ccc   IF ((COEF1.EQ.0.).OR.(COEF2.EQ.0.)) THEN
Ccc     MATLIN-00
ccc     CALL ERROR('MATLIN-00: Zero coefficient COEF1 or COEF2.')
ccc   ENDIF
C     Reading the header file of the input matrix M1:
      CALL RMATH(LU1,FILE1,FILED1,M1,M2,NSPAR1,SYM1,FORM1,NEL1)
      ISYM1=ISYM(SYM1)
C     Reading the properties of the output file:
      SYM3=SYM1
      CALL RSEP3T('MATFORM',FORM3,'formatted')
      CALL LOWER(FORM3)
      CALL RSEP3I('SPARSE',ISPAR3,0)
C
      IF ((FILE2.EQ.' ').OR.(COEF2.EQ.0.)) THEN
C       Matrix M2 is not given, matrix M1 is to be multiplied by COEF1
C       and written as output matrix M3:
        OA=1
C       Reading input matrix A:
        IF (NSPAR1.GE.0) THEN
          NA=2*NEL1
        ELSE
          NA=NEL1
        ENDIF
        IF (NA.GT.MRAM) THEN
C         MATLIN-04
          WRITE(TXTERR,'(A,I9,A)')
     *    'MATLIN-04: Array RAM too small,',NA-MRAM,' units missing.'
          CALL ERROR(TXTERR)
        ENDIF
        CALL RMATD(LU1,FILED1,NSPAR1,FORM1,NEL1,IRAM(OA),RAM(OA))
        ISTORA=0
        OC=OA
        NC=NA
        NSPARC=NSPAR1
        NELC=NEL1
        SYMC=SYM1
        ISYMC=ISYM1
        ISTORC=ISTORA
        IF (COEF1.NE.1.) THEN
          IF (NSPARC.GE.0) THEN
            DO 2, I1=OC+NELC,OC+NELC+NELC-1
              RAM(I1)=RAM(I1)*COEF1
   2        CONTINUE
          ELSE
            DO 4, I1=OC,OC+NELC-1
              RAM(I1)=RAM(I1)*COEF1
   4        CONTINUE
          ENDIF
        ENDIF
        GOTO 100
      ENDIF
C
C     Matrix M2 is given:
      CALL RMATH(LU1,FILE2,FILED2,M1B,M2B,NSPAR2,SYM2,FORM2,NEL2)
      ISYM2=ISYM(SYM2)
      IF ((M1.NE.M1B).OR.(M2.NE.M2B).OR.(ISYM1.NE.ISYM2)) THEN
C       MATLIN-05
        CALL ERROR('MATLIN-05: Wrong input matrices.')
C       The matrices M1 and M2 must have the same number of rows
C       and the same number of columns, and they must be of the
C       same symmetry.
      ENDIF
      IF ((NSPAR1.GE.0).AND.(NSPAR2.LT.0)) THEN
C       If first matrix is sparse and the second not, switching them:
        CALL SWITCH(FILE1,FILED1,NSPAR1,SYM1,ISYM1,FORM1,NEL1,COEF1,
     *              FILE2,FILED2,NSPAR2,SYM2,ISYM2,FORM2,NEL2,COEF2)
      ENDIF
C     Reading input matrix A:
      IF (NSPAR1.GE.0) THEN
        NA=2*NEL1
        OA=1+2*(NEL1+NEL2)
      ELSE
        NA=NEL1
        OA=1
      ENDIF
      IF (OA+NA-1.GT.MRAM) THEN
C       MATLIN-06
        WRITE(TXTERR,'(A,I9,A)')
     *  'MATLIN-06: Array RAM too small,',OA+NA-1-MRAM,' units missing.'
        CALL ERROR(TXTERR)
      ENDIF
      CALL RMATD(LU1,FILED1,NSPAR1,FORM1,NEL1,IRAM(OA),RAM(OA))
      ISTORA=0
C     Reading input matrix B:
      IF (NSPAR2.GE.0) THEN
        NB=2*NEL2
      ELSE
        NB=NEL2
      ENDIF
      OB=OA+NA
      IF (OB+NB-1.GT.MRAM) THEN
C       MATLIN-07
        WRITE(TXTERR,'(A,I9,A)')
     *  'MATLIN-07: Array RAM too small,',OB+NB-1-MRAM,' units missing.'
        CALL ERROR(TXTERR)
      ENDIF
      CALL RMATD(LU1,FILED2,NSPAR2,FORM2,NEL2,IRAM(OB),RAM(OB))
C
C     Summation:
      WRITE(*,'(A)') '+MATLIN: Calculating the linear combination ...  '
C
      IF (NSPAR1.GE.0) THEN
C       Both matrices are sparse:
        I1=1
        I2=1
        I3=0
C       Loop for matrix indices:
  20    CONTINUE
          IF     (IRAM(OA-1+I1).LT.IRAM(OB-1+I2)) THEN
            I3=I3+1
            IRAM(I3)=IRAM(OA-1+I1)
            I1=I1+1
            IF (I1.LE.NEL1) GOTO 20
          ELSEIF (IRAM(OA-1+I1).EQ.IRAM(OB-1+I2)) THEN
            IF ((RAM(OA-1+NEL1+I1)*COEF1+RAM(OB-1+NEL2+I2)*COEF2).NE.0.)
     *        THEN
              I3=I3+1
              IRAM(I3)=IRAM(OA-1+I1)
            ENDIF
            I1=I1+1
            I2=I2+1
            IF (I1.LE.NEL1) GOTO 20
            IF (I2.LE.NEL2) GOTO 20
          ELSEIF (IRAM(OA-1+I1).GT.IRAM(OB-1+I2)) THEN
            I3=I3+1
            IRAM(I3)=IRAM(OB-1+I2)
            I2=I2+1
            IF (I2.LE.NEL2) GOTO 20
          ENDIF
        DO 22, I4=I1,NEL1
          I3=I3+1
          IRAM(I3)=IRAM(OA-1+I4)
  22    CONTINUE
        DO 24, I4=I2,NEL2
          I3=I3+1
          IRAM(I3)=IRAM(OB-1+I4)
  24    CONTINUE
C       Matrix indices of matrix C are recorded.
        NELC=I3
        I1=1
        I2=1
        I3=0
C       Loop for matrix values:
  30    CONTINUE
          IF     (IRAM(OA-1+I1).LT.IRAM(OB-1+I2)) THEN
            I3=I3+1
            RAM(NELC+I3)=RAM(OA-1+NEL1+I1)*COEF1
            I1=I1+1
            IF (I1.LE.NEL1) GOTO 30
          ELSEIF (IRAM(OA-1+I1).EQ.IRAM(OB-1+I2)) THEN
            IF ((RAM(OA-1+NEL1+I1)*COEF1+RAM(OB-1+NEL2+I2)*COEF2).NE.0.)
     *        THEN
              I3=I3+1
              RAM(NELC+I3)=
     *          RAM(OA-1+NEL1+I1)*COEF1+RAM(OB-1+NEL2+I2)*COEF2
            ENDIF
            I1=I1+1
            I2=I2+1
            IF (I1.LE.NEL1) GOTO 30
            IF (I2.LE.NEL2) GOTO 30
          ELSEIF (IRAM(OA-1+I1).GT.IRAM(OB-1+I2)) THEN
            I3=I3+1
            RAM(NELC+I3)=RAM(OB-1+NEL2+I2)*COEF2
            I2=I2+1
            IF (I2.LE.NEL2) GOTO 30
          ENDIF
        DO 32, I4=I1,NEL1
          I3=I3+1
          RAM(NELC+I3)=RAM(OA-1+NEL1+I4)*COEF1
  32    CONTINUE
        DO 34, I4=I2,NEL2
          I3=I3+1
          RAM(NELC+I3)=RAM(OB-1+NEL2+I4)*COEF2
  34    CONTINUE
C
        OC=1
        NC=2*I3
        SYMC=SYM1
        ISYMC=ISYM1
        ISTORC=0
        NSPARC=NELMAT(M1,M2,ISYMC)-NELC
      ELSE
C       First matrix is non-sparse:
        IF (COEF1.NE.1.) THEN
          DO 40, I1=OA,OA+NEL1-1
            RAM(I1)=RAM(I1)*COEF1
   40     CONTINUE
        ENDIF
        IF (NSPAR2.GE.0) THEN
          DO 42, I1=1,NEL2
            IR=MOD((IRAM(OB-1+I1)-1),M1)+1
            IC=(IRAM(OB-1+I1)-1)/M1+1
            IF     (ISYM2.EQ.1) THEN
              I2=IR
            ELSEIF (ISYM2.EQ.2) THEN
              I2=(IC-1)*IC/2+IR
            ELSEIF (ISYM2.EQ.3) THEN
              I2=(IC-1)*(IC-2)/2+IR
            ELSE
              I2=(IC-1)*M1+IR
            ENDIF
            RAM(I2)=RAM(I2)+RAM(OB-1+NEL2+I1)*COEF2
   42     CONTINUE
        ELSE
          DO 50, I1=OA,OA+NEL1-1
            RAM(I1)=RAM(I1)+RAM(OB-1+I1)*COEF2
   50     CONTINUE
        ENDIF
        OC=1
        NELC=NEL1
        NC=NA
        SYMC=SYM1
        ISYMC=ISYM1
        ISTORC=0
        NSPARC=NSPAR1
      ENDIF
C
 100  CONTINUE
C
      IF ((ISPAR3.EQ.1).AND.(NSPARC.LT.0)) THEN
C       Changing non-sparse matrix to sparse matrix:
        CALL GSMATR(M1,M2,SYMC,NSPARC,NELC,RAM,IRAM,MRAM,OC,NC,ISTORC)
      ELSEIF ((ISPAR3.EQ.-1).AND.(NSPARC.GE.0)) THEN
C       Changing sparse matrix to non-sparse matrix:
        CALL SGMATR(M1,M2,SYMC,NSPARC,NELC,RAM,IRAM,MRAM,OC,NC,ISTORC)
      ELSEIF (ISPAR3.EQ.0) THEN
C       Automatic selection of the sparseness:
        IF (NSPARC.LT.0) THEN
C         Matrix is non-sparse, it will be changed to sparse if its
C         sparseness is 0.5 or more:
          RSPAR3=0.5
          CALL GSMAT(M1,M2,SYMC,NSPARC,NELC,RAM,IRAM,MRAM,
     *                  OC,NC,ISTORC,RSPAR3)
        ELSE
C         Matrix is sparse.
          RSPAR3=FLOAT(NSPARC)/FLOAT(NSPARC+NELC)
          IF (RSPAR3.LT.0.5) THEN
C           Matrix is sparse, sparseness is less than 0.5, thus changing
C           the sparse matrix to non-sparse matrix:
            CALL SGMATR(M1,M2,SYMC,NSPARC,NELC,RAM,IRAM,MRAM,
     *                  OC,NC,ISTORC)
          ENDIF
        ENDIF
      ENDIF
C
C     Writing output matrix C:
      CALL WMATH(LU1,FILE3,FILED3,M1,M2,NSPARC,SYM3,FORM3)
      CALL WMATD(LU1,FILED3,NSPARC,FORM3,NELC,IRAM(OC),RAM(OC))
C
      WRITE(*,'(A)') '+MATLIN: Done.                 '
C
      STOP
      END
C
C=======================================================================
C
      SUBROUTINE SWITCH(FILE1,FILED1,NSPAR1,SYM1,ISYM1,FORM1,NEL1,COEF1,
     *                  FILE2,FILED2,NSPAR2,SYM2,ISYM2,FORM2,NEL2,COEF2)
      CHARACTER*(*) FILE1,FILED1,FILE2,FILED2,SYM1,SYM2,FORM1,FORM2
      CHARACTER*80 FILE3,FILED3,SYM3,FORM3
      INTEGER NSPAR1,ISYM1,NEL1,NSPAR2,ISYM2,NEL2,NSPAR3,ISYM3,NEL3
      REAL COEF1,COEF2,COEF3
C-----------------------------------------------------------------------
       FILE3 =  FILE1
      FILED3 = FILED1
      NSPAR3 = NSPAR1
        SYM3 =   SYM1
       ISYM3 =  ISYM1
       FORM3 =  FORM1
        NEL3 =   NEL1
       COEF3 =  COEF1
       FILE1 =  FILE2
      FILED1 = FILED2
      NSPAR1 = NSPAR2
        SYM1 =   SYM2
       ISYM1 =  ISYM2
       FORM1 =  FORM2
        NEL1 =   NEL2
       COEF1 =  COEF2
       FILE2 =  FILE3
      FILED2 = FILED3
      NSPAR2 = NSPAR3
        SYM2 =   SYM3
       ISYM2 =  ISYM3
       FORM2 =  FORM3
        NEL2 =   NEL3
       COEF2 =  COEF3
      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 'indexi.for'
C     indexi.for
C
C=======================================================================
C