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: 5.50
C Date: 2002, May 21
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
C             matrices SM1, GM1, and DM1.
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 the input matrix.
C             No default, 'SM1' must be specified and cannot be blank.
C     GM1='string' ... Name of the file containing general
C             matrix of eigenvectors of matrix SM1 (output).
C             Default: GM1=' ' (the matrix is not written).
C     DM1='string' ... Name of the file containing 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     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 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: SM ... input symmetric matrix
C            N  ... number of rows (and columns) of the matrix
C            NN ... N*(N+1)/2
C     Output: 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: SM ... symmetric matrix
C            N  ... number of rows (and columns) of the matrix
C            NN ... N*(N+1)/2
C     Output: 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 '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