C
C Program MATFUN to compute diagonal matrix M2 as a function M2=FUN(M1)
C of diagonal matrix M1.
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             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
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
      UNDEF=UARRAY()
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
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