C
C Subroutine file 'mat.for' with subroutines dealing with matrices. C C Version: 6.20 C Date: 2008, June 12 C C....................................................................... C C This file consists of the following external procedures: C WMATH.. Subroutine designed to write a matrix header file. C WMATH C RMATH.. Subroutine designed to read a matrix header file. C RMATH C WMATD.. Subroutine designed to write a matrix data file. C WMATD C RMATD.. Subroutine designed to read a matrix data file. C RMATD C WMATR.. Auxiliary subroutine to WMATD, designed to write the array C of matrix elements into a file. C WMATR C TMATR...Subroutine to transpose a matrix. C TMATR C SMATRE..Subroutine to extend sparse symmetric or antisymmetric C matrix into sparse general matrix. C SMATRE C GSMATR..Subroutine to change normal (not sparse) matrix into C the sparse matrix. C GSMATR C SGMATR..Subroutine to change sparse matrix into the normal C (not sparse) matrix. C SGMATR C GSMAT...Subroutine to change normal (not sparse) matrix into C the sparse matrix if the matrix is sparser than given C limit. C GSMAT C VELEM...Real function to find the value of an element of a matrix. C NELMAT..Integer function to calculate total number of elements C of a non-sparse matrix. C NELMAT C NSPMAT..Integer function to calculate number of zero elements C of a non-sparse matrix. C NSPMAT C ISYM ...Integer function to assign the integer corresponding C to the symmetry. C ISYM C MSHIFT..Subroutine designed to shift the matrix in memory. C MSHIFT C Subroutines to enable functionality of some programs which C do not use above subroutines directly: C OMAT... Subroutine designed to set the form of the file for C reading or writing a matrix, and to open the file. C OMAT C WMAT... Subroutine designed to write a given matrix into the given C file. C WMAT C RMAT... Subroutine designed to read a matrix from the given file. C RMAT C C======================================================================= C C Storage of the matrices in the memory: C The matrices are stored in real array RAM of common block RAMC C defined in file ram.inc. C To store integers in this array, integer array IRAM coinciding C with real array RAM is also defined. C Normal (not sparse) matrices: C Normal (dense) matrices are identified by the value of C SPARSE=' ' in the matrix header file. For normal matrices C all values of the matrix are stored columnwise, starting C from the first value of the first column to the last C stored value of the first column, then the second column, C and so on to the value of the last stored element of the C last column. C If the matrix is other than general, only N elements C corresponding to the given symmetry are stored, i.e. for C a symmetric matrix, just elements from the first row to C the diagonal are stored for each column (N=M1*(M1+1)/2); C for a skew matrix, just elements above the diagonal C are stored for each column (N=M1*(M1-1)/2); for a diagonal C matrix, just diagonal elements are stored (N=M1). C Sparse matrices: C Sparse matrices are identified by the value of C SPARSE='CSC'. NELEM is the number of nonzero elements C of the matrix, calculating only the elements corresponding C to the given symmetry, see above. C The first M2 storage locations contain the addresses C of the beginnings of all columns of the matrix. C The nonzero elements of column J of the matrix are thus C stored in the positions from (I)RAM(IRAM(IMAT+J-1)) C to (I)RAM(IRAM(IMAT+J)-1), where IMAT is the address C of the first storage location of the matrix. C The (M2+1)th storage location contains the address in the C array just after the last stored element. This address C may be understood as the address of the beginning of a C fictitious column number M2+1. C Next 2*NELEM storage locations contain the nonzero C elements of the matrix stored columnwise. Each element C occupies two storage locations, the first one containing C index of the row of the element and written in the integer C array, the second one containing the value of the element C and written in the real array. In this form the matrix C occupies 1+M2+2*NELEM storage locations. C C For the description of the forms of the disk files with matrices C and of the parameters describing the matrices in the matrix header C files refer to file forms.htm. C C======================================================================= C C C SUBROUTINE WMATH(LU,FILEH,FILED,M1,M2,SPARSE,NELEM,SYMM,FORM) INTEGER LU,M1,M2,NELEM CHARACTER*(*) FILEH,FILED,SPARSE,SYMM,FORM C C Subroutine designed to write the matrix header file. C C Input: C LU... Logical unit number of the matrix header output file. C The output matrix header file is opened, written, and C closed. C FILEH.. String containing the name of the matrix header file. C The matrix header file has the form of the SEP parameter C file. C If FILEH=' ', no action is done. C FILED.. String containing the name of the matrix data file. C If FILED=' ', default filename is determined and returned. C The default matrix data filename is constructed from the C matrix header filename by either replacing its last C character by '@' if the filename has a three-character C extension, or by adding '@' in other cases. C M1... Number of rows of the matrix. C M2... Number of columns of the matrix. C SPARSE... Sparseness of the matrix data file: C SPARSE=' ': Matrix is dense (not sparse). C SPARSE='CSC' (case insensitive): Matrix is a sparse matrix C stored in the compressed column format. C NELEM...Number of nonzero elements of the sparse matrix, C calculating only the elements corresponding to the C symmetry of the matrix. C For SPARSE=' ', NELEM has no meaning and is not used. C SYMM... Symmetry of the matrix data file: C SYMM=' ': General matrix. C SYMM='sym': Symmetric matrix. C SYMM='skew': Skew matrix. C SYMM='diag': Diagonal matrix. C FORM... Form of the matrix data file: C FORM='FORMATTED': Formatted file. C FORM='UNFORMATTED': Unformatted file. C If FORM=' ', default FORM is determined according to the C SEP parameter FORMM and returned. Default for FORMM is C FORMM='FORMATTED'. C C Output: C FILED.. If FILED=' ' on input, default filename is determined C and returned. C Otherwise, the input value of FILED is returned. C FORM... If FORM=' ' on input, default FORM for writing C is determined, converted to uppercase, and returned. C C Input SEP parameter: C FORMM='string'... Form of the matrix data file: C FORMM='FORMATTED': Formatted file. C FORMM='UNFORMATTED': Unformatted file. C Default: FORMM='FORMATTED' C C Subroutines and external functions required: EXTERNAL LENGTH,UPPER,WSEP3I,WSEP3T,RSEP3T INTEGER LENGTH C LENGTH,UPPER ... File 'length.for'. C WSEP3I,WSEP3T,RSEP3T... File 'sep.for'. C C Coded by Ludek Klimes and Petr Bulant C C----------------------------------------------------------------------- C INTEGER L CHARACTER*72 TXTERR C IF(FILEH.EQ.' ') THEN RETURN END IF C C Opening the matrix header file: OPEN(LU,FILE=FILEH,ERR=100) C C Default data file name: IF(FILED.EQ.' ') THEN L=LENGTH(FILEH) IF (L.GE.4) THEN IF (FILEH(L-3:L-3).EQ.'.') THEN L=L-1 ENDIF ENDIF L=L+1 FILED=FILEH FILED(L:L)='@' END IF C C Default form of the data file: IF(FORM.EQ.' ') THEN CALL RSEP3T('FORMM',FORM,'FORMATTED') CALL UPPER(FORM) END IF C C Writing the values of the parameters describing the output matrix: CALL WSEP3T(LU,'IN',FILED) CALL WSEP3I(LU,'M1',M1) CALL WSEP3I(LU,'M2',M2) CALL WSEP3T(LU,'SPARSE',SPARSE) IF(SPARSE.NE.' ') THEN CALL WSEP3I(LU,'NELEM',NELEM) END IF CALL WSEP3T(LU,'SYMMETRY',SYMM) CALL WSEP3T(LU,'FORM',FORM) C CLOSE(LU) RETURN 100 CONTINUE C MAT-01 WRITE(TXTERR,'(A,A,A)') 'MAT-01: Error when opening file ''', *FILEH(1:MIN0(LENGTH(FILEH),37)),'''.' CALL ERROR(TXTERR) END C C======================================================================= C C C SUBROUTINE RMATH(LU,FILEH,FILED,M1,M2,SPARSE,NELEM,SYMM,FORM) INTEGER LU,M1,M2,NELEM CHARACTER*(*) FILEH,FILED,SPARSE,SYMM,FORM C C Subroutine designed to read the matrix header file. C C Input: C LU... Logical unit number of the input matrix header file. C The input matrix header file is opened, read, and C closed. C FILEH.. String containing the name of the matrix header file. C The matrix header file has the form of the SEP parameter C file. C If FILEH=' ', no file is read and defaults are returned. C The input parameters are not altered. C C Output: C FILED.. String containing the name of the matrix data file. C The default matrix data filename is constructed from the C matrix header filename by either replacing its last C character by '@' if the filename has a three-character C extension, or by adding '@' in other cases. C M1... Number of rows of the matrix. Default: M1=1 C M2... Number of columns of the matrix. Default: M2=1 C SPARSE... Sparseness of the matrix data file: C SPARSE=' ': Matrix is dense (not sparse). C SPARSE='CSC': Matrix is a sparse matrix stored in the C compressed column format. C Default: SPARSE=' ' C NELEM...Number of elements of the matrix stored in the matrix C data file. I.e. either the number of nonzero elements C of the sparse matrix, or the number of elements C of dense matrix. C SYMM... Symmetry of the matrix data file: C SYMM=' ': General matrix. C SYMM='sym': Symmetric matrix. C SYMM='skew': Skew matrix. C SYMM='diag': Diagonal matrix. C Default: SYMM=' ' C FORM... Form of the matrix data file: C FORM='FORMATTED': Formatted file. C FORM='UNFORMATTED': Unformatted file. C Default: FORM='FORMATTED' C Note: the values of SPARSE and FORM are converted to uppercase C on output, and the value of SYMM to lowercase, in order to C simplify comparing of the strings in the calling program. C C Subroutines and external functions required: EXTERNAL LENGTH,LOWER,UPPER,SSEP,RSEP1,RSEP3I,RSEP3T,ISYM,NELMAT INTEGER LENGTH,ISYM,NELMAT C LENGTH,LOWER,UPPER ... File 'length.for'. C SSEP,RSEP1,RSEP3I,RSEP3T... File 'sep.for'. C C Coded by Ludek Klimes and Petr Bulant C C----------------------------------------------------------------------- C CHARACTER*80 FILEDD INTEGER L,ISYMM C Indices of the sets of SEP parameters INTEGER ISEP,IOLD,ITMP SAVE ISEP DATA ISEP/0/ C C Switching to a new parameter set and reading the matrix header: CALL SSEP(ISEP,IOLD) CALL RSEP1(LU,FILEH) C C Reading the values of the parameters describing the input matrix: L=LENGTH(FILEH) IF (L.GE.4) THEN IF (FILEH(L-3:L-3).EQ.'.') THEN L=L-1 ENDIF ENDIF L=L+1 FILEDD=FILEH FILEDD(L:L)='@' CALL RSEP3T('IN',FILED,FILEDD) CALL RSEP3I('M1',M1,1) CALL RSEP3I('M2',M2,1) CALL RSEP3T('SPARSE',SPARSE,' ') CALL RSEP3T('SYMMETRY',SYMM,' ') CALL RSEP3T('FORM',FORM,'FORMATTED') C C Converting input strings to lowercase: CALL LOWER(SYMM) CALL UPPER(FORM) CALL UPPER(SPARSE) C C Number NELEM of stored matrix elements: ISYMM=ISYM(SYMM) NELEM=NELMAT(M1,M2,ISYMM) L=NELEM IF(SPARSE.NE.' ') THEN CALL RSEP3I('NELEM',NELEM,L) ENDIF C ISEP=-ISEP CALL SSEP(ISEP,ITMP) CALL SSEP(IOLD,ISEP) RETURN END C C======================================================================= C C C SUBROUTINE WMATD(LU,FILED,M1,M2,SPARSE,NELEM,FORM,IMAT) CHARACTER*(*) FILED,FORM,SPARSE INTEGER LU,M1,M2,NELEM,IMAT C C Subroutine designed to write a given matrix into the matrix data file. C The matrix to be written is assumed to be stored in common block RAMC C of file 'ram.inc'. C C Input: C LU... Logical unit number to be used for the output. C FILED...Destination filename. If not blank, the file will be C opened and closed. If blank, the file is assumed to be C already open, and will not be closed in this subroutine. C M1... Number of rows of the matrix. C Used only for SPARSE='CSC'. C M2... Number of columns of the matrix. C Used only for SPARSE='CSC'. C SPARSE... Sparseness of the matrix data file: C SPARSE=' ': Matrix is dense (not sparse). C SPARSE='CSC' (case insensitive): Matrix is a sparse matrix C stored in the compressed column format. C NELEM.. Number of elements of the matrix to be written. C FORM... Form of the matrix data file: C FORM='FORMATTED': Formatted file. C FORM='UNFORMATTED': Unformatted file. C IMAT... Address of the first storage location of the matrix C in array RAM (or IRAM) of common block RAMC. C No output. C C The description of storage of matrices in the memory may be found C above. C For description of storage of matrices in disk files refer to file C forms.htm. C C Subroutines and external functions required: EXTERNAL LENGTH,UPPER,WARRAI,WMATR INTEGER LENGTH C LENGTH,UPPER... File 'length.for'. C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc INTEGER IRAM(MRAM) EQUIVALENCE (RAM,IRAM) C C Coded by Ludek Klimes and Petr Bulant C C----------------------------------------------------------------------- C C Local storage locations: INTEGER O,I1,I2 CHARACTER*72 TXTERR CHARACTER*11 FORMU,SPARSU C O... Position of the first matrix element in (I)RAM. C C....................................................................... C C Form of the file with the matrix, opening the file: FORMU=FORM CALL UPPER(FORMU) IF(FILED.NE.' ') THEN WRITE(*,'(2A)') '+Writing: ',FILED(1:MIN0(LENGTH(FILED),70)) OPEN(LU,FILE=FILED,FORM=FORMU,ERR=100) END IF C C Writing the matrix: O=IMAT SPARSU=SPARSE CALL UPPER(SPARSU) IF(SPARSU.EQ.'CSC') THEN C Rewriting the array of pointers into the disk form: I2=(IMAT-1)+(M2+1)+1 DO 10, I1=IMAT,IMAT+M2 IRAM(I1)=(IRAM(I1)-I2)/2 10 CONTINUE C Writing the pointers: CALL WARRAI(LU,' ',FORMU,.FALSE.,0,.FALSE.,0,M2,IRAM(IMAT+1)) C Rewriting the array of pointers into the memory form: DO 20, I1=IMAT,IMAT+M2 IRAM(I1)=IRAM(I1)*2+I2 20 CONTINUE O=IMAT+M2+1 END IF CALL WMATR(LU,M1,M2,SPARSE,NELEM,FORMU,O) C C Closing output file: IF(FILED.NE.' ') THEN CLOSE(LU) WRITE(*,'(1A)') * '+ ' END IF RETURN 100 CONTINUE C MAT-02 WRITE(TXTERR,'(A,A,A)') 'MAT-02: Error when opening file ''', *FILED(1:MIN0(LENGTH(FILED),37)),'''.' CALL ERROR(TXTERR) END C C======================================================================= C C C SUBROUTINE RMATD(LU,FILED,M2,SPARSE,NELEM,FORM,IMAT) CHARACTER*(*) FILED,FORM,SPARSE INTEGER LU,M2,NELEM,IMAT C C Subroutine designed to read a matrix from the matrix data file C and to store it in common block RAMC of file 'ram.inc'. C C Input: C LU... Logical unit number to be used for the input. C FILED...Name of the matrix data file. If not blank, the file will C be opened and closed. If blank, the file is assumed to be C already open, and will not be closed in this subroutine. C M2... Number of columns of the matrix. C Used only for SPARSE='CSC'. C SPARSE... Sparseness of the matrix data file: C SPARSE=' ': Matrix is dense (not sparse). C SPARSE='CSC' (case insensitive): Matrix is a sparse matrix C stored in the compressed column format. C NELEM.. Number of elements of the matrix to be read. C FORM... Form of the matrix data file: C FORM='FORMATTED': Formatted file. C FORM='UNFORMATTED': Unformatted file. C IMAT... Address of the first storage location in array RAM C (or IRAM) of common block RAMC where to store the matrix. C C No output. C C Subroutines and external functions required: EXTERNAL LENGTH,UPPER,RARRAI,RARRAY INTEGER LENGTH C LENGTH,UPPER... File 'length.for'. C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc INTEGER IRAM(MRAM) EQUIVALENCE (RAM,IRAM) C C Coded by Ludek Klimes and Petr Bulant C C----------------------------------------------------------------------- C C Local storage locations: INTEGER O,I,I1,I2 CHARACTER*72 TXTERR CHARACTER*11 FORMU,SPARSU C O... Position of the first matrix element in (I)RAM. C C....................................................................... C C Opening the file with matrix data: FORMU=FORM CALL UPPER(FORMU) IF(FILED.NE.' ') THEN WRITE(*,'(2A)') '+Reading: ',FILED(1:MIN0(LENGTH(FILED),70)) OPEN(LU,FILE=FILED,FORM=FORMU,STATUS='OLD',ERR=100) END IF C C Reading the matrix: O=IMAT SPARSU=SPARSE CALL UPPER(SPARSU) IF(SPARSU.EQ.'CSC') THEN C Sparse matrix: C Reading the pointers: IRAM(O)=0 CALL RARRAI(LU,' ',FORMU,.FALSE.,0,M2,IRAM(O+1)) C Rewriting the array of pointers into the memory form: I2=(IMAT-1)+(M2+1)+1 DO 20, I1=IMAT,IMAT+M2 IRAM(I1)=IRAM(I1)*2+I2 20 CONTINUE C Reading the indices and values: O=IMAT+M2+1 IF(FORMU.EQ.'FORMATTED') THEN READ(LU,*) (IRAM(I),RAM(I+1),I=O,O+2*NELEM-1,2) ELSE READ(LU) (IRAM(I),RAM(I+1),I=O,O+2*NELEM-1,2) END IF ELSE CALL RARRAY(LU,' ',FORMU,.FALSE.,0.,NELEM,RAM(O)) END IF C C Closing input file: IF(FILED.NE.' ') THEN CLOSE(LU) WRITE(*,'(1A)') * '+ ' END IF RETURN 100 CONTINUE C MAT-03 WRITE(TXTERR,'(A,A,A)') 'MAT-03: Error when opening file ''', *FILED(1:MIN0(LENGTH(FILED),37)),'''.' CALL ERROR(TXTERR) END C C======================================================================= C C C SUBROUTINE WMATR(LU,M1,M2,SPARSE,NELEM,FORM,IMAT) CHARACTER*(*) SPARSE,FORM INTEGER LU,M1,M2,NELEM,IMAT C C Subroutine designed to write a given array of matrix elements into the C matrix data file. C C Input: C LU... Logical unit number to be used for the output. C The file should be already opened, and is not closed after C writing. C M1... Number of rows of the matrix. C Used only for SPARSE='CSC'. C M2... Number of columns of the matrix. C Used only for SPARSE='CSC'. C SPARSE... Sparseness of the matrix data file: C SPARSE=' ': Matrix is dense (not sparse). C SPARSE='CSC' (case insensitive): Matrix is a sparse matrix C stored in the compressed column format. C NELEM.. Number of matrix elements to be written. C FORM... Form of the matrix data file in lowercase: C FORM='FORMATTED': Formatted file. C FORM='UNFORMATTED': Unformatted file. C IMAT... Address of the first storage location in array RAM C (or IRAM) of common block RAMC where the matrix data C are stored (i.e. first location after pointers for C sparse matrix). C C No output. C C C Input SEP parameter: C NUMLINM=positive integer... Number of the numbers to be written C to each line of the output file. For sparse matrices C number of pairs of row indices and values of matrix C elements to be written to each line. C NUMLINM must be less than 100 (99 at most). C Default: NUMLINM=5 C C Subroutines and external functions required: EXTERNAL UPPER C UPPER... File 'length.for'. C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc INTEGER IRAM(MRAM) EQUIVALENCE (RAM,IRAM) C C Coded by Ludek Klimes and Petr Bulant C C----------------------------------------------------------------------- C C Local storage locations: INTEGER I CHARACTER*20 FORMAT CHARACTER*11 FORMU,SPARSU C FORMAT..String containing the output format. C INTEGER NUMLIN SAVE NUMLIN DATA NUMLIN/-1/ C C....................................................................... C C Setting output format: IF(NUMLIN.EQ.-1) THEN CALL RSEP3I('NUMLINM',NUMLIN,5) ENDIF FORMAT='(00(E13.7,1X))' FORMAT(3:3)=CHAR(ICHAR('0')+MOD(NUMLIN,10)) FORMAT(2:2)=CHAR(ICHAR('0')+ NUMLIN/10 ) C C Writing the array: FORMU=FORM CALL UPPER(FORMU) SPARSU=SPARSE CALL UPPER(SPARSU) IF(SPARSU.EQ.'CSC') THEN C Sparse matrix: IF(FORMU.EQ.'FORMATTED') THEN C FORMAT='(00(I0,1X,E13.7,1X))' FORMAT(5:20)='I0,1X,E13.7,1X))' IF (M1.GT.999999999) THEN C MAT-04 CALL ERROR * ('MAT-04: Insufficient space in format specification') ELSE IF (M1.GT.99999999) THEN FORMAT(6:6)='9' ELSE IF (M1.GT.9999999) THEN FORMAT(6:6)='8' ELSE IF (M1.GT.999999) THEN FORMAT(6:6)='7' ELSE IF (M1.GT.99999) THEN FORMAT(6:6)='6' ELSE IF (M1.GT.9999) THEN FORMAT(6:6)='5' ELSE IF (M1.GT.999) THEN FORMAT(6:6)='4' ELSE IF (M1.GT.99) THEN FORMAT(6:6)='3' ELSE IF (M1.GT.9) THEN FORMAT(6:6)='2' ELSE FORMAT(6:6)='1' END IF WRITE(LU,FORMAT) * (IRAM(I),RAM(I+1),I=IMAT,IMAT+2*NELEM-1,2) ELSE WRITE(LU) * (IRAM(I),RAM(I+1),I=IMAT,IMAT+2*NELEM-1,2) END IF ELSE C Dense (not sparse) matrix: IF(FORMU.EQ.'FORMATTED') THEN WRITE(LU,FORMAT) (RAM(I),I=IMAT,IMAT+NELEM-1) ELSE WRITE(LU) (RAM(I),I=IMAT,IMAT+NELEM-1) END IF END IF C RETURN END C======================================================================= C C C SUBROUTINE TMATR(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA) INTEGER M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA C C Subroutine designed to transpose 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 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 NSPAR.LT.0: Matrix is stored element by element. C NSPAR.GE.0: Matrix is stored as a sparse matrix, C NSPAR is the number of zero matrix elements. C NELEM...Number of elements of the matrix stored in array RAM. C MTMP .. Dimension of auxiliary arrays IATMP and ATMP. 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 For NSPAR.LT.0 and general matrix (symmetry=' '), C MAA-MIA+1 must be at least 2*NELEM. For NSPAR.LT.0 and C other symmetries, MAA-MIA+1 may equal NELEM. C For NSPAR.GE.0, MAA-MIA+1 must be at least NA+NAout. C IA ... Address of the first storage location in array (I)RAM C used for the matrix. C NA ... Number of storage locations for the input matrix. C C Output: C M1 ... Number of rows of the transposed matrix. C M2 ... Number of columns of the transposed matrix. C IA ... If possible, IAout equals IAin. For sparse matrix stored C at the end of (I)RAM and NAout.GT.NAin, IAout is set to C MAA-NAout+1. C NA ... Number of storage locations for the output matrix. C C Coded by Petr Bulant C C----------------------------------------------------------------------- C C Subroutines and external functions required: EXTERNAL ERROR 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,IAOLD,IANEW,IB,NB CHARACTER*72 TXTERR C C....................................................................... C IF (ISYM.EQ.1.OR.ISYM.EQ.2) THEN C Symmetric or diagonal matrix RETURN ENDIF IF (ISYM.EQ.3) THEN C Antisymmetric matrix IF (NSPAR.GE.0) THEN C Sparse matrix DO 1, I1=IA+M2+1+1,IA-1+NA,2 RAM(I1)=-RAM(I1) 1 CONTINUE ELSE C Normal (not sparse) matrix DO 10, I1=IA,IA-1+NA RAM(I1)=-RAM(I1) 10 CONTINUE ENDIF RETURN ENDIF C General matrix: IF (NSPAR.GE.0) THEN C Transposing sparse matrix: NB=NA-M2+M1 IF (NA+NB.GT.MAA-MIA+1) THEN C MAT-05 WRITE(TXTERR,'(A,I9,A)') * 'MAT-05: Array RAM too small,',NA+NB-(MAA-MIA+1), * ' units missing.' CALL ERROR(TXTERR) ENDIF IAOLD=IA IF (MAA-(IA+NA-1).LT.NB) THEN IANEW=MAA-(NA+NB)+1 CALL MSHIFT(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,IANEW) ENDIF IB=IA+NA C Computing lengths of columns of transp(A): DO 20, I1=0,M1 IRAM(IB+I1)=0 20 CONTINUE DO 24, I1=1,M2 DO 22, I2=IRAM(IA-1+I1),IRAM(IA+I1)-2,2 I3=IB+IRAM(I2) IRAM(I3)=IRAM(I3)+1 22 CONTINUE 24 CONTINUE C Computing pointers of transp(A) from lengths: IRAM(IB)=IB+M1+1 DO 26, I1=1,M1 IRAM(IB+I1)=IRAM(IB+I1-1)+2*IRAM(IB+I1) 26 CONTINUE C Creating transp(A): DO 30, I1=1,M2 DO 28, I2=IRAM(IA-1+I1),IRAM(IA+I1)-2,2 I3=IRAM(I2) I4=IRAM(IB-1+I3) IRAM(I4)=I1 RAM(I4+1)=RAM(I2+1) IRAM(IB-1+I3)=I4+2 28 CONTINUE 30 CONTINUE DO 32, I1=M1,1,-1 IRAM(IB+I1)=IRAM(IB-1+I1) 32 CONTINUE IRAM(IB)=IB+M1+1 C Shifting transposed matrix to the original position: IA=IB NA=NB I1=M1 M1=M2 M2=I1 IF (IAOLD-1+NA.LE.MAA) THEN CALL MSHIFT(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,IAOLD) ENDIF ELSE C Transposing normal (not sparse) matrix: IF (2*NELEM.GT.MAA-MIA+1) THEN C MAT-06 WRITE(TXTERR,'(A,I9,A)') * 'MAT-06: Array RAM too small,',2*NELEM-(MAA-MIA+1), * ' units missing.' CALL ERROR(TXTERR) ENDIF IAOLD=IA IF (MAA-(IA+NELEM-1).LT.NELEM) THEN IANEW=MAA-2*NELEM+1 CALL MSHIFT(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,IANEW) ENDIF I3=IA-1 I4=IA-1+NELEM DO 36, I2=1,M2 DO 34, I1=1,M1 RAM(I4+(I1-1)*M2+I2)=RAM(I3+(I2-1)*M1+I1) 34 CONTINUE 36 CONTINUE IA=I4+1 I1=M1 M1=M2 M2=I1 CALL MSHIFT(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,IAOLD) ENDIF C RETURN END C C======================================================================= C C C SUBROUTINE SMATRE(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA) INTEGER M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA C C Subroutine designed to extend sparse symmetric or antisymmetric C matrix to the sparse general matrix. C Sparse matrix to be extended must be stored in form "as on a disk", C i.e. matrix indices in IARRAY(1) to IARRAY(NELEM) and values of matrix C elements in ARRAY(NELEM+1) to ARRAY(NELEM+NELEM). C C Input: C M1... Number of rows of the matrix. C M2... Number of columns of the 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. Must be NSPAR.GE.0 on input. 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 At least NA+2*NNE storage locations must be available, C where NNE is the number of nonzero elements of the C extended matrix. C IA ... Address of the first storage location in array IRAM C used for the matrix. C NA ... Number of storage locations for the input matrix. C C Output: C ISYM ...ISYM=4, matrix is extended to general matrix. C NSPAR...Sparseness of the extended matrix (NSPAR.GE.0). C NELEM...Number of elements of the extended matrix. C IA ... Address of the first storage location in array IRAM C used for the extended matrix. C If IA-1+M2+1+2*NNE.LE.MAA, IA is not changed. C Otherwise IA is set to MAA-(M2+1+2*NNE)+1. C NA ... Number of storage locations for the extended matrix. C C Coded by Petr Bulant C C----------------------------------------------------------------------- C C Subroutines and external functions required: EXTERNAL ERROR,NELMAT,MSHIFT,VELEM INTEGER NELMAT REAL VELEM 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,NNE,IANEW,I2END REAL SYMMUL,VEL CHARACTER*72 TXTERR C C....................................................................... C IF ((NSPAR.LT.0).OR.((ISYM.NE.2).AND.(ISYM.NE.3))) THEN C MAT-07 CALL ERROR('MAT-07: Wrong invocation of SMATRE.') ENDIF C Calculating the new number of elements in the extended matrix: NNE=NELEM DO 3, I2=1,M2 DO 2, I3=IRAM(IA-1+I2),IRAM(IA+I2)-2,2 I1=IRAM(I3) IF (I1.NE.I2) NNE=NNE+1 2 CONTINUE 3 CONTINUE IF (M2+1+2*NNE.GT.MAA-MIA+1) THEN C MAT-08 WRITE(TXTERR,'(A,I9,A)') 'MAT-08: Array RAM too small,', * M2+1+2*NNE-(MAA-MIA+1),' units missing.' CALL ERROR(TXTERR) ENDIF C Multiplication switch for symmetric or antisymmetric matrix: IF (ISYM.EQ.2) THEN SYMMUL=1. ELSE SYMMUL=-1. ENDIF C Checking the available space for the extended matrix, C changing IA and shifting the matrix: IF (IA-1+M2+1+2*NNE.GT.MAA) THEN IANEW=MAA-(M2+1+2*NNE)+1 CALL MSHIFT(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,IANEW) ENDIF C Constructing the extended matrix: I4=IA+M2+1+2*NNE DO 14, I2=M2,1 C Address of the old end of the column I2: I2END=IRAM(IA-1+I2+1)-2 C Updating the address of the beginning of column I2+1: IRAM(IA-1+I2+1)=I4 C Recording the values under the diagonal: DO 12, I1=M1,I2+1,-1 VEL=SYMMUL*VELEM(M1,M2,ISYM,NSPAR,IA,I2,I1) IF (VEL.NE.0.) THEN I4=I4-2 IRAM(I4)=I1 RAM(I4+1)=VEL ENDIF 12 CONTINUE C Copying the values on and above the diagonal: DO 13, I3=I2END,IRAM(IA-1+I2),-2 I4=I4-2 IRAM(I4)=IRAM(I3) RAM(I4+1)=RAM(I3+1) 13 CONTINUE 14 CONTINUE C ISYM=4 NSPAR=NELMAT(M1,M2,ISYM)-NNE NELEM=NNE NA=M2+1+2*NNE RETURN END C C======================================================================= C C C SUBROUTINE GSMATR(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA) INTEGER M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA C C Subroutine designed to change normal (not sparse) matrix to the C sparse matrix. The sparse matrix is written to the end C of arrays I(RAM). C C Input: C M1... Number of rows of the matrix. C M2... Number of columns of the matrix. C ISYM... Index of symmetry of the matrix. C SYM='diag' ... ISYM=1 ... diagonal matrix C SYM='sym' ... ISYM=2 ... symmetric matrix C SYM='skew' ... ISYM=3 ... skew matrix C SYM=' ' ... ISYM=4 ... general matrix C NSPAR...Sparseness of the matrix. Must be NSPAR.LT.0 on input. 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 MAA-MIA should be at least M2+1+NAout, where NAout C stays for the value of NA on output, but higher value C of MAA-MIA is recommended. MAA-MIA equal to C NAin+M2+1+NAout is always sufficient. 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 C Output: C NSPAR...Sparseness of the matrix. NSPAR.GE.0 on output. C NELEM...Number of nonzero elements of the matrix in arrays (I)RAM. C IA ... New address of the first storage location in array RAM C used for the matrix. IA equals MAA-NA+1 on output. C NA ... New number of storage locations for the sparse matrix. C C Coded by Petr Bulant C C----------------------------------------------------------------------- C C Subroutines and external functions required: EXTERNAL ERROR,NSPMAT INTEGER NSPMAT 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,I11,I12,J1,J2,J3,NN,NSPARN,IANEW CHARACTER*72 TXTERR C C....................................................................... C IF ((NSPAR.GE.0).OR.(NELEM.NE.NA)) THEN C MAT-09 CALL ERROR('MAT-09: Wrong invocation of GSMATR.') ENDIF C C Number NSPARN of zero and NN of nonzero elements: NSPARN=NSPMAT(NA,RAM(IA)) NN=NA-NSPARN C C Moving the matrix to the optimum position: IF (IA.NE.MIA+M2+1) THEN IANEW=MAX0(MIA+M2+1,(MAA-(2*NN+NA)+1)) CALL MSHIFT(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,IANEW) ENDIF C IF ((IA+2*NN).GT.(MAA)) THEN C MAT-10 I1=IA+2*NN-MAA WRITE(TXTERR,'(A,I9,A)') * 'MAT-10: Array RAM too small,',I1,' units missing.' CALL ERROR(TXTERR) ENDIF C Moving the nonzero elements of the matrix to the end of RAM, C storing the corresponding row indices of the elements, C storing the column pointers: J1=MAA J2=IA-1+NA J3=IA-1-(M2+1) C J1 points where to write the next element of the sparse matrix C J2 I3 points to the next element of the dense matrix C J3 points before the array of pointers C Pointer for column M2+1: IRAM(J3+M2+1)=MAA+1 C Loop over indices of columns: DO 20, I2=M2,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 C Loop over indices of rows of the column I2: DO 18, I1=I11,I12,-1 IF (RAM(J2).NE.0.) THEN IF (J1-1.LT.J2) THEN C MAT-11 CALL ERROR('MAT-11: Array RAM too small') C An attemp to change normal matrix to sparse matrix C "in place" failed, more memory is required. ENDIF C Moving the value, recording the row index: RAM(J1)=RAM(J2) IRAM(J1-1)=I1 J1=J1-2 ENDIF J2=J2-1 18 CONTINUE C Writing the pointer for column I2: IRAM(J3+I2)=J1+1 20 CONTINUE IF (J2.NE.IA-1) THEN C MAT-12 CALL ERROR('MAT-12: Disorder in GSMATR') C This error should not appear. ENDIF C Updating IA: IA=J1+1-(M2+1) C Moving the array of pointers to the proper position: IF (IA.NE.J3+1) THEN DO 30, I1=M2+1,1,-1 IRAM(IA-1+I1)=IRAM(J3+I1) 30 CONTINUE ENDIF C Recording the numbers corresponding to the sparse matrix: NA=2*NN+M2+1 NSPAR=NSPARN NELEM=NN C RETURN END C C======================================================================= C C C SUBROUTINE SGMATR(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA) INTEGER M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA C C Subroutine designed to change sparse matrix to the normal (not sparse) C matrix. C C Input: C M1... Number of rows of the matrix. C M2... Number of columns of the 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 There must be always enough memory for columns 1 to K C of sparse matrix plus K to M2 of dense matrix. C (MAA-MIA+1).GE.(NAin+NAout) is always sufficient. C IA ... Address of the first storage location in array IRAM C used for the matrix. C NA ... Number of storage locations for the input matrix. C C Output: C NSPAR...Sparseness of the matrix. NSPAR=-1 on output. C NELEM...Number of elements of the matrix stored in array RAM. C IA ... New address of the first storage location in array RAM. C If IAin-1+NAout.LE.MAA, IAout equals IAin. Otherwise C IAout equals maximum of (MIA,MAA-NAout-NAin+1). C NA ... New number of storage locations for the matrix. C C Coded by Petr Bulant C C----------------------------------------------------------------------- C C Subroutines and external functions required: EXTERNAL ERROR,NELMAT,MSHIFT 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 I1,I2,I3,J1,J2,NNE,IANEW,IAOLD CHARACTER*72 TXTERR C C....................................................................... C IF (NSPAR.LT.0) RETURN C C Rewriting the matrix from the sparse form C to the normal (not sparse) form: C Number of elements of the dense matrix: NNE=NELMAT(M1,M2,ISYM) C Checking the memory, shifting the matrix if necessary: IAOLD=IA IF (IA-1+NA+NNE.GT.MAA) THEN IF (IA.GT.MIA) THEN IANEW=MAX0(MIA,MAA-NNE-NA+1) CALL MSHIFT(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,IANEW) ENDIF ENDIF C J2 points to the current position in the dense matrix: J2=MIN0(MAA,IA-1+NA+NNE) DO 20, I2=M2,1,-1 C J1 points before the beginning of column I2 in the dense matrix: IF (ISYM.EQ.1) THEN J1=J2-1 ELSEIF (ISYM.EQ.2) THEN J1=J2-I2 ELSEIF (ISYM.EQ.3) THEN J1=J2-I2+1 ELSE J1=J2-M1 ENDIF IF (J1.LT.IRAM(IA+I2)-1) THEN C MAT-13 WRITE(TXTERR,'(A,I9,A)') * 'MAT-13: Array RAM too small,',IRAM(IA+I2)-1-J1, * ' units missing.' CALL ERROR(TXTERR) ENDIF DO 8, I3=J1+1,J2 RAM(I3)=0. 8 CONTINUE C Rewriting column I2: DO 10, I3=IRAM(IA+I2)-1,IRAM(IA+I2-1)+1,-2 I1=IRAM(I3-1) IF (ISYM.EQ.1) THEN RAM(J1+1)=RAM(I3) ELSE RAM(J1+I1)=RAM(I3) ENDIF 10 CONTINUE J2=J1 20 CONTINUE IA=J2+1 NA=NNE NELEM=NA NSPAR=-1 C If possible, shifting the matrix to the original IA: IF (IA.NE.IAOLD) THEN IF (IAOLD-1+NA.LE.MAA) THEN CALL MSHIFT(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,IAOLD) ENDIF ENDIF C RETURN END C C======================================================================= C C C SUBROUTINE GSMAT(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,RSPAR) INTEGER M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA REAL RSPAR C C Subroutine designed to calculate the number NSPAR of zero elements C of a normal (not sparse) matrix, and to change the matrix into the C sparse matrix if NSPAR/NELEM.GE.RSPAR and the number of storage C locations required for the sparse matrix is lower than NA. C If the matrix is changed to sparse, it is written to the end C of arrays I(RAM). C C Input: C M1... Number of rows of the matrix. C M2... Number of columns of the 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. Must be NSPAR.LT.0 on input. 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 RSPAR...Minimum rate of sparseness to change the matrix C into the sparse matrix. C C Output: C For NSPAR.LT.0 on output, only the value of RSPAR is C calculated, and there are no other changes on output: C RSPAR...Rate of sparseness of the matrix. (Number of zero elements C of the matrix divided by the number of all elements.) C Otherwise: C NSPAR...Number of zero elements of the matrix. NSPAR.GE.0. C NELEM...Number of nonzero elements of the matrix. C IA ... New address of the first storage location in array IRAM. C NA ... New number of storage locations for the matrix. C RSPAR...Rate of sparseness of the matrix. (Number of zero elements C of the matrix divided by the number of all elements.) C C Coded by Petr Bulant C C----------------------------------------------------------------------- C C Subroutines and external functions required: EXTERNAL ERROR,NSPMAT,GSMATR,MSHIFT INTEGER NSPMAT 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 NN REAL RSP0 C C....................................................................... C IF ((NSPAR.GE.0).OR.(NELEM.NE.NA)) THEN C MAT-14 CALL ERROR('MAT-14: Wrong invocation of GSMAT.') ENDIF C Number of zero elements: NN=NSPMAT(NA,RAM(IA)) RSP0=RSPAR RSPAR=FLOAT(NN)/FLOAT(NA) IF ((RSPAR.LT.RSP0).OR.(M2+1+2*(NA-NN).GT.NA).OR. * (M2+1+NA.GT.MAA-MIA)) THEN C Matrix will stay non-sparse. RETURN ENDIF C Changing the matrix into the sparse matrix: CALL GSMATR(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA) C RETURN END C C======================================================================= C C C REAL FUNCTION VELEM(M1,M2,ISYM,NSPARS,IA,IROW,ICOL) INTEGER M1,M2,ISYM,NSPARS,IA,IROW,ICOL C C Subroutine designed to calculate Value of ELEMent 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 IROW .. Number of the row of the matrix element. C ICOL ...Number of the column of the matrix element. C For the detailed description of storage of matrices in the memory C refer to above. C C Output: C VELEM...Value of the matrix element. C C Coded by Petr Bulant 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, see C the detailed description of storage of matrices in C the memory given above. C C Local storage location: INTEGER IRO,ICO,II1,II2,II3,II21,IR1,IR2,IR3 REAL RMUL C....................................................................... C IF (NSPARS.GE.0) THEN C Sparse matrix, searching for VELEM by halving intervals: VELEM=0. C For diagonal of antisymmetric matrix: IF ((ISYM.EQ.3).AND.(IROW.EQ.ICOL)) RETURN RMUL=1. IRO=IROW ICO=ICOL IF (((ISYM.EQ.2).OR.(ISYM.EQ.3)).AND.(IROW.GT.ICOL)) THEN C (Anti)symmetric matrix: Aij=RMUL*Aji IRO=ICOL ICO=IROW IF (ISYM.EQ.3) RMUL=-1. ENDIF C Searching between beginning and end of column ICO. C II1 and II2 are the adresses of the two matrix elements between C which we search, IR1 and IR2 are their indices of rows. II1=IRAM(IA+ICO-1) II2=IRAM(IA+ICO)-2 IF (II2.GE.II1) THEN IR1=IRAM(II1) IR2=IRAM(II2) IF (IR1.EQ.IRO) THEN VELEM=RMUL*RAM(II1+1) RETURN ENDIF IF (IR2.EQ.IRO) THEN VELEM=RMUL*RAM(II2+1) RETURN ENDIF IF ((IR1.LT.IRO).AND.(IRO.LT.IR2)) THEN C IRO may be between IR1 and IR2, halving the interval: 10 CONTINUE II21=II2-II1 IF (II21.GT.2) THEN C II3 and IR3 is the half of the interval: II3=II1+(II21/4)*2 IR3=IRAM(II3) IF (IR3.EQ.IRO) THEN VELEM=RMUL*RAM(II3+1) RETURN ENDIF IF (IR3.LT.IRO) THEN II1=II3 IR1=IR3 GOTO 10 ELSEIF (IRO.LT.IR3) THEN II2=II3 IR2=IR3 GOTO 10 ENDIF ENDIF ENDIF ENDIF ELSE C Non-sparse matrix: IF (ISYM.EQ.2) THEN C 'sym' IF (IROW.LE.ICOL) THEN VELEM=RAM(IA-1+(ICOL-1)*ICOL/2+IROW) ELSE VELEM=RAM(IA-1+(IROW-1)*IROW/2+ICOL) ENDIF ELSEIF (ISYM.EQ.3) THEN C 'skew' IF (IROW.LT.ICOL) THEN VELEM=RAM(IA-1+(ICOL-1)*(ICOL-2)/2+IROW) ELSEIF (IROW.EQ.ICOL) THEN VELEM=0. ELSE VELEM=-RAM(IA-1+(IROW-1)*(IROW-2)/2+ICOL) ENDIF ELSEIF (ISYM.EQ.1) THEN C 'diag' IF (IROW.EQ.ICOL) THEN VELEM=RAM(IA-1+ICOL) ELSE VELEM=0. ENDIF ELSE C ' ' VELEM=RAM(IA-1+(ICOL-1)*M1+IROW) ENDIF ENDIF RETURN END C C======================================================================= C C C INTEGER FUNCTION NELMAT(M1,M2,ISYM) INTEGER M1,M2,ISYM C C Integer function to calculate number of elements of non-sparse 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 C Output: C NELMAT. Number of elements of M1*M2 non-sparse matrix C of given symmetry. Note that for indices corresponding to C symmetries 'diag', 'sym' and 'skew' number M2 is used C to calculate NELMAT. C C Coded by Petr Bulant C C----------------------------------------------------------------------- C IF (ISYM.EQ.4) THEN C General matrix NELMAT=M1*M2 ELSEIF (ISYM.EQ.2) THEN C Symmetric matrix NELMAT=M2*(M2+1)/2 ELSEIF (ISYM.EQ.3) THEN C Skew matrix NELMAT=M2*(M2-1)/2 ELSEIF (ISYM.EQ.1) THEN C Diagonal matrix NELMAT=M2 ELSE C MAT-15 CALL ERROR('MAT-15: Wrong index of matrix symmetry.') C Input argument ISYM should be equal to 1, 2, 3 or 4. ENDIF C RETURN END C C======================================================================= C C C INTEGER FUNCTION NSPMAT(NELEM,ARRAY) INTEGER NELEM REAL ARRAY(*) C C Integer function to calculate number of zero elements of non-sparse c matrix. C C Input: C NELEM...Number of elements of the matrix. C ARRAY...Elements of the matrix. C C Output: C NSPMAT..Number of zero elements of the matrix. C C Coded by Petr Bulant C INTEGER I C----------------------------------------------------------------------- C NSPMAT=0 DO 10, I=1,NELEM IF (ARRAY(I).EQ.0.) NSPMAT=NSPMAT+1 10 CONTINUE C RETURN END C C C======================================================================= C C C INTEGER FUNCTION ISYM(SYM) CHARACTER*(*) SYM C C Integer function to assign the integer corresponding to the symmetry. C C Input: C SYM ... Symmetry of a matrix. C C Output: C ISYM... Integer number corresponding to the symmetry SYM. C SYM='diag' ... ISYM=1 C SYM='sym' ... ISYM=2 C SYM='skew' ... ISYM=3 C SYM=' ' ... ISYM=4 C C Coded by Petr Bulant C C----------------------------------------------------------------------- C IF (SYM.EQ.'diag') THEN ISYM=1 ELSEIF (SYM.EQ.'sym' ) THEN ISYM=2 ELSEIF (SYM.EQ.'skew') THEN ISYM=3 ELSEIF (SYM.EQ.' ' ) THEN ISYM=4 ELSE C MAT-16 CALL ERROR('MAT-16: Wrong matrix symmetry.') C Input argument SYM should be equal C to one of strings ' ', 'sym', 'skew', 'diag'. ENDIF C RETURN END C C======================================================================= C C C SUBROUTINE MSHIFT(M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,IANEW) INTEGER M1,M2,ISYM,NSPAR,NELEM,MIA,MAA,IA,NA,IANEW C C Subroutine designed to shift the matrix in arrays (I)RAM C to the new position IANEW. C C Input: C M1... Number of rows of the matrix. C M2... Number of columns of the 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 arrays (I)RAM. C MIA,MAA ... Minimum and maximum address of arrays (I)RAM C available for the subroutine. C Only the storage locations between IA and IA+NA-1, c and between IANEW and IANEW+NA-1 are altered. C IA ... Address of the first storage location in arrays (I)RAM C used for the matrix. C NA ... Number of storage locations for the matrix. C IANEW.. Address of the first storage location in arrays (I)RAM C where the matrix is to be shifted. C C Output: C IA ... IA is changed to the value of IANEW. 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, see C the detailed description of storage of matrices in C the memory given above. C C Local storage location: INTEGER II,ISHIFT C C....................................................................... C IF (IANEW.EQ.IA) RETURN C IF ((IANEW.LT.MIA).OR.(IANEW+NA-1.GT.MAA).OR. * (IA.LT.MIA).OR.(IA+NA-1.GT.MAA)) THEN C MAT-17 CALL ERROR('MAT-17: Wrong invocation of MSHIFT.') C Storage locations IA to IA+NA-1 and C IANEW to IANEW+NA-1 must fit into the available memory. ENDIF C IF (NSPAR.LT.0) THEN C Non-sparse matrix: IF (IANEW.LT.IA) THEN DO 10, II=0,NA-1 RAM(IANEW+II)=RAM(IA+II) 10 CONTINUE ELSE DO 11, II=NA-1,0,-1 RAM(IANEW+II)=RAM(IA+II) 11 CONTINUE ENDIF ELSE C Sparse matrix, form CSC: ISHIFT=IANEW-IA IF (IANEW.LT.IA) THEN DO 30, II=0,M2 IRAM(IANEW+II)=IRAM(IA+II)+ISHIFT 30 CONTINUE DO 31, II=0,NELEM-1 IRAM(IANEW+M2+1+II*2) =IRAM(IA+M2+1+II*2) RAM(IANEW+M2+1+II*2+1)= RAM(IA+M2+1+II*2+1) 31 CONTINUE ELSE DO 32, II=NELEM-1,0,-1 RAM(IANEW+M2+1+II*2+1)= RAM(IA+M2+1+II*2+1) IRAM(IANEW+M2+1+II*2) =IRAM(IA+M2+1+II*2) 32 CONTINUE DO 33, II=M2,0,-1 IRAM(IANEW+II)=IRAM(IA+II)+ISHIFT 33 CONTINUE ENDIF ENDIF C IA=IANEW RETURN END C C======================================================================= C C C SUBROUTINE OMAT(LU,FILE,IRW,FORMM) CHARACTER*(*) FILE,FORMM INTEGER LU,IRW C C Subroutine for backward compatibility. C C Coded by Petr Bulant C C----------------------------------------------------------------------- RETURN END C C======================================================================= C C C SUBROUTINE WMAT(LU,FILE,M1,M2,OUT) CHARACTER*(*) FILE INTEGER LU,M1,M2 REAL OUT(*) C C Subroutine for backward compatibility. WMATH and WMATD should be used C instead of using WMAT !!! C Subroutine designed to write a given matrix into the file. C C Input: C LU... Logical unit number to be used for the output. C FILE... Destination header file name. Must not be blank. C M1... Number of rows of the given matrix. C M2... M2=0 for a symmetric matrix, C M2=-1 for a diagonal matrix, C M2=number of columns for a general matrix. C OUT... Components of the given matrix stored columnwise. C For a symmetric matrix, just components from the first row C to the diagonal are stored for each column, i.e., array C OUT has M1*(M1+1)/2 matrix components. C For a diagonal matrix, just M1 diagonal components are C stored. C C No output. C C Coded by Ludek Klimes and Petr Bulant C C----------------------------------------------------------------------- C C Local storage locations: C EXTERNAL ISYM,NELMAT,WMATD,WMATH,ERROR INTEGER ISYM,NELMAT CHARACTER*13 FORMM INTEGER I1,I2,NELEM C C FORMM ..Form of the files with matrices. C C....................................................................... C CHARACTER*80 FILED CHARACTER*13 FORMAT CHARACTER*4 SYMM C IF (FILE.EQ.' ') THEN C MAT-18 CALL ERROR('MAT-18: Matrix header file not given') C In this version of WMAT, the matrix header file C must be specified and sequential writing of matrix data file C is not allowed. ENDIF FORMM=' ' FILED=' ' IF(M2.EQ.0) THEN C Symmetric matrix I2=M1 SYMM='sym' ELSEIF(M2.EQ.-1) THEN C Diagonal matrix I2=M1 SYMM='diag' ELSE C General matrix I2=M2 SYMM=' ' END IF NELEM=NELMAT(M1,I2,ISYM(SYMM)) CALL WMATH(LU,FILE,FILED,M1,I2,' ',NELEM,SYMM,FORMM) C FORMAT='(5(G14.7,1X))' C Writing the matrix: OPEN(LU,FILE=FILED,FORM=FORMM) IF(M2.EQ.0) THEN C Symmetric matrix IF (FORMM.EQ.'FORMATTED') THEN DO 11 I2=1,M1 WRITE(LU,FORMAT) (OUT(I1),I1=I2*(I2-1)/2+1,I2*(I2+1)/2) 11 CONTINUE ELSE WRITE(LU) (OUT(I1),I1=1,M1*(M1+1)/2) ENDIF ELSE C Diagonal or general matrix IF (FORMM.EQ.'FORMATTED') THEN DO 12 I2=M1,M1*IABS(M2),M1 WRITE(LU,FORMAT) (OUT(I1),I1=I2-M1+1,I2) 12 CONTINUE ELSE WRITE(LU) (OUT(I1),I1=1,M1*IABS(M2)) ENDIF END IF CLOSE(LU) RETURN END C C======================================================================= C C C SUBROUTINE RMAT(LU,FILE,M1,M2,ARRAY) CHARACTER*(*) FILE INTEGER LU,M1,M2 REAL ARRAY(*) C C Subroutine for backward compatibility. RMATH and RMATD should be used C instead of using RMAT !!! C Subroutine designed to read a matrix from the file. C In this version of RMAT sequential reading of matrix data file C is not allowed. C C Input: C LU... Logical unit number to be used for the input. C FILE... Destination header file name. Must not be blank. C M1... Number of rows of the matrix. C M2... M2=0 for a symmetric matrix, C M2=1 for a diagonal matrix, C M2=number of columns for a general matrix. C C Output: C ARRAY...Components of the given matrix stored columnwise. C For a symmetric matrix, just components from the first row C to the diagonal are stored for each column, i.e., array C ARRAY has M1*(M1+1)/2 matrix components. C For a diagonal matrix, just M1 diagonal components are C stored. C C Coded by Ludek Klimes and Petr Bulant C C----------------------------------------------------------------------- C EXTERNAL RMATD,RMATH,ERROR C Local storage location: CHARACTER*13 FORMM,SPARSE CHARACTER*80 FILED CHARACTER*4 SYMM INTEGER I,M1R,M2R,NELEM C IF (FILE.EQ.' ') THEN C MAT-19 CALL ERROR('MAT-19: Matrix header file not given') C In this version of RMAT, the matrix header file C must be specified and sequential reading of matrix data file C is not allowed. ENDIF CALL RMATH(LU,FILE,FILED,M1R,M2R,SPARSE,NELEM,SYMM,FORMM) IF ((M1.NE.M1R).OR.(SPARSE.NE.' ').OR.(SYMM.EQ.'skew').OR. * ((SYMM.EQ.'sym' ).AND.(M2.NE.0)).OR. * ((SYMM.EQ.' ' ).AND.(M2.NE.M2R)).OR. * ((SYMM.EQ.'diag').AND.(M2.NE.1))) THEN C MAT-20 CALL ERROR('MAT-20: Unexpected values in matrix header file.') C Some of the values of the matrix header file differs from C the values estimated using input parameters M1 and M2. C For example, this version of RMAT cannot deal with sparse C matrices. ENDIF C OPEN(LU,FILE=FILED,FORM=FORMM) C Reading the matrix: IF(M2.LE.0) THEN C Symmetric matrix IF (FORMM.EQ.'FORMATTED') THEN READ(LU,*) (ARRAY(I),I=1,M1*(M1+1)/2) ELSE READ(LU) (ARRAY(I),I=1,M1*(M1+1)/2) ENDIF ELSE C Diagonal or general matrix IF (FORMM.EQ.'FORMATTED') THEN READ(LU,*) (ARRAY(I),I=1,M1*M2) ELSE READ(LU) (ARRAY(I),I=1,M1*M2) ENDIF END IF CLOSE(LU) C RETURN END C C======================================================================= C