C
C Program SMEIGEN to read a symetric matrix SM1 and to compute general C matrix GM1 of its eigenvectors and diagonal matrix DM1 of its C eigenvalues. C C Version: 6.70 C Date: 2012, November 7 C C Coded by Petr Bulant C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, 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 C matrices SM1, GM1, and DM1. 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 GM1='string'... Name of the header file of the general C matrix of eigenvectors of matrix SM1 (output). C Default: GM1=' ' (the matrix is not written). C DM1='string'... Name of the header file of the diagonal C matrix of eigenvalues of matrix SM1 (output). C Default: DM1=' ' (the matrix is not written). C If both GM1 and DM1 equal ' ', the program is stopped. C Recent version of the program cannot deal with sparse matrices. C For general description of the files with matrices refer to file C forms.htm. C Form of the output files with matrices: C FORMM='string'... Form of the output files with matrices. Allowed C values are FORMM='formatted' and FORMM='unformatted'. C C----------------------------------------------------------------------- C Subroutines and external functions required: EXTERNAL SMEIG,SMBLO,EIGEN,ERROR,RSEP1,RSEP3T,RMAT,WMAT C SMEIG,SMBLO... This file. C EIGEN... File eigennr.for. 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 INTEGER IRAM(MRAM) EQUIVALENCE (IRAM,RAM) C EXTERNAL IND,INDG INTEGER IND,INDG C CHARACTER*80 FILSEP,FILE1,FILE2,FILE3 INTEGER M1,N,NN,NB,LU1,I1,I2,I3,I4,IBMI,IBMA,J2,J3 PARAMETER (LU1=1) C C----------------------------------------------------------------------- C C Reading the name of the file with the input data: WRITE(*,'(A)') '+SMEIGEN: 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 SMEIGEN-01 CALL ERROR('SMEIGEN-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 N=M1 NN=M1*(M1+1)/2 MAXRAM=MRAM-2*N IF (2*NN+2*N*N+N+N+1+N+N.GT.MRAM) THEN C SMEIGEN-02 CALL ERROR('SMEIGEN-02: Small dimension MRAM of array RAM') C If possible, enlarge the dimension MRAM of array RAM in include C file ram.inc. END IF C C Reading the names of the files with the matrices: CALL RSEP3T('SM1',FILE1,' ') CALL RSEP3T('GM1',FILE2,' ') CALL RSEP3T('DM1',FILE3,' ') IF (FILE1.EQ.' ') THEN C SMEIGEN-03 CALL ERROR('SMEIGEN-03: Input file with matrix SM1 not given') ENDIF IF ((FILE2.EQ.' ').AND.(FILE3.EQ.' ')) THEN C SMEIGEN-04 CALL * ERROR('SMEIGEN-04: None of the output files GM1 or DM1 given') ENDIF C C Reading input matrix: CALL RMAT(LU1,FILE1,M1,0,RAM(MAXRAM-NN+1)) C WRITE(*,'(A)') '+SMEIGEN: Working... ' C C Identification of blocks: CALL SMBLO(RAM(MAXRAM-NN+1),N,NN,NB,IRAM(MAXRAM-NN-N)) C C Preparing array for results: DO 5, I1=MAXRAM-NN-N-1-N*N-N+1,MAXRAM-NN-N-1 RAM(I1)=0. 5 CONTINUE C C Computing the eigenvalues and the eigenvectors: DO 20, I1=1,NB C Current block: IBMI=IRAM(MAXRAM-NN-N+I1-1)+1 IBMA=IRAM(MAXRAM-NN-N+I1) C Moving the block to RAM(1): I4=0 DO 10, I2=IBMI,IBMA DO 8, I3=IBMI,I2 I4=I4+1 RAM(I4)=RAM(MAXRAM-NN+IND(I3,I2)) 8 CONTINUE 10 CONTINUE J2=IBMA-IBMI+1 J3=I2*(I2+1)/2 C Computing the eigenvalues and the eigenvectors: CALL SMEIG(RAM,J2,J3) C Moving the eigenvectors to RAM(MAXRAM-NN-N-1-N*N+1): I4=0 DO 14, I2=IBMI,IBMA DO 12, I3=IBMI,IBMA I4=I4+1 RAM(MAXRAM-NN-N-1-N*N+INDG(I3,I2,N))=RAM(J3+I4) 12 CONTINUE 14 CONTINUE C Moving the eigenvalues to RAM(MAXRAM-NN-N-1-N*N-N+1): DO 16, I4=1,J2 RAM(MAXRAM-NN-N-1-N*N-N+IBMI-1+I4)=RAM(J3+J2*J2+I4) 16 CONTINUE 20 CONTINUE C C Writing output matrix GM1: IF (FILE2.NE.' ') THEN CALL WMAT(LU1,FILE2,M1,M1,RAM(MAXRAM-NN-N-1-N*N+1)) ENDIF C Writing output matrix DM1: IF (FILE3.NE.' ') THEN CALL WMAT(LU1,FILE3,M1,1,RAM(MAXRAM-NN-N-1-N*N-N+1)) ENDIF WRITE(*,'(A)') '+SMEIGEN: Done. ' C STOP END C C======================================================================= C SUBROUTINE SMEIG(SM,N,NN) C C======================================================================= C Computes the eigenvectors and eigenvalues of symmetric matrix SM. INTEGER N,NN REAL SM(NN+N*N+N) C Input: C SM... Input symmetric matrix. C N... Number of rows (and columns) of the matrix. C NN... N*(N+1)/2 C Output: C SM(NN+1)... Eigenvectors of the input matrix. C SM(NN+N*N+1)... Eigenvalues of the input matrix. INTEGER I C----------------------------------------------------------------------- CALL EIGEN(SM,SM(NN+1),N,0) DO 21 I=1,N SM(NN+N*N+I)=SM(I*(I+1)/2) 21 CONTINUE RETURN END C C======================================================================= C SUBROUTINE SMBLO(SM,N,NN,NB,IB) C C======================================================================= C Subroutine to find blocks of a symmetric matrix INTEGER N,NN,NB,IB(0:N) REAL SM(NN) C Input: C SM... Symmetric matrix. C N... Number of rows (and columns) of the matrix. C NN... N*(N+1)/2 C Output: C NB... Number of blocks of the matrix. C IB... Array with numbers of lines (and rows), at which C the blocks finish. INTEGER I1,I2 EXTERNAL IND INTEGER IND C----------------------------------------------------------------------- C Initiating first block: NB=1 IB(0)=0 IB(1)=1 I1=0 10 CONTINUE C Loop for lines of the matrix: I1=I1+1 DO 20, I2=N,IB(NB)+1,-1 C Loop for rows of the line I1: IF (SM(IND(I1,I2)).NE.0.) THEN C The block is larger: IB(NB)=I2 GOTO 21 ENDIF 20 CONTINUE 21 CONTINUE IF (IB(NB).EQ.N) THEN C End of the search for blocks. RETURN ENDIF IF (IB(NB).EQ.I1) THEN C Block finished, continuing with the next block. NB=NB+1 IB(NB)=I1+1 ENDIF GOTO 10 END C======================================================================= INTEGER FUNCTION IND(I,J) INTEGER I,J C I... Index of a line. C J... Index of a row. IND=J*(J-1)/2+I RETURN END C======================================================================= INTEGER FUNCTION INDG(I,J,K) INTEGER I,J,K C I... Index of a line. C J... Index of a row. C K... Number of lines. INDG=K*(J-1)+I RETURN END C C======================================================================= C INCLUDE 'forms.for' C forms.for INCLUDE 'error.for' C error.for INCLUDE 'length.for' C length.for INCLUDE 'mat.for' C mat.for INCLUDE 'indexi.for' C indexi.for. INCLUDE 'sep.for' C sep.for INCLUDE 'eigennr.for' C eigennr.for INCLUDE 'indexx.for' C indexx.for of Numerical Recipes INCLUDE 'tred2.for' C tred2.for of Numerical Recipes INCLUDE 'tqli.for' C tqli.for of Numerical Recipes INCLUDE 'pythag.for' C pythag.for of Numerical Recipes C C======================================================================= C