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