C 22 DEC 06 - MWS - USE IMPLICIT STATEMENTS TO COVER BLAS NAME HACK
C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C  1 JUN 05 - DGF - GLDIAG: OPTION FORCING ALL DIAGS TO MASTER (FOR FMO)
C  3 FEB 97 - MWS - MAKE ESTPI A SUBROUTINE, AS IT MODIFIES ARG ANORM.
C 16 JAN 97 - STE - EINVIT: AVOID ACCIDENTAL ROW EXCHANGES
C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
C 31 MAR 94 - MWS - ADD A VARIABLE TO END OF MACHSW COMMON
C 26 JUN 93 - MWS - ETRED3: ADD RETURN FOR SPECIAL CASE N=1
C  4 JAN 92 - TLW - MAKE WRITES PARALLEL;ADD COMMON PAR
C 30 AUG 91 - MWS - JACDIA: LIMIT ITERATIONS, USE EPSLON IN TEST.
C 14 JUL 91 - MWS - JACOBI DIAGONALIZATION ALLOWS FOR LDVEC.NE.N
C 29 JAN 91 - TLW - GLDIAG: CHANGED COMMON DIAGSW TO MACHSW
C 29 OCT 90 - STE - FIX JACDIA UNDEFINED VARIABLE BUG
C 14 SEP 90 - MK  - NEW JACOBI DIAGONALIZATION (KDIAG=3)
C 27 MAR 88 - MWS - ALLOW FOR VECTOR ROUTINE IN GLDIAG
C 11 AUG 87 - MWS - SANITIZE CONSTANTS IN EQLRAT
C 15 FEB 87 - STE - FIX EINVIT SUB-MATRIX LOOP LIMIT
C                   SCRATCH ARRAYS ARE N*8 REAL AND N INTEGER
C  8 DEC 86 - STE - USE PERF INDEX FROM ESTPI1 TO JUDGE EINVIT FAILURE
C 30 NOV 86 - STE - DELETE LIGENB, MAKE EVVRSP DEFAULT
C                   (GIVEIS FAILS ON CRAY FOR BENCHMC AND BENCHCI)
C  7 JUL 86 - JAB - SANITIZE FLOATING POINT CONSTANTS
C 11 OCT 85 - STE - LIGENB,TQL2: USE DROT,DSWAP; TINVTB: SCALE VECTOR
C                   BEFORE NORMALIZING; GENERIC FUNCTIONS
C 24 FEB 84 - STE - INITIALIZE INDEX ARRAY FOR LIGENB IN GLDIAG
C  1 DEC 83 - STE - CHANGE MACHEP FROM 2**-54 TO 2**-50
C 28 SEP 82 - MWS - CONVERT TO IBM
C
C*MODULE EIGEN   *DECK EINVIT
      SUBROUTINE EINVIT(NM,N,D,E,E2,M,W,IND,Z,IERR,RV1,RV2,RV3,RV4,RV6)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C*
C*    AUTHORS-
C*       THIS IS A MODIFICATION OF TINVIT FROM EISPACK EDITION 3
C*       DATED AUGUST 1983.
C*       TINVIT IS A TRANSLATION OF THE INVERSE ITERATION TECHNIQUE
C*       IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON.
C*       HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971).
C*       THIS VERSION IS BY S. T. ELBERT (AMES LABORATORY-USDOE)
C*
C*    PURPOSE -
C*       THIS ROUTINE FINDS THOSE EIGENVECTORS OF A TRIDIAGONAL
C*       SYMMETRIC MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES.
C*
C*    METHOD -
C*       INVERSE ITERATION.
C*
C*    ON ENTRY -
C*       NM     - INTEGER
C*                MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C*                ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
C*                DIMENSION STATEMENT.
C*       N      - INTEGER
C*       D      - W.P. REAL (N)
C*                CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
C*       E      - W.P. REAL (N)
C*                CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
C*                IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
C*       E2     - W.P. REAL (N)
C*                CONTAINS THE SQUARES OF CORRESPONDING ELEMENTS OF E,
C*                WITH ZEROS CORRESPONDING TO NEGLIGIBLE ELEMENTS OF E.
C*                E(I) IS CONSIDERED NEGLIGIBLE IF IT IS NOT LARGER THAN
C*                THE PRODUCT OF THE RELATIVE MACHINE PRECISION AND THE
C*                SUM OF THE MAGNITUDES OF D(I) AND D(I-1).  E2(1) MUST
C*                CONTAIN 0.0 IF THE EIGENVALUES ARE IN ASCENDING ORDER,
C*                OR 2.0 IF THE EIGENVALUES ARE IN DESCENDING ORDER.
C*                IF TQLRAT, BISECT, TRIDIB, OR IMTQLV
C*                HAS BEEN USED TO FIND THE EIGENVALUES, THEIR
C*                OUTPUT E2 ARRAY IS EXACTLY WHAT IS EXPECTED HERE.
C*       M      - INTEGER
C*                THE NUMBER OF SPECIFIED EIGENVECTORS.
C*       W      - W.P. REAL (M)
C*                CONTAINS THE M EIGENVALUES IN ASCENDING
C*                OR DESCENDING ORDER.
C*       IND    - INTEGER (M)
C*                CONTAINS IN FIRST M POSITIONS THE SUBMATRIX INDICES
C*                ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W --
C*                1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX
C*                FROM THE TOP, 2 FOR THOSE BELONGING TO THE SECOND
C*                SUBMATRIX, ETC.
C*       IERR   - INTEGER (LOGICAL UNIT NUMBER)
C*                LOGICAL UNIT FOR ERROR MESSAGES
C*
C*    ON EXIT -
C*       ALL INPUT ARRAYS ARE UNALTERED.
C*       Z      - W.P. REAL (NM,M)
C*                CONTAINS THE ASSOCIATED SET OF ORTHONORMAL
C*                EIGENVECTORS. ANY VECTOR WHICH WHICH FAILS TO CONVERGE
C*                IS LEFT AS IS (BUT NORMALIZED) WHEN ITERATING STOPPED.
C*       IERR   - INTEGER
C*                SET TO
C*                ZERO    FOR NORMAL RETURN,
C*                -R      IF THE EIGENVECTOR CORRESPONDING TO THE R-TH
C*                        EIGENVALUE FAILS TO CONVERGE IN 5 ITERATIONS.
C*                        (ONLY LAST FAILURE TO CONVERGE IS REPORTED)
C*
C*       RV1, RV2, RV3, RV4, AND RV6 ARE TEMPORARY STORAGE ARRAYS.
C*
C*       RV1    - W.P. REAL (N)
C*                DIAGONAL ELEMENTS OF U FROM LU DECOMPOSITION
C*       RV2    - W.P. REAL (N)
C*                SUPER(1)-DIAGONAL ELEMENTS OF U FROM LU DECOMPOSITION
C*       RV3    - W.P. REAL (N)
C*                SUPER(2)-DIAGONAL ELEMENTS OF U FROM LU DECOMPOSITION
C*       RV4    - W.P. REAL (N)
C*                ELEMENTS DEFINING L IN LU DECOMPOSITION
C*       RV6    - W.P. REAL (N)
C*                APPROXIMATE EIGENVECTOR
C*
C*    DIFFERENCES FROM EISPACK 3 -
C*       EPS3 IS SCALED BY  EPSCAL  (ENHANCES CONVERGENCE, BUT
C*          LOWERS ACCURACY)!
C*       ONE MORE ITERATION (MINIMUM 2) IS PERFORMED AFTER CONVERGENCE
C*          (ENHANCES ACCURACY)!
C*       REPLACE LOOP WITH PYTHAG WITH SINGLE CALL TO DNRM2!
C*       IF NOT CONVERGED, USE PERFORMANCE INDEX TO DECIDE ON ERROR
C*          VALUE SETTING, BUT DO NOT STOP!
C*       L.U. FOR ERROR MESSAGES PASSED THROUGH IERR
C*       USE PARAMETER STATEMENTS AND GENERIC INTRINSIC FUNCTIONS
C*       USE LEVEL 1 BLAS
C*       USE IF-THEN-ELSE TO CLARIFY LOGIC
C*       LOOP OVER SUBSPACES MADE INTO DO LOOP.
C*       LOOP OVER INVERSE ITERATIONS MADE INTO DO LOOP
C*       ZERO ONLY REQUIRED PORTIONS OF OUTPUT VECTOR
C*
C*    NOTE -
C*       QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
C*       B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
C*
C
      LOGICAL CONVGD,GOPARR,DSKWRK,MASWRK
C
      INTEGER GROUP,I,IERR,ITS,J,JJ,M,N,NM,P,Q,R,S,SUBMAT,TAG
      INTEGER IND(M)
C
      DOUBLE PRECISION D(N),E(N),E2(N),W(M),Z(NM,M)
      DOUBLE PRECISION RV1(N),RV2(N),RV3(N),RV4(N),RV6(N)
      DOUBLE PRECISION ANORM,EPS2,EPS3,EPS4,NORM,ORDER,RHO,U,UK,V
      DOUBLE PRECISION X0,X1,XU
      DOUBLE PRECISION EPSCAL,GRPTOL,HUNDRD,ONE,TEN,ZERO
      DOUBLE PRECISION EPSLON
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO = 0.0D+00, ONE = 1.0D+00, GRPTOL = 0.001D+00)
      PARAMETER (EPSCAL = 0.5D+00, HUNDRD = 100.0D+00, TEN = 10.0D+00)
C
  001 FORMAT(' EIGENVECTOR ROUTINE EINVIT DID NOT CONVERGE FOR VECTOR'
     *      ,I5,'.  NORM =',1P,E10.2,' PERFORMANCE INDEX =',E10.2/
     *      ' (AN ERROR HALT WILL OCCUR IF THE PI IS GREATER THAN 100)')
C
C-----------------------------------------------------------------------
C
      LUEMSG = IERR
      IERR = 0
      X0 = ZERO
      UK = ZERO
      NORM = ZERO
      EPS2 = ZERO
      EPS3 = ZERO
      EPS4 = ZERO
      GROUP = 0
      TAG = 0
      ORDER = ONE - E2(1)
      Q = 0
      DO 930 SUBMAT = 1, N
         P = Q + 1
C
C        .......... ESTABLISH AND PROCESS NEXT SUBMATRIX ..........
C
         DO 120 Q = P, N-1
            IF (E2(Q+1) .EQ. ZERO) GO TO 140
  120    CONTINUE
         Q = N
C
C        .......... FIND VECTORS BY INVERSE ITERATION ..........
C
  140    CONTINUE
         TAG = TAG + 1
         ANORM = ZERO
         S = 0
C
         DO 920 R = 1, M
            IF (IND(R) .NE. TAG) GO TO 920
            ITS = 1
            X1 = W(R)
            IF (S .NE. 0) GO TO 510
C
C        .......... CHECK FOR ISOLATED ROOT ..........
C
            XU = ONE
            IF (P .EQ. Q) THEN
               RV6(P) = ONE
               CONVGD = .TRUE.
               GO TO 860
C
            END IF
            NORM = ABS(D(P))
            DO 500 I = P+1, Q
               NORM = MAX( NORM, ABS(D(I)) + ABS(E(I)) )
  500       CONTINUE
C
C        .......... EPS2 IS THE CRITERION FOR GROUPING,
C                   EPS3 REPLACES ZERO PIVOTS AND EQUAL
C                   ROOTS ARE MODIFIED BY EPS3,
C                   EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW .........
C
            EPS2 = GRPTOL * NORM
            EPS3 = EPSCAL * EPSLON(NORM)
            UK = Q - P + 1
            EPS4 = UK * EPS3
            UK = EPS4 / SQRT(UK)
            S = P
            GROUP = 0
            GO TO 520
C
C        .......... LOOK FOR CLOSE OR COINCIDENT ROOTS ..........
C
  510       IF (ABS(X1-X0) .GE. EPS2) THEN
C
C                 ROOTS ARE SEPERATE
C
               GROUP = 0
            ELSE
C
C                 ROOTS ARE CLOSE
C
               GROUP = GROUP + 1
               IF (ORDER * (X1 - X0) .LE. EPS3) X1 = X0 + ORDER * EPS3
            END IF
C
C        .......... ELIMINATION WITH INTERCHANGES AND
C                   INITIALIZATION OF VECTOR ..........
C
  520       CONTINUE
C
            U = D(P) - X1
            V = E(P+1)
            RV6(P) = UK
            DO 550 I = P+1, Q
               RV6(I) = UK
               IF (ABS(E(I)) .GT. ABS(U)) THEN
C
C                 EXCHANGE ROWS BEFORE ELIMINATION
C
C                  *** WARNING -- A DIVIDE CHECK MAY OCCUR HERE IF
C                      E2 ARRAY HAS NOT BEEN SPECIFIED CORRECTLY .......
C
                  XU = U / E(I)
                  RV4(I) = XU
                  RV1(I-1) = E(I)
                  RV2(I-1) = D(I) - X1
                  RV3(I-1) = E(I+1)
                  U = V - XU * RV2(I-1)
                  V = -XU * RV3(I-1)
C
               ELSE
C
C                    STRAIGHT ELIMINATION
C
                  XU = E(I) / U
                  RV4(I) = XU
                  RV1(I-1) = U
                  RV2(I-1) = V
                  RV3(I-1) = ZERO
                  U = D(I) - X1 - XU * V
                  V = E(I+1)
               END IF
  550       CONTINUE
C
            IF (ABS(U) .LE. EPS3) U = EPS3
            RV1(Q) = U
            RV2(Q) = ZERO
            RV3(Q) = ZERO
C
C              DO INVERSE ITERATIONS
C
            CONVGD = .FALSE.
            DO 800 ITS = 1, 5
               IF (ITS .EQ. 1) GO TO 600
C
C                    .......... FORWARD SUBSTITUTION ..........
C
                  IF (NORM .EQ. ZERO) THEN
                     RV6(S) = EPS4
                     S = S + 1
                     IF (S .GT. Q) S = P
                  ELSE
                     XU = EPS4 / NORM
                     CALL DSCAL(Q-P+1, XU, RV6(P), 1)
                  END IF
C
C                     ... ELIMINATION OPERATIONS ON NEXT VECTOR
C
                  DO 590 I = P+1, Q
                     U = RV6(I)
C
C                         IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE
C                         WAS PERFORMED EARLIER IN THE
C                         TRIANGULARIZATION PROCESS ..........
C
C-STE-               IF (RV1(I-1) .EQ. E(I)) THEN
                     IF (RV1(I-1) .EQ. E(I) .AND. RV4(I) .NE. ONE) THEN
                        U = RV6(I-1)
                        RV6(I-1) = RV6(I)
                     ELSE
                        U = RV6(I)
                     END IF
                     RV6(I) = U - RV4(I) * RV6(I-1)
  590             CONTINUE
  600          CONTINUE
C
C           .......... BACK SUBSTITUTION
C
               RV6(Q) = RV6(Q) / RV1(Q)
               V = U
               U = RV6(Q)
               NORM = ABS(U)
               DO 620 I = Q-1, P, -1
                  RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I)
                  V = U
                  U = RV6(I)
                  NORM = NORM + ABS(U)
  620          CONTINUE
               IF (GROUP .EQ. 0) GO TO 700
C
C                 ....... ORTHOGONALIZE WITH RESPECT TO PREVIOUS
C                         MEMBERS OF GROUP ..........
C
                  J = R
                  DO 680 JJ = 1, GROUP
  630                J = J - 1
                     IF (IND(J) .NE. TAG) GO TO 630
                     CALL DAXPY(Q-P+1, -DDOT(Q-P+1,RV6(P),1,Z(P,J),1),
     *                          Z(P,J),1,RV6(P),1)
  680             CONTINUE
                  NORM = DASUM(Q-P+1, RV6(P), 1)
  700          CONTINUE
C
               IF (CONVGD) GO TO 840
               IF (NORM .GE. ONE) CONVGD = .TRUE.
  800       CONTINUE
C
C        .......... NORMALIZE SO THAT SUM OF SQUARES IS
C                   1 AND EXPAND TO FULL ORDER ..........
C
  840       CONTINUE
C
            XU = ONE / DNRM2(Q-P+1,RV6(P),1)
C
  860       CONTINUE
            DO 870 I = 1, P-1
               Z(I,R) = ZERO
  870       CONTINUE
            DO 890 I = P,Q
               Z(I,R) = RV6(I) * XU
  890       CONTINUE
            DO 900 I = Q+1, N
               Z(I,R) = ZERO
  900       CONTINUE
C
            IF (.NOT.CONVGD) THEN
               CALL ESTPI(RHO,Q-P+1,X1,D(P),E(P),Z(P,R),ANORM)
               IF (RHO .GE. TEN .AND. LUEMSG .GT. 0 .AND. MASWRK)
     *             WRITE(LUEMSG,001) R,NORM,RHO
C
C               *** SET ERROR -- NON-CONVERGED EIGENVECTOR ..........
C
               IF (RHO .GT. HUNDRD) IERR = -R
            END IF
C
            X0 = X1
  920    CONTINUE
C
         IF (Q .EQ. N) GO TO 940
  930 CONTINUE
  940 CONTINUE
      RETURN
      END
C*MODULE EIGEN   *DECK ELAUM
      SUBROUTINE ELAU(HINV,L,D,A,E)
C
      DOUBLE PRECISION A(*)
      DOUBLE PRECISION D(L)
      DOUBLE PRECISION E(L)
      DOUBLE PRECISION F
      DOUBLE PRECISION G
      DOUBLE PRECISION HALF
      DOUBLE PRECISION HH
      DOUBLE PRECISION HINV
      DOUBLE PRECISION ZERO
C
      PARAMETER (ZERO = 0.0D+00, HALF = 0.5D+00)
C
      JL = L
      E(1) = A(1) * D(1)
      JK = 2
      DO 210 J = 2, JL
         F = D(J)
         G = ZERO
         JM1 = J - 1
C
         DO 200 K = 1, JM1
            G = G + A(JK) * D(K)
            E(K) = E(K) + A(JK) * F
            JK = JK + 1
  200    CONTINUE
C
         E(J) = G + A(JK) * F
         JK = JK + 1
  210 CONTINUE
C
C        .......... FORM P ..........
C
      F = ZERO
      DO 245 J = 1, L
         E(J) = E(J) * HINV
         F = F + E(J) * D(J)
  245 CONTINUE
C
C     .......... FORM Q ..........
C
      HH = F * HALF * HINV
      DO 250 J = 1, L
  250 E(J) = E(J) - HH * D(J)
C
      RETURN
      END
C*MODULE EIGEN   *DECK EPSLON
      DOUBLE PRECISION FUNCTION EPSLON (X)
C*
C*    AUTHORS -
C*       THIS ROUTINE WAS TAKEN FROM EISPACK EDITION 3 DATED 4/6/83
C*       THIS VERSION IS BY S. T. ELBERT, AMES LABORATORY-USDOE NOV 1986
C*
C*    PURPOSE -
C*       ESTIMATE UNIT ROUNDOFF IN QUANTITIES OF SIZE X.
C*
C*    ON ENTRY -
C*       X      - WORKING PRECISION REAL
C*                VALUES TO FIND EPSLON FOR
C*
C*    ON EXIT -
C*       EPSLON - WORKING PRECISION REAL
C*                SMALLEST POSITIVE VALUE SUCH THAT X+EPSLON .NE. ZERO
C*
C*    QUALIFICATIONS -
C*       THIS ROUTINE SHOULD PERFORM PROPERLY ON ALL SYSTEMS
C*       SATISFYING THE FOLLOWING TWO ASSUMPTIONS,
C*          1.  THE BASE USED IN REPRESENTING FLOATING POINT
C*              NUMBERS IS NOT A POWER OF THREE.
C*          2.  THE QUANTITY  A  IN STATEMENT 10 IS REPRESENTED TO
C*              THE ACCURACY USED IN FLOATING POINT VARIABLES
C*              THAT ARE STORED IN MEMORY.
C*       THE STATEMENT NUMBER 10 AND THE GO TO 10 ARE INTENDED TO
C*       FORCE OPTIMIZING COMPILERS TO GENERATE CODE SATISFYING
C*       ASSUMPTION 2.
C*       UNDER THESE ASSUMPTIONS, IT SHOULD BE TRUE THAT,
C*              A  IS NOT EXACTLY EQUAL TO FOUR-THIRDS,
C*              B  HAS A ZERO FOR ITS LAST BIT OR DIGIT,
C*              C  IS NOT EXACTLY EQUAL TO ONE,
C*              EPS  MEASURES THE SEPARATION OF 1.0 FROM
C*                   THE NEXT LARGER FLOATING POINT NUMBER.
C*       THE DEVELOPERS OF EISPACK WOULD APPRECIATE BEING INFORMED
C*       ABOUT ANY SYSTEMS WHERE THESE ASSUMPTIONS DO NOT HOLD.
C*
C*    DIFFERENCES FROM EISPACK 3 -
C*       USE IS MADE OF PARAMETER STATEMENTS AND INTRINSIC FUNCTIONS
C*       --NO EXECUTEABLE CODE CHANGES--
C*
C*    NOTE -
C*       QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
C*       B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
C
      DOUBLE PRECISION A,B,C,EPS,X
      DOUBLE PRECISION ZERO, ONE, THREE, FOUR
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, THREE=3.0D+00, FOUR=4.0D+00)
C
C-----------------------------------------------------------------------
C
      A = FOUR/THREE
   10 B = A - ONE
      C = B + B + B
      EPS = ABS(C - ONE)
      IF (EPS .EQ. ZERO) GO TO 10
      EPSLON = EPS*ABS(X)
      RETURN
      END
C*MODULE EIGEN   *DECK EQLRAT
      SUBROUTINE EQLRAT(N,DIAG,E,E2IN,D,IND,IERR,E2)
C*
C*    AUTHORS -
C*       THIS IS A MODIFICATION OF ROUTINE EQLRAT FROM EISPACK EDITION 3
C*       DATED AUGUST 1983.
C*       TQLRAT IS A TRANSLATION OF THE ALGOL PROCEDURE TQLRAT,
C*       ALGORITHM 464, COMM. ACM 16, 689(1973) BY REINSCH.
C*       THIS VERSION IS BY S. T. ELBERT (AMES LABORATORY-USDOE)
C*
C*    PURPOSE -
C*       THIS ROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC
C*       TRIDIAGONAL MATRIX
C*
C*    METHOD -
C*       RATIONAL QL
C*
C*    ON ENTRY -
C*       N      - INTEGER
C*                THE ORDER OF THE MATRIX.
C*       D      - W.P. REAL (N)
C*                CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
C*       E2     - W.P. REAL (N)
C*                CONTAINS THE SQUARES OF THE SUBDIAGONAL ELEMENTS OF
C*                THE INPUT MATRIX IN ITS LAST N-1 POSITIONS.
C*                E2(1) IS ARBITRARY.
C*
C*     ON EXIT -
C*       D      - W.P. REAL (N)
C*                CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
C*                ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND
C*                ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE
C*                THE SMALLEST EIGENVALUES.
C*       E2     - W.P. REAL (N)
C*                DESTROYED.
C*       IERR   - INTEGER
C*                SET TO
C*                ZERO       FOR NORMAL RETURN,
C*                J          IF THE J-TH EIGENVALUE HAS NOT BEEN
C*                           DETERMINED AFTER 30 ITERATIONS.
C*
C*    DIFFERENCES FROM EISPACK 3 -
C*       G=G+B INSTEAD OF IF(G.EQ.0) G=B ; B=B/4
C*       F77 BACKWARD LOOPS INSTEAD OF F66 CONSTRUCT
C*       GENERIC INTRINSIC FUNCTIONS
C*       ARRARY  IND  ADDED FOR USE BY EINVIT
C*
C*    NOTE -
C*       QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
C*       B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
C
      INTEGER I,J,L,M,N,II,L1,IERR
      INTEGER IND(N)
C
      DOUBLE PRECISION D(N),E(N),E2(N),DIAG(N),E2IN(N)
      DOUBLE PRECISION B,C,F,G,H,P,R,S,T,EPSLON
      DOUBLE PRECISION SCALE,ZERO,ONE
C
      PARAMETER (ZERO = 0.0D+00, SCALE= 1.0D+00/64.0D+00, ONE = 1.0D+00)
C
C-----------------------------------------------------------------------
      IERR = 0
      D(1)=DIAG(1)
      IND(1) = 1
      K = 0
      ITAG = 0
      IF (N .EQ. 1) GO TO 1001
C
      DO 100 I = 2, N
         D(I)=DIAG(I)
  100 E2(I-1) = E2IN(I)
C
      F = ZERO
      T = ZERO
      B = EPSLON(ONE)
      C = B *B
      B = B * SCALE
      E2(N) = ZERO
C
      DO 290 L = 1, N
         H = ABS(D(L)) + ABS(E(L))
         IF (T .GE. H) GO TO 105
            T = H
            B = EPSLON(T)
            C = B * B
            B = B * SCALE
  105    CONTINUE
C     .......... LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT ..........
         M = L - 1
  110    M = M + 1
         IF (E2(M) .GT. C) GO TO 110
C     .......... E2(N) IS ALWAYS ZERO, SO THERE IS AN EXIT
C                FROM THE LOOP ..........
C
         IF (M .LE. K) GO TO 125
            IF (M .NE. N) E2IN(M+1) = ZERO
            K = M
            ITAG = ITAG + 1
  125    CONTINUE
         IF (M .EQ. L) GO TO 210
C
C           ITERATE
C
         DO 205 J = 1, 30
C              .......... FORM SHIFT ..........
            L1 = L + 1
            S = SQRT(E2(L))
            G = D(L)
            P = (D(L1) - G) / (2.0D+00 * S)
            R = SQRT(P*P+1.0D+00)
            D(L) = S / (P + SIGN(R,P))
            H = G - D(L)
C
            DO 140 I = L1, N
  140       D(I) = D(I) - H
C
            F = F + H
C              .......... RATIONAL QL TRANSFORMATION ..........
            G = D(M) + B
            H = G
            S = ZERO
            DO 200 I = M-1,L,-1
               P = G * H
               R = P + E2(I)
               E2(I+1) = S * R
               S = E2(I) / R
               D(I+1) = H + S * (H + D(I))
               G = D(I) - E2(I) / G   + B
               H = G * P / R
  200       CONTINUE
C
            E2(L) = S * G
            D(L) = H
C              .......... GUARD AGAINST UNDERFLOW IN CONVERGENCE TEST
            IF (H .EQ. ZERO) GO TO 210
            IF (ABS(E2(L)) .LE. ABS(C/H)) GO TO 210
            E2(L) = H * E2(L)
            IF (E2(L) .EQ. ZERO) GO TO 210
  205    CONTINUE
C     .......... SET ERROR -- NO CONVERGENCE TO AN
C                EIGENVALUE AFTER 30 ITERATIONS ..........
      IERR = L
      GO TO 1001
C
C           CONVERGED
C
  210    P = D(L) + F
C           .......... ORDER EIGENVALUES ..........
         I = 1
         IF (L .EQ. 1) GO TO 250
            IF (P .LT. D(1)) GO TO 230
               I = L
C           .......... LOOP TO FIND ORDERED POSITION
  220          I = I - 1
               IF (P .LT. D(I)) GO TO 220
C
               I = I + 1
               IF (I .EQ. L) GO TO 250
  230       CONTINUE
            DO 240 II = L, I+1, -1
               D(II) = D(II-1)
               IND(II) = IND(II-1)
  240       CONTINUE
C
  250    CONTINUE
         D(I) = P
         IND(I) = ITAG
  290 CONTINUE
C
 1001 RETURN
      END
C*MODULE EIGEN   *DECK ESTPI
      SUBROUTINE ESTPI(ESTPI1,N,EVAL,D,E,X,ANORM)
C*
C*    AUTHOR -
C*       STEPHEN T. ELBERT (AMES LABORATORY-USDOE) DATE: 5 DEC 1986
C*
C*    PURPOSE -
C*       EVALUATE SYMMETRIC TRIDIAGONAL MATRIX PERFORMANCE INDEX
C*       *        *         *                  *           *
C*       FOR 1 EIGENVECTOR
C*           *
C*
C*    METHOD -
C*       THIS ROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX A*X-X*EVAL
C*       WHERE  A  IS A SYMMETRIC TRIDIAGONAL MATRIX STORED
C*       IN THE DIAGONAL (D) AND SUB-DIAGONAL (E) VECTORS, EVAL IS THE
C*       EIGENVALUE OF AN EIGENVECTOR OF  A,  NAMELY  X.
C*       THIS NORM IS SCALED BY MACHINE ACCURACY FOR THE PROBLEM SIZE.
C*       ALL NORMS APPEARING IN THE COMMENTS BELOW ARE 1-NORMS.
C*
C*    ON ENTRY -
C*       N      - INTEGER
C*                THE ORDER OF THE MATRIX  A.
C*       EVAL   - W.P. REAL
C*                THE EIGENVALUE CORRESPONDING TO VECTOR  X.
C*       D      - W.P. REAL (N)
C*                THE DIAGONAL VECTOR OF  A.
C*       E      - W.P. REAL (N)
C*                THE SUB-DIAGONAL VECTOR OF  A.
C*       X      - W.P. REAL (N)
C*                AN EIGENVECTOR OF  A.
C*       ANORM  - W.P. REAL
C*                THE NORM OF  A  IF IT HAS BEEN PREVIOUSLY COMPUTED.
C*
C*    ON EXIT -
C*       ANORM  - W.P. REAL
C*                THE NORM OF  A, COMPUTED IF INITIALLY ZERO.
C*       ESTPI1 - W.P. REAL
C*          !!A*X-X*EVAL!! / (EPSLON(10*N)*!!A!!*!!X!!);
C*          WHERE EPSLON(X) IS THE SMALLEST NUMBER SUCH THAT
C*             X + EPSLON(X) .NE. X
C*
C*          ESTPI1 .LT. 1 == SATISFACTORY PERFORMANCE
C*                 .GE. 1 AND .LE. 100 == MARGINAL PERFORMANCE
C*                 .GT. 100 == POOR PERFORMANCE
C*          (SEE LECT. NOTES IN COMP. SCI. VOL.6 PP 124-125)
C
      DOUBLE PRECISION ESTPI1,ANORM,EVAL,RNORM,SIZE,XNORM
      DOUBLE PRECISION D(N), E(N), X(N)
      DOUBLE PRECISION EPSLON, ONE, ZERO
C
      PARAMETER (ZERO = 0.0D+00, ONE = 1.0D+00)
C
C-----------------------------------------------------------------------
C
      ESTPI1 = ZERO
      IF( N .LE. 1 ) RETURN
      SIZE = 10 * N
      IF (ANORM .EQ. ZERO) THEN
C
C              COMPUTE NORM OF  A
C
         ANORM = MAX( ABS(D(1)) + ABS(E(2))
     *               ,ABS(D(N)) + ABS(E(N)))
         DO 110 I = 2, N-1
            ANORM = MAX( ANORM, ABS(E(I))+ABS(D(I))+ABS(E(I+1)))
  110    CONTINUE
         IF(ANORM .EQ. ZERO) ANORM = ONE
      END IF
C
C           COMPUTE NORMS OF RESIDUAL AND EIGENVECTOR
C
      XNORM = ABS(X(1)) + ABS(X(N))
      RNORM = ABS( (D(1)-EVAL)*X(1) + E(2)*X(2))
     *       +ABS( (D(N)-EVAL)*X(N) + E(N)*X(N-1))
      DO 120 I = 2, N-1
         XNORM = XNORM + ABS(X(I))
         RNORM = RNORM + ABS(E(I)*X(I-1) + (D(I)-EVAL)*X(I)
     *                       + E(I+1)*X(I+1))
  120 CONTINUE
C
      ESTPI1 = RNORM / (EPSLON(SIZE)*ANORM*XNORM)
      RETURN
      END
C*MODULE EIGEN   *DECK ETRBK3
      SUBROUTINE ETRBK3(NM,N,NV,A,M,Z)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C*
C*    AUTHORS-
C*       THIS IS A MODIFICATION OF ROUTINE TRBAK3 FROM EISPACK EDITION 3
C*       DATED AUGUST 1983.
C*       EISPACK TRBAK3 IS A TRANSLATION OF THE ALGOL PROCEDURE TRBAK3,
C*       NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
C*       HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
C*       THIS VERSION IS BY S. T. ELBERT (AMES LABORATORY-USDOE)
C*
C*    PURPOSE -
C*       THIS ROUTINE FORMS THE EIGENVECTORS OF A REAL SYMMETRIC
C*       MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
C*       SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY  ETRED3.
C*
C*    METHOD -
C*       THE CALCULATION IS CARRIED OUT BY FORMING THE MATRIX PRODUCT
C*          Q*Z
C*       WHERE  Q  IS A PRODUCT OF THE ORTHOGONAL SYMMETRIC MATRICES
C*                Q = PROD(I)[1 - U(I)*.TRANSPOSE.U(I)*H(I)]
C*       U  IS THE AUGMENTED SUB-DIAGONAL ROWS OF  A  AND
C*       Z  IS THE SET OF EIGENVECTORS OF THE TRIDIAGONAL
C*       MATRIX  F  WHICH WAS FORMED FROM THE ORIGINAL SYMMETRIC
C*       MATRIX  C  BY THE SIMILARITY TRANSFORMATION
C*                F = Q(TRANSPOSE) C Q
C*       NOTE THAT ETRBK3 PRESERVES VECTOR EUCLIDEAN NORMS.
C*
C*
C*    COMPLEXITY -
C*       M*N**2
C*
C*    ON ENTRY-
C*       NM     - INTEGER
C*                MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C*                ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
C*                DIMENSION STATEMENT.
C*       N      - INTEGER
C*                THE ORDER OF THE MATRIX  A.
C*       NV     - INTEGER
C*                MUST BE SET TO THE DIMENSION OF THE ARRAY  A  AS
C*                DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT.
C*       A      - W.P. REAL (NV)
C*                CONTAINS INFORMATION ABOUT THE ORTHOGONAL
C*                TRANSFORMATIONS USED IN THE REDUCTION BY  ETRED3  IN
C*                ITS FIRST  NV = N*(N+1)/2 POSITIONS.
C*       M      - INTEGER
C*                THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.
C*       Z      - W.P REAL (NM,M)
C*                CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
C*                IN ITS FIRST M COLUMNS.
C*
C*    ON EXIT-
C*       Z      - W.P. REAL (NM,M)
C*                CONTAINS THE TRANSFORMED EIGENVECTORS
C*                IN ITS FIRST M COLUMNS.
C*
C*    DIFFERENCES WITH EISPACK 3 -
C*       THE TWO INNER LOOPS ARE REPLACED BY DDOT AND DAXPY.
C*       MULTIPLICATION USED INSTEAD OF DIVISION TO FIND S.
C*       OUTER LOOP RANGE CHANGED FROM 2,N TO 3,N.
C*       ADDRESS POINTERS FOR  A  SIMPLIFIED.
C*
C*    NOTE -
C*       QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
C*       B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
C
      INTEGER I,II,IM1,IZ,J,M,N,NM,NV
C
      DOUBLE PRECISION A(NV),Z(NM,M)
      DOUBLE PRECISION H,S,ZERO
C
      PARAMETER (ZERO = 0.0D+00)
C
C-----------------------------------------------------------------------
C
      IF (M .EQ. 0) RETURN
      IF (N .LE. 2) RETURN
C
      II=3
      DO 140 I = 3, N
         IZ=II+1
         II=II+I
         H = A(II)
         IF (H .EQ. ZERO) GO TO 140
            IM1 = I - 1
            DO 130 J = 1, M
               S = -( DDOT(IM1,A(IZ),1,Z(1,J),1) * H) * H
               CALL DAXPY(IM1,S,A(IZ),1,Z(1,J),1)
  130       CONTINUE
  140 CONTINUE
      RETURN
      END
C*MODULE EIGEN   *DECK ETRED3
      SUBROUTINE ETRED3(N,NV,A,D,E,E2)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C*
C*    AUTHORS -
C*       THIS IS A MODIFICATION OF ROUTINE TRED3 FROM EISPACK EDITION 3
C*       DATED AUGUST 1983.
C*       EISPACK TRED3 IS A TRANSLATION OF THE ALGOL PROCEDURE TRED3,
C*       NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
C*       HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
C*       THIS VERSION IS BY S. T. ELBERT, AMES LABORATORY-USDOE JUN 1986
C*
C*    PURPOSE -
C*       THIS ROUTINE REDUCES A REAL SYMMETRIC (PACKED) MATRIX, STORED
C*       AS A ONE-DIMENSIONAL ARRAY, TO A SYMMETRIC TRIDIAGONAL MATRIX
C*       USING ORTHOGONAL SIMILARITY TRANSFORMATIONS, PRESERVING THE
C*       INFORMATION ABOUT THE TRANSFORMATIONS IN  A.
C*
C*    METHOD -
C*       THE TRIDIAGONAL REDUCTION IS PERFORMED IN THE FOLLOWING WAY.
C*       STARTING WITH J=N, THE ELEMENTS IN THE J-TH ROW TO THE
C*       LEFT OF THE DIAGONAL ARE FIRST SCALED, TO AVOID POSSIBLE
C*       UNDERFLOW IN THE TRANSFORMATION THAT MIGHT RESULT IN SEVERE
C*       DEPARTURE FROM ORTHOGONALITY.  THE SUM OF SQUARES  SIGMA  OF
C*       THESE SCALED ELEMENTS IS NEXT FORMED.  THEN, A VECTOR  U  AND
C*       A SCALAR
C*                      H = U(TRANSPOSE) * U / 2
C*       DEFINE A REFLECTION OPERATOR
C*                      P = I - U * U(TRANSPOSE) / H
C*       WHICH IS ORTHOGONAL AND SYMMETRIC AND FOR WHICH THE
C*       SIMILIARITY TRANSFORMATION  PAP  ELIMINATES THE ELEMENTS IN
C*       THE J-TH ROW OF  A  TO THE LEFT OF THE SUBDIAGONAL AND THE
C*       SYMMETRICAL ELEMENTS IN THE J-TH COLUMN.
C*
C*       THE NON-ZERO COMPONENTS OF  U  ARE THE ELEMENTS OF THE J-TH
C*       ROW TO THE LEFT OF THE DIAGONAL WITH THE LAST OF THEM
C*       AUGMENTED BY THE SQUARE ROOT OF  SIGMA  PREFIXED BY THE SIGN
C*       OF THE SUBDIAGONAL ELEMENT.  BY STORING THE TRANSFORMED SUB-
C*       DIAGONAL ELEMENT IN  E(J)  AND NOT OVERWRITING THE ROW
C*       ELEMENTS ELIMINATED IN THE TRANSFORMATION, FULL INFORMATION
C*       ABOUT  P  IS SAVE FOR LATER USE IN  ETRBK3.
C*
C*       THE TRANSFORMATION SETS  E2(J)  EQUAL TO  SIGMA  AND  E(J)
C*       EQUAL TO THE SQUARE ROOT OF  SIGMA  PREFIXED BY THE SIGN
C*       OF THE REPLACED SUBDIAGONAL ELEMENT.
C*
C*       THE ABOVE STEPS ARE REPEATED ON FURTHER ROWS OF THE
C*       TRANSFORMED  A  IN REVERSE ORDER UNTIL  A  IS REDUCED TO TRI-
C*       DIAGONAL FORM, THAT IS, REPEATED FOR  J = N-1,N-2,...,3.
C*
C*    COMPLEXITY -
C*       2/3 N**3
C*
C*    ON ENTRY-
C*       N      - INTEGER
C*                THE ORDER OF THE MATRIX.
C*       NV     - INTEGER
C*                MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A
C*                AS DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT
C*       A      - W.P. REAL (NV)
C*                CONTAINS THE LOWER TRIANGLE OF THE REAL SYMMETRIC
C*                INPUT MATRIX, STORED ROW-WISE AS A ONE-DIMENSIONAL
C*                ARRAY, IN ITS FIRST N*(N+1)/2 POSITIONS.
C*
C*    ON EXIT-
C*       A      - W.P. REAL (NV)
C*                CONTAINS INFORMATION ABOUT THE ORTHOGONAL
C*                TRANSFORMATIONS USED IN THE REDUCTION.
C*       D      - W.P. REAL (N)
C*                CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL
C*                MATRIX.
C*       E      - W.P. REAL (N)
C*                CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
C*                MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO
C*       E2     - W.P. REAL (N)
C*                CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF
C*                E. MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
C*
C*    DIFFERENCES FROM EISPACK 3 -
C*       OUTER LOOP CHANGED FROM II=1,N TO I=N,3,-1
C*       PARAMETER STATEMENT AND GENERIC INTRINSIC FUNCTIONS USED
C*       SCALE.NE.0 TEST NOW SPOTS TRI-DIAGONAL FORM
C*       VALUES LESS THAN EPSLON CLEARED TO ZERO
C*       USE BLAS(1)
C*       U NOT COPIED TO D, LEFT IN A
C*       E2 COMPUTED FROM E
C*       INNER LOOPS SPLIT INTO ROUTINES ELAU AND FREDA
C*       INVERSE OF H STORED INSTEAD OF H
C*
C*    NOTE -
C*       QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
C*       B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
C
      INTEGER I,IIA,IZ0,L,N,NV
C
      DOUBLE PRECISION A(NV),D(N),E(N),E2(N)
      DOUBLE PRECISION AIIMAX,F,G,H,HROOT,SCALE,SCALEI
      DOUBLE PRECISION ONE, ZERO
C
      PARAMETER (ZERO = 0.0D+00, ONE = 1.0D+00)
C
C-----------------------------------------------------------------------
C
      IF (N .LE. 2) GO TO 310
      IZ0 = (N*N+N)/2
      AIIMAX = ABS(A(IZ0))
      DO 300 I = N, 3, -1
         L = I - 1
         IIA = IZ0
         IZ0 = IZ0 - I
         AIIMAX = MAX(AIIMAX, ABS(A(IIA)))
         SCALE = DASUM(L, A(IZ0+1), 1)
         IF(SCALE .EQ. ABS(A(IIA-1)) .OR. AIIMAX+SCALE .EQ. AIIMAX) THEN
C
C           THIS ROW IS ALREADY IN TRI-DIAGONAL FORM
C
            D(I) = A(IIA)
            IF (AIIMAX+D(I) .EQ. AIIMAX) D(I) = ZERO
            E(I) = A(IIA-1)
            IF (AIIMAX+E(I) .EQ. AIIMAX) E(I) = ZERO
            E2(I) = E(I)*E(I)
            A(IIA) = ZERO
            GO TO 300
C
         END IF
C
         SCALEI = ONE / SCALE
         CALL DSCAL(L,SCALEI,A(IZ0+1),1)
         HROOT = DNRM2(L,A(IZ0+1),1)
C
         F = A(IZ0+L)
         G = -SIGN(HROOT,F)
         E(I) = SCALE * G
         E2(I) = E(I)*E(I)
         H = HROOT*HROOT - F * G
         A(IZ0+L) = F - G
         D(I) = A(IIA)
         A(IIA) = ONE / SQRT(H)
C           .......... FORM P THEN Q IN E(1:L) ..........
         CALL ELAU(ONE/H,L,A(IZ0+1),A,E)
C           .......... FORM REDUCED A ..........
         CALL FREDA(L,A(IZ0+1),A,E)
C
  300 CONTINUE
  310 CONTINUE
      E(1) = ZERO
      E2(1)= ZERO
      D(1) = A(1)
      IF(N.EQ.1) RETURN
C
      E(2) = A(2)
      E2(2)= A(2)*A(2)
      D(2) = A(3)
      RETURN
      END
C*MODULE EIGEN   *DECK EVVRSP
      SUBROUTINE EVVRSP(MSGFL,N,NVECT,LENA,NV,A,B,IND,ROOT,
     *                  VECT,IORDER,IERR)
C*
C*    AUTHOR:  S. T. ELBERT, AMES LABORATORY-USDOE, JUNE 1985
C*
C*    PURPOSE -
C*       FINDS   (ALL) EIGENVALUES    AND    (SOME OR ALL) EIGENVECTORS
C*                     *    *                                   *
C*       OF A REAL SYMMETRIC PACKED MATRIX.
C*            *    *         *
C*
C*    METHOD -
C*       THE METHOD AS PRESENTED IN THIS ROUTINE CONSISTS OF FOUR STEPS:
C*       FIRST, THE INPUT MATRIX IS REDUCED TO TRIDIAGONAL FORM BY THE
C*       HOUSEHOLDER TECHNIQUE (ORTHOGONAL SIMILARITY TRANSFORMATIONS).
C*       SECOND, THE ROOTS ARE LOCATED USING THE RATIONAL QL METHOD.
C*       THIRD, THE VECTORS OF THE TRIDIAGONAL FORM ARE EVALUATED BY THE
C*       INVERSE ITERATION TECHNIQUE.  VECTORS FOR DEGENERATE OR NEAR-
C*       DEGENERATE ROOTS ARE FORCED TO BE ORTHOGONAL.
C*       FOURTH, THE TRIDIAGONAL VECTORS ARE ROTATED TO VECTORS OF THE
C*       ORIGINAL ARRAY.
C*
C*       THESE ROUTINES ARE MODIFICATIONS OF THE EISPACK 3
C*       ROUTINES TRED3, TQLRAT, TINVIT AND TRBAK3
C*
C*       FOR FURTHER DETAILS, SEE EISPACK USERS GUIDE, B. T. SMITH
C*       ET AL, SPRINGER-VERLAG, LECTURE NOTES IN COMPUTER SCIENCE,
C*       VOL. 6, 2-ND EDITION, 1976.  ANOTHER GOOD REFERENCE IS
C*       THE SYMMETRIC EIGENVALUE PROBLEM BY B. N. PARLETT
C*       PUBLISHED BY PRENTICE-HALL, INC., ENGLEWOOD CLIFFS, N.J. (1980)
C*
C*    ON ENTRY -
C*       MSGFL  - INTEGER (LOGICAL UNIT NO.)
C*                FILE WHERE ERROR MESSAGES WILL BE PRINTED.
C*                IF MSGFL IS 0, ERROR MESSAGES WILL BE PRINTED ON LU 6.
C*                IF MSGFL IS NEGATIVE, NO ERROR MESSAGES PRINTED.
C*       N      - INTEGER
C*                ORDER OF MATRIX A.
C*       NVECT  - INTEGER
C*                NUMBER OF VECTORS DESIRED.  0 .LE. NVECT .LE. N.
C*       LENA   - INTEGER
C*                DIMENSION OF  A  IN CALLING ROUTINE.  MUST NOT BE LESS
C*                THAN (N*N+N)/2.
C*       NV     - INTEGER
C*                ROW DIMENSION OF VECT IN CALLING ROUTINE.   N .LE. NV.
C*       A      - WORKING PRECISION REAL (LENA)
C*                INPUT MATRIX, ROWS OF THE LOWER TRIANGLE PACKED INTO
C*                LINEAR ARRAY OF DIMENSION N*(N+1)/2.  THE PACKED ORDER
C*                IS A(1,1), A(2,1), A(2,2), A(3,1), A(3,2), ...
C*       B      - WORKING PRECISION REAL (N,8)
C*                SCRATCH ARRAY, 8*N ELEMENTS
C*       IND    - INTEGER (N)
C*                SCRATCH ARRAY OF LENGTH N.
C*       IORDER - INTEGER
C*                ROOT ORDERING FLAG.
C*                = 0, ROOTS WILL BE PUT IN ASCENDING ORDER.
C*                = 2, ROOTS WILL BE PUT IN DESCENDING ORDER.
C*
C*    ON EXIT -
C*       A      - DESTORYED.  NOW HOLDS REFLECTION OPERATORS.
C*       ROOT   - WORKING PRECISION REAL (N)
C*                ALL EIGENVALUES IN ASCENDING OR DESCENDING ORDER.
C*                  IF IORDER = 0, ROOT(1) .LE. ... .LE. ROOT(N)
C*                  IF IORDER = 2, ROOT(1) .GE. ... .GE. ROOT(N)
C*       VECT   - WORKING PRECISION REAL (NV,NVECT)
C*                EIGENVECTORS FOR ROOT(1), ..., ROOT(NVECT).
C*       IERR   - INTEGER
C*                = 0 IF NO ERROR DETECTED,
C*                = K IF ITERATION FOR K-TH EIGENVALUE FAILED,
C*                = -K IF ITERATION FOR K-TH EIGENVECTOR FAILED.
C*                (FAILURES SHOULD BE VERY RARE.  CONTACT C. MOLER.)
C*
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DOUBLE PRECISION A(LENA)
      DOUBLE PRECISION B(N,8)
      DOUBLE PRECISION ROOT(N)
      DOUBLE PRECISION T
      DOUBLE PRECISION VECT(NV,*)
C
      INTEGER IND(N)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
  900 FORMAT(26H0*** EVVRSP PARAMETERS ***/
     +       14H ***      N = ,I8,4H ***/
     +       14H ***  NVECT = ,I8,4H ***/
     +       14H ***   LENA = ,I8,4H ***/
     +       14H ***     NV = ,I8,4H ***/
     +       14H *** IORDER = ,I8,4H ***/
     +       14H ***   IERR = ,I8,4H ***)
  901 FORMAT(37H VALUE OF LENA IS LESS THAN (N*N+N)/2)
  902 FORMAT(39H EQLRAT HAS FAILED TO CONVERGE FOR ROOT,I5)
  903 FORMAT(18H NV IS LESS THAN N)
  904 FORMAT(41H EINVIT HAS FAILED TO CONVERGE FOR VECTOR,I5)
  905 FORMAT(51H VALUE OF IORDER MUST BE 0 (SMALLEST ROOT FIRST) OR
     *      ,23H 2 (LARGEST ROOT FIRST))
  906 FORMAT(' VALUE OF N IS LESS THAN OR EQUAL ZERO')
C
C-----------------------------------------------------------------------
C
      LMSGFL=MSGFL
      IF (MSGFL .EQ. 0) LMSGFL=6
      IERR = N - 1
      IF (N .LE. 0) GO TO 800
      IERR = N + 1
      IF ( (N*N+N)/2 .GT. LENA) GO TO 810
C
C        REDUCE REAL SYMMETRIC MATRIX A TO TRIDIAGONAL FORM
C
      CALL ETRED3(N,LENA,A,B(1,1),B(1,2),B(1,3))
C
C        FIND ALL EIGENVALUES OF TRIDIAGONAL MATRIX
C
      CALL EQLRAT(N,B(1,1),B(1,2),B(1,3),ROOT,IND,IERR,B(1,4))
      IF (IERR .NE. 0) GO TO 820
C
C         CHECK THE DESIRED ORDER OF THE EIGENVALUES
C
      B(1,3) = IORDER
      IF (IORDER .EQ. 0) GO TO 300
         IF (IORDER .NE. 2) GO TO 850
C
C         ORDER ROOTS IN DESCENDING ORDER (LARGEST FIRST)...
C        TURN ROOT AND IND ARRAYS END FOR END
C
         DO 210 I = 1, N/2
            J = N+1-I
            T = ROOT(I)
            ROOT(I) = ROOT(J)
            ROOT(J) = T
            L = IND(I)
            IND(I) = IND(J)
            IND(J) = L
  210    CONTINUE
C
C           FIND I AND J MARKING THE START AND END OF A SEQUENCE
C           OF DEGENERATE ROOTS
C
         I=0
  220    CONTINUE
            I = I+1
            IF (I .GT. N) GO TO 300
            DO 230 J=I,N
               IF (ROOT(J) .NE. ROOT(I)) GO TO 240
  230       CONTINUE
            J = N+1
  240       CONTINUE
            J = J-1
            IF (J .EQ. I) GO TO 220
C
C                    TURN AROUND IND BETWEEN I AND J
C
            JSV = J
            KLIM = (J-I+1)/2
            DO 250 K=1,KLIM
               L = IND(J)
               IND(J) = IND(I)
               IND(I) = L
               I = I+1
               J = J-1
  250       CONTINUE
            I = JSV
         GO TO 220
C
  300 CONTINUE
C
      IF (NVECT .LE. 0) RETURN
      IF (NV .LT. N) GO TO 830
C
C        FIND EIGENVECTORS OF TRI-DIAGONAL MATRIX VIA INVERSE ITERATION
C
      IERR = LMSGFL
      CALL EINVIT(NV,N,B(1,1),B(1,2),B(1,3),NVECT,ROOT,IND,
     +            VECT,IERR,B(1,4),B(1,5),B(1,6),B(1,7),B(1,8))
      IF (IERR .NE. 0) GO TO 840
C
C        FIND EIGENVECTORS OF SYMMETRIC MATRIX VIA BACK TRANSFORMATION
C
  400 CONTINUE
      CALL ETRBK3(NV,N,LENA,A,NVECT,VECT)
      RETURN
C
C        ERROR MESSAGE SECTION
C
  800 IF (LMSGFL .LT. 0) RETURN
      IF (MASWRK) WRITE(LMSGFL,906)
      GO TO 890
C
  810 IF (LMSGFL .LT. 0) RETURN
      IF (MASWRK) WRITE(LMSGFL,901)
      GO TO 890
C
  820 IF (LMSGFL .LT. 0) RETURN
      IF (MASWRK) WRITE(LMSGFL,902) IERR
      GO TO 890
C
  830 IF (LMSGFL .LT. 0) RETURN
      IF (MASWRK) WRITE(LMSGFL,903)
      GO TO 890
C
  840 CONTINUE
      IF ((LMSGFL .GT. 0).AND.MASWRK) WRITE(LMSGFL,904) -IERR
      GO TO 400
C
  850 IERR=-1
      IF (LMSGFL .LT. 0) RETURN
      IF (MASWRK) WRITE(LMSGFL,905)
      GO TO 890
C
  890 CONTINUE
      IF (MASWRK) WRITE(LMSGFL,900) N,NVECT,LENA,NV,IORDER,IERR
      RETURN
      END
C*MODULE EIGEN   *DECK FREDA
      SUBROUTINE FREDA(L,D,A,E)
C
      DOUBLE PRECISION A(*)
      DOUBLE PRECISION D(L)
      DOUBLE PRECISION E(L)
      DOUBLE PRECISION F
      DOUBLE PRECISION G
C
      JK = 1
C
C     .......... FORM REDUCED A ..........
C
      DO 280 J = 1, L
         F = D(J)
         G = E(J)
C
         DO 260 K = 1, J
            A(JK) = A(JK) - F * E(K) - G * D(K)
            JK = JK + 1
  260    CONTINUE
C
  280 CONTINUE
      RETURN
      END
C*MODULE EIGEN   *DECK GIVEIS
      SUBROUTINE GIVEIS(N,NVECT,NV,A,B,INDB,ROOT,VECT,IERR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(*),B(N,8),INDB(N),ROOT(N),VECT(NV,NVECT)
C
C     EISPACK-BASED SUBSTITUTE FOR QCPE ROUTINE GIVENS.
C     FINDS ALL EIGENVALUES AND SOME EIGENVECTORS OF A REAL SYMMETRIC
C     MATRIX.   AUTHOR.. C. MOLER AND D. SPANGLER, N.R.C.C., 4/1/79.
C
C     INPUT..
C     N     = ORDER OF MATRIX .
C     NVECT = NUMBER OF VECTORS DESIRED.  0 .LE. NVECT .LE. N .
C     NV    = LEADING DIMENSION OF VECT .
C     A     = INPUT MATRIX, COLUMNS OF THE UPPER TRIANGLE PACKED INTO
C             LINEAR ARRAY OF DIMENSION N*(N+1)/2 .
C     B     = SCRATCH ARRAY, 8*N ELEMENTS (NOTE THIS IS MORE THAN
C             PREVIOUS VERSIONS OF GIVENS.)
C    IND    = INDEX ARRAY OF N ELEMENTS
C
C     OUTPUT..
C     A       DESTROYED .
C     ROOT  = ALL EIGENVALUES, ROOT(1) .LE. ... .LE. ROOT(N) .
C             (FOR OTHER ORDERINGS, SEE BELOW.)
C     VECT  = EIGENVECTORS FOR ROOT(1),..., ROOT(NVECT) .
C     IERR  = 0 IF NO ERROR DETECTED,
C           = K IF ITERATION FOR K-TH EIGENVALUE FAILED,
C           = -K IF ITERATION FOR K-TH EIGENVECTOR FAILED.
C             (FAILURES SHOULD BE VERY RARE.  CONTACT MOLER.)
C
C     CALLS MODIFIED EISPACK ROUTINES TRED3B, IMTQLV, TINVTB, AND
C     TRBK3B.  THE ROUTINES TRED3B, TINVTB, AND TRBK3B.
C     THE ORIGINAL EISPACK ROUTINES TRED3, TINVIT, AND TRBAK3
C     WERE MODIFIED BY THE INTRODUCTION OF TWO ROUTINES FROM THE
C     BLAS LIBRARY - DDOT AND DAXPY.
C
C         IF TINVIT FAILS TO CONVERGE, TQL2 IS CALLED
C
C         SEE EISPACK USERS GUIDE, B. T. SMITH ET AL, SPRINGER-VERLAG
C     LECTURE NOTES IN COMPUTER SCIENCE, VOL. 6, 2-ND EDITION, 1976 .
C     NOTE THAT IMTQLV AND TINVTB HAVE INTERNAL MACHINE
C     DEPENDENT CONSTANTS.
C
      DATA ONE, ZERO /1.0D+00, 0.0D+00/
      CALL TRED3B(N,(N*N+N)/2,A,B(1,1),B(1,2),B(1,3))
      CALL IMTQLV(N,B(1,1),B(1,2),B(1,3),ROOT,INDB,IERR,B(1,4))
      IF (IERR .NE. 0) RETURN
C
C     TO REORDER ROOTS...
C     K = N/2
C     B(1,3) = 2.0D+00
C     DO 50 I = 1, K
C        J = N+1-I
C        T = ROOT(I)
C        ROOT(I) = ROOT(J)
C        ROOT(J) = T
C 50  CONTINUE
C
      IF (NVECT .LE. 0) RETURN
      CALL TINVTB(NV,N,B(1,1),B(1,2),B(1,3),NVECT,ROOT,INDB,VECT,IERR,
     +     B(1,4),B(1,5),B(1,6),B(1,7),B(1,8))
      IF (IERR .EQ. 0) GO TO 160
C
C      IF INVERSE ITERATION GIVES AN ERROR IN DETERMINING THE
C      EIGENVECTORS, TRY THE QL ALGORITHM IF ALL THE EIGENVECTORS
C      ARE DESIRED.
C
      IF (NVECT .NE. N) RETURN
      DO 120 I = 1, NVECT
      DO 100 J = 1, N
      VECT(I,J) = ZERO
  100 CONTINUE
      VECT(I,I) = ONE
  120 CONTINUE
      CALL TQL2 (NV,N,B(1,1),B(1,2),VECT,IERR)
      DO 140 I = 1, NVECT
      ROOT(I) = B(I,1)
  140 CONTINUE
      IF (IERR .NE. 0) RETURN
  160 CALL TRBK3B(NV,N,(N*N+N)/2,A,NVECT,VECT)
      RETURN
      END
C*MODULE EIGEN   *DECK GLDIAG
      SUBROUTINE GLDIAG(LDVECT,NVECT,N,H,WRK,EIG,VECTOR,IERR,IWRK)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION H(*),WRK(N,8),EIG(N),VECTOR(LDVECT,NVECT),IWRK(N)
C
      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(950)
      COMMON /MACHSW/ KDIAGG,ICORFL,IXDR
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
C     ----- GENERAL ROUTINE TO DIAGONALIZE A SYMMETRIC MATRIX -----
C     IF KDIAG = 0, USE A ROUTINE FROM THE VECTOR LIBRARY,
C                   IF AVAILABLE (SEE THE SUBROUTINE 'GLDIAG'
C                   IN VECTOR.SRC), OR EVVRSP OTHERWISE
C              = 1, USE EVVRSP
C              = 2, USE GIVEIS
C              = 3, USE JACOBI
C     100 ADDED TO KDIAG FORCES DIAGONALISING ONLY ON MASTER AND
C     BROADCASTING THE RESULTS ON SLAVES (USEFUL FOR MIXED NODE CLUSTERS
C     AND/OR FOR CLUSTERS WITH CHARGES ACRUED BY CPU TIME).
C
C           N      = DIMENSION (ORDER) OF MATRIX TO BE SOLVED
C           LDVECT = LEADING DIMENSION OF VECTOR
C           NVECT  = NUMBER OF VECTORS DESIRED
C           H      = MATRIX TO BE DIAGONALIZED
C           WRK    = N*8 W.P. REAL WORDS OF SCRATCH SPACE
C           EIG    = EIGENVALUES  (OUTPUT)
C           VECTOR = EIGENVECTORS (OUTPUT)
C           IERR   = ERROR FLAG (OUTPUT)
C           IWRK   = N INTEGER WORDS OF SCRATCH SPACE
C
C        TRADITIONAL RUNS WITH -KDIAG- WILL HAVE LDIAG=0 AND KDIAG SAME.
C        ADDING 100 TO KDIAG BEFORE CALLING THIS WILL RESULT IN THE
C        DIAGONALIZATION OCCURING ENTIRELY ON THE MASTER NODE.
C
      KDIAG=MOD(KDIAGG,100)
      LDIAG=KDIAGG/100
C
      IF(LDIAG.EQ.0.OR.MASWRK) THEN
         IERR = 0
C
C         ----- USE STEVE ELBERT'S ROUTINE -----
C
         IF(KDIAG.LE.1 .OR. KDIAG.GT.3) THEN
            LENH = (N*N+N)/2
            KORDER =0
            CALL EVVRSP(IW,N,NVECT,LENH,LDVECT,H,WRK,IWRK,EIG,VECTOR
     *                 ,KORDER,IERR)
         END IF
C
C         ----- USE MODIFIED EISPAK ROUTINE -----
C
         IF(KDIAG.EQ.2)
     *      CALL GIVEIS(N,NVECT,LDVECT,H,WRK,IWRK,EIG,VECTOR,IERR)
C
C         ----- USE JACOBI ROTATION ROUTINE -----
C
         IF(KDIAG.EQ.3) THEN
            IF(NVECT.EQ.N) THEN
               CALL JACDG(H,VECTOR,EIG,IWRK,WRK,LDVECT,N)
            ELSE
               IF (MASWRK) WRITE(IW,9000) N,NVECT,LDVECT
               CALL ABRT
            END IF
         END IF
      ENDIF
C
C        BROADCAST THE RESULTS
C
      IF(LDIAG.NE.0.AND.GOPARR) THEN
         CALL DDI_BCAST(2421,'I',IERR,1,MASTER)
         CALL DDI_BCAST(2421,'F',EIG,N,MASTER)
         CALL DDI_BCAST(2421,'F',VECTOR,LDVECT*NVECT,MASTER)
      ENDIF
C
      RETURN
C
 9000 FORMAT(1X,'IN -GLDIAG-, N,NVECT,LDVECT=',3I8/
     *       1X,'THE JACOBI CODE CANNOT COPE WITH N.NE.NVECT!'/
     *       1X,'SO THIS RUN DOES NOT PERMIT KDIAG=3.')
      END
C*MODULE EIGEN   *DECK IMTQLV
      SUBROUTINE IMTQLV(N,D,E,E2,W,IND,IERR,RV1)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER TAG
      DOUBLE PRECISION MACHEP
      DIMENSION D(N),E(N),E2(N),W(N),RV1(N),IND(N)
C
C     THIS ROUTINE IS A VARIANT OF  IMTQL1  WHICH IS A TRANSLATION OF
C     ALGOL PROCEDURE IMTQL1, NUM. MATH. 12, 377-383(1968) BY MARTIN AND
C     WILKINSON, AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971).
C
C     THIS ROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC TRIDIAGONAL
C     MATRIX BY THE IMPLICIT QL METHOD AND ASSOCIATES WITH THEM
C     THEIR CORRESPONDING SUBMATRIX INDICES.
C
C     ON INPUT-
C
C        N IS THE ORDER OF THE MATRIX,
C
C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
C
C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY,
C
C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
C          E2(1) IS ARBITRARY.
C
C     ON OUTPUT-
C
C        D AND E ARE UNALTERED,
C
C        ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED
C          AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE
C          MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES.
C          E2(1) IS ALSO SET TO ZERO,
C
C        W CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
C          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND
C          ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE
C          THE SMALLEST EIGENVALUES,
C
C        IND CONTAINS THE SUBMATRIX INDICES ASSOCIATED WITH THE
C          CORRESPONDING EIGENVALUES IN W -- 1 FOR EIGENVALUES
C          BELONGING TO THE FIRST SUBMATRIX FROM THE TOP,
C          2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.,
C
C        IERR IS SET TO
C          ZERO       FOR NORMAL RETURN,
C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
C                     DETERMINED AFTER 30 ITERATIONS,
C
C        RV1 IS A TEMPORARY STORAGE ARRAY.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
C     ------------------------------------------------------------------
C
C     ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
C
C                **********
      MACHEP = 2.0D+00**(-50)
C
      IERR = 0
      K = 0
      TAG = 0
C
      DO 100 I = 1, N
      W(I) = D(I)
      IF (I .NE. 1) RV1(I-1) = E(I)
  100 CONTINUE
C
      E2(1) = 0.0D+00
      RV1(N) = 0.0D+00
C
      DO 360 L = 1, N
      J = 0
C     ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT **********
  120 DO 140 M = L, N
      IF (M .EQ. N) GO TO 160
      IF (ABS(RV1(M)) .LE. MACHEP * (ABS(W(M)) + ABS(W(M+1)))) GO TO
     +     160
C     ********** GUARD AGAINST UNDERFLOWED ELEMENT OF E2 **********
      IF (E2(M+1) .EQ. 0.0D+00) GO TO 180
  140 CONTINUE
C
  160 IF (M .LE. K) GO TO 200
      IF (M .NE. N) E2(M+1) = 0.0D+00
  180 K = M
      TAG = TAG + 1
  200 P = W(L)
      IF (M .EQ. L) GO TO 280
      IF (J .EQ. 30) GO TO 380
      J = J + 1
C     ********** FORM SHIFT **********
      G = (W(L+1) - P) / (2.0D+00 * RV1(L))
      R = SQRT(G*G+1.0D+00)
      G = W(M) - P + RV1(L) / (G + SIGN(R,G))
      S = 1.0D+00
      C = 1.0D+00
      P = 0.0D+00
      MML = M - L
C     ********** FOR I=M-1 STEP -1 UNTIL L DO -- **********
      DO 260 II = 1, MML
      I = M - II
      F = S * RV1(I)
      B = C * RV1(I)
      IF (ABS(F) .LT. ABS(G)) GO TO 220
      C = G / F
      R = SQRT(C*C+1.0D+00)
      RV1(I+1) = F * R
      S = 1.0D+00 / R
      C = C * S
      GO TO 240
  220 S = F / G
      R = SQRT(S*S+1.0D+00)
      RV1(I+1) = G * R
      C = 1.0D+00 / R
      S = S * C
  240 G = W(I+1) - P
      R = (W(I) - G) * S + 2.0D+00 * C * B
      P = S * R
      W(I+1) = G + P
      G = C * R - B
  260 CONTINUE
C
      W(L) = W(L) - P
      RV1(L) = G
      RV1(M) = 0.0D+00
      GO TO 120
C     ********** ORDER EIGENVALUES **********
  280 IF (L .EQ. 1) GO TO 320
C     ********** FOR I=L STEP -1 UNTIL 2 DO -- **********
      DO 300 II = 2, L
      I = L + 2 - II
      IF (P .GE. W(I-1)) GO TO 340
      W(I) = W(I-1)
      IND(I) = IND(I-1)
  300 CONTINUE
C
  320 I = 1
  340 W(I) = P
      IND(I) = TAG
  360 CONTINUE
C
      GO TO 400
C     ********** SET ERROR -- NO CONVERGENCE TO AN
C                EIGENVALUE AFTER 30 ITERATIONS **********
  380 IERR = L
  400 RETURN
C     ********** LAST CARD OF IMTQLV **********
      END
C*MODULE EIGEN   *DECK JACDG
      SUBROUTINE JACDG(A,VEC,EIG,JBIG,BIG,LDVEC,N)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION A(*),VEC(LDVEC,N),EIG(N),JBIG(N),BIG(N)
C
      PARAMETER (ONE=1.0D+00)
C
C     ----- JACOBI DIAGONALIZATION OF SYMMETRIC MATRIX -----
C     SYMMETRIC MATRIX -A- OF DIMENSION -N- IS DESTROYED ON EXIT.
C     ALL EIGENVECTORS ARE FOUND, SO -VEC- MUST BE SQUARE,
C     UNLESS SOMEONE TAKES THE TROUBLE TO LOOK AT -NMAX- BELOW.
C     -BIG- AND -JBIG- ARE SCRATCH WORK ARRAYS.
C
      CALL VCLR(VEC,1,LDVEC*N)
      DO 20 I = 1,N
        VEC(I,I) = ONE
   20 CONTINUE
C
      NB1 = N
      NB2 = (NB1*NB1+NB1)/2
      NMIN = 1
      NMAX = NB1
C
      CALL JACDIA(A,VEC,NB1,NB2,LDVEC,NMIN,NMAX,BIG,JBIG)
C
      DO 30 I=1,N
        EIG(I) = A((I*I+I)/2)
   30 CONTINUE
C
      CALL JACORD(VEC,EIG,NB1,LDVEC)
      RETURN
      END
C*MODULE EIGEN   *DECK JACDIA
      SUBROUTINE JACDIA(F,VEC,NB1,NB2,LDVEC,NMIN,NMAX,BIG,JBIG)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      LOGICAL GOPARR,DSKWRK,MASWRK
      DIMENSION F(NB2),VEC(LDVEC,NB1),BIG(NB1),JBIG(NB1)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ROOT2=0.707106781186548D+00 )
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, D1050=1.05D+00,
     *           D1500=1.5D+00, D3875=3.875D+00,
     *           D0500=0.5D+00, D1375=1.375D+00, D0250=0.25D+00 )
      PARAMETER (C2=1.0D-12, C3=4.0D-16,
     *           C4=2.0D-16, C5=8.0D-09, C6=3.0D-06 )
C
C      F IS THE MATRIX TO BE DIAGONALIZED, F IS STORED TRIANGULAR
C      VEC IS THE ARRAY OF EIGENVECTORS, DIMENSION NB1*NB1
C      BIG AND JBIG ARE TEMPORARY SCRATCH AREAS OF DIMENSION NB1
C      THE ROTATIONS AMONG THE FIRST NMIN BASIS FUNCTIONS ARE NOT
C      ACCOUNTED FOR.
C      THE ROTATIONS AMONG THE LAST NB1-NMAX BASIS FUNCTIONS ARE NOT
C      ACCOUNTED FOR.
C
      IEAA=0
      IEAB=0
      TT=ZERO
      EPS = 64.0D+00*EPSLON(ONE)
C
C      LOOP OVER COLUMNS (K) OF TRIANGULAR MATRIX TO DETERMINE
C      LARGEST OFF-DIAGONAL ELEMENTS IN ROW(I).
C
      DO 20 I=1,NB1
         BIG(I)=ZERO
         JBIG(I)=0
         IF(I.LT.NMIN  .OR.  I.EQ.1) GO TO 20
         II = (I*I-I)/2
         J=MIN(I-1,NMAX)
         DO 10 K=1,J
            IF(ABS(BIG(I)).GE.ABS(F(II+K))) GO TO 10
            BIG(I)=F(II+K)
            JBIG(I)=K
   10    CONTINUE
   20 CONTINUE
C
C     ----- 2X2 JACOBI ITERATIONS BEGIN HERE -----
C
      MAXIT=MAX(NB2*20,500)
      ITER=0
   30 CONTINUE
      ITER=ITER+1
C
C      FIND SMALLEST DIAGONAL ELEMENT
C
      SD=D1050
      JJ=0
      DO 40 J=1,NB1
         JJ=JJ+J
         SD= MIN(SD,ABS(F(JJ)))
   40 CONTINUE
      TEST = MAX(EPS, C2*MAX(SD,C6))
C
C      FIND LARGEST OFF-DIAGONAL ELEMENT
C
      T=ZERO
      I1=MAX(2,NMIN)
      IB = I1
      DO 50 I=I1,NB1
         IF(T.GE.ABS(BIG(I))) GO TO 50
         T= ABS(BIG(I))
         IB=I
   50 CONTINUE
C
C      TEST FOR CONVERGENCE, THEN DETERMINE ROTATION.
C
      IF(T.LT.TEST) RETURN
C                   ******
C
      IF(ITER.GT.MAXIT) THEN
         IF (MASWRK) THEN
            WRITE(6,*) 'JACOBI DIAGONALIZATION FAILS, DIMENSION=',NB1
            WRITE(6,9020) ITER,T,TEST,SD
         ENDIF
         CALL ABRT
         STOP
      END IF
C
      IA=JBIG(IB)
      IAA=IA*(IA-1)/2
      IBB=IB*(IB-1)/2
      DIF=F(IAA+IA)-F(IBB+IB)
      IF(ABS(DIF).GT.C3*T) GO TO 70
      SX=ROOT2
      CX=ROOT2
      GO TO 110
   70 T2X2=BIG(IB)/DIF
      T2X25=T2X2*T2X2
      IF(T2X25 . GT . C4) GO TO 80
      CX=ONE
      SX=T2X2
      GO TO 110
   80 IF(T2X25 . GT . C5) GO TO 90
      SX=T2X2*(ONE-D1500*T2X25)
      CX=ONE-D0500*T2X25
      GO TO 110
   90 IF(T2X25 . GT . C6) GO TO 100
      CX=ONE+T2X25*(T2X25*D1375 - D0500)
      SX= T2X2*(ONE + T2X25*(T2X25*D3875 - D1500))
      GO TO 110
  100 T=D0250  / SQRT(D0250   + T2X25)
      CX= SQRT(D0500   + T)
      SX= SIGN( SQRT(D0500   - T),T2X2)
  110 IEAR=IAA+1
      IEBR=IBB+1
C
      DO 230 IR=1,NB1
         T=F(IEAR)*SX
         F(IEAR)=F(IEAR)*CX+F(IEBR)*SX
         F(IEBR)=T-F(IEBR)*CX
         IF(IR-IA) 220,120,130
  120    TT=F(IEBR)
         IEAA=IEAR
         IEAB=IEBR
         F(IEBR)=BIG(IB)
         IEAR=IEAR+IR-1
         IF(JBIG(IR)) 200,220,200
  130    T=F(IEAR)
         IT=IA
         IEAR=IEAR+IR-1
         IF(IR-IB) 180,150,160
  150    F(IEAA)=F(IEAA)*CX+F(IEAB)*SX
         F(IEAB)=TT*CX+F(IEBR)*SX
         F(IEBR)=TT*SX-F(IEBR)*CX
         IEBR=IEBR+IR-1
         GO TO 200
  160    IF(  ABS(T) . GE .  ABS(F(IEBR))) GO TO 170
         IF(IB.GT.NMAX) GO TO 170
         T=F(IEBR)
         IT=IB
  170    IEBR=IEBR+IR-1
  180    IF(  ABS(T) . LT .  ABS(BIG(IR))) GO TO 190
         BIG(IR) = T
         JBIG(IR) = IT
         GO TO 220
  190    IF(IA . NE . JBIG(IR) . AND . IB . NE . JBIG(IR))  GO TO 220
  200    KQ=IEAR-IR-IA+1
         BIG(IR)=ZERO
         IR1=MIN(IR-1,NMAX)
         DO 210 I=1,IR1
            K=KQ+I
            IF(ABS(BIG(IR)) . GE . ABS(F(K)))  GO TO 210
            BIG(IR) = F(K)
            JBIG(IR)=I
  210    CONTINUE
  220    IEAR=IEAR+1
  230    IEBR=IEBR+1
C
      DO 240 I=1,NB1
         T1=VEC(I,IA)*CX + VEC(I,IB)*SX
         T2=VEC(I,IA)*SX - VEC(I,IB)*CX
         VEC(I,IA)=T1
         VEC(I,IB)=T2
  240 CONTINUE
      GO TO 30
C
 9020 FORMAT(1X,'ITER=',I6,' T,TEST,SD=',1P,3E20.10)
      END
C*MODULE EIGEN   *DECK JACORD
      SUBROUTINE JACORD(VEC,EIG,N,LDVEC)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION VEC(LDVEC,N),EIG(N)
C
C     ---- SORT EIGENDATA INTO ASCENDING ORDER -----
C
      DO 290 I = 1, N
         JJ = I
         DO 270 J = I, N
            IF (EIG(J) .LT. EIG(JJ)) JJ = J
  270    CONTINUE
         IF (JJ .EQ. I) GO TO 290
         T = EIG(JJ)
         EIG(JJ) = EIG(I)
         EIG(I) = T
         DO 280 J = 1, N
            T = VEC(J,JJ)
            VEC(J,JJ) = VEC(J,I)
            VEC(J,I) = T
  280    CONTINUE
  290 CONTINUE
      RETURN
      END
C*MODULE EIGEN   *DECK TINVTB
      SUBROUTINE TINVTB(NM,N,D,E,E2,M,W,IND,Z,
     *                  IERR,RV1,RV2,RV3,RV4,RV6)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION D(N),E(N),E2(N),W(M),Z(NM,M),
     *          RV1(N),RV2(N),RV3(N),RV4(N),RV6(N),IND(M)
      DOUBLE PRECISION MACHEP,NORM
      INTEGER P,Q,R,S,TAG,GROUP
C     ------------------------------------------------------------------
C
C     THIS ROUTINE IS A TRANSLATION OF THE INVERSE ITERATION TECH-
C     NIQUE IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971).
C
C     THIS ROUTINE FINDS THOSE EIGENVECTORS OF A TRIDIAGONAL
C     SYMMETRIC MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES,
C     USING INVERSE ITERATION.
C
C     ON INPUT-
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
C          DIMENSION STATEMENT,
C
C        N IS THE ORDER OF THE MATRIX,
C
C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
C
C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY,
C
C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E,
C          WITH ZEROS CORRESPONDING TO NEGLIGIBLE ELEMENTS OF E.
C          E(I) IS CONSIDERED NEGLIGIBLE IF IT IS NOT LARGER THAN
C          THE PRODUCT OF THE RELATIVE MACHINE PRECISION AND THE SUM
C          OF THE MAGNITUDES OF D(I) AND D(I-1).  E2(1) MUST CONTAIN
C          0.0 IF THE EIGENVALUES ARE IN ASCENDING ORDER, OR 2.0
C          IF THE EIGENVALUES ARE IN DESCENDING ORDER.  IF  BISECT,
C          TRIDIB, OR  IMTQLV  HAS BEEN USED TO FIND THE EIGENVALUES,
C          THEIR OUTPUT E2 ARRAY IS EXACTLY WHAT IS EXPECTED HERE,
C
C        M IS THE NUMBER OF SPECIFIED EIGENVALUES,
C
C        W CONTAINS THE M EIGENVALUES IN ASCENDING OR DESCENDING ORDER,
C
C        IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES
C          ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W --
C          1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM
C          THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.
C
C     ON OUTPUT-
C
C        ALL INPUT ARRAYS ARE UNALTERED,
C
C        Z CONTAINS THE ASSOCIATED SET OF ORTHONORMAL EIGENVECTORS.
C          ANY VECTOR WHICH FAILS TO CONVERGE IS SET TO ZERO,
C
C        IERR IS SET TO
C          ZERO       FOR NORMAL RETURN,
C          -R         IF THE EIGENVECTOR CORRESPONDING TO THE R-TH
C                     EIGENVALUE FAILS TO CONVERGE IN 5 ITERATIONS,
C
C        RV1, RV2, RV3, RV4, AND RV6 ARE TEMPORARY STORAGE ARRAYS.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
C     ------------------------------------------------------------------
C
C     ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
C
C                **********
      MACHEP = 2.0D+00**(-50)
C
      IERR = 0
      IF (M .EQ. 0) GO TO 680
      TAG = 0
      ORDER = 1.0D+00 - E2(1)
      XU = 0.0D+00
      UK = 0.0D+00
      X0 = 0.0D+00
      U  = 0.0D+00
      EPS2 = 0.0D+00
      EPS3 = 0.0D+00
      EPS4 = 0.0D+00
      GROUP = 0
      Q = 0
C     ********** ESTABLISH AND PROCESS NEXT SUBMATRIX **********
  100 P = Q + 1
      IP = P + 1
C
      DO 120 Q = P, N
      IF (Q .EQ. N) GO TO 140
      IF (E2(Q+1) .EQ. 0.0D+00) GO TO 140
  120 CONTINUE
C     ********** FIND VECTORS BY INVERSE ITERATION **********
  140 TAG = TAG + 1
      IQMP = Q - P + 1
      S = 0
C
      DO 660 R = 1, M
      IF (IND(R) .NE. TAG) GO TO 660
      ITS = 1
      X1 = W(R)
      IF (S .NE. 0) GO TO 220
C     ********** CHECK FOR ISOLATED ROOT **********
      XU = 1.0D+00
      IF (P .NE. Q) GO TO 160
      RV6(P) = 1.0D+00
      GO TO 600
  160 NORM = ABS(D(P))
C
      DO 180 I = IP, Q
  180 NORM = NORM + ABS(D(I)) + ABS(E(I))
C     ********** EPS2 IS THE CRITERION FOR GROUPING,
C                EPS3 REPLACES ZERO PIVOTS AND EQUAL
C                ROOTS ARE MODIFIED BY EPS3,
C                EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW **********
      EPS2 = 1.0D-03 * NORM
      EPS3 = MACHEP * NORM
      UK = IQMP
      EPS4 = UK * EPS3
      UK = EPS4 / SQRT(UK)
      S = P
  200 GROUP = 0
      GO TO 240
C     ********** LOOK FOR CLOSE OR COINCIDENT ROOTS **********
  220 IF (ABS(X1-X0) .GE. EPS2) GO TO 200
      GROUP = GROUP + 1
      IF (ORDER * (X1 - X0) .LE. 0.0D+00) X1 = X0 + ORDER * EPS3
C     ********** ELIMINATION WITH INTERCHANGES AND
C                INITIALIZATION OF VECTOR **********
  240 V = 0.0D+00
C
      DO 300 I = P, Q
      RV6(I) = UK
      IF (I .EQ. P) GO TO 280
      IF (ABS(E(I)) .LT. ABS(U)) GO TO 260
C     ********** WARNING -- A DIVIDE CHECK MAY OCCUR HERE IF
C                E2 ARRAY HAS NOT BEEN SPECIFIED CORRECTLY **********
      XU = U / E(I)
      RV4(I) = XU
      RV1(I-1) = E(I)
      RV2(I-1) = D(I) - X1
      RV3(I-1) = 0.0D+00
      IF (I .NE. Q) RV3(I-1) = E(I+1)
      U = V - XU * RV2(I-1)
      V = -XU * RV3(I-1)
      GO TO 300
  260 XU = E(I) / U
      RV4(I) = XU
      RV1(I-1) = U
      RV2(I-1) = V
      RV3(I-1) = 0.0D+00
  280 U = D(I) - X1 - XU * V
      IF (I .NE. Q) V = E(I+1)
  300 CONTINUE
C
      IF (U .EQ. 0.0D+00) U = EPS3
      RV1(Q) = U
      RV2(Q) = 0.0D+00
      RV3(Q) = 0.0D+00
C     ********** BACK SUBSTITUTION
C                FOR I=Q STEP -1 UNTIL P DO -- **********
  320 DO 340 II = P, Q
      I = P + Q - II
      RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I)
      V = U
      U = RV6(I)
  340 CONTINUE
C     ********** ORTHOGONALIZE WITH RESPECT TO PREVIOUS
C                MEMBERS OF GROUP **********
      IF (GROUP .EQ. 0) GO TO 400
      J = R
C
      DO 380 JJ = 1, GROUP
  360 J = J - 1
      IF (IND(J) .NE. TAG) GO TO 360
      XU = DDOT(IQMP,RV6(P),1,Z(P,J),1)
C
      CALL DAXPY(IQMP,-XU,Z(P,J),1,RV6(P),1)
C
  380 CONTINUE
C
  400 NORM = 0.0D+00
C
      DO 420 I = P, Q
  420 NORM = NORM + ABS(RV6(I))
C
      IF (NORM .GE. 1.0D+00) GO TO 560
C     ********** FORWARD SUBSTITUTION **********
      IF (ITS .EQ. 5) GO TO 540
      IF (NORM .NE. 0.0D+00) GO TO 440
      RV6(S) = EPS4
      S = S + 1
      IF (S .GT. Q) S = P
      GO TO 480
  440 XU = EPS4 / NORM
C
      DO 460 I = P, Q
  460 RV6(I) = RV6(I) * XU
C     ********** ELIMINATION OPERATIONS ON NEXT VECTOR
C                ITERATE **********
  480 DO 520 I = IP, Q
      U = RV6(I)
C     ********** IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE
C                WAS PERFORMED EARLIER IN THE
C                TRIANGULARIZATION PROCESS **********
      IF (RV1(I-1) .NE. E(I)) GO TO 500
      U = RV6(I-1)
      RV6(I-1) = RV6(I)
  500 RV6(I) = U - RV4(I) * RV6(I-1)
  520 CONTINUE
C
      ITS = ITS + 1
      GO TO 320
C     ********** SET ERROR -- NON-CONVERGED EIGENVECTOR **********
  540 IERR = -R
      XU = 0.0D+00
      GO TO 600
C     ********** NORMALIZE SO THAT SUM OF SQUARES IS
C                1 AND EXPAND TO FULL ORDER **********
  560 U = 0.0D+00
C
      DO 580 I = P, Q
      RV6(I) = RV6(I) / NORM
  580 U = U + RV6(I)**2
C
      XU = 1.0D+00 / SQRT(U)
C
  600 DO 620 I = 1, N
  620 Z(I,R) = 0.0D+00
C
      DO 640 I = P, Q
  640 Z(I,R) = RV6(I) * XU
C
      X0 = X1
  660 CONTINUE
C
      IF (Q .LT. N) GO TO 100
  680 RETURN
C     ********** LAST CARD OF TINVIT **********
      END
C*MODULE EIGEN   *DECK TQL2
C
C     ------------------------------------------------------------------
C
      SUBROUTINE TQL2(NM,N,D,E,Z,IERR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION MACHEP
      DIMENSION D(N),E(N),Z(NM,N)
C
C     THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL2,
C     NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND
C     WILKINSON.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971).
C
C     THIS ROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
C     OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE QL METHOD.
C     THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO
C     BE FOUND IF  TRED2  HAS BEEN USED TO REDUCE THIS
C     FULL MATRIX TO TRIDIAGONAL FORM.
C
C     ON INPUT-
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
C          DIMENSION STATEMENT,
C
C        N IS THE ORDER OF THE MATRIX,
C
C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
C
C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY,
C
C        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE
C          REDUCTION BY  TRED2, IF PERFORMED.  IF THE EIGENVECTORS
C          OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN
C          THE IDENTITY MATRIX.
C
C      ON OUTPUT-
C
C        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
C          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT
C          UNORDERED FOR INDICES 1,2,...,IERR-1,
C
C        E HAS BEEN DESTROYED,
C
C        Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC
C          TRIDIAGONAL (OR FULL) MATRIX.  IF AN ERROR EXIT IS MADE,
C          Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED
C          EIGENVALUES,
C
C        IERR IS SET TO
C          ZERO       FOR NORMAL RETURN,
C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
C                     DETERMINED AFTER 30 ITERATIONS.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
C     ------------------------------------------------------------------
C
C     ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
C
C                **********
      MACHEP = 2.0D+00**(-50)
C
      IERR = 0
      IF (N .EQ. 1) GO TO 400
C
      DO 100 I = 2, N
  100 E(I-1) = E(I)
C
      F = 0.0D+00
      B = 0.0D+00
      E(N) = 0.0D+00
C
      DO 300 L = 1, N
      J = 0
      H = MACHEP * (ABS(D(L)) + ABS(E(L)))
      IF (B .LT. H) B = H
C     ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT **********
      DO 120 M = L, N
      IF (ABS(E(M)) .LE. B) GO TO 140
C     ********** E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
C                THROUGH THE BOTTOM OF THE LOOP **********
  120 CONTINUE
C
  140 IF (M .EQ. L) GO TO 280
  160 IF (J .EQ. 30) GO TO 380
      J = J + 1
C     ********** FORM SHIFT **********
      L1 = L + 1
      G = D(L)
      P = (D(L1) - G) / (2.0D+00 * E(L))
      R = SQRT(P*P+1.0D+00)
      D(L) = E(L) / (P + SIGN(R,P))
      H = G - D(L)
C
      DO 180 I = L1, N
  180 D(I) = D(I) - H
C
      F = F + H
C     ********** QL TRANSFORMATION **********
      P = D(M)
      C = 1.0D+00
      S = 0.0D+00
      MML = M - L
C     ********** FOR I=M-1 STEP -1 UNTIL L DO -- **********
      DO 260 II = 1, MML
      I = M - II
      G = C * E(I)
      H = C * P
      IF (ABS(P) .LT. ABS(E(I))) GO TO 200
      C = E(I) / P
      R = SQRT(C*C+1.0D+00)
      E(I+1) = S * P * R
      S = C / R
      C = 1.0D+00 / R
      GO TO 220
  200 C = P / E(I)
      R = SQRT(C*C+1.0D+00)
      E(I+1) = S * E(I) * R
      S = 1.0D+00 / R
      C = C * S
  220 P = C * D(I) - S * G
      D(I+1) = H + S * (C * G + S * D(I))
C     ********** FORM VECTOR **********
      CALL DROT(N,Z(1,I+1),1,Z(1,I),1,C,S)
C
  260 CONTINUE
C
      E(L) = S * P
      D(L) = C * P
      IF (ABS(E(L)) .GT. B) GO TO 160
  280 D(L) = D(L) + F
  300 CONTINUE
C     ********** ORDER EIGENVALUES AND EIGENVECTORS **********
      DO 360 II = 2, N
      I = II - 1
      K = I
      P = D(I)
C
      DO 320 J = II, N
      IF (D(J) .GE. P) GO TO 320
      K = J
      P = D(J)
  320 CONTINUE
C
      IF (K .EQ. I) GO TO 360
      D(K) = D(I)
      D(I) = P
C
      CALL DSWAP(N,Z(1,I),1,Z(1,K),1)
C
  360 CONTINUE
C
      GO TO 400
C     ********** SET ERROR -- NO CONVERGENCE TO AN
C                EIGENVALUE AFTER 30 ITERATIONS **********
  380 IERR = L
  400 RETURN
C     ********** LAST CARD OF TQL2 **********
      END
C*MODULE EIGEN   *DECK TRBK3B
C
C     ------------------------------------------------------------------
C
      SUBROUTINE TRBK3B(NM,N,NV,A,M,Z)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(NV),Z(NM,M)
C
C     THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRBAK3,
C     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
C
C     THIS ROUTINE FORMS THE EIGENVECTORS OF A REAL SYMMETRIC
C     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
C     SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY  TRED3B.
C
C     ON INPUT-
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
C          DIMENSION STATEMENT,
C
C        N IS THE ORDER OF THE MATRIX,
C
C        NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A
C          AS DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT,
C
C        A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANSFORMATIONS
C          USED IN THE REDUCTION BY  TRED3B IN ITS FIRST
C          N*(N+1)/2 POSITIONS,
C
C        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED,
C
C        Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
C          IN ITS FIRST M COLUMNS.
C
C     ON OUTPUT-
C
C        Z CONTAINS THE TRANSFORMED EIGENVECTORS
C          IN ITS FIRST M COLUMNS.
C
C     NOTE THAT TRBAK3 PRESERVES VECTOR EUCLIDEAN NORMS.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
C     ------------------------------------------------------------------
C
      IF (M .EQ. 0) GO TO 140
      IF (N .EQ. 1) GO TO 140
C
      DO 120 I = 2, N
      L = I - 1
      IZ = (I * L) / 2
      IK = IZ + I
      H = A(IK)
      IF (H .EQ. 0.0D+00) GO TO 120
C
      DO 100 J = 1, M
      S = -DDOT(L,A(IZ+1),1,Z(1,J),1)
C
C     ********** DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW **********
      S = (S / H) / H
C
      CALL DAXPY(L,S,A(IZ+1),1,Z(1,J),1)
C
  100 CONTINUE
C
  120 CONTINUE
C
  140 RETURN
C     ********** LAST CARD OF TRBAK3 **********
      END
C*MODULE EIGEN   *DECK TRED3B
C
C     ------------------------------------------------------------------
C
      SUBROUTINE TRED3B(N,NV,A,D,E,E2)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(NV),D(N),E(N),E2(N)
C
C     THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED3,
C     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
C
C     THIS ROUTINE REDUCES A REAL SYMMETRIC MATRIX, STORED AS
C     A ONE-DIMENSIONAL ARRAY, TO A SYMMETRIC TRIDIAGONAL MATRIX
C     USING ORTHOGONAL SIMILARITY TRANSFORMATIONS.
C
C     ON INPUT-
C
C        N IS THE ORDER OF THE MATRIX,
C
C        NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A
C          AS DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT,
C
C        A CONTAINS THE LOWER TRIANGLE OF THE REAL SYMMETRIC
C          INPUT MATRIX, STORED ROW-WISE AS A ONE-DIMENSIONAL
C          ARRAY, IN ITS FIRST N*(N+1)/2 POSITIONS.
C
C     ON OUTPUT-
C
C        A CONTAINS INFORMATION ABOUT THE ORTHOGONAL
C          TRANSFORMATIONS USED IN THE REDUCTION,
C
C        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX,
C
C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO,
C
C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
C          E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
C     ------------------------------------------------------------------
C
C     ********** FOR I=N STEP -1 UNTIL 1 DO -- **********
      DO 300 II = 1, N
      I = N + 1 - II
      L = I - 1
      IZ = (I * L) / 2
      H = 0.0D+00
      SCALE = 0.0D+00
      IF (L .LT. 1) GO TO 120
C     ********** SCALE ROW (ALGOL TOL THEN NOT NEEDED) **********
      DO 100 K = 1, L
      IZ = IZ + 1
      D(K) = A(IZ)
      SCALE = SCALE + ABS(D(K))
  100 CONTINUE
C
      IF (SCALE .NE. 0.0D+00) GO TO 140
  120 E(I) = 0.0D+00
      E2(I) = 0.0D+00
      GO TO 280
C
  140 DO 160 K = 1, L
      D(K) = D(K) / SCALE
      H = H + D(K) * D(K)
  160 CONTINUE
C
      E2(I) = SCALE * SCALE * H
      F = D(L)
      G = -SIGN(SQRT(H),F)
      E(I) = SCALE * G
      H = H - F * G
      D(L) = F - G
      A(IZ) = SCALE * D(L)
      IF (L .EQ. 1) GO TO 280
      F = 0.0D+00
C
      JK = 1
      DO 220 J = 1, L
      JM1 = J - 1
      DT = D(J)
      G = 0.0D+00
C     ********** FORM ELEMENT OF A*U **********
      IF (JM1 .EQ. 0) GO TO 200
      DO 180 K = 1, JM1
      E(K) = E(K) + DT * A(JK)
      G = G + D(K) * A(JK)
      JK = JK + 1
  180 CONTINUE
  200 E(J) = G + A(JK) * DT
      JK = JK + 1
C     ********** FORM ELEMENT OF P **********
  220 CONTINUE
      F = 0.0D+00
      DO 240 J = 1, L
      E(J) = E(J) / H
      F = F + E(J) * D(J)
  240 CONTINUE
C
      HH = F / (H + H)
      JK = 0
C     ********** FORM REDUCED A **********
      DO 260 J = 1, L
      F = D(J)
      G = E(J) - HH * F
      E(J) = G
C
      DO 260 K = 1, J
      JK = JK + 1
      A(JK) = A(JK) - F * E(K) - G * D(K)
  260 CONTINUE
C
  280 D(I) = A(IZ+1)
      A(IZ+1) = SCALE * SQRT(H)
  300 CONTINUE
C
      RETURN
C     ********** LAST CARD OF TRED3 **********
      END
