C
C Program MATFUN to compute diagonal matrix M2 as a function M2=FUN(M1) C of diagonal matrix M1. C C Version: 7.40 C Date: 2017, May 19 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 diagonal. C No default, MATIN1 must be specified and cannot be blank. 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 Function to be applied to the matrix M1: C MATFUN='string'... String specifying the function: C 'abs'... absolute value C 'aint'... part of the number before the decimal point C 'anint'... nearest integer (note that output file with C matrix M2 is always written using real numbers, even C for functions 'aint' and 'anint') C 'sqrt'... square root C 'exp'... exponential function C 'alog'... natural logarithm C 'alog10'... logarithm of 10 C 'sin'... sinus C 'cos'... cosinus C 'tan'... tangent C 'asin'... arcus sinus C 'acos'... arcus cosinus C 'atan'... arcus tangent C 'sinh'... sinus hyperbolicus C 'cosh'... cosinus hyperbolicus C 'tanh'... tangent hyperbolicus C 'inv'... inverse function (M2=1./M1) C No default, MATFUN must be specified and cannot be blank. 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 Value of undefined quantities: C UNDEF=real... The value to be used for undefined real quantities. C Default: UNDEF=undefined value used in forms.for C C For detailed description of storage of matrices in the memory C refer to file mat.for. C C----------------------------------------------------------------------- C Subroutines and external functions required: EXTERNAL ERROR,RSEP1,RSEP3I,RSEP3T,UARRAY,LOWER,UPPER, *RMATH,RMATD,WMATH,WMATD,GSMAT,GSMATR,SGMATR,NELMAT,ISYM INTEGER NELMAT,ISYM REAL UARRAY C ERROR ... File error.for. C RSEP1,RSEP3I,RSEP3T... File sep.for. C LOWER,UPPER... File length.for. C UARRAY... File forms.for. C RMATH,RMATD,WMATH,WMATD,GSMAT,GSMATR,SGMATR,NELMAT,ISYM... C File mat.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,TFUN,SYM1,FORM1, * FORM3,SPARS1 INTEGER M1,M2,NEL1,IA,NA,ISYM1,NSPAR1,ISPAR3,LU1,I1,I11,I12,I13, * IFUN REAL RSPAR1,UNDEF CHARACTER*72 TXTERR PARAMETER (LU1=1) C C----------------------------------------------------------------------- C C Reading a name of the file with the input data: WRITE(*,'(A)') '+MATFUN: 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 MATFUN-01 CALL ERROR('MATFUN-01: SEP file not given.') ENDIF C UNDEF=UARRAY() C C Reading the names of matrices header files: CALL RSEP3T('MATIN1',FILE1,' ') IF (FILE1.EQ.' ') THEN C MATFUN-02 CALL ERROR('MATFUN-02: File MATIN1 not given.') ENDIF CALL RSEP3T('MATOUT',FILE3,' ') FILED3=' ' IF (FILE3.EQ.' ') THEN C MATFUN-03 CALL ERROR('MATFUN-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 MATFUN-08 CALL ERROR('MATFUN-08: Wrong format of input matrix.') C In this version only dense or sparse CSC matrix is allowed. ENDIF IF (SYM1.NE.'diag') THEN C MATFUN-04 CALL ERROR('MATFUN-04: Input matrix is not diagonal.') ENDIF C Reading the properties of the output file: CALL RSEP3T('FORMM',FORM3,'FORMATTED') CALL UPPER(FORM3) CALL RSEP3I('SPARSE',ISPAR3,0) C C Reading the function: CALL RSEP3T('MATFUN',TFUN,' ') IF (TFUN.EQ.' ') THEN C MATFUN-05 CALL ERROR('MATFUN-05: Function MATFUN not given.') ENDIF CALL LOWER(TFUN) C C Registration of the function: C (numbers of functions correspond to 'grdcal.for') IFUN=-1 IF (TFUN.EQ.'abs') THEN IFUN= 6 ELSEIF (TFUN.EQ.'aint') THEN IFUN= 7 ELSEIF (TFUN.EQ.'anint') THEN IFUN= 8 ELSEIF (TFUN.EQ.'sqrt') THEN IFUN=14 ELSEIF (TFUN.EQ.'exp') THEN IFUN=15 ELSEIF (TFUN.EQ.'alog') THEN IFUN=16 ELSEIF (TFUN.EQ.'alog10') THEN IFUN=17 ELSEIF (TFUN.EQ.'sin') THEN IFUN=18 ELSEIF (TFUN.EQ.'cos') THEN IFUN=19 ELSEIF (TFUN.EQ.'tan') THEN IFUN=20 ELSEIF (TFUN.EQ.'asin') THEN IFUN=21 ELSEIF (TFUN.EQ.'acos') THEN IFUN=22 ELSEIF (TFUN.EQ.'atan') THEN IFUN=23 ELSEIF (TFUN.EQ.'sinh') THEN IFUN=25 ELSEIF (TFUN.EQ.'cosh') THEN IFUN=26 ELSEIF (TFUN.EQ.'tanh') THEN IFUN=27 ELSEIF (TFUN.EQ.'inv') THEN IFUN=28 ELSE GOTO 201 END IF C C Reading input matrix A: IA=1 IF (NSPAR1.GE.0) THEN NA=M2+1+2*NEL1 ELSE NA=NEL1 ENDIF IF (NA+1.GT.MRAM) THEN C MATFUN-07 WRITE(TXTERR,'(A,I9,A)') * 'MATFUN-07: Array RAM too small,',NA+1-MRAM,' units missing.' CALL ERROR(TXTERR) ENDIF CALL RMATD(LU1,FILED1,M2,SPARS1,NEL1,FORM1,IA) C C Applying the function: IF (NSPAR1.GE.0) THEN I11=IRAM(IA)+1 I12=IRAM(IA+M2)-1 I13=2 ELSE I11=1 I12=NEL1 I13=1 ENDIF DO 200, I1=I11,I12,I13 GO TO (101,102,103,104,105,106,107,108,109,110, * 111,112,113,114,115,116,117,118,119,120, * 121,122,123,124,125,126,127,128) IFUN GO TO 201 101 CONTINUE c RAM(I1)=RAM(I1)+RNAME(I2) GO TO 201 102 CONTINUE c RAM(I1)=RAM(I1)-RNAME(I2) GO TO 201 103 CONTINUE c RAM(I1)=RAM(I1)*RNAME(I2) GO TO 201 104 CONTINUE c IF(RNAME(I2).EQ.0.) THEN c IF(RAM(I1).EQ.0.) THEN c RAM(I1)=0. c ELSE c RAM(I1)=UNDEF c END IF c ELSE c RAM(I1)=RAM(I1)/RNAME(I2) c END IF GO TO 201 105 CONTINUE c IF(RAM(I1).LT.0.) THEN c RAM(I1)=UNDEF c ELSE c RAM(I1)=RAM(I1)**RNAME(I2) c END IF GO TO 201 106 CONTINUE RAM(I1)=ABS(RAM(I1)) GO TO 200 107 CONTINUE RAM(I1)=AINT(RAM(I1)) GO TO 200 108 CONTINUE RAM(I1)=ANINT(RAM(I1)) GO TO 200 109 CONTINUE c IF(RNAME(I2).EQ.0.) THEN c RAM(I1)=UNDEF c ELSE c RAM(I1)=AMOD(RAM(I1),RNAME(I2)) c END IF GO TO 201 110 CONTINUE c RAM(I1)=SIGN(RAM(I1),RNAME(I2)) GO TO 201 111 CONTINUE c RAM(I1)=DIM(RAM(I1),RNAME(I2)) GO TO 201 112 CONTINUE c RAM(I1)=AMAX1(RAM(I1),RNAME(I2)) GO TO 201 113 CONTINUE c RAM(I1)=AMIN1(RAM(I1),RNAME(I2)) GO TO 201 114 CONTINUE IF(RAM(I1).LT.0.) THEN RAM(I1)=UNDEF ELSE RAM(I1)=SQRT(RAM(I1)) END IF GO TO 200 115 CONTINUE RAM(I1)=EXP(RAM(I1)) GO TO 200 116 CONTINUE IF(RAM(I1).LE.0.) THEN RAM(I1)=UNDEF ELSE RAM(I1)=ALOG(RAM(I1)) END IF GO TO 200 117 CONTINUE IF(RAM(I1).LE.0.) THEN RAM(I1)=UNDEF ELSE RAM(I1)=ALOG10(RAM(I1)) END IF GO TO 200 118 CONTINUE RAM(I1)=SIN(RAM(I1)) GO TO 200 119 CONTINUE RAM(I1)=COS(RAM(I1)) GO TO 200 120 CONTINUE RAM(I1)=TAN(RAM(I1)) GO TO 200 121 CONTINUE IF(ABS(RAM(I1)).GT.1.) THEN RAM(I1)=UNDEF ELSE RAM(I1)=ASIN(RAM(I1)) END IF GO TO 200 122 CONTINUE IF(ABS(RAM(I1)).GT.1.) THEN RAM(I1)=UNDEF ELSE RAM(I1)=ACOS(RAM(I1)) END IF GO TO 200 123 CONTINUE IF(ABS(RAM(I1)).GT.1.) THEN RAM(I1)=UNDEF ELSE RAM(I1)=ATAN(RAM(I1)) END IF GO TO 200 124 CONTINUE c IF(RAM(I1).EQ.0..AND.RNAME(I2).EQ.0.) THEN c RAM(I1)=UNDEF c ELSE c RAM(I1)=ATAN2(RAM(I1),RNAME(I2)) c END IF GO TO 201 125 CONTINUE RAM(I1)=SINH(RAM(I1)) GO TO 200 126 CONTINUE RAM(I1)=COSH(RAM(I1)) GO TO 200 127 CONTINUE RAM(I1)=TANH(RAM(I1)) GO TO 200 128 CONTINUE IF(RAM(I1).EQ.0.) THEN RAM(I1)=UNDEF ELSE RAM(I1)=1./RAM(I1) END IF GO TO 200 200 CONTINUE GOTO 202 201 CONTINUE C MATFUN-06 CALL ERROR ('MATFUN-06: Unsupported function') C The function given by input parameter MATFUN is not coded. 202 CONTINUE C IF ((ISPAR3.EQ.1).AND.(NSPAR1.LT.0)) THEN C Changing non-sparse matrix to sparse matrix: CALL GSMATR(M1,M2,ISYM1,NSPAR1,NEL1,1,MRAM,IA,NA) ELSEIF ((ISPAR3.EQ.-1).AND.(NSPAR1.GE.0)) THEN C Changing sparse matrix to non-sparse matrix: CALL SGMATR(M1,M2,ISYM1,NSPAR1,NEL1,1,MRAM,IA,NA) ELSEIF (ISPAR3.EQ.0) THEN C Automatic selection of the sparseness: IF (NSPAR1.LT.0) THEN 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) ELSE C Matrix is sparse. RSPAR1=FLOAT(NSPAR1)/FLOAT(NSPAR1+NEL1) IF (RSPAR1.LT.0.5) THEN C Matrix is sparse, sparseness is less than 0.5, thus changing C the sparse matrix to non-sparse matrix: CALL SGMATR(M1,M2,ISYM1,NSPAR1,NEL1,1,MRAM,IA,NA) ENDIF ENDIF 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)') '+MATFUN: Done. ' C STOP 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