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.70 C Date: 2012, November 7 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 For MATIN2=' ' (no input matrix M2), matrix M1 multiplied C by COEF1 will be written on output. This option may be C used to multiply matrix M1 by a constant, or to change C its form. C Default: MATIN2=' ' C MATOUT='string'... Name of the header file of the output matrix. 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... File sep.for. C UPPER... 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