--- /dev/null
+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)
+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, ESTPI1, DASUM, DDOT, DNRM2
+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
+ IF (RV1(I-1) .EQ. E(I)) 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
+ RHO = ESTPI1(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 ESTPI1
+ DOUBLE PRECISION FUNCTION 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 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)
+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,DDOT,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)
+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 DASUM, DNRM2
+ 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(400)
+ COMMON /MACHSW/ KDIAG,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
+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 = EIGENVECTORS (OUTPUT)
+C VECTOR = EIGENVECTORS (OUTPUT)
+C IERR = ERROR FLAG (OUTPUT)
+C IWRK = N INTEGER WORDS OF SCRATCH SPACE
+C
+ 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
+ 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
projects / unres.git / blobdiff