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