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.50 C Date: 2000, October 20 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 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 FILE1 CHARACTER*72 TXTERR 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 I1=M1M2+M2+M1M1-MRAM WRITE(TXTERR,'(A,I9,A)') * 'GMDMGMT-02: Array RAM too small,',I1,' units missing.' CALL ERROR(TXTERR) 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