C 
C Program SMSMSM to compute product SM3=SM1*SM2*SM1 of symmetric
C matrices SM1 and SM2.
C
C Version: 5.50
C Date: 2001, May 12
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 (and columns) of symmetric
C             matrices SM1, SM2 and SM3.
C             Default: M1=' ' means that the number is 1.
C Filenames of the files with the matrices:
C     SM1='string'... Name of the file containing matrix SM1 (input).
C             No default, SM1 must be specified and cannot be blank.
C     SM2='string'... Name of the file containing matrix SM2 (input).
C             No default, SM2 must be specified and cannot be blank.
C     SM3='string'... Name of the file containing symmetric
C             matrix SM3=SM1*SM2*SM1 (output).
C             No default, SM3 must be specified and cannot be blank.
C     For general description of the files with matrices refer to file
C     forms.htm.
C Form of the files with matrices:
C     FORMM='string' ... Form of the files with matrices. Allowed values
C             are FORMM='formatted' and FORMM='unformatted'. If the form
C             differs for input and for output files, FORMMR and FORMMW
C             should be used instead of FORMM.
C             Default: FORMM='formatted'
C     FORMMR='string' ... Form of the files with matrices to be read.
C             Default: FORMMR=FORMM
C     FORMMW='string' ... Form of the files with matrices to be written.
C             Default: FORMMW=FORMM
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 FILSEP,FILE1,FILE2,FILE3
      INTEGER LU1,M1,I0,I1,I2,J0,J1,J2
      PARAMETER (LU1=1)
      REAL AUX
C
C-----------------------------------------------------------------------
C
C     Reading a name of the file with the input data:
      WRITE(*,'(A)') '+SMSMSM: 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       SMSMSM-01
        CALL ERROR('SMSMSM-01: SEP file not given')
      ENDIF
C
C     Reading the dimension 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
C
C     Memory for 3 general matrices GM1, GM2 and GM3:
      IF (3*M1*M1.GT.MRAM) THEN
C       SMSMSM-02
        CALL ERROR('SMSMSM-02: Small dimension MRAM of array RAM')
      END IF
C
C     Reading the names of the files with the matrices:
      CALL RSEP3T('SM1',FILE1,' ')
      CALL RSEP3T('SM2',FILE2,' ')
      CALL RSEP3T('SM3',FILE3,' ')
      IF (FILE1.EQ.' ') THEN
C       SMSMSM-03
        CALL ERROR('SMSMSM-03: Input file with matrix SM1 not given.')
      ENDIF
      IF (FILE2.EQ.' ') THEN
C       SMSMSM-04
        CALL ERROR('SMSMSM-04: Input file with matrix SM2 not given.')
      ENDIF
      IF (FILE3.EQ.' ') THEN
C       SMSMSM-05
        CALL ERROR('SMSMSM-05: Output file with matrix SM3 not given.')
      ENDIF
C
C     Reading input matrices:
C     Reading SM1
      CALL RMAT(LU1,FILE1,M1,0,RAM(2*M1*M1+1))
C     Storing SM1 (RAM index J0) as GM1 (RAM indices J1 and J2)
      J0=2*M1*M1
      DO 12 I2=1,M1
        J1=M1*(I2-1)
        J2=I2-M1
        DO 11 I1=1,I2
          J0=J0+1
          J1=J1+1
          J2=J2+M1
          AUX=RAM(J0)
          RAM(J1)=AUX
          RAM(J2)=AUX
   11   CONTINUE
   12 CONTINUE
C     Reading SM2
      CALL RMAT(LU1,FILE2,M1,0,RAM(2*M1*M1+1))
C     Storing SM2 (RAM index J0) as GM2 (RAM indices J1 and J2)
      J0=2*M1*M1
      DO 22 I2=1,M1
        J1=M1*M1+M1*(I2-1)
        J2=M1*M1+I2-M1
        DO 21 I1=1,I2
          J0=J0+1
          J1=J1+1
          J2=J2+M1
          AUX=RAM(J0)
          RAM(J1)=AUX
          RAM(J2)=AUX
   21   CONTINUE
   22 CONTINUE
      WRITE(*,'(A)') '+SMSMSM: Working...            '
C
C     Multiplication:
C     Storing GM2*GM1 as GM3
      DO 32 I2=1,M1
        DO 31 I1=1,M1
          J1=M1*M1+M1*(I1-1)
          J2=      M1*(I2-1)
          AUX=0.
          DO 30 I0=1,M1
            J1=J1+1
            J2=J2+1
            AUX=AUX+RAM(J1)*RAM(J2)
   30     CONTINUE
          RAM(2*M1*M1+M1*(I2-1)+I1)=AUX
   31   CONTINUE
   32 CONTINUE
C     Storing SM3=GM1*GM3 in place of GM2
      DO 42 I2=1,M1
        DO 41 I1=1,I2
          J1=        M1*(I1-1)
          J2=2*M1*M1+M1*(I2-1)
          AUX=0.
          DO 40 I0=1,M1
            J1=J1+1
            J2=J2+1
            AUX=AUX+RAM(J1)*RAM(J2)
   40     CONTINUE
          RAM(M1*M1+I2*(I2-1)/2+I1)=AUX
   41   CONTINUE
   42 CONTINUE
C
C     Writing output matrix SM3:
      CALL WMAT(LU1,FILE3,M1,0,RAM(M1*M1+1))
      WRITE(*,'(A)') '+SMSMSM: 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
C
C=======================================================================
C