C
C Program MATFUN to compute diagonal matrix M2 as a function M2=FUN(M1)
C of diagonal matrix M1.
C
C Version: 6.10
C Date: 2007, June 8
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 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 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     MATFORM='string' ... Form of the output file with the matrix.
C             Possible values of MATFORM are:
C               MATFORM='formatted'   ... formatted file
C               MATFORM='unformatted' ... unformatted file
C             Default: MATFORM='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               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     NUMLIN=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,RMATH,RMATD,WMATH,WMATD,
     *GSMAT,GSMATR,SGMATR
      REAL UARRAY
C     ERROR ... File error.for.
C     RSEP1,RSEP3I,RSEP3T ...
C             File sep.for.
C     UARRAY ... File forms.for.
C     RMATH,RMATD,WMATH,WMATD,GSMAT,GSMATR,SGMATR ...
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
      INTEGER M1,M2,NEL1,OA,NA,NSPAR1,ISPAR3,ISTORA,LU1,I1,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,NSPAR1,SYM1,FORM1,NEL1)
      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('MATFORM',FORM3,'formatted')
      CALL LOWER(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')
      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
C       MATFUN-06
        CALL ERROR ('MATFUN-06: Unknown function')
      END IF
C
C     Reading input matrix A:
      OA=1
      IF (NSPAR1.GE.0) THEN
        NA=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,NSPAR1,FORM1,NEL1,IRAM(OA),RAM(OA))
      ISTORA=0
      IF (NSPAR1.GE.0) THEN
C       Changing A to non-sparse:
        CALL SGMATR(M1,M2,SYM1,NSPAR1,NEL1,RAM,IRAM,MRAM,
     *              OA,NA,ISTORA)
      ENDIF
C
C     Applying the function:
      DO 200, I1=OA,OA+NA-1
        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 199
  101   CONTINUE
c         RAM(I1)=RAM(I1)+RNAME(I2)
          GO TO 199
  102   CONTINUE
c         RAM(I1)=RAM(I1)-RNAME(I2)
          GO TO 199
  103   CONTINUE
c         RAM(I1)=RAM(I1)*RNAME(I2)
          GO TO 199
  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 199
  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 199
  106   CONTINUE
          RAM(I1)=ABS(RAM(I1))
          GO TO 199
  107   CONTINUE
          RAM(I1)=AINT(RAM(I1))
          GO TO 199
  108   CONTINUE
          RAM(I1)=ANINT(RAM(I1))
          GO TO 199
  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 199
  110   CONTINUE
c         RAM(I1)=SIGN(RAM(I1),RNAME(I2))
          GO TO 199
  111   CONTINUE
c         RAM(I1)=DIM(RAM(I1),RNAME(I2))
          GO TO 199
  112   CONTINUE
c         RAM(I1)=AMAX1(RAM(I1),RNAME(I2))
          GO TO 199
  113   CONTINUE
c         RAM(I1)=AMIN1(RAM(I1),RNAME(I2))
          GO TO 199
  114   CONTINUE
          IF(RAM(I1).LT.0.) THEN
            RAM(I1)=UNDEF
          ELSE
            RAM(I1)=SQRT(RAM(I1))
          END IF
          GO TO 199
  115   CONTINUE
          RAM(I1)=EXP(RAM(I1))
          GO TO 199
  116   CONTINUE
          IF(RAM(I1).LE.0.) THEN
            RAM(I1)=UNDEF
          ELSE
            RAM(I1)=ALOG(RAM(I1))
          END IF
          GO TO 199
  117   CONTINUE
          IF(RAM(I1).LE.0.) THEN
            RAM(I1)=UNDEF
          ELSE
            RAM(I1)=ALOG10(RAM(I1))
          END IF
          GO TO 199
  118   CONTINUE
          RAM(I1)=SIN(RAM(I1))
          GO TO 199
  119   CONTINUE
          RAM(I1)=COS(RAM(I1))
          GO TO 199
  120   CONTINUE
          RAM(I1)=TAN(RAM(I1))
          GO TO 199
  121   CONTINUE
          IF(ABS(RAM(I1)).GT.1.) THEN
            RAM(I1)=UNDEF
          ELSE
            RAM(I1)=ASIN(RAM(I1))
          END IF
          GO TO 199
  122   CONTINUE
          IF(ABS(RAM(I1)).GT.1.) THEN
            RAM(I1)=UNDEF
          ELSE
            RAM(I1)=ACOS(RAM(I1))
          END IF
          GO TO 199
  123   CONTINUE
          IF(ABS(RAM(I1)).GT.1.) THEN
            RAM(I1)=UNDEF
          ELSE
            RAM(I1)=ATAN(RAM(I1))
          END IF
          GO TO 199
  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 199
  125   CONTINUE
          RAM(I1)=SINH(RAM(I1))
          GO TO 199
  126   CONTINUE
          RAM(I1)=COSH(RAM(I1))
          GO TO 199
  127   CONTINUE
          RAM(I1)=TANH(RAM(I1))
          GO TO 199
  128   CONTINUE
          IF(RAM(I1).EQ.0.) THEN
            RAM(I1)=UNDEF
          ELSE
            RAM(I1)=1./RAM(I1)
          END IF
          GO TO 199
  199   CONTINUE
  200 CONTINUE
C
      IF (ISPAR3.EQ.1) THEN
C       Changing non-sparse matrix to sparse matrix:
        CALL GSMATR(M1,M2,SYM1,NSPAR1,NEL1,RAM,IRAM,MRAM,OA,NA,ISTORA)
      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,SYM1,NSPAR1,NEL1,RAM,IRAM,MRAM,
     *                  OA,NA,ISTORA,RSPAR1)
      ENDIF
C
C     Writing output matrix:
      CALL WMATH(LU1,FILE3,FILED3,M1,M2,NSPAR1,SYM1,FORM3)
      CALL WMATD(LU1,FILED3,NSPAR1,FORM3,NEL1,IRAM(OA),RAM(OA))
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
      INCLUDE 'indexi.for'
C     indexi.for
C
C=======================================================================
C