C
C Program MATINV to compute symmetric matrix M2, which is inverse to C input symmetric matrix M1. C C Version: 6.20 C Date: 2008, June 12 C C Coded by Petr Bulant C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C C....................................................................... C C Description of data files: C C Input data read from the standard input device (*): C The data are read by the list directed input (free format) and C consist of a single string 'SEP': C 'SEP'...String in apostrophes containing the name of the input C SEP parameter or history file with the input data. C No default, 'SEP' must be specified and cannot be blank. C C C Input data file 'SEP': C File 'SEP' has the form of the SEP C parameter file. The parameters, which do not differ from their C defaults, need not be specified in file 'SEP'. C Filenames of the header files of the matrices: C MATIN1='string' .. Name of the header file of the input matrix M1. C Matrix M1 must be symmetric. C No default, MATIN1 must be specified and cannot be blank. C MATOUT='string' . Name of the header file of the output matrix M2. 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 M2: 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 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 SMBLO,IND,ERROR,RSEP1,RSEP3I,RSEP3T,UPPER, * RMATH,RMATD,WMATH,WMATD,GSMAT,GSMATR,SGMATR,NELMAT,ISYM,SINV INTEGER IND,NELMAT,ISYM C SMBLO,IND ... This file. C ERROR ... File error.for. C RSEP1,RSEP3I,RSEP3T ... C File sep.for. C UPPER ... File C length.for. C RMATH,RMATD,WMATH,WMATD,GSMAT,GSMATR,SGMATR,NELMAT,ISYM ... C File mat.for. C SINV ... File sinv.for. C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C INTEGER IRAM(MRAM) EQUIVALENCE (IRAM,RAM) C CHARACTER*80 FILSEP,FILE1,FILE3,FILED1,FILED3,SYM1,FORM1, * FORM3,SPARS1 INTEGER M1,M2,NEL1,IA,NA,NSPAR1,ISPAR3,LU1,I1,I2,I3,I4, * N,NN,IANEW,NB,IBMI,IBMA,IER,ISYM1 REAL RSPAR1,EPS CHARACTER*72 TXTERR PARAMETER (LU1=1) C C----------------------------------------------------------------------- C EPS=0.000001 C C Reading a name of the file with the input data: WRITE(*,'(A)') '+MATINV: 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 MATINV-01 CALL ERROR('MATINV-01: SEP file not given.') ENDIF C C Reading the names of matrices header files: CALL RSEP3T('MATIN1',FILE1,' ') IF (FILE1.EQ.' ') THEN C MATINV-02 CALL ERROR('MATINV-02: File MATIN1 not given.') ENDIF CALL RSEP3T('MATOUT',FILE3,' ') FILED3=' ' IF (FILE3.EQ.' ') THEN C MATINV-03 CALL ERROR('MATINV-03: File MATOUT not given.') ENDIF C Reading the header file of the input matrix: 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 MATINV-06 CALL ERROR('MATINV-06: Wrong sparseness of the input matrix.') C In this version only dense or sparse CSC matrix is allowed. ENDIF IF (SYM1.NE.'sym') THEN C MATINV-07 CALL ERROR('MATINV-07: Input matrix is not symmetric.') ENDIF C Reading the properties of the output file: CALL RSEP3T('FORMM',FORM3,'FORMATTED') CALL UPPER(FORM3) CALL RSEP3I('SPARSE',ISPAR3,0) C N=M1 NN=M1*(M1+1)/2 IF (3*NN+N+1.GT.MRAM) THEN C MATINV-04 WRITE(TXTERR,'(A,I9,A)') * 'MATINV-04: Array RAM too small,',3*NN+N+1-MRAM,' units missing' CALL ERROR(TXTERR) END IF C C Reading input matrix: IF (NSPAR1.GE.0) THEN IA=1 NA=2*NEL1 ELSE IA=MRAM-NN+1 NA=NEL1 ENDIF CALL RMATD(LU1,FILED1,M2,SPARS1,NEL1,FORM1,IA) IF (NSPAR1.GE.0) THEN C Changing input matrix to non-sparse: CALL SGMATR(M1,M2,ISYM1,NSPAR1,NEL1,1,MRAM,IA,NA) IANEW=MRAM-NN+1 CALL MSHIFT(M1,M2,ISYM1,NSPAR1,NEL1,1,MRAM,IA,NA,IANEW) ENDIF C WRITE(*,'(A)') '+MATINV: Working... ' C C Identification of blocks: CALL SMBLO(RAM(MRAM-NN+1),N,NN,NB,IRAM(MRAM-NN-N)) C C Preparing array for results: DO 5, I1=MRAM-2*NN-N,MRAM-NN-N-1 RAM(I1)=0. 5 CONTINUE C C Computing the inverse matrix: DO 20, I1=1,NB C Current block: IBMI=IRAM(MRAM-NN-N+I1-1)+1 IBMA=IRAM(MRAM-NN-N+I1) C Moving the block to RAM(1): I4=0 DO 10, I2=IBMI,IBMA DO 8, I3=IBMI,I2 I4=I4+1 RAM(I4)=RAM(MRAM-NN+IND(I3,I2)) 8 CONTINUE 10 CONTINUE I2=IBMA-IBMI+1 I3=I2*(I2+1)/2 C Inverting matrix: CALL SINV(RAM,I2,EPS,IER) IF(IER.NE.0) THEN C MATINV-05 WRITE(TXTERR,'(A,I5,A)') * 'MATINV-05: Error',IER,' in subroutine SINV' CALL ERROR(TXTERR) ENDIF C Moving the block to RAM(MRAM-NN-N-1-NN+1): I4=0 DO 14, I2=IBMI,IBMA DO 12, I3=IBMI,I2 I4=I4+1 RAM(MRAM-2*NN-N-1+IND(I3,I2))=RAM(I4) 12 CONTINUE 14 CONTINUE 20 CONTINUE C C IA=MRAM-2*NN-N IF (ISPAR3.EQ.1) THEN C Changing non-sparse matrix to sparse matrix: CALL GSMATR(M1,M2,ISYM1,NSPAR1,NEL1,1,MRAM,IA,NA) ELSEIF (ISPAR3.EQ.0) THEN C Automatic selection of the sparseness: C Matrix is non-sparse, it will be changed to sparse if its C sparseness is 0.5 or more: RSPAR1=0.5 CALL GSMAT(M1,M2,ISYM1,NSPAR1,NEL1,1,MRAM,IA,NA,RSPAR1) ENDIF C C Writing output matrix: SPARS1=' ' IF (NSPAR1.GE.0) SPARS1='CSC' CALL WMATH(LU1,FILE3,FILED3,M1,M2,SPARS1,NEL1,SYM1,FORM3) CALL WMATD(LU1,FILED3,M1,M2,SPARS1,NEL1,FORM3,IA) C WRITE(*,'(A)') '+MATINV: Done. ' C STOP END C C======================================================================= C SUBROUTINE SMBLO(SM,N,NN,NB,IB) C C======================================================================= C Subroutine to find blocks of a symmetric matrix INTEGER N,NN,NB,IB(0:N) REAL SM(NN) C Input: SM ... symmetric matrix C N ... number of rows (and columns) of the matrix C NN ... N*(N+1)/2 C Output: NB ... number of blocks of the matrix C IB ... array with numbers of lines (and rows), at which C the blocks finish INTEGER I1,I2 EXTERNAL IND INTEGER IND C----------------------------------------------------------------------- C Initiating first block: NB=1 IB(0)=0 IB(1)=1 I1=0 10 CONTINUE C Loop for lines of the matrix: I1=I1+1 DO 20, I2=N,IB(NB)+1,-1 C Loop for rows of the line I1: IF (SM(IND(I1,I2)).NE.0.) THEN C The block is larger: IB(NB)=I2 GOTO 21 ENDIF 20 CONTINUE 21 CONTINUE IF (IB(NB).EQ.N) THEN C End of the search for blocks. RETURN ENDIF IF (IB(NB).EQ.I1) THEN C Block finished, continuing with the next block. NB=NB+1 IB(NB)=I1+1 ENDIF GOTO 10 END C======================================================================= INTEGER FUNCTION IND(I,J) INTEGER I,J C I ... Index of a line. C J ... Index of a row. IND=J*(J-1)/2+I 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 'sinv.for' C sinv.for INCLUDE 'mfsd.for' C mfsd.for C C======================================================================= C