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.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             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     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 Coefficients to multiply the input matrices:
C     COEF1=real ... coefficient to multiply matrix M1. Must not be 0.
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     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 ERROR,RSEP1,RSEP3I,RSEP3T,UPPER,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     UPPER ...
C     File length.for.
C     RMATH,RMATD,WMATH,WMATD,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,FORM1,FORM2,FORM3
      CHARACTER*80 SPARS1,SPARS2,SPARSC
      INTEGER M1,M2,M1B,M2B,NEL1,NEL2,NELC,IA,NA,
     *IB,NB,IC,NC,NSPAR1,NSPAR2,ISPAR3,NSPARC,ISYM1,ISYM2,ISYMC,
     *LU1,I1,I2,I3,I4,IROW,ICOL,KK
      REAL RSPAR3,COEF1,COEF2,AAA
      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.)
      IF (COEF1.EQ.0.) THEN
C       MATLIN-10
        CALL ERROR('MATLIN-10: COEF1 equal zero.')
      ENDIF
C     Reading the header file of the input matrix M1:
      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       MATLIN-08
        CALL ERROR('MATLIN-08: Wrong format of input matrix.')
C       In this version only dense or sparse CSC matrix is allowed.
      ENDIF
C     Reading the properties of the output file:
      SYM3=SYM1
      CALL RSEP3T('FORMM',FORM3,'FORMATTED')
      CALL UPPER(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:
        IA=1
C       Reading input matrix A:
        IF (NSPAR1.GE.0) THEN
          NA=M2+1+2*NEL1
        ELSE
          NA=NEL1
        ENDIF
        IF (IA-1+NA.GT.MRAM) THEN
C         MATLIN-04
          WRITE(TXTERR,'(A,I9,A)') 'MATLIN-04: Array RAM too small,',
     *    IA-1+NA-MRAM,' units missing.'
          CALL ERROR(TXTERR)
        ENDIF
        CALL RMATD(LU1,FILED1,M2,SPARS1,NEL1,FORM1,IA)
        IC=IA
        NC=NA
        NSPARC=NSPAR1
        NELC=NEL1
        ISYMC=ISYM1
        IF (COEF1.NE.1.) THEN
          IF (NSPARC.GE.0) THEN
            DO 2, I1=IRAM(IC)+1,IRAM(IC+M2)-1,2
              RAM(I1)=RAM(I1)*COEF1
   2        CONTINUE
          ELSE
            DO 4, I1=IC,IC+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,SPARS2,NEL2,SYM2,FORM2)
      ISYM2=ISYM(SYM2)
      IF (SPARS2.EQ.' ') THEN
        NSPAR2=-1
      ELSEIF (SPARS2.EQ.'CSC') THEN
        NSPAR2=NELMAT(M1B,M2B,ISYM2)-NEL2
      ELSE
C       MATLIN-09
        CALL ERROR('MATLIN-09: Wrong format of input sparse matrix.')
C       In this version only dense or sparse CSC matrix is allowed.
      ENDIF
      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,SPARS1,NSPAR1,SYM1,ISYM1,FORM1,NEL1,COEF1,
     *       FILE2,FILED2,SPARS2,NSPAR2,SYM2,ISYM2,FORM2,NEL2,COEF2)
      ENDIF
C     Reading input matrix A:
      IF (NSPAR1.GE.0) THEN
        NA=M2+1+2*NEL1
        IA=M2+1+2*(NEL1+NEL2)
      ELSE
        NA=NEL1
        IA=1
      ENDIF
      IF (IA+NA-1.GT.MRAM) THEN
C       MATLIN-06
        WRITE(TXTERR,'(A,I9,A)')
     *  'MATLIN-06: Array RAM too small,',IA+NA-1-MRAM,' units missing.'
        CALL ERROR(TXTERR)
      ENDIF
      CALL RMATD(LU1,FILED1,M2,SPARS1,NEL1,FORM1,IA)
C     Reading input matrix B:
      IF (NSPAR2.GE.0) THEN
        NB=M2+1+2*NEL2
      ELSE
        NB=NEL2
      ENDIF
      IB=IA+NA
      IF (IB+NB-1.GT.MRAM) THEN
C       MATLIN-07
        WRITE(TXTERR,'(A,I9,A)')
     *  'MATLIN-07: Array RAM too small,',IB+NB-1-MRAM,' units missing.'
        CALL ERROR(TXTERR)
      ENDIF
      CALL RMATD(LU1,FILED2,M2,SPARS2,NEL2,FORM2,IB)
C
C     Summation:
      WRITE(*,'(A)') '+MATLIN: Calculating the linear combination ...  '
C
      IF (NSPAR1.GE.0) THEN
C       Both matrices are sparse:
        IC=1
        I3=IC+M2+1-2
        IRAM(IC)=I3+2
C       Loop over columns:
        DO 28, KK=1,M2
          I1=IRAM(IA-1+KK)
          I2=IRAM(IB-1+KK)
          IF ((I1.LT.IRAM(IA+KK)).AND.(I2.LT.IRAM(IB+KK))) THEN
  20        CONTINUE
              IF     (IRAM(I1).LT.IRAM(I2)) THEN
                I3=I3+2
                IRAM(I3)=IRAM(I1)
                RAM(I3+1)=RAM(I1+1)*COEF1
                I1=I1+2
                IF (I1.LT.IRAM(IA+KK)) GOTO 20
              ELSEIF (IRAM(I1).EQ.IRAM(I2)) THEN
                AAA=RAM(I1+1)*COEF1+RAM(I2+1)*COEF2
                IF (AAA.NE.0.)
     *            THEN
                  I3=I3+2
                  IRAM(I3)=IRAM(I1)
                  RAM(I3+1)=AAA
                ENDIF
                I1=I1+2
                I2=I2+2
                IF (I1.LT.IRAM(IA+KK)) GOTO 20
                IF (I2.LT.IRAM(IB+KK)) GOTO 20
              ELSEIF (IRAM(I1).GT.IRAM(I2)) THEN
                I3=I3+2
                IRAM(I3)=IRAM(I2)
                RAM(I3+1)=RAM(I2+1)*COEF2
                I2=I2+2
                IF (I2.LT.IRAM(IB+KK)) GOTO 20
              ENDIF
          ENDIF
          DO 22, I4=I1,IRAM(IA+KK)-2
            I3=I3+2
            IRAM(I3)=IRAM(I4)
            RAM(I3+1)=RAM(I4+1)*COEF1
  22      CONTINUE
          DO 24, I4=I2,IRAM(IB+KK)-2
            I3=I3+2
            IRAM(I3)=IRAM(I4)
            RAM(I3+1)=RAM(I4+1)*COEF2
  24      CONTINUE
          IRAM(IC+KK)=I3+2
  28    CONTINUE
        NELC=(IRAM(IC+M2)-IRAM(IC))/2
        NC=M2+1+2*NELC
        ISYMC=ISYM1
        NSPARC=NELMAT(M1,M2,ISYMC)-NELC
      ELSE
C       First matrix is non-sparse:
        IF (COEF1.NE.1.) THEN
          DO 40, I1=IA,IA+NEL1-1
            RAM(I1)=RAM(I1)*COEF1
   40     CONTINUE
        ENDIF
        IF (NSPAR2.GE.0) THEN
          DO 44, ICOL=1,M2
            DO 42, I1=IRAM(IB-1+ICOL),IRAM(IB+ICOL)-2,2
              IROW=IRAM(I1)
              IF     (ISYM2.EQ.1) THEN
                I2=IA-1+IROW
              ELSEIF (ISYM2.EQ.2) THEN
                I2=IA-1+(ICOL-1)*ICOL/2+IROW
              ELSEIF (ISYM2.EQ.3) THEN
                I2=IA-1+(ICOL-1)*(ICOL-2)/2+IROW
              ELSE
                I2=IA-1+(ICOL-1)*M1+IROW
              ENDIF
              RAM(I2)=RAM(I2)+RAM(I1+1)*COEF2
   42       CONTINUE
   44     CONTINUE
        ELSE
          DO 50, I1=1,NEL1
            RAM(IA-1+I1)=RAM(IA-1+I1)+RAM(IB-1+I1)*COEF2
   50     CONTINUE
        ENDIF
        IC=1
        NELC=NEL1
        NC=NA
        ISYMC=ISYM1
        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,ISYMC,NSPARC,NELC,1,MRAM,IC,NC)
      ELSEIF ((ISPAR3.EQ.-1).AND.(NSPARC.GE.0)) THEN
C       Changing sparse matrix to non-sparse matrix:
        CALL SGMATR(M1,M2,ISYMC,NSPARC,NELC,1,MRAM,IC,NC)
      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,ISYMC,NSPARC,NELC,1,MRAM,IC,NC,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,ISYMC,NSPARC,NELC,1,MRAM,IC,NC)
          ENDIF
        ENDIF
      ENDIF
C
C     Writing output matrix C:
      SPARSC=' '
      IF (NSPARC.GE.0) SPARSC='CSC'
      CALL WMATH(LU1,FILE3,FILED3,M1,M2,SPARSC,NELC,SYM3,FORM3)
      CALL WMATD(LU1,FILED3,M1,M2,SPARSC,NELC,FORM3,IC)
C
      WRITE(*,'(A)') '+MATLIN: Done.                 '
C
      STOP
      END
C
C=======================================================================
C
      SUBROUTINE SWITCH(
     *           FILE1,FILED1,SPARS1,NSPAR1,SYM1,ISYM1,FORM1,NEL1,COEF1,
     *           FILE2,FILED2,SPARS2,NSPAR2,SYM2,ISYM2,FORM2,NEL2,COEF2)
      CHARACTER*(*) FILE1,FILED1,SPARS1,FILE2,FILED2,SPARS2,SYM1,SYM2,
     *              FORM1,FORM2
      CHARACTER*80 FILE3,FILED3,SPARS3,SYM3,FORM3
      INTEGER NSPAR1,ISYM1,NEL1,NSPAR2,ISYM2,NEL2,NSPAR3,ISYM3,NEL3
      REAL COEF1,COEF2,COEF3
C-----------------------------------------------------------------------
       FILE3 =  FILE1
      FILED3 = FILED1
      SPARS3 = SPARS1
      NSPAR3 = NSPAR1
        SYM3 =   SYM1
       ISYM3 =  ISYM1
       FORM3 =  FORM1
        NEL3 =   NEL1
       COEF3 =  COEF1
       FILE1 =  FILE2
      FILED1 = FILED2
      SPARS1 = SPARS2
      NSPAR1 = NSPAR2
        SYM1 =   SYM2
       ISYM1 =  ISYM2
       FORM1 =  FORM2
        NEL1 =   NEL2
       COEF1 =  COEF2
       FILE2 =  FILE3
      FILED2 = FILED3
      SPARS2 = SPARS3
      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
C
C=======================================================================
C