C 
C Subroutine to compute eigenvalues and eigenvectors of a real
C symmetric matrix, using subroutines from Numerical Recipes. Form
C of the subroutine is the same as the form of subroutine 'EIGEN'
C from the IBM Scientific Subroutine Package.
C
C Version: 5.50
C Date: 2001, June 4
C
C Coded by Petr Bulant
C     bulant@seis.karlov.mff.cuni.cz
C
C-----------------------------------------------------------------------
C Usage:
C     CALL EIGEN(A,R,N,MV)
C
C Input:
C   A ... Original matrix (symmetric), destroyed in computation.
C   N ... Order of matrices A and R.
C   MV .. Input code:
C         0 ... Compute eigenvalues and eigenvectors.
C         1 ... Compute eigenvalues only (currently not possible).
C Output:
C   A ... Eigenvalues of the input matrix are developed in diagonal of
C         matrix A in descending order.
C   R ... Resultant matrix of eigenvectors (stored columnwise,
C         in same sequence as eigenvalues).
C
C Remarks:
C   Original matrix A must be real symmetric (storage mode=1).
C   Matrix A cannot be in the same location as matrix R.
C   Value of MAXRAM must be set by calling program, MRAM-MAXRAM.GE.2*N.
C
C-----------------------------------------------------------------------
      SUBROUTINE EIGEN(A,R,N,MV)
      INTEGER N,MV
      REAL A(N*(N+1)/2+1),R(N,N)
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
C
      REAL WORK(MRAM)
      INTEGER I1,I2,I3,INDX(MRAM)
      EQUIVALENCE (WORK,RAM)
      EQUIVALENCE (INDX,IRAM)
C Subroutines and external functions required:
      EXTERNAL ERROR,TRED2,TQLI,INDEXX
C     ERROR ... File error.for.
C     TRED2 ... File tred2.for.
C     TQLI ... File tqli.for.
C     INDEXX ... File indexx.for.
C-----------------------------------------------------------------------
      IF (N.GT.(MRAM-MAXRAM)/2) THEN
C       EIGENNR-03
        CALL ERROR('EIGENNR-03: Small dimension of auxiliary arrays.')
C       Value of MAXRAM must be set by the calling program,
C       MRAM-MAXRAM must be greater or equal 2*N.
      ENDIF
C
      IF (MV.EQ.0) THEN
C       Computation of eigenvalues and eigenvectors.
C
C       Preparing matrix R:
        I3=0
        DO 12, I1=1,N
          DO 10, I2=1,I1
            I3=I3+1
            R(I1,I2)=A(I3)
            R(I2,I1)=A(I3)
   10     CONTINUE
   12   CONTINUE
C
C       Preparing a tridiagonal matrix from R:
        CALL TRED2(R,N,N,A,A(N+1))
C
C       Computing eigenvalues and eigenvectors of matrix R:
        CALL TQLI(A,A(N+1),N,N,R)
C
C       Sorting the eigenvalues into descending order:
        CALL INDEXX(N,A,INDX(MAXRAM+N+1))
C       Reorganizing eigenvalues:
        DO 14, I2=1,N
C         Store the element to WORK:
          WORK(MAXRAM+I2)=A(I2)
   14   CONTINUE
        DO 16, I2=1,N
C         Copy it back in the rearranged order:
          A(I2)=WORK(MAXRAM+INDX(MAXRAM+N+N+1-I2))
   16   CONTINUE
C       Reorganizing eigenvectors:
C       Loop over lines:
        DO 25, I1=1,N
C         For each element of the line:
          DO 18, I2=1,N
C           Store the element to WORK:
            WORK(MAXRAM+I2)=R(I1,I2)
   18     CONTINUE
          DO 20, I2=1,N
C           Copy it back in the rearranged order:
            R(I1,I2)=WORK(MAXRAM+INDX(MAXRAM+N+N+1-I2))
   20     CONTINUE
   25   CONTINUE
C
C       Writing eigenvalues to the diagonal of A:
        DO 27, I1=N,2,-1
          A((I1*(I1+1))/2)=A(I1)
   27   CONTINUE
      ELSEIF (MV.EQ.1) THEN
C       Computation of eigenvalues only.
C       EIGENNR-01
        CALL ERROR('EIGENNR-01: MV=1 currently not supported.')
      ELSE
C       EIGENNR-02
        CALL ERROR('EIGENNR-02: Wrong value of MV.')
      ENDIF
      RETURN
      END