C 
C Program GMDMGMT to compute product SM1=GM1*DM1*GM1T of general matrix
C GM1, diagonal matrix DM1 and transposed matrix GM1.  Resulting matrix
C SM1 is symmetric.
C
C Version: 5.40
C Date: 2000, February 21
C
C Coded by Petr Bulant
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: bulant@seis.karlov.mff.cuni.cz
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 Data specifying dimensions of the matrices:
C     M1='string'... Name of the file containing a single integer number
C             specifying the number of rows of general matrix GM1 and
C             rows and columns of symmetric matrix SM1.
C             Default: M1=' ' means that the number is 1.
C     M2='string'... Name of the file containing a single integer number
C             specifying the number of columns of matrix GM1 and rows
C             and columns of diagonal matrix DM1.
C             Default: M2=' ' means that the number is 1.
C Filenames of the files with the matrices:
C     GM1='string' ... Name of the input file containing matrix GM1.
C             No default, 'GM1' must be specified and cannot be blank.
C     DM1='string' ... Name of the input file containing matrix DM1.
C             If DM1 is blank (default), a unit matrix is used in
C             place of DM1.
C             Default: 'DM1'=' ' (means that DM1 is unit matrix).
C     SM1='string' ... Name of the output file to contain matrix SM1.
C             No default, 'SM1' must be specified and cannot be blank.
C     For general description of the files with matrices refer to file
C     forms.htm.
C
C-----------------------------------------------------------------------
C Subroutines and external functions required:
      EXTERNAL ERROR,RSEP1,RSEP3T,RMAT,WMAT
C     ERROR ... File error.for.
C     RSEP1,RSEP3T ... File sep.for.
C     RMAT,WMAT ... File forms.for.
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      CHARACTER*80 FILE1
      INTEGER M1,M2,M1M2,M1M1,LU1,I1,I2,I3,I
      PARAMETER (LU1=1)
C-----------------------------------------------------------------------
C
C     Reading name of SEP file with input data:
      WRITE(*,'(A)') '+GMDMGMT: Enter input filename: '
      FILE1=' '
      READ(*,*) FILE1
C
C     Reading all data from the SEP file into the memory:
      IF (FILE1.NE.' ') THEN
        CALL RSEP1(LU1,FILE1)
      ELSE
C       GMDMGMT-01
        CALL ERROR('GMDMGMT-01: SEP file not given')
      ENDIF
C
C     Reading the dimensions of the matrices:
      CALL RSEP3T('M1',FILE1,' ')
      IF (FILE1.EQ.' ') THEN
        M1=1
      ELSE
        OPEN(LU1,FILE=FILE1,STATUS='OLD')
        READ(LU1,*) M1
        CLOSE(LU1)
      ENDIF
      CALL RSEP3T('M2',FILE1,' ')
      IF (FILE1.EQ.' ') THEN
        M2=1
      ELSE
        OPEN(LU1,FILE=FILE1,STATUS='OLD')
        READ(LU1,*) M2
        CLOSE(LU1)
      ENDIF
      M1M2=M1*M2
      M1M1=M1*(M1+1)/2
C
      IF (M1M2+M2+M1M1.GT.MRAM) THEN
C       GMDMGMT-02
        CALL ERROR('GMDMGMT-02: Small dimension MRAM of array RAM')
      END IF
C
C     Reading input matrices:
      CALL RSEP3T('GM1',FILE1,' ')
      IF (FILE1.EQ.' ') THEN
C       GMDMGMT-03
        CALL ERROR('GMDMGMT-03: Input file with matrix GM1 not given.')
      ENDIF
      CALL RMAT(LU1,FILE1,M1,M2,RAM)
      CALL RSEP3T('DM1',FILE1,' ')
      IF (FILE1.EQ.' ') THEN
        DO 5 I=M1M2+1,M1M2+M2
          RAM(I)=1.
   5    CONTINUE
      ELSE
        CALL RMAT(LU1,FILE1,M2,1,RAM(M1M2+1))
      ENDIF
      CALL RSEP3T('SM1',FILE1,' ')
      IF (FILE1.EQ.' ') THEN
C       GMDMGMT-05
        CALL ERROR('GMDMGMT-05: Output file with matrix SM1 not given.')
      ENDIF
C
C     Multiplication:
      WRITE(*,'(A)') '+GMDMGMT: Calculating...        '
      DO 10 I=M1M2+M2+1,M1M2+M2+M1M1
        RAM(I)=0.
  10  CONTINUE
      DO 13 I3=1,M2
        I=M1M2+M2
        DO 12 I2=M1*(I3-1)+1,M1*I3
          AUX=RAM(M1M2+I3)*RAM(I2)
          DO 11 I1=M1*(I3-1)+1,I2
            I=I+1
            RAM(I)=RAM(I)+RAM(I1)*AUX
  11      CONTINUE
  12    CONTINUE
  13  CONTINUE
C
C     Writing output matrix SM1:
      CALL WMAT(LU1,FILE1,M1,0,RAM(M1M2+M2+1))
C
      WRITE(*,'(A)') '+GMDMGMT: Done.                 '
      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
C
C=======================================================================
C