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.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 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 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 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 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 in the CSC format C SPARSE=0... automatic selection of the sparseness: C matrix with half or more zero elements will be sparse C in the CSC format, otherwise non-sparse matrix C SPARSE=-1... normal (not sparse) matrix C Default: SPARSE=0 (automatic selection of the sparseness) C For general 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 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----------------------------------------------------------------------- C Subroutines and external functions required: EXTERNAL ERROR,RSEP1,RSEP3I,RSEP3T,RMATH,RMATD,WMATH,WMATD, *TMATR,SMATRE,GSMAT,GSMATR,SGMATR,ISYM,NELMAT,MSHIFT,VELEM, *MMILOC,MMIRVE,OMULS,OMULNS,LOWER,UPPER INTEGER ISYM,NELMAT REAL VELEM C ERROR... File error.for. C RSEP1,RSEP3I,RSEP3T,SSEP...File sep.for. C RMATH,RMATD,WMATH,WMATD,TMATR,SMATRE,GSMAT,GSMATR,SGMATR, C ISYM,NELMAT,MSHIFT,VELEM... File mat.for. C MMILOC,MMIRVE,OMULS,OMULNS,GSPART... This file. C LOWER,UPPER... File length.for. 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... File forms.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,FORM1,FORM2,FORM3 CHARACTER*80 SPARS1,SPARS2,SPARSC INTEGER MATT1,MATT2,M1,M2,M2B,M3,NEL1,NEL2,NELC,IA,IA1,IAN,NA, *IB,IB1,IBN,NB,IC,NC,MC,ICNEW,NCTMP,NCCOL,NTMP1,NTMP2, *NSPAR1,NSPAR2,ISPAR3,NSPARC,ISYM1,ISYM2,ISYM3,ISYMC,LU1, *I1,II,JJ,KK,II1,II2,II3,IIA,JJ1,JJ2,JJ3,JJA,IIKK,ICTMP1,JJKK,IIJJ REAL AA,BB,CIK,RSPAR1,RSPAR2,RSPAR3,BBB LOGICAL LSPARC CHARACTER*72 TXTERR PARAMETER (LU1=1) C C----------------------------------------------------------------------- C MINRAM=1 MAXRAM=MRAM 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,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 MATMUL-15 CALL ERROR('MATMUL-15: Wrong value of SPARSE of input matrix') C Only SPARSE=' ' or SPARSE='CSC' is allowed. ENDIF IF (FILE2.NE.' ') THEN CALL RMATH(LU1,FILE2,FILED2,M2B,M3,SPARS2,NEL2,SYM2,FORM2) ISYM2=ISYM(SYM2) IF (SPARS2.EQ.' ') THEN NSPAR2=-1 ELSEIF (SPARS2.EQ.'CSC') THEN NSPAR2=NELMAT(M2B,M3,ISYM2)-NEL2 ELSE C MATMUL-16 CALL ERROR('MATMUL-16: Wrong value of SPARSE of input matrix') C Only SPARSE=' ' or SPARSE='CSC' is allowed. ENDIF 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('FORMM',FORM3,'FORMATTED') CALL UPPER(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=M2+1+2*NEL1 ELSE NA=NEL1 ENDIF IA=MAXRAM-NA+1 IF (NA+1.GT.MAXRAM) THEN C MATMUL-04 WRITE(TXTERR,'(A,I9,A)') * 'MATMUL-04: Array RAM too small,',NA+1-MAXRAM,' units missing.' CALL ERROR(TXTERR) C This is the memory needed to read the first matrix M1. ENDIF CALL RMATD(LU1,FILED1,M2,SPARS1,NEL1,FORM1,IA) 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,ISYM1,NSPAR1,NEL1,MINRAM,MAXRAM,IA,NA,RSPAR1) ELSE RSPAR1=FLOAT(NSPAR1)/FLOAT(NSPAR1+NEL1) ENDIF C Transposing input matrix A: IF (MATT1.EQ.1) THEN CALL TMATR(M1,M2,ISYM1,NSPAR1,NEL1,MINRAM,MAXRAM,IA,NA) CALL MSHIFT(M1,M2,ISYM1,NSPAR1,NEL1,MINRAM,MAXRAM,IA,NA, * MAXRAM-NA+1) ENDIF C Extending sparse matrix A: CCC IF ((NSPAR1.GE.0).AND.((SYM1.EQ.'sym').OR.(SYM1.EQ.'skew'))) THEN CCC CALL SMATRE CCC * (M1,M2,SYM1,NSPAR1,NEL1,IRAM(IA),RAM(IA),MAXRAM,NA,ISTORA) CCC ISYM1=ISYM(SYM1) CCC ENDIF C IF (FILE2.EQ.' ') THEN C Matrix A is to be written as output matrix C: M3=M2 IF (ISYM3.EQ.ISYM1) THEN C No change of matrix symmetry: IC=IA NC=NA NSPARC=NSPAR1 NELC=NEL1 ISYMC=ISYM1 C Copying matrix A to the position of C: ICNEW=MINRAM CALL MSHIFT(M1,M3,ISYMC,NSPARC,NELC,MINRAM,MAXRAM,IC,NC,ICNEW) 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 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 IC=MINRAM MC=IA-1 NSPARC=-1 NELC=NC ISYMC=ISYM3 C NCTMP=0 ICTMP1=IC-1 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 (IC-1+NCTMP+NCCOL.GT.MC) THEN C Erasing unneeded columns 1 to KK-1 of A: IF (NSPAR1.LT.0) THEN C IA stays unchanged, MC is redefined: IF (ISYM1.EQ.1) THEN MC=IA+KK-1-1 ELSEIF (ISYM1.EQ.2) THEN MC=IA+(KK-1)*KK/2-1 ELSEIF (ISYM1.EQ.3) THEN MC=IA+(KK-2)*(KK-1)/2-1 ELSE MC=IA+M1*(KK-1)-1 ENDIF ELSE C IA must be changed, addresses of the ends of the columns C must be moved: IAN=IRAM(IA+KK-1)-M2-1 DO 12, I1=M2,0,-1 IF (I1.GE.KK-1) THEN IRAM(IAN+I1)=IRAM(IA+I1) ELSE IRAM(IAN+I1)=IRAM(IA+KK-1) ENDIF 12 CONTINUE IA=IAN MC=IA-1 ENDIF IF (IC-1+NCTMP+NCCOL.GT.MC) THEN C Changing the written part of C to sparse matrix: CALL GSPART(M1,M3,ISYMC,NSPARC,NELC,MINRAM,IA-1, * IC,NCTMP,KK-1,ICTMP1) ENDIF IF (IC-1+NCTMP+NCCOL.GT.MC) THEN C MATMUL-06 WRITE(TXTERR,'(A,I9,A)') * 'MATMUL-06: Array RAM too small,',IC-1+NCTMP+NCCOL-MC, * ' units missing.' CALL ERROR(TXTERR) ENDIF ENDIF C Initiating the values in column KK of C: DO 11, I1=IC-1+NCTMP+1,IC-1+NCTMP+NCCOL RAM(I1)=0. 11 CONTINUE C Preparing loop over column KK of A: CALL MMILOC(M1,M2,ISYM1,NSPAR1,IA,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,IA,KK,IIA,II1,II3,II,AA) IF ((ISYMC.EQ.1).AND.(II.EQ.KK)) THEN IIKK=ICTMP1+II ELSEIF ((ISYMC.EQ.2).AND.(II.LE.KK)) THEN IIKK=ICTMP1+(KK-1)*KK/2+II ELSEIF ((ISYMC.EQ.3).AND.(II.LT.KK)) THEN IIKK=ICTMP1+(KK-2)*(KK-1)/2+II ELSEIF (ISYMC.EQ.4) THEN IIKK=ICTMP1+(KK-1)*M1+II ELSE GOTO 13 ENDIF C Storing the element: RAM(IIKK)=AA IF (((ISYM1.EQ.2).OR.(ISYM1.EQ.3)).AND.(II.NE.KK).AND. * (ISYMC.EQ.4)) THEN C A is (anti)symmetric, storing also Aki=(-)Aik: IIKK=ICTMP1+(II-1)*M1+KK IF (ISYM1.EQ.3) AA=-AA RAM(IIKK)=AA ENDIF 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,MINRAM,MAXRAM, * IC,NCTMP,M3,ICTMP1) 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=M3+1+2*NEL2 ELSE NB=NEL2 ENDIF IB=IA-NB IF (NB+1.GT.MAXRAM-MINRAM-NA) THEN C MATMUL-09 WRITE(TXTERR,'(A,I9,A)') * 'MATMUL-09: Array RAM too small,',NB+1-(MAXRAM-MINRAM-NA), * ' units missing.' CALL ERROR(TXTERR) ENDIF CALL RMATD(LU1,FILED2,M3,SPARS2,NEL2,FORM2,IB) 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,ISYM2,NSPAR2,NEL2,MINRAM,IA-1,IB,NB,RSPAR2) ELSE RSPAR2=FLOAT(NSPAR2)/FLOAT(NSPAR2+NEL2) ENDIF C Transposing input matrix B: IF (MATT2.EQ.1) THEN CALL TMATR(M2B,M3,ISYM2,NSPAR2,NEL2,MINRAM,IA-1,IB,NB) CALL MSHIFT(M2B,M3,ISYM2,NSPAR2,NEL2,MINRAM,IA-1,IB,NB,IA-NB) ENDIF C Extending sparse matrix B: IF ((NSPAR1.GE.0).AND.((SYM1.EQ.'sym').OR.(SYM1.EQ.'skew')).AND. * ((SYM2.EQ.'sym').OR.(SYM2.EQ.'skew'))) THEN CALL SMATRE(M2,M3,ISYM2,NSPAR2,NEL2,MINRAM,IA-1,IB,NB) SYM2=' ' 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 IC for the case that B will be sparse: CALL OMULS(M1,ISYM3,M2,M3,ISYM2,NSPAR2,NEL2,MINRAM,IA-1, * IB,NB,IC) IAN=MINRAM+IA-IC C Shifting IC if A is to be changed to sparse: IF (NSPAR1.LT.0) THEN IC=IC + NA - (2*NINT((1.-RSPAR1)*FLOAT(NA))+M2+1) ENDIF IF (IC.GE.MINRAM) THEN C Matrices may be (and will be) changed to sparse matrices C and shifted to optimum position for multiplication: IF (NSPAR2.LT.0) THEN C Changing B to sparse: CALL GSMATR(M2,M3,ISYM2,NSPAR2,NEL2,MINRAM,IA-1,IB,NB) ENDIF C Shifting B to the beginning of RAM: CALL MSHIFT(M2,M3,ISYM2,NSPAR2,NEL2,MINRAM,IA-1,IB,NB,MINRAM) IF (NSPAR1.LT.0) THEN C Changing A to sparse: CALL GSMATR(M1,M2,ISYM1,NSPAR1,NEL1,IB+NB,MAXRAM,IA,NA) ENDIF C Shifting A to the optimum position for multiplication: CALL MSHIFT(M1,M2,ISYM1,NSPAR1,NEL1,IB+NB,MAXRAM,IA,NA,IAN) C Shifting B to the optimum position for multiplication: CALL MSHIFT(M2,M3,ISYM2,NSPAR2,NEL2,MINRAM,IA-1,IB,NB,IA-NB) ENDIF ELSE C Matrices are to be multiplied as non-sparse matrices: C Checking the memory needed for conversion of matrices C into non-sparse form: NTMP1=0 NTMP2=0 IF (NSPAR1.GE.0) NTMP1=M1 IF (NSPAR2.GE.0) NTMP2=M2 IF ((NELMAT(M1,M2,ISYM1)+NELMAT(M2,M3,ISYM2)+MAX0(NTMP1,NTMP2)) * .LE.MAXRAM) THEN C There is enough RAM for both matrices non-sparse plus the C auxiliary array to change them to non-sparse. C Calculating IC for the case that B will be non-sparse: CALL OMULNS(M1,ISYM3,M2,M3,ISYM2,NSPAR2,NEL2,MINRAM,IA-1, * IB,NB,IC) IAN=MINRAM+IA-IC C Shifting IC if A is to be changed to non-sparse: IF (NSPAR1.GE.0) THEN IC=IC + NA - NELMAT(M1,M2,ISYM1) ENDIF C Shifting B to the beginning of RAM: CALL MSHIFT(M2,M3,ISYM2,NSPAR2,NEL2,MINRAM,MAXRAM, * IB,NB,MINRAM) IF (NSPAR2.GE.0) THEN C Changing B to non-sparse: CALL SGMATR(M2,M3,ISYM2,NSPAR2,NEL2,MINRAM,IA-1,IB,NB) ENDIF IF (NSPAR1.GE.0) THEN C Shifting A behind B: CALL MSHIFT(M1,M2,ISYM1,NSPAR1,NEL1,IB+NB,MAXRAM, * IA,NA,IB+NB) C Changing A to non-sparse: CALL SGMATR(M1,M2,ISYM1,NSPAR1,NEL1,IB+NB,MAXRAM,IA,NA) ENDIF C Shifting A to the optimum position for multiplication: IF (IC.LT.MINRAM) IAN=MAXRAM-NA+1 CALL MSHIFT(M1,M2,ISYM1,NSPAR1,NEL1,IB+NB,MAXRAM,IA,NA,IAN) C Shifting B to the optimum position for multiplication: CALL MSHIFT(M2,M3,ISYM2,NSPAR2,NEL2,MINRAM,IA-1,IB,NB,IA-NB) ENDIF ENDIF C C Multiplication: WRITE(*,'(A)') '+MATMUL: Multiplying... ' C IF ((NSPAR1.LT.0).AND.(NSPAR2.LT.0).AND. ccc * ((M1.LT.IB).AND.(M1*M3.LT.IA-M2))) THEN * (IC.GE.MINRAM)) THEN C Multiplication for non-sparse matrices: IC=MINRAM NC=M1*M3 NSPARC=-1 NELC=NC ISYMC=ISYM3 IIKK=IC-1 II1=1 II2=M1 IA1=IA-1 IB1=IB-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(IA1+JJ) ELSE GOTO 25 ENDIF ELSEIF (ISYM1.EQ.2) THEN IF (II.LE.JJ) THEN AA=RAM(IA1+(JJ-1)*JJ/2+II) ELSE AA=RAM(IA1+(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(IA1+(JJ-1)*(JJ-2)/2+II) ELSE AA=-RAM(IA1+(II-1)*(II-2)/2+JJ) ENDIF ELSE AA=RAM(IA1+(JJ-1)*M1+II) ENDIF C Element of the matrix B: IF (ISYM2.EQ.1) THEN IF (JJ.EQ.KK) THEN BB=RAM(IB1+KK) ELSE GOTO 25 ENDIF ELSEIF (ISYM2.EQ.2) THEN IF (JJ.LE.KK) THEN BB=RAM(IB1+(KK-1)*KK/2+JJ) ELSE BB=RAM(IB1+(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(IB1+(KK-1)*(KK-2)/2+JJ) ELSE BB=-RAM(IB1+(JJ-1)*(JJ-2)/2+KK) ENDIF ELSE BB=RAM(IB1+(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 IC is negative (no space for whole matrix C): IC=MINRAM LSPARC=.TRUE. 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=IB-1 NSPARC=-1 NELC=NC ISYMC=ISYM3 C NCTMP=0 ICTMP1=IC-1 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 (IC-1+NCTMP+NCCOL.GT.MC) THEN C Erasing unneeded columns 1 to KK-1 of B: IF (NSPAR2.LT.0) THEN C IB stays unchanged, MC is redefined: IF (ISYM2.EQ.1) THEN MC=IB+KK-1-1 ELSEIF (ISYM2.EQ.2) THEN MC=IB+(KK-1)*KK/2-1 ELSEIF (ISYM2.EQ.3) THEN MC=IB+(KK-2)*(KK-1)/2-1 ELSE MC=IB+M2*(KK-1)-1 ENDIF ELSE C IB must be changed, addresses of the ends of the columns C must be moved: IBN=IRAM(IB-1+KK)-(M3+1) DO 28, I1=M3,0,-1 IF (I1.GE.KK-1) THEN IRAM(IBN+I1)=IRAM(IB+I1) ELSE IRAM(IBN+I1)=IRAM(IB-1+KK) ENDIF 28 CONTINUE IB=IBN MC=IB-1 ENDIF IF ((IC-1+NCTMP+NCCOL.GT.MC).AND.(LSPARC)) THEN C Changing the calculated part of C to sparse matrix: CALL GSPART(M1,M3,ISYMC,NSPARC,NELC,MINRAM,IB-1, * IC,NCTMP,KK-1,ICTMP1) ENDIF IF (IC-1+NCTMP+NCCOL.GT.MC) THEN C MATMUL-10 WRITE(TXTERR,'(A,I9,A)') * 'MATMUL-10: Array RAM too small,',IC-1+NCTMP+NCCOL-MC, * ' units missing.' CALL ERROR(TXTERR) ENDIF ENDIF C Initiating the values in column KK of C: DO 29, I1=IC-1+NCTMP+1,IC-1+NCTMP+NCCOL RAM(I1)=0. 29 CONTINUE C Preparing loop over column KK of B: CALL MMILOC(M2,M3,ISYM2,NSPAR2,IB,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,IB,KK,JJA,JJ1,JJ3,JJ,BB) C Preparing loop over column JJ of A: CALL MMILOC(M1,M2,ISYM1,NSPAR1,IA,JJ,II1,II2,II3) C Loop over nonzero elements of column JJ of A: DO 32, IIA=II1,II2,II3 C Computing index II of row and value AA of element of A: CALL MMIRVE(M1,M2,ISYM1,NSPAR1,IA,JJ,IIA,II1,II3,II,AA) IF ((ISYMC.EQ.1).AND.(II.EQ.KK)) THEN IIKK=ICTMP1+II ELSEIF ((ISYMC.EQ.2).AND.(II.LE.KK)) THEN IIKK=ICTMP1+(KK-1)*KK/2+II ELSEIF ((ISYMC.EQ.3).AND.(II.LT.KK)) THEN IIKK=ICTMP1+(KK-2)*(KK-1)/2+II ELSEIF (ISYMC.EQ.4) THEN IIKK=ICTMP1+(KK-1)*M1+II ELSE GOTO 30 ENDIF C Multiplying Aij*Bjk: RAM(IIKK)=RAM(IIKK) + AA*BB 30 CONTINUE IF (((ISYM1.EQ.2).OR.(ISYM1.EQ.3)).AND.(II.NE.JJ)) THEN C A is (anti)symmetric, multiplying also Aji=(-)Aij: IF ((ISYMC.EQ.1).AND.(JJ.EQ.KK)) THEN JJKK=ICTMP1+JJ ELSEIF ((ISYMC.EQ.2).AND.(JJ.LE.KK)) THEN JJKK=ICTMP1+(KK-1)*KK/2+JJ ELSEIF ((ISYMC.EQ.3).AND.(JJ.LT.KK)) THEN JJKK=ICTMP1+(KK-2)*(KK-1)/2+JJ ELSEIF (ISYMC.EQ.4) THEN JJKK=ICTMP1+(KK-1)*M1+JJ ELSE GOTO 32 ENDIF IF (ISYM1.EQ.3) AA=-AA BBB=VELEM(M2,M3,ISYM2,NSPAR2,IB,II,KK) C Multiplying Aji*Bik: RAM(JJKK)=RAM(JJKK) + AA*BBB ENDIF 32 CONTINUE IF (((ISYM2.EQ.2).OR.(ISYM2.EQ.3)).AND.(JJ.NE.KK)) THEN C B is (anti)symmetric, multiplying also Bkj=(-)Bjk: IF (ISYM2.EQ.3) BB=-BB C In this case matrix C may not be changed to sparse: LSPARC=.FALSE. C Preparing loop over column KK of A: CALL MMILOC(M1,M2,ISYM1,NSPAR1,IA,KK,II1,II2,II3) C Loop over nonzero elements stored for column KK of A: DO 34, IIA=II1,II2,II3 C Computing index II of row and value AA of element of A: CALL MMIRVE(M1,M2,ISYM1,NSPAR1,IA,KK,IIA,II1,II3,II,AA) IF ((ISYMC.EQ.1).AND.(II.EQ.JJ)) THEN IIJJ=ICTMP1+II ELSEIF ((ISYMC.EQ.2).AND.(II.LE.JJ)) THEN IIJJ=ICTMP1+(JJ-1)*JJ/2+II ELSEIF ((ISYMC.EQ.3).AND.(II.LT.JJ)) THEN IIJJ=ICTMP1+(JJ-2)*(JJ-1)/2+II ELSEIF (ISYMC.EQ.4) THEN IIJJ=ICTMP1+(JJ-1)*M1+II ELSE GOTO 34 ENDIF C Multiplying Aik*Bkj: RAM(IIJJ)=RAM(IIJJ) + AA*BB 34 CONTINUE C If A is (anti)symmetric dense, multiplying the rest C of column KK of A under the diagonal. C (Here B is (anti)symmetric and thus A cannot be C (anti)symmetric sparse.) IF ((ISYM1.EQ.2).OR.(ISYM1.EQ.3)) THEN DO 36, II=KK+1,M1 IF ((ISYMC.EQ.1).AND.(II.EQ.JJ)) THEN IIJJ=ICTMP1+II ELSEIF ((ISYMC.EQ.2).AND.(II.LE.JJ)) THEN IIJJ=ICTMP1+(JJ-1)*JJ/2+II ELSEIF ((ISYMC.EQ.3).AND.(II.LT.JJ)) THEN IIJJ=ICTMP1+(JJ-2)*(JJ-1)/2+II ELSEIF (ISYMC.EQ.4) THEN IIJJ=ICTMP1+(JJ-1)*M1+II ELSE GOTO 36 ENDIF IF (ISYM1.EQ.2) THEN C Aik=Aki AA=RAM(IA-1+(II-1)*II/2+KK) C Multiplying Aik*Bkj: RAM(IIJJ)=RAM(IIJJ) + AA*BB ELSE C Aik=-Aki AA=-RAM(IA-1+(II-1)*(II-2)/2+KK) C Multiplying Aik*Bkj: RAM(IIJJ)=RAM(IIJJ) + AA*BB ENDIF 36 CONTINUE ENDIF ENDIF 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,MINRAM,MAXRAM, * IC,NCTMP,M3,ICTMP1) 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,ISYMC,NSPARC,NELC,MINRAM,MAXRAM,IC,NC) ELSEIF ((ISPAR3.EQ.-1).AND.(NSPARC.GE.0)) THEN C Changing sparse matrix to non-sparse matrix: CALL SGMATR(M1,M3,ISYMC,NSPARC,NELC,MINRAM,MAXRAM,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,M3,ISYMC,NSPARC,NELC,MINRAM,MAXRAM,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,M3,ISYMC,NSPARC,NELC,MINRAM,MAXRAM,IC,NC) ENDIF ENDIF ENDIF C C Writing output matrix C: SPARSC=' ' IF (NSPARC.GE.0) SPARSC='CSC' CALL WMATH(LU1,FILE3,FILED3,M1,M3,SPARSC,NELC,SYM3,FORM3) CALL WMATD(LU1,FILED3,M1,M3,SPARSC,NELC,FORM3,IC) C WRITE(*,'(A)') '+MATMUL: Done. ' C STOP END C C======================================================================= C C C SUBROUTINE MMILOC(M1,M2,ISYM,NSPARS,IA,ICOLUM, * IELEM1,IELEM2,IELEM3) INTEGER M1,M2,ISYM,NSPARS,IA,ICOLUM,IELEM1,IELEM2,IELEM3 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 IA... 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 C ICOLUM in array ARRAY. C IELEM2... Address of the value of the last element of column C ICOLUM in array ARRAY. C IELEM3... Step between two consecutive values of the matrix C elements in array ARRAY. C C Coded by Petr Bulant C C----------------------------------------------------------------------- C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C INTEGER IRAM(MRAM) EQUIVALENCE (IRAM,RAM) C IRAM,RAM... Indices of elements and elements of the matrix. C....................................................................... C C Preparing loop over column ICOLUM of the matrix: IF (NSPARS.GE.0) THEN IELEM1=IRAM(IA-1+ICOLUM)+1 IELEM2=IRAM(IA+ICOLUM)-1 IELEM3=2 ELSE IF (ISYM.EQ.1) THEN IELEM1=IA+ICOLUM-1 IELEM2=IELEM1 IELEM3=1 ELSEIF (ISYM.EQ.2) THEN IELEM1=IA+(ICOLUM-1)*ICOLUM/2 IELEM2=IELEM1+M1-1 IELEM3=1 ELSEIF (ISYM.EQ.3) THEN IELEM1=IA+(ICOLUM-1)*(ICOLUM-2)/2 IELEM2=IELEM1+M1-1 IELEM3=1 ELSE IELEM1=IA+(ICOLUM-1)*M1 IELEM2=IELEM1+M1-1 IELEM3=1 ENDIF ENDIF RETURN END C C======================================================================= C C C SUBROUTINE MMIRVE(M1,M2,ISYM,NSPARS,IA,ICOLUM, * IELEMA,IELEM1,IELEM3, * IROW,VELEM) INTEGER M1,M2,ISYM,NSPARS,IA,ICOLUM,IELEMA,IELEM1,IELEM3,IROW REAL VELEM 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 IA... Address of the first storage location of the matrix. C ICOLUM... Number of the column under consideration. C IELEMA... Address of the value of the current element of column C ICOLUM in array RAM. C IELEM1... Address of the value of the first element of column C ICOLUM in array RAM. C IELEM3... Step between two consecutive values of the matrix C elements in array RAM. 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 C----------------------------------------------------------------------- C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C INTEGER IRAM(MRAM) EQUIVALENCE (IRAM,RAM) C IRAM,RAM... Indices of elements and elements of the matrix. C....................................................................... C C Calculating index IROW of row of the matrix: IF (NSPARS.GE.0) THEN IROW=IRAM(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=RAM(IELEMA) ELSE IF (ISYM.EQ.2) THEN C 'sym' IF (IROW.LE.ICOLUM) THEN VELEM=RAM(IELEMA) ELSE VELEM=RAM(IA-1+(IROW-1)*IROW/2+ICOLUM) ENDIF ELSEIF (ISYM.EQ.3) THEN C 'skew' IF (IROW.LT.ICOLUM) THEN VELEM=RAM(IELEMA) ELSEIF (IROW.EQ.ICOLUM) THEN VELEM=0. ELSE VELEM=-RAM(IA-1+(IROW-1)*(IROW-2)/2+ICOLUM) ENDIF ELSE C 'diag' or ' ' VELEM=RAM(IELEMA) ENDIF ENDIF RETURN END C C======================================================================= C C C SUBROUTINE OMULS(M1C,ISYMC,M1,M2,ISYM,NSPAR,NELEM, * MIA,MAA,IA,NA,IC) INTEGER M1C,ISYMC,M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,IC C C Subroutine to estimate optimum origin IC 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 IC 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 MIA,MAA... Minimum and maximum address of arrays (I)RAM C available for the subroutine. C IA... Address of the first storage location of matrix B. C NA... Number of storage locations of matrix B. C For detailed description of storage of matrices in the memory C refer to file mat.for. C C Output: C IC... Address of the origin where to write matrix C. C C Coded by Petr Bulant C C----------------------------------------------------------------------- C C Subroutines and external functions required: EXTERNAL ERROR,NELMAT INTEGER NELMAT C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C INTEGER IRAM(MRAM) EQUIVALENCE (IRAM,RAM) C IRAM,RAM... Indices of elements and elements of the matrix, see C the detailed description of storage of matrices in C the memory given above. C C Local storage location: INTEGER II,KK,JJ,NN,JB,II1,II2,ITMP C C....................................................................... 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: ITMP=IA-M2-1 IF (ITMP.LE.MIA) THEN C There is no space for ends of rows, thus IC is put to -1: IC=-1 RETURN ENDIF JB=IA+NA-1 JJ=JB IRAM(ITMP+M2)=JJ+1 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 (RAM(JB).NE.0.) NN=NN+1 JB=JB-1 10 CONTINUE JJ=JJ-2*NN IRAM(ITMP+KK-1)=JJ+1 12 CONTINUE ELSE ITMP=IA ENDIF IC=ITMP DO 20, KK=1,M2 IC=MIN0(IC,IRAM(ITMP-1+KK)-M2-1-NELMAT(M1C,KK,ISYMC)) 20 CONTINUE RETURN END C C======================================================================= C C C SUBROUTINE OMULNS(M1C,ISYMC,M1,M2,ISYM,NSPAR,NELEM, * MIA,MAA,IA,NA,IC) INTEGER M1C,ISYMC,M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,IC C C Subroutine to estimate optimum origin IC where to write matrix C, C assuming that matrix B given to OMULNS 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 IC 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 MIA,MAA... Minimum and maximum address of arrays (I)RAM C available for the subroutine. C IA... Address of the first storage location of matrix B. C NA... Number of storage locations of matrix B. C For detailed description of storage of matrices in the memory C refer to file mat.for. C C Output: C IC... Origin where to write matrix C. C C Coded by Petr Bulant C C----------------------------------------------------------------------- C C Subroutines and external functions required: EXTERNAL NELMAT INTEGER NELMAT C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C INTEGER IRAM(MRAM) EQUIVALENCE (IRAM,RAM) C IRAM,RAM... Indices of elements and elements of the matrix, see C the detailed description of storage of matrices in C the memory given above. C C Local storage location: INTEGER KK,ITMP C C....................................................................... C IF (NSPAR.LT.0) THEN ITMP=IA ELSE ITMP=IA+NA-NELMAT(M1,M2,ISYM) ENDIF IC=ITMP DO 20, KK=1,M2 IC=MIN0(IC,ITMP+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,MIA,MAA,IA,NA,ICOL,ITMP) INTEGER M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,ICOL,ITMP 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 RAM, and this part is changed to sparse matrix and written C to RAM and IRAM together with M2+1 addresses of beginnings C of columns. For second and other invocations, it is assumed that C IRAM and RAM contain M2+1 addresses of beginnings 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+1 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 SYM='diag' ... ISYM=1 C SYM='sym' ... ISYM=2 C SYM='skew' ... ISYM=3 C SYM=' ' ... ISYM=4 C NSPAR... Sparseness of the matrix. C NELEM... Number of elements of the matrix stored in array RAM. C MIA,MAA... Minimum and maximum address of arrays (I)RAM C available for the subroutine. Entire arrays (I)RAM C from MIA to MAA may be used for temporary storage. C IA... Address of the first storage location in array RAM C used for the matrix. C NA... Number of storage locations for the input 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 NA... New number of storage locations for the M1*ICOL matrix. C ITMP... Fictitious position in RAM 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 http://sw3d.cz/staff/bulant.htm C C----------------------------------------------------------------------- C C Subroutines and external functions required: EXTERNAL ERROR,NSPMAT,NELMAT INTEGER NSPMAT,NELMAT C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C INTEGER IRAM(MRAM) EQUIVALENCE (IRAM,RAM) C IRAM,RAM... Indices of elements and elements of the matrix, see C the detailed description of storage of matrices in C the memory given above. C C Local storage location: INTEGER I1,I2,I3,I4,I11,I12,IAR,IAS,NAR,NSPA,NELE,ICOL0,ISHIFT CHARACTER*72 TXTERR SAVE ICOL0 DATA ICOL0/0/ C C....................................................................... C C Beginning IAR and number of elements NAR of the non-sparse part C of RAM: IF (ICOL0.EQ.0) THEN IAR=IA NAR=NA ELSE IAR=IRAM(IA+M2) NAR=NA-(IAR-IA) ENDIF C Number NSPA of zero and NELE of nonzero elements: NSPA=NSPMAT(NAR,RAM(IAR)) NELE=NAR-NSPA IF (((ICOL0.EQ.0).AND.(IAR+M2+1+2*NELE.GT.MAA)).OR. * ((ICOL0.NE.0).AND.(IAR+ 2*NELE.GT.MAA))) THEN C MATMUL-12 IF (ICOL0.EQ.0) THEN I1=IAR+M2+1+2*NELE - MAA ELSE I1=IAR+ 2*NELE - MAA ENDIF WRITE(TXTERR,'(A,I9,A)') * 'MATMUL-12: Array RAM too small,',I1,' units missing.' CALL ERROR(TXTERR) ENDIF C Moving the dense part of the matrix to the end of available RAM: DO 10, I1=0,NAR-1 RAM(MAA-I1)=RAM(IAR+NAR-1-I1) 10 CONTINUE IAR=MAA-(NAR-1) C Beginning IAS where to write the new sparse part of the matrix, C shifting index ISHIFT for case of ICOL0.eq.0: IF (ICOL0.EQ.0) THEN IAS=IA+ICOL+1 ISHIFT=M2+1-(ICOL+1) IRAM(IA)=IAS+ISHIFT ELSE IAS=IRAM(IA+M2) ISHIFT=0 ENDIF C Moving the nonzero elements of the matrix, C storing their row indices, recording column indices: I3=IAR I4=IAS C I3 points to the dense array, C I4 points to the new part of the sparse array DO 20, I2=ICOL0+1,ICOL IF (ISYM.EQ.1) THEN I12=I2 I11=I2 ELSEIF (ISYM.EQ.2) THEN I12=I2 I11=1 ELSEIF (ISYM.EQ.3) THEN I12=I2-1 I11=1 ELSE I12=M1 I11=1 ENDIF DO 18, I1=I11,I12 IF (RAM(I3).NE.0.) THEN IF (I4+1.LE.I3) THEN C MATMUL-13 CALL ERROR('MATMUL-13: Array RAM too small') C The output matrix is being formed in the memory, C but there was not enough space in RAM for it. C Thus the already formed part of the matrix is being C converted to the sparse matrix. There is enough space for C the sparse matrix, but the conversion is not possible. C May be that only several more memory locations will solve C the problem, try to enlarge MRAM, if possible. ENDIF IRAM(I4)=I1 RAM(I4+1)=RAM(I3) I4=I4+2 ENDIF I3=I3+1 18 CONTINUE IRAM(IA-1+I2+1)=I4+ISHIFT 20 CONTINUE C C I4 now points just after the last nonzero element IF (ISHIFT.NE.0) THEN C Shifting the sparse part of the matrix: IF (I4+ISHIFT.GT.MAA) THEN C MATMUL-14 WRITE(TXTERR,'(A,I9,A)') 'MATMUL-14: Array RAM too small,', * I4+ISHIFT-MAA,' units missing.' CALL ERROR(TXTERR) ENDIF DO 30, I1=I4-2,IAS,-2 RAM(I1+1+ISHIFT)=RAM(I1+1) IRAM(I1+ISHIFT)=IRAM(I1) 30 CONTINUE ENDIF C C Recording pointers for columns ICOL+2 to M2+1: DO 32, I1=IA-1+ICOL+2,IA-1+M2+1 IRAM(I1)=IRAM(IA-1+ICOL+1) 32 CONTINUE C C ICOL0=ICOL NA=IRAM(IA+M2)-IA NELEM=(NA-M2-1)/2 I1=NELMAT(M1,ICOL,ISYM) NSPAR=I1-NELEM ITMP=IRAM(IA+M2)-1-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 C C======================================================================= C