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