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.40
C Date: 2000, February 21
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
C-----------------------------------------------------------------------
      SUBROUTINE EIGEN(A,R,N,MV)
      INTEGER N,MV
      REAL A(N*(N+1)/2+1),R(N,N)
      INTEGER M
      PARAMETER (M=1000)
      REAL WORK(M)
      INTEGER I1,I2,I3,INDX(M)
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.M) THEN
C       EIGENNR-03
        CALL ERROR('EIGENNR-03: Small dimension M of auxiliary arrays.')
      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)
C       Reorganizing eigenvalues:
        DO 14, I2=1,N
C         Store the element to WORK:
          WORK(I2)=A(I2)
   14   CONTINUE
        DO 16, I2=1,N
C         Copy it back in the rearranged order:
          A(I2)=WORK(INDX(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(I2)=R(I1,I2)
   18     CONTINUE
          DO 20, I2=1,N
C           Copy it back in the rearranged order:
            R(I1,I2)=WORK(INDX(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