X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-NEWSC%2Feigen.f;fp=source%2Funres%2Fsrc_MD-NEWSC%2Feigen.f;h=e4088eed6b3ab8cf9551eaa1f8556bb69c8300cb;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git diff --git a/source/unres/src_MD-NEWSC/eigen.f b/source/unres/src_MD-NEWSC/eigen.f new file mode 100644 index 0000000..e4088ee --- /dev/null +++ b/source/unres/src_MD-NEWSC/eigen.f @@ -0,0 +1,2351 @@ +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