C
C Program MATMUL to compute product M3=M1*M2 of two matrices M1 and M2. C Matrices M1 and M2 may be transposed before the multiplication. C C Version: 6.10 C Date: 2007, June 13 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 transpose C matrix M1, or to change the form of the matrix M1. 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 SYMMETRY='string' ... Symmetry of the matrix. For SYMMETRY='diag', C only the diagonal components of the matrix are calculated C and written to the output file, other components are not C calculated. Simmilarly for other symmetries. C Possible values of SYMMETRY are: C SYMMETRY='diag' ... diagonal matrix C SYMMETRY='sym' ... symmetric matrix C SYMMETRY='skew' ... antisymmetric matrix C SYMMETRY=' ' ... general matrix C Default: C For matrix M2 not specified (MATIN2=' '): C SYMMETRY equal to symmetry of matrix M1 C For both input matrices M1 and M2 diagonal: C SYMMETRY='diag' (output matrix M3 will be diagonal) C In other cases: C SYMMETRY=' ' (output matrix M3 will be general) 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 Switches identifying transposition of the input matrices: C MATT1=integer ... Identifies whether the matrix M1 is to be C transposed before the multiplication: C MATT1=0 ... no transposition C MATT1=1 ... matrix M1 is to be transposed C Default: MATT1=0 (no transposition of M1) C MATT2=integer ... Same as MATT1, but for matrix M2. C Default: MATT2=0 (no transposition of M2) 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,RMATH,RMATD,WMATH,WMATD, *TMATR,SMATRE,SMATR,GSMAT,GSMATR,SGMATR,ISYM,NELMAT,MSHIFT, *MMILOC,MMIRVE,OMULS,OMULNS INTEGER ISYM,NELMAT C ERROR ... File error.for. C RSEP1,RSEP3I,RSEP3T,SSEP ... C File sep.for. C RMATH,RMATD,WMATH,WMATD ... C TMATR,SMATRE,SMATR,GSMAT,GSMATR,SGMATR,ISYM,NELMAT,MSHIFT ... C File forms.for. C MMILOC,MMIRVE,OMULS,OMULNS ... This file. C Subroutines and external functions referred by above subroutines: C LENGTH ... File length.for. C WSEP3I,WSEP3T ... File sep.for. C RARRAI,RARRAY,WARRAI,WARRAY ... C File forms.for. C INDEXI ... File indexi.for. 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,SYM3I,SYMC,FORM1,FORM2,FORM3 INTEGER MATT1,MATT2,M1,M2,M2B,M3,NEL1,NEL2,NELC,OA,OA1,OAN,NA, *OB,OB1,OBN,NB,OC,NC,MC,NCTMP,NCCOL,OTMP,MTMP, *NSPAR1,NSPAR2,ISPAR3,NSPARC,ISYM1,ISYM2,ISYM3,ISYMC, *ISTORA,ISTORB,ISTORC,LU1, *I1,II,JJ,KK,II1,II2,II3,IIA,JJ1,JJ2,JJ3,JJA,IIKK,OCTMP1 REAL AA,BB,CIK,RSPAR1,RSPAR2,RSPAR3 CHARACTER*72 TXTERR PARAMETER (LU1=1) C C----------------------------------------------------------------------- C C Reading a name of the file with the input data: WRITE(*,'(A)') '+MATMUL: 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 MATMUL-01 CALL ERROR('MATMUL-01: SEP file not given.') ENDIF C C Reading the names of matrices header files: CALL RSEP3T('MATIN1',FILE1,' ') IF (FILE1.EQ.' ') THEN C MATMUL-02 CALL ERROR('MATMUL-02: File MATIN1 not given.') ENDIF CALL RSEP3T('MATIN2',FILE2,' ') CALL RSEP3T('MATOUT',FILE3,' ') FILED3=' ' IF (FILE3.EQ.' ') THEN C MATMUL-03 CALL ERROR('MATMUL-03: File MATOUT not given.') ENDIF C Reading the header files of the input matrices: CALL RMATH(LU1,FILE1,FILED1,M1,M2,NSPAR1,SYM1,FORM1,NEL1) ISYM1=ISYM(SYM1) IF (FILE2.NE.' ') THEN CALL RMATH(LU1,FILE2,FILED2,M2B,M3,NSPAR2,SYM2,FORM2,NEL2) ISYM2=ISYM(SYM2) ENDIF C Reading the properties of the output file: IF (FILE2.EQ.' ') THEN SYM3I=SYM1 ELSEIF ((SYM1.EQ.'diag').AND.(SYM2.EQ.'diag')) THEN SYM3I='diag' ELSE SYM3I=' ' ENDIF CALL RSEP3T('SYMMETRY',SYM3,SYM3I) CALL LOWER(SYM3) ISYM3=ISYM(SYM3) CALL RSEP3T('MATFORM',FORM3,'formatted') CALL LOWER(FORM3) CALL RSEP3I('SPARSE',ISPAR3,0) C Reading the transposition switches: CALL RSEP3I('MATT1',MATT1,0) CALL RSEP3I('MATT2',MATT2,0) C C Reading input matrix A: IF (NSPAR1.GE.0) THEN NA=2*NEL1 OA=1 OTMP=NA+1 MTMP=MRAM-NA ELSE NA=NEL1 OA=MRAM-NA+1 OTMP=1 MTMP=OA-1 ENDIF IF (NA+1.GT.MRAM) THEN C MATMUL-04 WRITE(TXTERR,'(A,I9,A)') * 'MATMUL-04: Array RAM too small,',NA+1-MRAM,' units missing.' CALL ERROR(TXTERR) ENDIF CALL RMATD(LU1,FILED1,NSPAR1,FORM1,NEL1,IRAM(OA),RAM(OA)) ISTORA=0 C Checking whether the non-sparse matrix is to be changed to the C sparse matrix, changing the matrix: IF (NSPAR1.LT.0) THEN RSPAR1=0.5 CALL GSMAT(M1,M2,SYM1,NSPAR1,NEL1,RAM,IRAM,MRAM, * OA,NA,ISTORA,RSPAR1) IF (NSPAR1.GE.0) THEN OTMP=NA+1 MTMP=MRAM-NA ENDIF ELSE RSPAR1=FLOAT(NSPAR1)/FLOAT(NSPAR1+NEL1) ENDIF C Transposing input matrix A: IF (MATT1.EQ.1) THEN CALL TMATR(M1,M2,SYM1,NSPAR1,NEL1,IRAM(OA),RAM(OA), * IRAM(OTMP),RAM(OTMP),MTMP) ENDIF C Extending sparse matrix A: IF ((NSPAR1.GE.0).AND.((SYM1.EQ.'sym').OR.(SYM1.EQ.'skew'))) THEN CALL SMATRE * (M1,M2,SYM1,NSPAR1,NEL1,IRAM(OA),RAM(OA),MRAM,NA,ISTORA) ISYM1=ISYM(SYM1) ENDIF C Reorganizing sparse matrix A in the memory: IF (NSPAR1.GE.0) THEN CALL SMATR(M1,M2,NSPAR1,NEL1,IRAM,RAM,MRAM, * OA,NA,ISTORA) ENDIF C IF (FILE2.EQ.' ') THEN C Matrix A is to be written as output matrix C: M3=M2 IF (SYM3.EQ.SYM1) THEN C No change of matrix symmetry: OC=OA NC=NA NSPARC=NSPAR1 NELC=NEL1 SYMC=SYM1 ISYMC=ISYM1 ISTORC=ISTORA C Copying matrix A to the position of C: IF (NSPARC.GE.0) THEN DO 2, I1=1,M3+1 IRAM(I1)=IRAM(OC-1+I1)-OC+1 2 CONTINUE DO 3, I1=1,NC-M3-1-1,2 IRAM(1+M3+I1)=IRAM(OC+M3+I1) RAM(1+M3+I1+1)=RAM(OC+M3+I1+1) 3 CONTINUE ELSE DO 4, I1=1,NC RAM(I1)=RAM(OC-1+I1) 4 CONTINUE ENDIF OC=1 ELSE C Changing matrix symmetry: IF ((ISYM3.LE.3).AND.(M1.NE.M3)) THEN C MATMUL-05 CALL ERROR('MATMUL-05: Wrong symmetry of output matrix.') C If the output matrix should have symmetry 'diag', 'sym' or C 'skew', it must have the same number of rows and columns. ENDIF OC=1 IF (ISYM3.EQ.2) THEN C Symmetric matrix NC=M1*(M1+1)/2 ELSEIF (ISYM3.EQ.3) THEN C Antisymmetric matrix NC=M1*(M1-1)/2 ELSEIF (ISYM3.EQ.1) THEN C Diagonal matrix NC=M1 ELSE C General matrix NC=M1*M3 ENDIF MC=OA-1 NSPARC=-1 NELC=NC SYMC=SYM3 ISYMC=ISYM3 ISTORC=0 C NCTMP=0 OCTMP1=0 DO 15, KK=1,M3 C Preparing sufficient number of the storage locations: C Number NCCOL of storage locations for column KK: IF (ISYMC.EQ.1) THEN NCCOL=1 ELSEIF (ISYMC.EQ.2) THEN NCCOL=KK ELSEIF (ISYMC.EQ.3) THEN NCCOL=KK-1 ELSE NCCOL=M1 ENDIF IF (NCTMP+NCCOL.GT.MC) THEN C Erasing unneeded columns 1 to KK-1 of A: IF (NSPAR1.LT.0) THEN C OA stays unchanged, MC is redefined: IF (ISYM1.EQ.1) THEN MC=OA+KK-1-1 ELSEIF (ISYM1.EQ.2) THEN MC=OA+(KK-1)*KK/2-1 ELSEIF (ISYM1.EQ.3) THEN MC=OA+(KK-2)*(KK-1)/2-1 ELSE MC=OA+M1*(KK-1)-1 ENDIF ELSE C OA must be changed, addresses of the ends of the columns C must be moved: OAN=IRAM(OA+KK-1)-M2 DO 12, I1=M2,0,-1 IF (I1.GE.KK) THEN IRAM(OAN+I1)=IRAM(OA+I1) ELSE IRAM(OAN+I1)=IRAM(OA+KK-1) ENDIF 12 CONTINUE OA=OAN MC=OA-1 ENDIF IF (NCTMP+NCCOL.GT.MC) THEN C Changing the written part of C to sparse matrix: CALL GSPART(M1,M3,ISYMC,NSPARC,NELC,RAM,IRAM,OA-1, * OC,NCTMP,ISTORC,KK-1,OCTMP1) ENDIF IF (NCTMP+NCCOL.GT.MC) THEN C MATMUL-06 WRITE(TXTERR,'(A,I9,A)') * 'MATMUL-06: Array RAM too small,',NCTMP+NCCOL-MC, * ' units missing.' CALL ERROR(TXTERR) ENDIF ENDIF C Initiating the values in column KK of C: DO 11, I1=NCTMP+1,NCTMP+NCCOL RAM(I1)=0. 11 CONTINUE C Preparing loop over column KK of A: CALL MMILOC(M1,M2,ISYM1,NSPAR1,RAM,IRAM,OA,KK,II1,II2,II3) C Loop over nonzero elements of column KK of A: DO 13, IIA=II1,II2,II3 C Computing index II of row and value AA of element of A: CALL MMIRVE(M1,M2,ISYM1,NSPAR1,RAM,IRAM,OA,KK,IIA,II1,II3, * II,AA) IF ((ISYMC.EQ.1).AND.(II.EQ.KK)) THEN IIKK=OCTMP1+II ELSEIF ((ISYMC.EQ.2).AND.(II.LE.KK)) THEN IIKK=OCTMP1+(KK-1)*KK/2+II ELSEIF ((ISYMC.EQ.3).AND.(II.LT.KK)) THEN IIKK=OCTMP1+(KK-2)*(KK-1)/2+II ELSEIF (ISYMC.EQ.4) THEN IIKK=OCTMP1+(KK-1)*M1+II ELSE GOTO 13 ENDIF C Storing the element: RAM(IIKK)=AA 13 CONTINUE NCTMP=NCTMP+NCCOL 15 CONTINUE IF (NSPARC.GE.0) THEN C Changing the remaining part of C to sparse matrix: CALL GSPART(M1,M3,ISYMC,NSPARC,NELC,RAM,IRAM,OA-1, * OC,NCTMP,ISTORC,M3,OCTMP1) ENDIF ENDIF GOTO 100 ENDIF C C Reading input matrix B: IF (((MATT2.EQ.0).AND.(M2B.NE.M2)).OR. * ((MATT2.EQ.1).AND.(M3 .NE.M2))) THEN C MATMUL-07 CALL ERROR('MATMUL-07: Mismatch in dimensions of the matrices.') ENDIF IF ((ISYM3.LE.3).AND.(((MATT2.EQ.0).AND.(M1.NE.M3)).OR. * ((MATT2.EQ.1).AND.(M2.NE.M3)))) THEN C MATMUL-08 CALL ERROR('MATMUL-08: Wrong symmetry of output matrix.') C If the output matrix should have symmetry 'diag', 'sym' or C 'skew', it must have the same number of rows and columns. ENDIF IF (NSPAR2.GE.0) THEN NB=2*NEL2 OB=1 OTMP=NB+1 MTMP=MRAM-NA-NB ELSE NB=NEL2 OB=MRAM-NA-NB+1 OTMP=1 MTMP=OB-1 ENDIF IF (NB+1.GT.MRAM-NA) THEN C MATMUL-09 WRITE(TXTERR,'(A,I9,A)') * 'MATMUL-09: Array RAM too small,',NB+1-MRAM+NA,' units missing.' CALL ERROR(TXTERR) ENDIF CALL RMATD(LU1,FILED2,NSPAR2,FORM2,NEL2,IRAM(OB),RAM(OB)) ISTORB=0 C Checking whether the non-sparse matrix is to be changed to the C sparse matrix, changing the matrix: IF (NSPAR2.LT.0) THEN RSPAR2=0.5 CALL GSMAT(M2B,M3,SYM2,NSPAR2,NEL2,RAM,IRAM,OA-1, * OB,NB,ISTORB,RSPAR2) IF (NSPAR2.GE.0) THEN OTMP=NB+1 MTMP=MRAM-NA-NB ENDIF ELSE RSPAR2=FLOAT(NSPAR2)/FLOAT(NSPAR2+NEL2) ENDIF C Transposing input matrix B: IF (MATT2.EQ.1) THEN CALL TMATR(M2B,M3,SYM2,NSPAR2,NEL2,IRAM(OB),RAM(OB), * IRAM(OTMP),RAM(OTMP),MTMP) ENDIF C Extending sparse matrix B: IF ((NSPAR2.GE.0).AND.((SYM2.EQ.'sym').OR.(SYM2.EQ.'skew'))) THEN CALL SMATRE * (M2,M3,SYM2,NSPAR2,NEL2,IRAM(OB),RAM(OB),OA-1,NB,ISTORB) ISYM2=ISYM(SYM2) ENDIF C Reorganizing sparse matrix B in the memory: IF (NSPAR2.GE.0) THEN CALL SMATR(M2,M3,NSPAR2,NEL2,IRAM,RAM,OA-1, * OB,NB,ISTORB) ENDIF C C Optimizations for the speed of calculation: IF (((1.-RSPAR1)*(1.-RSPAR2).LE.0.5).AND. * (.NOT.((NSPAR1.LT.0).AND.(ISYM1.EQ.2.OR.ISYM1.EQ.3))).AND. * (.NOT.((NSPAR2.LT.0).AND.(ISYM2.EQ.2.OR.ISYM2.EQ.3)))) THEN C Matrices are to be multiplied as sparse matrices: C Checking if there is enough of memory for multiplication: C Calculating OC for the case that B will be sparse: CALL OMULS(M1,ISYM3,M2,M3,ISYM2,NSPAR2,NEL2,IRAM,RAM,OA-1, * OB,NB,ISTORB,OC) OAN=OA-OC+1 C Shifting OC if A is to be changed to sparse: IF (NSPAR1.LT.0) THEN OC=OC + NA - (2*NINT((1.-RSPAR1)*FLOAT(NA))+M2+1) ENDIF IF (OC.GT.0) THEN C Matrices may be (and will be) changed to sparse matrices C and shifted to optimum position for multiplication: C Shifting B to the beginning of RAM: CALL MSHIFT(M2,M3,ISYM2,NSPAR2,NEL2,IRAM,RAM,MRAM, * OB,NB,ISTORB,1) IF (NSPAR2.LT.0) THEN C Changing B to sparse, form "as in the memory": ISTORB=1 CALL GSMATR(M2,M3,SYM2,NSPAR2,NEL2,RAM,IRAM,OA-1, * OB,NB,ISTORB) ENDIF IF (NSPAR1.LT.0) THEN C Shifting A behind B: CALL MSHIFT(M1,M2,ISYM1,NSPAR1,NEL1,IRAM,RAM,MRAM, * OA,NA,ISTORA,NB+1) C Changing A to sparse, form "as in the memory": ISTORA=1 CALL GSMATR(M1,M2,SYM1,NSPAR1,NEL1,RAM,IRAM,MRAM, * OA,NA,ISTORA) ENDIF C Shifting A to the optimum position for multiplication: CALL MSHIFT(M1,M2,ISYM1,NSPAR1,NEL1,IRAM,RAM,MRAM, * OA,NA,ISTORA,OAN) C Shifting B to the optimum position for multiplication: CALL MSHIFT(M2,M3,ISYM2,NSPAR2,NEL2,IRAM,RAM,MRAM, * OB,NB,ISTORB,OAN-NB) ENDIF ELSE C Matrices are to be multiplied as non-sparse matrices: C Calculating OC for the case that B will be non-sparse: CALL OMULNS(M1,ISYM3,M2,M3,ISYM2,NSPAR2,NEL2,IRAM,RAM,OA-1, * OB,NB,ISTORB,OC) OAN=OA-OC+1 C Shifting OC if A is to be changed to non-sparse: IF (NSPAR1.GE.0) THEN OC=OC + NA - NELMAT(M1,M2,ISYM1) ENDIF C Shifting B to the beginning of RAM: CALL MSHIFT(M2,M3,ISYM2,NSPAR2,NEL2,IRAM,RAM,MRAM, * OB,NB,ISTORB,1) IF (NSPAR2.GE.0) THEN C Changing B to non-sparse: CALL SGMATR(M2,M3,SYM2,NSPAR2,NEL2,RAM,IRAM,OA-1, * OB,NB,ISTORB) ENDIF IF (NSPAR1.GE.0) THEN C Shifting A behind B: CALL MSHIFT(M1,M2,ISYM1,NSPAR1,NEL1,IRAM,RAM,MRAM, * OA,NA,ISTORA,NB+1) C Changing A to non-sparse: CALL SGMATR(M1,M2,SYM1,NSPAR1,NEL1,RAM,IRAM,MRAM, * OA,NA,ISTORA) ENDIF C Shifting A to the optimum position for multiplication: IF (OC.LE.0) OAN=MRAM-NA+1 CALL MSHIFT(M1,M2,ISYM1,NSPAR1,NEL1,IRAM,RAM,MRAM, * OA,NA,ISTORA,OAN) C Shifting B to the optimum position for multiplication: CALL MSHIFT(M2,M3,ISYM2,NSPAR2,NEL2,IRAM,RAM,MRAM, * OB,NB,ISTORB,OAN-NB) ENDIF C C Multiplication: WRITE(*,'(A)') '+MATMUL: Multiplying... ' C * IF ((SYM2.EQ.'diag').AND. * * (((SYM1.EQ.'DIAG').OR.(SYM1.EQ.' ')).AND.(SYM1.EQ.SYM3)))THEN *C Matrix B is diagonal, A is diagonal or general, *C thus C will be A multiplied by the diagonal of B: * OC=OA * NC=NA * NSPARC=NSPAR1 * NELC=NEL1 * SYMC=SYM1 * ISYMC=ISYM1 * ISTORC=ISTORA *C Loop over diagonal of B: * DO 10, KK=1,M2 * WRITE(*,'(2(A,1I9))') '+MATMUL: Multiplying...',KK,' /',M2 *C Calculating the value BB of diagonal element of B: * IF (NSPAR2.GE.0) THEN * IF (IRAM(OB+KK).EQ.IRAM(OB+KK-1)) THEN * BB=0. * ELSE * BB=RAM(IRAM(OB+KK)) * ENDIF * ELSE * BB=RAM(OB-1+KK) * ENDIF *C Preparing loop over column KK of A: * IF (NSPAR1.GE.0) THEN * II1=IRAM(OA+KK-1)+2 * II2=IRAM(OA+KK) * II3=2 * ELSE * IF (SYM1.EQ.'diag') THEN * II1=OA+KK-1 * II2=II1 * II3=1 * ELSEIF (SYM1.EQ.'sym') THEN * II1=OA+(KK-1)*KK/2 * II2=II1+KK-1 * II3=1 * ELSEIF (SYM1.EQ.'skew') THEN * II1=OA+(KK-1)*(KK-2)/2 * II2=II1+KK-2 * II3=1 * ELSE * II1=OA+(KK-1)*M1 * II2=II1+M1-1 * II3=1 * ENDIF * ENDIF *C Loop over column KK of A: * DO 8, II=II1,II2,II3 *C Multiplying: * RAM(II)=RAM(II)*BB * 8 CONTINUE * 10 CONTINUE *C Shifting matrix C: * IF (NSPARC.GE.0) THEN * DO 12, I1=1,M3+1 * IRAM(I1)=IRAM(OC-1+I1)-OC+1 * 12 CONTINUE * DO 13, I1=1,NC-M3-1-1,2 * IRAM(1+M3+I1)=IRAM(OC+M3+I1) * RAM(1+M3+I1+1)=RAM(OC+M3+I1+1) * 13 CONTINUE * ELSE * DO 14, I1=1,NC * RAM(I1)=RAM(OC-1+I1) * 14 CONTINUE * ENDIF * OC=1 * ELSEIF ((SYM1.EQ.'diag').AND. * * (((SYM2.EQ.'DIAG').OR.(SYM2.EQ.' ')).AND.(SYM2.EQ.SYM3)))THEN *C Matrix A is diagonal, B is diagonal or general, *C thus C will be B multiplied by the diagonal of A: * OC=OB * NC=NB * NELC=NEL2 * NSPARC=NSPAR2 * SYMC=SYM2 * ISYMC=ISYM2 * ISTORC=ISTORB *C Loop over all columns of B: * DO 20, KK=1,M3 * WRITE(*,'(2(A,1I9))') '+MATMUL: Multiplying...',KK,' /',M3 *C Preparing loop over column KK of B: * IF (NSPAR2.GE.0) THEN * JJ1=IRAM(OB+KK-1)+2 * JJ2=IRAM(OB+KK) * JJ3=2 * ELSE * IF (SYM2.EQ.'sym') THEN * JJ1=OB+(KK-1)*KK/2 * JJ2=JJ1+KK-1 * JJ3=1 * ELSEIF (SYM2.EQ.'skew') THEN * JJ1=OB+(KK-1)*(KK-2)/2 * JJ2=JJ1+KK-2 * JJ3=1 * ELSE * JJ1=OB+(KK-1)*M2 * JJ2=JJ1+M2-1 * JJ3=1 * ENDIF * ENDIF *C Loop over column KK of B: * DO 18, JJ=JJ1,JJ2,JJ3 *C Calculating index II of row of B: * IF (NSPAR2.GE.0) THEN * II=IRAM(JJ-1) * ELSE * II=(JJ-JJ1+JJ3)/JJ3 * ENDIF *C Calculating the value AA of diagonal element of A: * IF (NSPAR1.GE.0) THEN * IF (IRAM(OA+II).EQ.IRAM(OA+II-1)) THEN * AA=0. * ELSE * AA=RAM(IRAM(OA+II)) * ENDIF * ELSE * AA=RAM(OA-1+II) * ENDIF *C Multiplying: * RAM(JJ)=RAM(JJ)*AA * 18 CONTINUE * 20 CONTINUE *C Shifting matrix C: * IF (NSPARC.GE.0) THEN * DO 22, I1=1,M3+1 * IRAM(I1)=IRAM(OC+I1-1)-OC+1 * 22 CONTINUE * DO 23, I1=1,NC-M3-1-1,2 * IRAM(1+M3+I1)=IRAM(OC+M3+I1) * RAM(1+M3+I1+1)=RAM(OC+M3+I1+1) * 23 CONTINUE * ELSE * DO 24, I1=1,NC * RAM(I1)=RAM(OC+I1-1) * 24 CONTINUE * ENDIF * OC=1 * ELSEIF ((SYM1.EQ.' ').AND.(SYM2.EQ.' ').AND. IF ((NSPAR1.LT.0).AND.(NSPAR2.LT.0).AND. ccc * ((M1.LT.OB).AND.(M1*M3.LT.OA-M2))) THEN * (OC.GT.0)) THEN C Multiplication for non-sparse matrices: OC=1 NC=M1*M3 NSPARC=-1 NELC=NC SYMC=SYM3 ISYMC=ISYM3 ISTORC=0 IIKK=0 II1=1 II2=M1 OA1=OA-1 OB1=OB-1 C Loop over columns of C: DO 27, KK=1,M3 WRITE(*,'(2(A,1I9))') '+MATMUL: Multiplying...',KK,' /',M3 C Loop over lines of column KK of matrix C: IF (ISYMC.EQ.1) THEN II1=KK II2=KK ELSEIF (ISYMC.EQ.2) THEN II2=KK ELSEIF (ISYMC.EQ.3) THEN II2=KK-1 ENDIF DO 26, II=II1,II2 IIKK=IIKK+1 CIK=0. DO 25, JJ=1,M2 C Element of the matrix A: IF (ISYM1.EQ.1) THEN IF (II.EQ.JJ) THEN AA=RAM(OA1+JJ) ELSE GOTO 25 ENDIF ELSEIF (ISYM1.EQ.2) THEN IF (II.LE.JJ) THEN AA=RAM(OA1+(JJ-1)*JJ/2+II) ELSE AA=RAM(OA1+(II-1)*II/2+JJ) ENDIF ELSEIF (ISYM1.EQ.3) THEN IF (II.EQ.JJ) THEN GOTO 25 ELSEIF (II.LT.JJ) THEN AA= RAM(OA1+(JJ-1)*(JJ-2)/2+II) ELSE AA=-RAM(OA1+(II-1)*(II-2)/2+JJ) ENDIF ELSE AA=RAM(OA1+(JJ-1)*M1+II) ENDIF C Element of the matrix B: IF (ISYM2.EQ.1) THEN IF (JJ.EQ.KK) THEN BB=RAM(OB1+KK) ELSE GOTO 25 ENDIF ELSEIF (ISYM2.EQ.2) THEN IF (JJ.LE.KK) THEN BB=RAM(OB1+(KK-1)*KK/2+JJ) ELSE BB=RAM(OB1+(JJ-1)*JJ/2+KK) ENDIF ELSEIF (ISYM2.EQ.3) THEN IF (JJ.EQ.KK) THEN GOTO 25 ELSEIF (JJ.LT.KK) THEN BB= RAM(OB1+(KK-1)*(KK-2)/2+JJ) ELSE BB=-RAM(OB1+(JJ-1)*(JJ-2)/2+KK) ENDIF ELSE BB=RAM(OB1+(KK-1)*M2+JJ) ENDIF CIK=CIK+AA*BB 25 CONTINUE RAM(IIKK)=CIK 26 CONTINUE 27 CONTINUE ELSE C Multiplication for the case that at least one matrix is sparse, C or OC is negative (no space for whole matrix C): OC=1 IF (SYM3.EQ.'sym') THEN C Symmetric matrix NC=M1*(M1+1)/2 ELSEIF (SYM3.EQ.'skew') THEN C Antisymmetric matrix NC=M1*(M1-1)/2 ELSEIF (SYM3.EQ.'diag') THEN C Diagonal matrix NC=M1 ELSE C General matrix NC=M1*M3 ENDIF MC=OB-1 NSPARC=-1 NELC=NC SYMC=SYM3 ISYMC=ISYM3 ISTORC=0 C NCTMP=0 OCTMP1=0 DO 50, KK=1,M3 WRITE(*,'(2(A,1I9))') '+MATMUL: Multiplying...',KK,' /',M3 C Preparing sufficient number of the storage locations: C Number NCCOL of storage locations for column KK: IF (ISYMC.EQ.1) THEN NCCOL=1 ELSEIF (ISYMC.EQ.2) THEN NCCOL=KK ELSEIF (ISYMC.EQ.3) THEN NCCOL=KK-1 ELSE NCCOL=M1 ENDIF IF (NCTMP+NCCOL.GT.MC) THEN C Erasing unneeded columns 1 to KK-1 of B: IF (NSPAR2.LT.0) THEN C OB stays unchanged, MC is redefined: IF (ISYM2.EQ.1) THEN MC=OB+KK-1-1 ELSEIF (ISYM2.EQ.2) THEN MC=OB+(KK-1)*KK/2-1 ELSEIF (ISYM2.EQ.3) THEN MC=OB+(KK-2)*(KK-1)/2-1 ELSE MC=OB+M2*(KK-1)-1 ENDIF ELSE C OB must be changed, addresses of the ends of the columns C must be moved: OBN=IRAM(OB+KK-1)-M3 DO 28, I1=M3,0,-1 IF (I1.GE.KK) THEN IRAM(OBN+I1)=IRAM(OB+I1) ELSE IRAM(OBN+I1)=IRAM(OB+KK-1) ENDIF 28 CONTINUE OB=OBN MC=OB-1 ENDIF IF (NCTMP+NCCOL.GT.MC) THEN C Changing the calculated part of C to sparse matrix: CALL GSPART(M1,M3,ISYMC,NSPARC,NELC,RAM,IRAM,OB-1, * OC,NCTMP,ISTORC,KK-1,OCTMP1) ENDIF IF (NCTMP+NCCOL.GT.MC) THEN C MATMUL-10 WRITE(TXTERR,'(A,I9,A)') * 'MATMUL-10: Array RAM too small,',NCTMP+NCCOL-MC, * ' units missing.' CALL ERROR(TXTERR) ENDIF ENDIF C Initiating the values in column KK of C: DO 29, I1=NCTMP+1,NCTMP+NCCOL RAM(I1)=0. 29 CONTINUE C Preparing loop over column KK of B: CALL MMILOC(M2,M3,ISYM2,NSPAR2,RAM,IRAM,OB,KK,JJ1,JJ2,JJ3) C Loop over nonzero elements of column KK of B: DO 40, JJA=JJ1,JJ2,JJ3 C Computing index JJ of row and value BB of the element of B: CALL MMIRVE(M2,M3,ISYM2,NSPAR2,RAM,IRAM,OB,KK,JJA,JJ1,JJ3, * JJ,BB) C Preparing loop over column JJ of A: CALL MMILOC(M1,M2,ISYM1,NSPAR1,RAM,IRAM,OA,JJ,II1,II2,II3) C Loop over nonzero elements of column JJ of A: DO 30, IIA=II1,II2,II3 C Computing index II of row and value AA of element of A: CALL MMIRVE(M1,M2,ISYM1,NSPAR1,RAM,IRAM,OA,JJ,IIA,II1,II3, * II,AA) IF ((ISYMC.EQ.1).AND.(II.EQ.KK)) THEN IIKK=OCTMP1+II ELSEIF ((ISYMC.EQ.2).AND.(II.LE.KK)) THEN IIKK=OCTMP1+(KK-1)*KK/2+II ELSEIF ((ISYMC.EQ.3).AND.(II.LT.KK)) THEN IIKK=OCTMP1+(KK-2)*(KK-1)/2+II ELSEIF (ISYMC.EQ.4) THEN IIKK=OCTMP1+(KK-1)*M1+II ELSE GOTO 30 ENDIF C Multiplying: RAM(IIKK)=RAM(IIKK) + AA*BB 30 CONTINUE 40 CONTINUE NCTMP=NCTMP+NCCOL 50 CONTINUE IF (NSPARC.GE.0) THEN C Changing the remaining part of C to sparse matrix: CALL GSPART(M1,M3,ISYMC,NSPARC,NELC,RAM,IRAM,OB-1, * OC,NCTMP,ISTORC,M3,OCTMP1) ENDIF 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,M3,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,M3,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,M3,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,M3,SYMC,NSPARC,NELC,RAM,IRAM,MRAM, * OC,NC,ISTORC) ENDIF ENDIF ENDIF C IF ((NSPARC.GE.0).AND.(ISTORC.EQ.1)) THEN C Reorganizing sparse matrix in the memory: CALL SMATR(M1,M3,NSPARC,NELC,IRAM,RAM,MRAM, * OC,NC,ISTORC) ENDIF C C C Writing output matrix C: CALL WMATH(LU1,FILE3,FILED3,M1,M3,NSPARC,SYM3,FORM3) CALL WMATD(LU1,FILED3,NSPARC,FORM3,NELC,IRAM(OC),RAM(OC)) C WRITE(*,'(A)') '+MATMUL: Done. ' C STOP END C C======================================================================= C C C SUBROUTINE MMILOC(M1,M2,ISYM,NSPARS,ARRAY,IARRAY,OARRAY,ICOLUM, * IELEM1,IELEM2,IELEM3) INTEGER M1,M2,ISYM,NSPARS,OARRAY,ICOLUM,IELEM1,IELEM2,IELEM3 REAL ARRAY(*) INTEGER IARRAY(*) C C Subroutine designed to Initiate Loop Over Column ICOLUM of a matrix. C C Input: C M1 ... Number of rows of the matrix. C M2 ... Number of columns of the matrix. C ISYM... Index of the symmetry of the matrix. C NSPARS . Sparseness of the matrix. C ARRAY,IARRAY ... Array with indices and values of the matrix. C Sparse matrix must be stored in the form "as in memory". C OARRAY..Address of the first storage location of the matrix. C ICOLUM..Number of the column under considertation. C For detailed description of storage of matrices in the memory C refer to file mat.for. C C Output: C IELEM1..Address of the value of the first element of column ICOLUM C in array ARRAY. C IELEM2..Address of the value of the last element of column ICOLUM C in array ARRAY. C IELEM3..Step between two consecutive values of the matrix elements C in array ARRAY. C C Coded by Petr Bulant C E-mail: bulant@seis.karlov.mff.cuni.cz C C----------------------------------------------------------------------- C....................................................................... C C Preparing loop over column ICOLUM of the matrix: IF (NSPARS.GE.0) THEN IELEM1=IARRAY(OARRAY+ICOLUM-1)+2 IELEM2=IARRAY(OARRAY+ICOLUM) IELEM3=2 ELSE IF (ISYM.EQ.1) THEN IELEM1=OARRAY+ICOLUM-1 IELEM2=IELEM1 IELEM3=1 ELSEIF (ISYM.EQ.2) THEN IELEM1=OARRAY+(ICOLUM-1)*ICOLUM/2 IELEM2=IELEM1+M1-1 IELEM3=1 ELSEIF (ISYM.EQ.3) THEN IELEM1=OARRAY+(ICOLUM-1)*(ICOLUM-2)/2 IELEM2=IELEM1+M1-1 IELEM3=1 ELSE IELEM1=OARRAY+(ICOLUM-1)*M1 IELEM2=IELEM1+M1-1 IELEM3=1 ENDIF ENDIF RETURN END C C======================================================================= C C C SUBROUTINE MMIRVE(M1,M2,ISYM,NSPARS,ARRAY,IARRAY,OARRAY,ICOLUM, * IELEMA,IELEM1,IELEM3, * IROW,VELEM) INTEGER M1,M2,ISYM,NSPARS,OARRAY,ICOLUM,IELEMA,IELEM1,IELEM3,IROW REAL ARRAY(*),VELEM INTEGER IARRAY(*) C C Subroutine designed to calculate Index of Row and Value of Element C of a matrix. C C Input: C M1 ... Number of rows of the matrix. C M2 ... Number of columns of the matrix. C ISYM... Index of the symmetry of the matrix. C NSPARS. Sparseness of the matrix. C ARRAY,IARRAY ... Array with indices and values of the matrix. C Sparse matrix must be stored in the form "as in memory". C OARRAY..Address of the first storage location of the matrix. C ICOLUM..Number of the column under considertation. C IELEMA..Address of the value of the current element of column C ICOLUM in array ARRAY. C IELEM1..Address of the value of the first element of column ICOLUM C in array ARRAY. C IELEM3..Step between two consecutive values of the matrix elements C in array ARRAY. C For detailed description of storage of matrices in the memory C refer to file mat.for. C C Output: C IROW .. Number of the row corresponding to the matrix element C under considertation. C VELEM...Value of the considered matrix element. C C Coded by Petr Bulant C E-mail: bulant@seis.karlov.mff.cuni.cz C C----------------------------------------------------------------------- C....................................................................... C C Calculating index IROW of row of the matrix: IF (NSPARS.GE.0) THEN IROW=IARRAY(IELEMA-1) ELSE IF (ISYM.EQ.1) THEN IROW=ICOLUM ELSE IROW=(IELEMA-IELEM1+IELEM3)/IELEM3 ENDIF ENDIF C Calculating the value VELEM of element of the matrix: IF (NSPARS.GE.0) THEN VELEM=ARRAY(IELEMA) ELSE IF (ISYM.EQ.2) THEN C 'sym' IF (IROW.LE.ICOLUM) THEN VELEM=ARRAY(IELEMA) ELSE VELEM=ARRAY(OARRAY-1+(IROW-1)*IROW/2+ICOLUM) ENDIF ELSEIF (ISYM.EQ.3) THEN C 'skew' IF (IROW.LT.ICOLUM) THEN VELEM=ARRAY(IELEMA) ELSEIF (IROW.EQ.ICOLUM) THEN VELEM=0. ELSE VELEM=-ARRAY(OARRAY-1+(IROW-1)*(IROW-2)/2+ICOLUM) ENDIF ELSE C 'diag' or ' ' VELEM=ARRAY(IELEMA) ENDIF ENDIF RETURN END C C======================================================================= C C C SUBROUTINE OMULS(M1C,ISYMC,M1,M2,ISYM,NSPAR,NELEM, * IARRAY,ARRAY,MARRAY,OARRAY,NARRAY,ISTOR,OC) INTEGER M1C,ISYMC,M1,M2,ISYM,NSPAR,NELEM,MARRAY,OARRAY,NARRAY, * ISTOR,OC REAL ARRAY(*) INTEGER IARRAY(*) C C Subroutine to estimate optimum origin OC where to write matrix C, C assuming that matrix B given to OMUL will be stored as sparse matrix, C and will be stored so that it will terminate at the same position C where it terminates now. C OC is estimated as minimum over KK=1,M2 of End(KK-1)-M2-AllKKC, C where End(KK-1) is the address of end of row KK-1 of matrix B if C matrix B were stored as sparse matrix in the form "as in the memory", C and AllKKC is the number of elements of columns 1 to KK of matrix C. C C Input: C M1C ... Number of rows of matrix C. C ISYMC.. Index of symmetry of matrix C. C M1... Number of rows of matrix B. C M2... Number of columns of matrix B. C ISYM .. Index of symmetry of matrix B. C NSPAR...Sparseness of matrix B. C NELEM...Number of elements of matrix B. C IARRAY..Components of matrix B. C ARRAY...Components of matrix B. C MARRAY..Dimension of array ARRAY. C OARRAY..First storage location of matrix B. C NARRAY..Number of storage locations for matrix B. C ISTOR...Form of matrix B. C For detailed description of the parameters describing the matrices C refer to file forms.htm. C Detailed description of storage of matrices in the memory C is given above. C C Output: C OC ... Origin where to write matrix C. C C Coded by Petr Bulant C E-mail: bulant@seis.karlov.mff.cuni.cz C C----------------------------------------------------------------------- C C Subroutines and external functions required: EXTERNAL ERROR,NELMAT INTEGER NELMAT C C Local storage location: INTEGER II,KK,JJ,NN,IB,II1,II2,OTMP C C....................................................................... C IF (NSPAR.GE.0.AND.ISTOR.EQ.0) THEN C MATMUL-11 CALL ERROR('MATMUL-11: Wrong ISTOR.') C This option is not yet coded, sparse matrix must be in the form C "as in the memory". ENDIF C IF (NSPAR.LT.0) THEN C Calculating addresses of ends of columns of B if B were stored C as sparse matrix ending at the same position where it ends now: OTMP=OARRAY-M2-1 IF (OTMP.LE.0) THEN C There is no space for ends of rows, thus OC is put to -1: OC=-1 RETURN ENDIF IB=OARRAY+NARRAY-1 JJ=IB IARRAY(OTMP+M2)=JJ II1=1 II2=M1 C Loop over columns of B: DO 12, KK=M2,1,-1 C Loop over lines of B: IF (ISYM.EQ.1) THEN II1=KK II2=KK ELSEIF (ISYM.EQ.2) THEN II2=KK ELSEIF (ISYM.EQ.3) THEN II2=KK-1 ENDIF C Number NN of nonzero elements of column KK: NN=0 DO 10, II=II1,II2 IF (ARRAY(IB).NE.0.) NN=NN+1 IB=IB-1 10 CONTINUE JJ=JJ-2*NN IARRAY(OTMP+KK-1)=JJ 12 CONTINUE ELSE OTMP=OARRAY ENDIF OC=OTMP DO 20, KK=1,M2 OC=MIN0(OC,IARRAY(OTMP+KK-1)-M2-NELMAT(M1C,KK,ISYMC)) 20 CONTINUE RETURN END C C======================================================================= C C C SUBROUTINE OMULNS(M1C,ISYMC,M1,M2,ISYM,NSPAR,NELEM, * IARRAY,ARRAY,MARRAY,OARRAY,NARRAY,ISTOR,OC) INTEGER M1C,ISYMC,M1,M2,ISYM,NSPAR,NELEM,MARRAY,OARRAY,NARRAY, * ISTOR,OC REAL ARRAY(*) INTEGER IARRAY(*) C C Subroutine to estimate optimum origin OC where to write matrix C, C assuming that matrix B given to OMUL will be stored as non-sparse C matrix, and will be stored so that it will terminate at the same C position where it terminates now. C OC is estimated as minimum over KK=1,M2 of End(KK-1)-AllKKC+1, C where End(KK-1) is the address of end of row KK-1 of matrix B if C matrix B were stored as non-sparse matrix, C and AllKKC is the number of elements of columns 1 to KK of matrix C. C C Input: C M1C ... Number of rows of matrix C. C ISYMC.. Index of symmetry of matrix C. C M1... Number of rows of matrix B. C M2... Number of columns of matrix B. C ISYM .. Index of symmetry of matrix B. C NSPAR...Sparseness of matrix B. C NELEM...Number of elements of matrix B. C IARRAY..Components of matrix B. C ARRAY...Components of matrix B. C MARRAY..Dimension of array ARRAY. C OARRAY..First storage location of matrix B. C NARRAY..Number of storage locations for matrix B. C ISTOR...Form of matrix B. C For detailed description of the parameters describing the matrices C refer to file forms.htm. C Detailed description of storage of matrices in the memory C is given above. C C Output: C OC ... Origin where to write matrix C. C C Coded by Petr Bulant C E-mail: bulant@seis.karlov.mff.cuni.cz C C----------------------------------------------------------------------- C C Subroutines and external functions required: EXTERNAL NELMAT INTEGER NELMAT C C Local storage location: INTEGER KK,OTMP C C....................................................................... C IF (NSPAR.LT.0) THEN OTMP=OARRAY ELSE OTMP=OARRAY+NARRAY-NELMAT(M1,M2,ISYM) ENDIF OC=OTMP DO 20, KK=1,M2 OC=MIN0(OC,OTMP+NELMAT(M1,KK-1,ISYM)-NELMAT(M1C,KK,ISYMC)) 20 CONTINUE RETURN END C C======================================================================= C C C SUBROUTINE GSPART(M1,M2,ISYM,NSPAR,NELEM,ARRAY,IARRAY,MARRAY, * OARRAY,NARRAY,ISTOR,ICOL,OTMP) INTEGER M1,M2,ISYM,NSPAR,NELEM,OARRAY,NARRAY,MARRAY,ISTOR, * ICOL,OTMP REAL ARRAY(*) INTEGER IARRAY(*) C C Subroutine designed to consecutively change parts of non-sparse matrix C to the sparse matrix. For the first invocation, it is assumed on input C that non-sparse part of the matrix for columns 1 to ICOL is stored C in ARRAY, and this part is changed to sparse matrix, form "as in the C memory" and written to ARRAY and IARRAY together with M2+1 addresses C of ends of columns. For second and other invocations, it is assumed C that IARRAY and ARRAY contain M2+1 addresses of ends of columns of the C sparse part of the matrix, followed by the sparse part of the matrix C (columns 1 to ICOL0), and finally by non-sparse part of the matrix C (columns ICOL0+1 to ICOL). The non-sparse part is converted to the C sparse and written after the already stored sparse part, and the C addressses of columns ICOL0+1 to M2 are updated. C C Input: C M1... Number of rows of the full matrix. C M2... Number of columns of the full matrix. C ISYM ...Index of symmetry of the matrix. C ARRAY...Components of the given matrix. C IARRAY..Components of the given matrix. C MARRAY..Dimension of array ARRAY. C OARRAY..Address of the first storage location in array ARRAY C used for the matrix. C NARRAY..Number of storage locations for the matrix. C ICOL ...Index of the last calculated column of the matrix. C C Output: C NSPAR...Sparseness of the M1*ICOL matrix. C NELEM...Number of elements of the M1*ICOL matrix. C NARRAY..New number of storage locations for the M1*ICOL matrix. C IARRAY..Components of the M1*ICOL matrix. C ARRAY...Components of the M1*ICOL matrix. C ISTOR...Form of the matrix. C ISTOR=1 ... form "as in the memory" C OTMP ...Fictitious position in ARRAY just before the first storage C location of the matrix if the matrix would stay non-sparse C and would terminate at the same location where it C terminates after its conversion to the sparse matrix. C C Coded by Petr Bulant C E-mail: bulant@seis.karlov.mff.cuni.cz C C----------------------------------------------------------------------- C C Subroutines and external functions required: EXTERNAL ERROR,NSPMAT,NELMAT INTEGER NSPMAT,NELMAT C C Local storage location: INTEGER I1,I2,I3,I4,I20,I11,I12,J1,OAR,NAR,NSPA,NELE,ICOL0 CHARACTER*72 TXTERR SAVE ICOL0 DATA ICOL0/0/ C C....................................................................... C C Beginning OAR and number of elements NAR of the non-sparse part C of ARRAY: IF (ICOL0.EQ.0) THEN OAR=OARRAY NAR=NARRAY ELSE OAR=IARRAY(OARRAY+M2)+1 NAR=NARRAY-(OAR-OARRAY) ENDIF C Number NSPA of zero and NELE of nonzero elements: NSPA=NSPMAT(NAR,ARRAY(OAR)) NELE=NAR-NSPA IF (((ICOL0.EQ.0).AND.(OAR+M2+1+2*NELE.GT.MARRAY)).OR. * ((ICOL0.NE.0).AND.(OAR+ 2*NELE.GT.MARRAY))) THEN C MATMUL-12 IF (ICOL0.EQ.0) THEN I1=OAR+M2+1+2*NELE - MARRAY ELSE I1=OAR+ 2*NELE - MARRAY ENDIF WRITE(TXTERR,'(A,I9,A)') * 'MATMUL-12: Array RAM too small,',I1,' units missing.' CALL ERROR(TXTERR) ENDIF C Moving the nonzero elements of the matrix to the end of ARRAY, C storing the corresponding matrix indices of the elements: J1=MARRAY I3=OARRAY-1+NARRAY DO 20, I2=ICOL,ICOL0+1,-1 IF (ISYM.EQ.1) THEN I11=I2 I12=I2 ELSEIF (ISYM.EQ.2) THEN I11=I2 I12=1 ELSEIF (ISYM.EQ.3) THEN I11=I2-1 I12=1 ELSE I11=M1 I12=1 ENDIF DO 18, I1=I11,I12,-1 IF (ARRAY(I3).NE.0.) THEN IF (J1-1.LT.I3) THEN C MATMUL-13 WRITE(TXTERR,'(A,I9,A)') * 'MATMUL-13: Array RAM too small,',I3-J1+1, * ' units missing.' CALL ERROR(TXTERR) ENDIF ARRAY(J1)=ARRAY(I3) IARRAY(J1-1)=(I2-1)*M1+I1 J1=J1-2 ENDIF I3=I3-1 18 CONTINUE 20 CONTINUE IF (I3.NE.OAR-1) THEN C MATMUL-14 WRITE(TXTERR,'(A,I9,A)') * 'MATMUL-14: Array RAM too small,',OAR-1-I3,' units missing.' CALL ERROR(TXTERR) ENDIF C J1 now points just before the index of the first nonzero element C ISTOR=1 IF (ICOL0.EQ.0) OAR=OARRAY+M2+1 C Rearranging the matrix to the form 'as in the memory': C Moving the indices and the values to the new positions: DO 42, I3=1,NELE ARRAY(OAR-1+2*I3)=ARRAY(J1+2*I3) IARRAY(OAR-1+2*I3-1)=IARRAY(J1+2*I3-1) 42 CONTINUE C Calculating end addresses of the columns, C changing matrix indices of the elements to the row numbers: I20=M2+1 DO 44, I3=NELE,1,-1 I1=MOD(IARRAY(OAR-1+2*I3-1),M1) IF (I1.EQ.0) I1=M1 I2=(IARRAY(OAR-1+2*I3-1)-1)/M1+1 IARRAY(OAR-1+2*I3-1)=I1 IF (I2.NE.I20) THEN DO 43, I4=I20-1,I2,-1 IARRAY(OARRAY+I4)=OAR-1+2*I3 43 CONTINUE I20=I2 ENDIF 44 CONTINUE DO 45, I4=I20-1,ICOL0,-1 IARRAY(OARRAY+I4)=OAR-1 45 CONTINUE C ICOL0=ICOL NARRAY=IARRAY(OARRAY+M2)-OARRAY+1 NELEM=(NARRAY-M2-1)/2 I1=NELMAT(M1,ICOL,ISYM) NSPAR=I1-NELEM OTMP=IARRAY(OARRAY+M2)-I1 C 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