+++ /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