C                                                            <HTML><PRE>
C Program TRSMSM to compute trace tr(SM1*SM2) of the product of two
C symmetric matrices SM1 and SM2
C
C Version: 6.20
C Date: 2008, June 12
C
C Coded by Petr Bulant and Ludek Klimes
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C
C.......................................................................
C                                                    <A NAME="data"></A>
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                                                     <A NAME="SEP"></A>
C Input data file 'SEP':
C     File 'SEP' has the form of the <A HREF="../forms/sep.htm">SEP</A>
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 and SM2.
C             Default: M1=' ' means that the number is 1.
C Filenames of the matrix header files:
C     SM1='string' ... Name of the header file of the input symmetric
C             matrix SM1.
C             No default, SM1 must be specified and cannot be blank.
C     SM2='string'... Name of the  header file of the input symmetric
C             matrix SM2.
C             No default, SM2 must be specified and cannot be blank.
C     TRACE='string'... Name of the header file of the output
C             general matrix containing a single value, trace
C             tr(SM1*SM2) of the product of two symmetric matrices
C             SM1 and SM2.
C             No default, TRACE must be specified and cannot be blank.
C     Recent version of the program cannot deal with sparse matrices.
C     For general description of the files with matrices refer to file
C     <A HREF="../forms/forms.htm#MAT">forms.htm</A>.
C Form of the output file with matrix TRACE:
C     FORMM='string' ... Form of the output files with matrices. Allowed
C             values are FORMM='formatted' and FORMM='unformatted'.
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL ERROR,RSEP1,RSEP3T,RMAT,WMAT
C     ERROR ... File <A HREF="error.for" TYPE="text/html">error.for</A>.
C     RSEP1,RSEP3T ... File
C     <A HREF="sep.for" TYPE="text/html">sep.for</A>.
C     RMAT,WMAT ... File <A HREF="mat.for" TYPE="text/html">mat.for</A>.
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     <A HREF="ram.inc" TYPE="text/html">ram.inc</A>
C
      CHARACTER*80 FILSEP,FILE1,FILE2,FILE3
      INTEGER LU1,M1,NN,I1,I2,I
      PARAMETER (LU1=1)
      REAL TRACE
C
C-----------------------------------------------------------------------
C
C     Reading a name of the file with the input data:
      FILSEP=' '
      WRITE(*,'(A)') '+TRSMSM: Enter input filename: '
      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       <A NAME="error01"></A><FONT COLOR="RED">TRSMSM-01</FONT>
        CALL ERROR('TRSMSM-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
      NN=M1*(M1+1)/2
C
      IF (2*NN.GT.MRAM) THEN
C       <A NAME="error02"></A><FONT COLOR="RED">TRSMSM-02</FONT>
        CALL ERROR('TRSMSM-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('TRACE',FILE3,' ')
      IF (FILE1.EQ.' ') THEN
C       <A NAME="error03"></A><FONT COLOR="RED">TRSMSM-03</FONT>
        CALL ERROR('TRSMSM-03: Input file with matrix SM1 not given')
      ENDIF
      IF (FILE2.EQ.' ') THEN
C       <A NAME="error04"></A><FONT COLOR="RED">TRSMSM-04</FONT>
        CALL ERROR('TRSMSM-04: Input file with matrix SM2 not given')
      ENDIF
      IF (FILE3.EQ.' ') THEN
C       <A NAME="error05"></A><FONT COLOR="RED">TRSMSM-05</FONT>
        CALL ERROR('TRSMSM-05: Output filename TRACE not specified')
      ENDIF
C
C     Reading input matrices:
      CALL RMAT(LU1,FILE1,M1,0,RAM)
      CALL RMAT(LU1,FILE2,M1,0,RAM(NN+1))
C
      WRITE(*,'(A)') '+TRSMSM: Working...            '
C
C     Calculating the trace:
      TRACE=0.
      I=0
C     Loop over columns:
      DO 12, I2=1,M1
C       Loop over rows:
        DO 11, I1=1,I2-1
          I=I+1
          TRACE=TRACE+2.*RAM(I)*RAM(NN+I)
  11    CONTINUE
        I=I+1
        TRACE=TRACE+RAM(I)*RAM(NN+I)
  12  CONTINUE
C
C     Writing output file TRACE:
      RAM(1)=TRACE
      IF (FILE3.NE.' ') THEN
        CALL WMAT(LU1,FILE3,1,1,RAM)
      ENDIF
      WRITE(*,'(A)') '+TRSMSM: Done.                 '
C
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     <A HREF="../forms/error.for" TYPE="text/html">error.for</A>
      INCLUDE 'sep.for'
C     <A HREF="../forms/sep.for" TYPE="text/html">sep.for</A>
      INCLUDE 'forms.for'
C     <A HREF="../forms/forms.for" TYPE="text/html">forms.for</A>
      INCLUDE 'length.for'
C     <A HREF="../forms/length.for" TYPE="text/html">length.for</A>
      INCLUDE 'mat.for'
C     <A HREF="../forms/mat.for" TYPE="text/html">mat.for</A>
C
C=======================================================================
C                                                                 </PRE>
