1 C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
2 C 31 MAR 94 - MWS - ADD A VARIABLE TO END OF MACHSW COMMON
3 C 26 JUN 93 - MWS - ETRED3: ADD RETURN FOR SPECIAL CASE N=1
4 C 4 JAN 92 - TLW - MAKE WRITES PARALLEL;ADD COMMON PAR
5 C 30 AUG 91 - MWS - JACDIA: LIMIT ITERATIONS, USE EPSLON IN TEST.
6 C 14 JUL 91 - MWS - JACOBI DIAGONALIZATION ALLOWS FOR LDVEC.NE.N
7 C 29 JAN 91 - TLW - GLDIAG: CHANGED COMMON DIAGSW TO MACHSW
8 C 29 OCT 90 - STE - FIX JACDIA UNDEFINED VARIABLE BUG
9 C 14 SEP 90 - MK - NEW JACOBI DIAGONALIZATION (KDIAG=3)
10 C 27 MAR 88 - MWS - ALLOW FOR VECTOR ROUTINE IN GLDIAG
11 C 11 AUG 87 - MWS - SANITIZE CONSTANTS IN EQLRAT
12 C 15 FEB 87 - STE - FIX EINVIT SUB-MATRIX LOOP LIMIT
13 C SCRATCH ARRAYS ARE N*8 REAL AND N INTEGER
14 C 8 DEC 86 - STE - USE PERF INDEX FROM ESTPI1 TO JUDGE EINVIT FAILURE
15 C 30 NOV 86 - STE - DELETE LIGENB, MAKE EVVRSP DEFAULT
16 C (GIVEIS FAILS ON CRAY FOR BENCHMC AND BENCHCI)
17 C 7 JUL 86 - JAB - SANITIZE FLOATING POINT CONSTANTS
18 C 11 OCT 85 - STE - LIGENB,TQL2: USE DROT,DSWAP; TINVTB: SCALE VECTOR
19 C BEFORE NORMALIZING; GENERIC FUNCTIONS
20 C 24 FEB 84 - STE - INITIALIZE INDEX ARRAY FOR LIGENB IN GLDIAG
21 C 1 DEC 83 - STE - CHANGE MACHEP FROM 2**-54 TO 2**-50
22 C 28 SEP 82 - MWS - CONVERT TO IBM
24 C*MODULE EIGEN *DECK EINVIT
25 SUBROUTINE EINVIT(NM,N,D,E,E2,M,W,IND,Z,IERR,RV1,RV2,RV3,RV4,RV6)
28 C* THIS IS A MODIFICATION OF TINVIT FROM EISPACK EDITION 3
30 C* TINVIT IS A TRANSLATION OF THE INVERSE ITERATION TECHNIQUE
31 C* IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON.
32 C* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971).
33 C* THIS VERSION IS BY S. T. ELBERT (AMES LABORATORY-USDOE)
36 C* THIS ROUTINE FINDS THOSE EIGENVECTORS OF A TRIDIAGONAL
37 C* SYMMETRIC MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES.
44 C* MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
45 C* ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
46 C* DIMENSION STATEMENT.
49 C* CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
51 C* CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
52 C* IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY.
54 C* CONTAINS THE SQUARES OF CORRESPONDING ELEMENTS OF E,
55 C* WITH ZEROS CORRESPONDING TO NEGLIGIBLE ELEMENTS OF E.
56 C* E(I) IS CONSIDERED NEGLIGIBLE IF IT IS NOT LARGER THAN
57 C* THE PRODUCT OF THE RELATIVE MACHINE PRECISION AND THE
58 C* SUM OF THE MAGNITUDES OF D(I) AND D(I-1). E2(1) MUST
59 C* CONTAIN 0.0 IF THE EIGENVALUES ARE IN ASCENDING ORDER,
60 C* OR 2.0 IF THE EIGENVALUES ARE IN DESCENDING ORDER.
61 C* IF TQLRAT, BISECT, TRIDIB, OR IMTQLV
62 C* HAS BEEN USED TO FIND THE EIGENVALUES, THEIR
63 C* OUTPUT E2 ARRAY IS EXACTLY WHAT IS EXPECTED HERE.
65 C* THE NUMBER OF SPECIFIED EIGENVECTORS.
67 C* CONTAINS THE M EIGENVALUES IN ASCENDING
68 C* OR DESCENDING ORDER.
70 C* CONTAINS IN FIRST M POSITIONS THE SUBMATRIX INDICES
71 C* ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W --
72 C* 1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX
73 C* FROM THE TOP, 2 FOR THOSE BELONGING TO THE SECOND
75 C* IERR - INTEGER (LOGICAL UNIT NUMBER)
76 C* LOGICAL UNIT FOR ERROR MESSAGES
79 C* ALL INPUT ARRAYS ARE UNALTERED.
80 C* Z - W.P. REAL (NM,M)
81 C* CONTAINS THE ASSOCIATED SET OF ORTHONORMAL
82 C* EIGENVECTORS. ANY VECTOR WHICH WHICH FAILS TO CONVERGE
83 C* IS LEFT AS IS (BUT NORMALIZED) WHEN ITERATING STOPPED.
86 C* ZERO FOR NORMAL RETURN,
87 C* -R IF THE EIGENVECTOR CORRESPONDING TO THE R-TH
88 C* EIGENVALUE FAILS TO CONVERGE IN 5 ITERATIONS.
89 C* (ONLY LAST FAILURE TO CONVERGE IS REPORTED)
91 C* RV1, RV2, RV3, RV4, AND RV6 ARE TEMPORARY STORAGE ARRAYS.
93 C* RV1 - W.P. REAL (N)
94 C* DIAGONAL ELEMENTS OF U FROM LU DECOMPOSITION
95 C* RV2 - W.P. REAL (N)
96 C* SUPER(1)-DIAGONAL ELEMENTS OF U FROM LU DECOMPOSITION
97 C* RV3 - W.P. REAL (N)
98 C* SUPER(2)-DIAGONAL ELEMENTS OF U FROM LU DECOMPOSITION
99 C* RV4 - W.P. REAL (N)
100 C* ELEMENTS DEFINING L IN LU DECOMPOSITION
101 C* RV6 - W.P. REAL (N)
102 C* APPROXIMATE EIGENVECTOR
104 C* DIFFERENCES FROM EISPACK 3 -
105 C* EPS3 IS SCALED BY EPSCAL (ENHANCES CONVERGENCE, BUT
107 C* ONE MORE ITERATION (MINIMUM 2) IS PERFORMED AFTER CONVERGENCE
108 C* (ENHANCES ACCURACY)!
109 C* REPLACE LOOP WITH PYTHAG WITH SINGLE CALL TO DNRM2!
110 C* IF NOT CONVERGED, USE PERFORMANCE INDEX TO DECIDE ON ERROR
111 C* VALUE SETTING, BUT DO NOT STOP!
112 C* L.U. FOR ERROR MESSAGES PASSED THROUGH IERR
113 C* USE PARAMETER STATEMENTS AND GENERIC INTRINSIC FUNCTIONS
115 C* USE IF-THEN-ELSE TO CLARIFY LOGIC
116 C* LOOP OVER SUBSPACES MADE INTO DO LOOP.
117 C* LOOP OVER INVERSE ITERATIONS MADE INTO DO LOOP
118 C* ZERO ONLY REQUIRED PORTIONS OF OUTPUT VECTOR
121 C* QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
122 C* B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
125 LOGICAL CONVGD,GOPARR,DSKWRK,MASWRK
127 INTEGER GROUP,I,IERR,ITS,J,JJ,M,N,NM,P,Q,R,S,SUBMAT,TAG
130 DOUBLE PRECISION D(N),E(N),E2(N),W(M),Z(NM,M)
131 DOUBLE PRECISION RV1(N),RV2(N),RV3(N),RV4(N),RV6(N)
132 DOUBLE PRECISION ANORM,EPS2,EPS3,EPS4,NORM,ORDER,RHO,U,UK,V
133 DOUBLE PRECISION X0,X1,XU
134 DOUBLE PRECISION EPSCAL,GRPTOL,HUNDRD,ONE,TEN,ZERO
135 DOUBLE PRECISION EPSLON, ESTPI1, DASUM, DDOT, DNRM2
137 COMMON /PAR / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
139 PARAMETER (ZERO = 0.0D+00, ONE = 1.0D+00, GRPTOL = 0.001D+00)
140 PARAMETER (EPSCAL = 0.5D+00, HUNDRD = 100.0D+00, TEN = 10.0D+00)
142 001 FORMAT(' EIGENVECTOR ROUTINE EINVIT DID NOT CONVERGE FOR VECTOR'
143 * ,I5,'. NORM =',1P,E10.2,' PERFORMANCE INDEX =',E10.2/
144 * ' (AN ERROR HALT WILL OCCUR IF THE PI IS GREATER THAN 100)')
146 C-----------------------------------------------------------------------
163 C .......... ESTABLISH AND PROCESS NEXT SUBMATRIX ..........
166 IF (E2(Q+1) .EQ. ZERO) GO TO 140
170 C .......... FIND VECTORS BY INVERSE ITERATION ..........
178 IF (IND(R) .NE. TAG) GO TO 920
181 IF (S .NE. 0) GO TO 510
183 C .......... CHECK FOR ISOLATED ROOT ..........
194 NORM = MAX( NORM, ABS(D(I)) + ABS(E(I)) )
197 C .......... EPS2 IS THE CRITERION FOR GROUPING,
198 C EPS3 REPLACES ZERO PIVOTS AND EQUAL
199 C ROOTS ARE MODIFIED BY EPS3,
200 C EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW .........
203 EPS3 = EPSCAL * EPSLON(NORM)
211 C .......... LOOK FOR CLOSE OR COINCIDENT ROOTS ..........
213 510 IF (ABS(X1-X0) .GE. EPS2) THEN
223 IF (ORDER * (X1 - X0) .LE. EPS3) X1 = X0 + ORDER * EPS3
226 C .......... ELIMINATION WITH INTERCHANGES AND
227 C INITIALIZATION OF VECTOR ..........
236 IF (ABS(E(I)) .GT. ABS(U)) THEN
238 C EXCHANGE ROWS BEFORE ELIMINATION
240 C *** WARNING -- A DIVIDE CHECK MAY OCCUR HERE IF
241 C E2 ARRAY HAS NOT BEEN SPECIFIED CORRECTLY .......
248 U = V - XU * RV2(I-1)
253 C STRAIGHT ELIMINATION
260 U = D(I) - X1 - XU * V
265 IF (ABS(U) .LE. EPS3) U = EPS3
270 C DO INVERSE ITERATIONS
274 IF (ITS .EQ. 1) GO TO 600
276 C .......... FORWARD SUBSTITUTION ..........
278 IF (NORM .EQ. ZERO) THEN
284 CALL DSCAL (Q-P+1, XU, RV6(P), 1)
287 C ... ELIMINATION OPERATIONS ON NEXT VECTOR
292 C IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE
293 C WAS PERFORMED EARLIER IN THE
294 C TRIANGULARIZATION PROCESS ..........
296 IF (RV1(I-1) .EQ. E(I)) THEN
302 RV6(I) = U - RV4(I) * RV6(I-1)
306 C .......... BACK SUBSTITUTION
308 RV6(Q) = RV6(Q) / RV1(Q)
312 DO 620 I = Q-1, P, -1
313 RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I)
318 IF (GROUP .EQ. 0) GO TO 700
320 C ....... ORTHOGONALIZE WITH RESPECT TO PREVIOUS
321 C MEMBERS OF GROUP ..........
326 IF (IND(J) .NE. TAG) GO TO 630
327 CALL DAXPY(Q-P+1, -DDOT(Q-P+1,RV6(P),1,Z(P,J),1),
330 NORM = DASUM(Q-P+1, RV6(P), 1)
333 IF (CONVGD) GO TO 840
334 IF (NORM .GE. ONE) CONVGD = .TRUE.
337 C .......... NORMALIZE SO THAT SUM OF SQUARES IS
338 C 1 AND EXPAND TO FULL ORDER ..........
342 XU = ONE / DNRM2(Q-P+1,RV6(P),1)
355 IF (.NOT.CONVGD) THEN
356 RHO = ESTPI1(Q-P+1,X1,D(P),E(P),Z(P,R),ANORM)
357 IF (RHO .GE. TEN .AND. LUEMSG .GT. 0 .AND. MASWRK)
358 * WRITE(LUEMSG,001) R,NORM,RHO
360 C *** SET ERROR -- NON-CONVERGED EIGENVECTOR ..........
362 IF (RHO .GT. HUNDRD) IERR = -R
368 IF (Q .EQ. N) GO TO 940
373 C*MODULE EIGEN *DECK ELAUM
374 SUBROUTINE ELAU(HINV,L,D,A,E)
376 DOUBLE PRECISION A(*)
377 DOUBLE PRECISION D(L)
378 DOUBLE PRECISION E(L)
381 DOUBLE PRECISION HALF
383 DOUBLE PRECISION HINV
384 DOUBLE PRECISION ZERO
386 PARAMETER (ZERO = 0.0D+00, HALF = 0.5D+00)
398 E(K) = E(K) + A(JK) * F
406 C .......... FORM P ..........
414 C .......... FORM Q ..........
418 250 E(J) = E(J) - HH * D(J)
422 C*MODULE EIGEN *DECK EPSLON
423 DOUBLE PRECISION FUNCTION EPSLON (X)
426 C* THIS ROUTINE WAS TAKEN FROM EISPACK EDITION 3 DATED 4/6/83
427 C* THIS VERSION IS BY S. T. ELBERT, AMES LABORATORY-USDOE NOV 1986
430 C* ESTIMATE UNIT ROUNDOFF IN QUANTITIES OF SIZE X.
433 C* X - WORKING PRECISION REAL
434 C* VALUES TO FIND EPSLON FOR
437 C* EPSLON - WORKING PRECISION REAL
438 C* SMALLEST POSITIVE VALUE SUCH THAT X+EPSLON .NE. ZERO
441 C* THIS ROUTINE SHOULD PERFORM PROPERLY ON ALL SYSTEMS
442 C* SATISFYING THE FOLLOWING TWO ASSUMPTIONS,
443 C* 1. THE BASE USED IN REPRESENTING FLOATING POINT
444 C* NUMBERS IS NOT A POWER OF THREE.
445 C* 2. THE QUANTITY A IN STATEMENT 10 IS REPRESENTED TO
446 C* THE ACCURACY USED IN FLOATING POINT VARIABLES
447 C* THAT ARE STORED IN MEMORY.
448 C* THE STATEMENT NUMBER 10 AND THE GO TO 10 ARE INTENDED TO
449 C* FORCE OPTIMIZING COMPILERS TO GENERATE CODE SATISFYING
451 C* UNDER THESE ASSUMPTIONS, IT SHOULD BE TRUE THAT,
452 C* A IS NOT EXACTLY EQUAL TO FOUR-THIRDS,
453 C* B HAS A ZERO FOR ITS LAST BIT OR DIGIT,
454 C* C IS NOT EXACTLY EQUAL TO ONE,
455 C* EPS MEASURES THE SEPARATION OF 1.0 FROM
456 C* THE NEXT LARGER FLOATING POINT NUMBER.
457 C* THE DEVELOPERS OF EISPACK WOULD APPRECIATE BEING INFORMED
458 C* ABOUT ANY SYSTEMS WHERE THESE ASSUMPTIONS DO NOT HOLD.
460 C* DIFFERENCES FROM EISPACK 3 -
461 C* USE IS MADE OF PARAMETER STATEMENTS AND INTRINSIC FUNCTIONS
462 C* --NO EXECUTEABLE CODE CHANGES--
465 C* QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
466 C* B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
468 DOUBLE PRECISION A,B,C,EPS,X
469 DOUBLE PRECISION ZERO, ONE, THREE, FOUR
471 PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, THREE=3.0D+00, FOUR=4.0D+00)
473 C-----------------------------------------------------------------------
479 IF (EPS .EQ. ZERO) GO TO 10
483 C*MODULE EIGEN *DECK EQLRAT
484 SUBROUTINE EQLRAT(N,DIAG,E,E2IN,D,IND,IERR,E2)
487 C* THIS IS A MODIFICATION OF ROUTINE EQLRAT FROM EISPACK EDITION 3
488 C* DATED AUGUST 1983.
489 C* TQLRAT IS A TRANSLATION OF THE ALGOL PROCEDURE TQLRAT,
490 C* ALGORITHM 464, COMM. ACM 16, 689(1973) BY REINSCH.
491 C* THIS VERSION IS BY S. T. ELBERT (AMES LABORATORY-USDOE)
494 C* THIS ROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC
495 C* TRIDIAGONAL MATRIX
502 C* THE ORDER OF THE MATRIX.
504 C* CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
505 C* E2 - W.P. REAL (N)
506 C* CONTAINS THE SQUARES OF THE SUBDIAGONAL ELEMENTS OF
507 C* THE INPUT MATRIX IN ITS LAST N-1 POSITIONS.
508 C* E2(1) IS ARBITRARY.
512 C* CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN
513 C* ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND
514 C* ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE
515 C* THE SMALLEST EIGENVALUES.
516 C* E2 - W.P. REAL (N)
520 C* ZERO FOR NORMAL RETURN,
521 C* J IF THE J-TH EIGENVALUE HAS NOT BEEN
522 C* DETERMINED AFTER 30 ITERATIONS.
524 C* DIFFERENCES FROM EISPACK 3 -
525 C* G=G+B INSTEAD OF IF(G.EQ.0) G=B ; B=B/4
526 C* F77 BACKWARD LOOPS INSTEAD OF F66 CONSTRUCT
527 C* GENERIC INTRINSIC FUNCTIONS
528 C* ARRARY IND ADDED FOR USE BY EINVIT
531 C* QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
532 C* B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
534 INTEGER I,J,L,M,N,II,L1,IERR
537 DOUBLE PRECISION D(N),E(N),E2(N),DIAG(N),E2IN(N)
538 DOUBLE PRECISION B,C,F,G,H,P,R,S,T,EPSLON
539 DOUBLE PRECISION SCALE,ZERO,ONE
541 PARAMETER (ZERO = 0.0D+00, SCALE= 1.0D+00/64.0D+00, ONE = 1.0D+00)
543 C-----------------------------------------------------------------------
549 IF (N .EQ. 1) GO TO 1001
553 100 E2(I-1) = E2IN(I)
563 H = ABS(D(L)) + ABS(E(L))
564 IF (T .GE. H) GO TO 105
570 C .......... LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT ..........
573 IF (E2(M) .GT. C) GO TO 110
574 C .......... E2(N) IS ALWAYS ZERO, SO THERE IS AN EXIT
575 C FROM THE LOOP ..........
577 IF (M .LE. K) GO TO 125
578 IF (M .NE. N) E2IN(M+1) = ZERO
582 IF (M .EQ. L) GO TO 210
587 C .......... FORM SHIFT ..........
591 P = (D(L1) - G) / (2.0D+00 * S)
592 R = SQRT(P*P+1.0D+00)
593 D(L) = S / (P + SIGN(R,P))
600 C .......... RATIONAL QL TRANSFORMATION ..........
609 D(I+1) = H + S * (H + D(I))
610 G = D(I) - E2(I) / G + B
616 C .......... GUARD AGAINST UNDERFLOW IN CONVERGENCE TEST
617 IF (H .EQ. ZERO) GO TO 210
618 IF (ABS(E2(L)) .LE. ABS(C/H)) GO TO 210
620 IF (E2(L) .EQ. ZERO) GO TO 210
622 C .......... SET ERROR -- NO CONVERGENCE TO AN
623 C EIGENVALUE AFTER 30 ITERATIONS ..........
630 C .......... ORDER EIGENVALUES ..........
632 IF (L .EQ. 1) GO TO 250
633 IF (P .LT. D(1)) GO TO 230
635 C .......... LOOP TO FIND ORDERED POSITION
637 IF (P .LT. D(I)) GO TO 220
640 IF (I .EQ. L) GO TO 250
642 DO 240 II = L, I+1, -1
654 C*MODULE EIGEN *DECK ESTPI1
655 DOUBLE PRECISION FUNCTION ESTPI1 (N,EVAL,D,E,X,ANORM)
658 C* STEPHEN T. ELBERT (AMES LABORATORY-USDOE) DATE: 5 DEC 1986
661 C* EVALUATE SYMMETRIC TRIDIAGONAL MATRIX PERFORMANCE INDEX
667 C* THIS ROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX A*X-X*EVAL
668 C* WHERE A IS A SYMMETRIC TRIDIAGONAL MATRIX STORED
669 C* IN THE DIAGONAL (D) AND SUB-DIAGONAL (E) VECTORS, EVAL IS THE
670 C* EIGENVALUE OF AN EIGENVECTOR OF A, NAMELY X.
671 C* THIS NORM IS SCALED BY MACHINE ACCURACY FOR THE PROBLEM SIZE.
672 C* ALL NORMS APPEARING IN THE COMMENTS BELOW ARE 1-NORMS.
676 C* THE ORDER OF THE MATRIX A.
678 C* THE EIGENVALUE CORRESPONDING TO VECTOR X.
680 C* THE DIAGONAL VECTOR OF A.
682 C* THE SUB-DIAGONAL VECTOR OF A.
684 C* AN EIGENVECTOR OF A.
686 C* THE NORM OF A IF IT HAS BEEN PREVIOUSLY COMPUTED.
690 C* THE NORM OF A, COMPUTED IF INITIALLY ZERO.
691 C* ESTPI1 - W.P. REAL
692 C* !!A*X-X*EVAL!! / (EPSLON(10*N)*!!A!!*!!X!!);
693 C* WHERE EPSLON(X) IS THE SMALLEST NUMBER SUCH THAT
694 C* X + EPSLON(X) .NE. X
696 C* ESTPI1 .LT. 1 == SATISFACTORY PERFORMANCE
697 C* .GE. 1 AND .LE. 100 == MARGINAL PERFORMANCE
698 C* .GT. 100 == POOR PERFORMANCE
699 C* (SEE LECT. NOTES IN COMP. SCI. VOL.6 PP 124-125)
701 DOUBLE PRECISION ANORM,EVAL,RNORM,SIZE,XNORM
702 DOUBLE PRECISION D(N), E(N), X(N)
703 DOUBLE PRECISION EPSLON, ONE, ZERO
705 PARAMETER (ZERO = 0.0D+00, ONE = 1.0D+00)
707 C-----------------------------------------------------------------------
710 IF( N .LE. 1 ) RETURN
712 IF (ANORM .EQ. ZERO) THEN
716 ANORM = MAX( ABS(D(1)) + ABS(E(2))
717 * ,ABS(D(N)) + ABS(E(N)))
719 ANORM = MAX( ANORM, ABS(E(I))+ABS(D(I))+ABS(E(I+1)))
721 IF(ANORM .EQ. ZERO) ANORM = ONE
724 C COMPUTE NORMS OF RESIDUAL AND EIGENVECTOR
726 XNORM = ABS(X(1)) + ABS(X(N))
727 RNORM = ABS( (D(1)-EVAL)*X(1) + E(2)*X(2))
728 * +ABS( (D(N)-EVAL)*X(N) + E(N)*X(N-1))
730 XNORM = XNORM + ABS(X(I))
731 RNORM = RNORM + ABS(E(I)*X(I-1) + (D(I)-EVAL)*X(I)
735 ESTPI1 = RNORM / (EPSLON(SIZE)*ANORM*XNORM)
738 C*MODULE EIGEN *DECK ETRBK3
739 SUBROUTINE ETRBK3(NM,N,NV,A,M,Z)
742 C* THIS IS A MODIFICATION OF ROUTINE TRBAK3 FROM EISPACK EDITION 3
743 C* DATED AUGUST 1983.
744 C* EISPACK TRBAK3 IS A TRANSLATION OF THE ALGOL PROCEDURE TRBAK3,
745 C* NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
746 C* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
747 C* THIS VERSION IS BY S. T. ELBERT (AMES LABORATORY-USDOE)
750 C* THIS ROUTINE FORMS THE EIGENVECTORS OF A REAL SYMMETRIC
751 C* MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
752 C* SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY ETRED3.
755 C* THE CALCULATION IS CARRIED OUT BY FORMING THE MATRIX PRODUCT
757 C* WHERE Q IS A PRODUCT OF THE ORTHOGONAL SYMMETRIC MATRICES
758 C* Q = PROD(I)[1 - U(I)*.TRANSPOSE.U(I)*H(I)]
759 C* U IS THE AUGMENTED SUB-DIAGONAL ROWS OF A AND
760 C* Z IS THE SET OF EIGENVECTORS OF THE TRIDIAGONAL
761 C* MATRIX F WHICH WAS FORMED FROM THE ORIGINAL SYMMETRIC
762 C* MATRIX C BY THE SIMILARITY TRANSFORMATION
763 C* F = Q(TRANSPOSE) C Q
764 C* NOTE THAT ETRBK3 PRESERVES VECTOR EUCLIDEAN NORMS.
772 C* MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
773 C* ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
774 C* DIMENSION STATEMENT.
776 C* THE ORDER OF THE MATRIX A.
778 C* MUST BE SET TO THE DIMENSION OF THE ARRAY A AS
779 C* DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT.
780 C* A - W.P. REAL (NV)
781 C* CONTAINS INFORMATION ABOUT THE ORTHOGONAL
782 C* TRANSFORMATIONS USED IN THE REDUCTION BY ETRED3 IN
783 C* ITS FIRST NV = N*(N+1)/2 POSITIONS.
785 C* THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.
786 C* Z - W.P REAL (NM,M)
787 C* CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
788 C* IN ITS FIRST M COLUMNS.
791 C* Z - W.P. REAL (NM,M)
792 C* CONTAINS THE TRANSFORMED EIGENVECTORS
793 C* IN ITS FIRST M COLUMNS.
795 C* DIFFERENCES WITH EISPACK 3 -
796 C* THE TWO INNER LOOPS ARE REPLACED BY DDOT AND DAXPY.
797 C* MULTIPLICATION USED INSTEAD OF DIVISION TO FIND S.
798 C* OUTER LOOP RANGE CHANGED FROM 2,N TO 3,N.
799 C* ADDRESS POINTERS FOR A SIMPLIFIED.
802 C* QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
803 C* B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
805 INTEGER I,II,IM1,IZ,J,M,N,NM,NV
807 DOUBLE PRECISION A(NV),Z(NM,M)
808 DOUBLE PRECISION H,S,DDOT,ZERO
810 PARAMETER (ZERO = 0.0D+00)
812 C-----------------------------------------------------------------------
822 IF (H .EQ. ZERO) GO TO 140
825 S = -( DDOT(IM1,A(IZ),1,Z(1,J),1) * H) * H
826 CALL DAXPY(IM1,S,A(IZ),1,Z(1,J),1)
831 C*MODULE EIGEN *DECK ETRED3
832 SUBROUTINE ETRED3(N,NV,A,D,E,E2)
835 C* THIS IS A MODIFICATION OF ROUTINE TRED3 FROM EISPACK EDITION 3
836 C* DATED AUGUST 1983.
837 C* EISPACK TRED3 IS A TRANSLATION OF THE ALGOL PROCEDURE TRED3,
838 C* NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
839 C* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
840 C* THIS VERSION IS BY S. T. ELBERT, AMES LABORATORY-USDOE JUN 1986
843 C* THIS ROUTINE REDUCES A REAL SYMMETRIC (PACKED) MATRIX, STORED
844 C* AS A ONE-DIMENSIONAL ARRAY, TO A SYMMETRIC TRIDIAGONAL MATRIX
845 C* USING ORTHOGONAL SIMILARITY TRANSFORMATIONS, PRESERVING THE
846 C* INFORMATION ABOUT THE TRANSFORMATIONS IN A.
849 C* THE TRIDIAGONAL REDUCTION IS PERFORMED IN THE FOLLOWING WAY.
850 C* STARTING WITH J=N, THE ELEMENTS IN THE J-TH ROW TO THE
851 C* LEFT OF THE DIAGONAL ARE FIRST SCALED, TO AVOID POSSIBLE
852 C* UNDERFLOW IN THE TRANSFORMATION THAT MIGHT RESULT IN SEVERE
853 C* DEPARTURE FROM ORTHOGONALITY. THE SUM OF SQUARES SIGMA OF
854 C* THESE SCALED ELEMENTS IS NEXT FORMED. THEN, A VECTOR U AND
856 C* H = U(TRANSPOSE) * U / 2
857 C* DEFINE A REFLECTION OPERATOR
858 C* P = I - U * U(TRANSPOSE) / H
859 C* WHICH IS ORTHOGONAL AND SYMMETRIC AND FOR WHICH THE
860 C* SIMILIARITY TRANSFORMATION PAP ELIMINATES THE ELEMENTS IN
861 C* THE J-TH ROW OF A TO THE LEFT OF THE SUBDIAGONAL AND THE
862 C* SYMMETRICAL ELEMENTS IN THE J-TH COLUMN.
864 C* THE NON-ZERO COMPONENTS OF U ARE THE ELEMENTS OF THE J-TH
865 C* ROW TO THE LEFT OF THE DIAGONAL WITH THE LAST OF THEM
866 C* AUGMENTED BY THE SQUARE ROOT OF SIGMA PREFIXED BY THE SIGN
867 C* OF THE SUBDIAGONAL ELEMENT. BY STORING THE TRANSFORMED SUB-
868 C* DIAGONAL ELEMENT IN E(J) AND NOT OVERWRITING THE ROW
869 C* ELEMENTS ELIMINATED IN THE TRANSFORMATION, FULL INFORMATION
870 C* ABOUT P IS SAVE FOR LATER USE IN ETRBK3.
872 C* THE TRANSFORMATION SETS E2(J) EQUAL TO SIGMA AND E(J)
873 C* EQUAL TO THE SQUARE ROOT OF SIGMA PREFIXED BY THE SIGN
874 C* OF THE REPLACED SUBDIAGONAL ELEMENT.
876 C* THE ABOVE STEPS ARE REPEATED ON FURTHER ROWS OF THE
877 C* TRANSFORMED A IN REVERSE ORDER UNTIL A IS REDUCED TO TRI-
878 C* DIAGONAL FORM, THAT IS, REPEATED FOR J = N-1,N-2,...,3.
885 C* THE ORDER OF THE MATRIX.
887 C* MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A
888 C* AS DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT
889 C* A - W.P. REAL (NV)
890 C* CONTAINS THE LOWER TRIANGLE OF THE REAL SYMMETRIC
891 C* INPUT MATRIX, STORED ROW-WISE AS A ONE-DIMENSIONAL
892 C* ARRAY, IN ITS FIRST N*(N+1)/2 POSITIONS.
895 C* A - W.P. REAL (NV)
896 C* CONTAINS INFORMATION ABOUT THE ORTHOGONAL
897 C* TRANSFORMATIONS USED IN THE REDUCTION.
899 C* CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL
902 C* CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
903 C* MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO
904 C* E2 - W.P. REAL (N)
905 C* CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF
906 C* E. MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
908 C* DIFFERENCES FROM EISPACK 3 -
909 C* OUTER LOOP CHANGED FROM II=1,N TO I=N,3,-1
910 C* PARAMETER STATEMENT AND GENERIC INTRINSIC FUNCTIONS USED
911 C* SCALE.NE.0 TEST NOW SPOTS TRI-DIAGONAL FORM
912 C* VALUES LESS THAN EPSLON CLEARED TO ZERO
914 C* U NOT COPIED TO D, LEFT IN A
915 C* E2 COMPUTED FROM E
916 C* INNER LOOPS SPLIT INTO ROUTINES ELAU AND FREDA
917 C* INVERSE OF H STORED INSTEAD OF H
920 C* QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
921 C* B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
923 INTEGER I,IIA,IZ0,L,N,NV
925 DOUBLE PRECISION A(NV),D(N),E(N),E2(N)
926 DOUBLE PRECISION AIIMAX,F,G,H,HROOT,SCALE,SCALEI
927 DOUBLE PRECISION DASUM, DNRM2
928 DOUBLE PRECISION ONE, ZERO
930 PARAMETER (ZERO = 0.0D+00, ONE = 1.0D+00)
932 C-----------------------------------------------------------------------
934 IF (N .LE. 2) GO TO 310
941 AIIMAX = MAX(AIIMAX, ABS(A(IIA)))
942 SCALE = DASUM (L, A(IZ0+1), 1)
943 IF(SCALE .EQ. ABS(A(IIA-1)) .OR. AIIMAX+SCALE .EQ. AIIMAX) THEN
945 C THIS ROW IS ALREADY IN TRI-DIAGONAL FORM
948 IF (AIIMAX+D(I) .EQ. AIIMAX) D(I) = ZERO
950 IF (AIIMAX+E(I) .EQ. AIIMAX) E(I) = ZERO
958 CALL DSCAL(L,SCALEI,A(IZ0+1),1)
959 HROOT = DNRM2(L,A(IZ0+1),1)
965 H = HROOT*HROOT - F * G
968 A(IIA) = ONE / SQRT(H)
969 C .......... FORM P THEN Q IN E(1:L) ..........
970 CALL ELAU(ONE/H,L,A(IZ0+1),A,E)
971 C .......... FORM REDUCED A ..........
972 CALL FREDA(L,A(IZ0+1),A,E)
986 C*MODULE EIGEN *DECK EVVRSP
987 SUBROUTINE EVVRSP(MSGFL,N,NVECT,LENA,NV,A,B,IND,ROOT,
990 C* AUTHOR: S. T. ELBERT, AMES LABORATORY-USDOE, JUNE 1985
993 C* FINDS (ALL) EIGENVALUES AND (SOME OR ALL) EIGENVECTORS
995 C* OF A REAL SYMMETRIC PACKED MATRIX.
999 C* THE METHOD AS PRESENTED IN THIS ROUTINE CONSISTS OF FOUR STEPS:
1000 C* FIRST, THE INPUT MATRIX IS REDUCED TO TRIDIAGONAL FORM BY THE
1001 C* HOUSEHOLDER TECHNIQUE (ORTHOGONAL SIMILARITY TRANSFORMATIONS).
1002 C* SECOND, THE ROOTS ARE LOCATED USING THE RATIONAL QL METHOD.
1003 C* THIRD, THE VECTORS OF THE TRIDIAGONAL FORM ARE EVALUATED BY THE
1004 C* INVERSE ITERATION TECHNIQUE. VECTORS FOR DEGENERATE OR NEAR-
1005 C* DEGENERATE ROOTS ARE FORCED TO BE ORTHOGONAL.
1006 C* FOURTH, THE TRIDIAGONAL VECTORS ARE ROTATED TO VECTORS OF THE
1009 C* THESE ROUTINES ARE MODIFICATIONS OF THE EISPACK 3
1010 C* ROUTINES TRED3, TQLRAT, TINVIT AND TRBAK3
1012 C* FOR FURTHER DETAILS, SEE EISPACK USERS GUIDE, B. T. SMITH
1013 C* ET AL, SPRINGER-VERLAG, LECTURE NOTES IN COMPUTER SCIENCE,
1014 C* VOL. 6, 2-ND EDITION, 1976. ANOTHER GOOD REFERENCE IS
1015 C* THE SYMMETRIC EIGENVALUE PROBLEM BY B. N. PARLETT
1016 C* PUBLISHED BY PRENTICE-HALL, INC., ENGLEWOOD CLIFFS, N.J. (1980)
1019 C* MSGFL - INTEGER (LOGICAL UNIT NO.)
1020 C* FILE WHERE ERROR MESSAGES WILL BE PRINTED.
1021 C* IF MSGFL IS 0, ERROR MESSAGES WILL BE PRINTED ON LU 6.
1022 C* IF MSGFL IS NEGATIVE, NO ERROR MESSAGES PRINTED.
1024 C* ORDER OF MATRIX A.
1026 C* NUMBER OF VECTORS DESIRED. 0 .LE. NVECT .LE. N.
1028 C* DIMENSION OF A IN CALLING ROUTINE. MUST NOT BE LESS
1031 C* ROW DIMENSION OF VECT IN CALLING ROUTINE. N .LE. NV.
1032 C* A - WORKING PRECISION REAL (LENA)
1033 C* INPUT MATRIX, ROWS OF THE LOWER TRIANGLE PACKED INTO
1034 C* LINEAR ARRAY OF DIMENSION N*(N+1)/2. THE PACKED ORDER
1035 C* IS A(1,1), A(2,1), A(2,2), A(3,1), A(3,2), ...
1036 C* B - WORKING PRECISION REAL (N,8)
1037 C* SCRATCH ARRAY, 8*N ELEMENTS
1038 C* IND - INTEGER (N)
1039 C* SCRATCH ARRAY OF LENGTH N.
1041 C* ROOT ORDERING FLAG.
1042 C* = 0, ROOTS WILL BE PUT IN ASCENDING ORDER.
1043 C* = 2, ROOTS WILL BE PUT IN DESCENDING ORDER.
1046 C* A - DESTORYED. NOW HOLDS REFLECTION OPERATORS.
1047 C* ROOT - WORKING PRECISION REAL (N)
1048 C* ALL EIGENVALUES IN ASCENDING OR DESCENDING ORDER.
1049 C* IF IORDER = 0, ROOT(1) .LE. ... .LE. ROOT(N)
1050 C* IF IORDER = 2, ROOT(1) .GE. ... .GE. ROOT(N)
1051 C* VECT - WORKING PRECISION REAL (NV,NVECT)
1052 C* EIGENVECTORS FOR ROOT(1), ..., ROOT(NVECT).
1054 C* = 0 IF NO ERROR DETECTED,
1055 C* = K IF ITERATION FOR K-TH EIGENVALUE FAILED,
1056 C* = -K IF ITERATION FOR K-TH EIGENVECTOR FAILED.
1057 C* (FAILURES SHOULD BE VERY RARE. CONTACT C. MOLER.)
1060 LOGICAL GOPARR,DSKWRK,MASWRK
1062 DOUBLE PRECISION A(LENA)
1063 DOUBLE PRECISION B(N,8)
1064 DOUBLE PRECISION ROOT(N)
1066 DOUBLE PRECISION VECT(NV,*)
1070 COMMON /PAR / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
1072 900 FORMAT(26H0*** EVVRSP PARAMETERS ***/
1073 + 14H *** N = ,I8,4H ***/
1074 + 14H *** NVECT = ,I8,4H ***/
1075 + 14H *** LENA = ,I8,4H ***/
1076 + 14H *** NV = ,I8,4H ***/
1077 + 14H *** IORDER = ,I8,4H ***/
1078 + 14H *** IERR = ,I8,4H ***)
1079 901 FORMAT(37H VALUE OF LENA IS LESS THAN (N*N+N)/2)
1080 902 FORMAT(39H EQLRAT HAS FAILED TO CONVERGE FOR ROOT,I5)
1081 903 FORMAT(18H NV IS LESS THAN N)
1082 904 FORMAT(41H EINVIT HAS FAILED TO CONVERGE FOR VECTOR,I5)
1083 905 FORMAT(51H VALUE OF IORDER MUST BE 0 (SMALLEST ROOT FIRST) OR
1084 * ,23H 2 (LARGEST ROOT FIRST))
1085 906 FORMAT(' VALUE OF N IS LESS THAN OR EQUAL ZERO')
1087 C-----------------------------------------------------------------------
1090 IF (MSGFL .EQ. 0) LMSGFL=6
1092 IF (N .LE. 0) GO TO 800
1094 IF ( (N*N+N)/2 .GT. LENA) GO TO 810
1096 C REDUCE REAL SYMMETRIC MATRIX A TO TRIDIAGONAL FORM
1098 CALL ETRED3(N,LENA,A,B(1,1),B(1,2),B(1,3))
1100 C FIND ALL EIGENVALUES OF TRIDIAGONAL MATRIX
1102 CALL EQLRAT(N,B(1,1),B(1,2),B(1,3),ROOT,IND,IERR,B(1,4))
1103 IF (IERR .NE. 0) GO TO 820
1105 C CHECK THE DESIRED ORDER OF THE EIGENVALUES
1108 IF (IORDER .EQ. 0) GO TO 300
1109 IF (IORDER .NE. 2) GO TO 850
1111 C ORDER ROOTS IN DESCENDING ORDER (LARGEST FIRST)...
1112 C TURN ROOT AND IND ARRAYS END FOR END
1124 C FIND I AND J MARKING THE START AND END OF A SEQUENCE
1125 C OF DEGENERATE ROOTS
1130 IF (I .GT. N) GO TO 300
1132 IF (ROOT(J) .NE. ROOT(I)) GO TO 240
1137 IF (J .EQ. I) GO TO 220
1139 C TURN AROUND IND BETWEEN I AND J
1155 IF (NVECT .LE. 0) RETURN
1156 IF (NV .LT. N) GO TO 830
1158 C FIND EIGENVECTORS OF TRI-DIAGONAL MATRIX VIA INVERSE ITERATION
1161 CALL EINVIT(NV,N,B(1,1),B(1,2),B(1,3),NVECT,ROOT,IND,
1162 + VECT,IERR,B(1,4),B(1,5),B(1,6),B(1,7),B(1,8))
1163 IF (IERR .NE. 0) GO TO 840
1165 C FIND EIGENVECTORS OF SYMMETRIC MATRIX VIA BACK TRANSFORMATION
1168 CALL ETRBK3(NV,N,LENA,A,NVECT,VECT)
1171 C ERROR MESSAGE SECTION
1173 800 IF (LMSGFL .LT. 0) RETURN
1174 IF (MASWRK) WRITE(LMSGFL,906)
1177 810 IF (LMSGFL .LT. 0) RETURN
1178 IF (MASWRK) WRITE(LMSGFL,901)
1181 820 IF (LMSGFL .LT. 0) RETURN
1182 IF (MASWRK) WRITE(LMSGFL,902) IERR
1185 830 IF (LMSGFL .LT. 0) RETURN
1186 IF (MASWRK) WRITE(LMSGFL,903)
1190 IF ((LMSGFL .GT. 0).AND.MASWRK) WRITE(LMSGFL,904) -IERR
1194 IF (LMSGFL .LT. 0) RETURN
1195 IF (MASWRK) WRITE(LMSGFL,905)
1199 IF (MASWRK) WRITE(LMSGFL,900) N,NVECT,LENA,NV,IORDER,IERR
1202 C*MODULE EIGEN *DECK FREDA
1203 SUBROUTINE FREDA(L,D,A,E)
1205 DOUBLE PRECISION A(*)
1206 DOUBLE PRECISION D(L)
1207 DOUBLE PRECISION E(L)
1213 C .......... FORM REDUCED A ..........
1220 A(JK) = A(JK) - F * E(K) - G * D(K)
1227 C*MODULE EIGEN *DECK GIVEIS
1228 SUBROUTINE GIVEIS(N,NVECT,NV,A,B,INDB,ROOT,VECT,IERR)
1229 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1230 DIMENSION A(*),B(N,8),INDB(N),ROOT(N),VECT(NV,NVECT)
1232 C EISPACK-BASED SUBSTITUTE FOR QCPE ROUTINE GIVENS.
1233 C FINDS ALL EIGENVALUES AND SOME EIGENVECTORS OF A REAL SYMMETRIC
1234 C MATRIX. AUTHOR.. C. MOLER AND D. SPANGLER, N.R.C.C., 4/1/79.
1237 C N = ORDER OF MATRIX .
1238 C NVECT = NUMBER OF VECTORS DESIRED. 0 .LE. NVECT .LE. N .
1239 C NV = LEADING DIMENSION OF VECT .
1240 C A = INPUT MATRIX, COLUMNS OF THE UPPER TRIANGLE PACKED INTO
1241 C LINEAR ARRAY OF DIMENSION N*(N+1)/2 .
1242 C B = SCRATCH ARRAY, 8*N ELEMENTS (NOTE THIS IS MORE THAN
1243 C PREVIOUS VERSIONS OF GIVENS.)
1244 C IND = INDEX ARRAY OF N ELEMENTS
1248 C ROOT = ALL EIGENVALUES, ROOT(1) .LE. ... .LE. ROOT(N) .
1249 C (FOR OTHER ORDERINGS, SEE BELOW.)
1250 C VECT = EIGENVECTORS FOR ROOT(1),..., ROOT(NVECT) .
1251 C IERR = 0 IF NO ERROR DETECTED,
1252 C = K IF ITERATION FOR K-TH EIGENVALUE FAILED,
1253 C = -K IF ITERATION FOR K-TH EIGENVECTOR FAILED.
1254 C (FAILURES SHOULD BE VERY RARE. CONTACT MOLER.)
1256 C CALLS MODIFIED EISPACK ROUTINES TRED3B, IMTQLV, TINVTB, AND
1257 C TRBK3B. THE ROUTINES TRED3B, TINVTB, AND TRBK3B.
1258 C THE ORIGINAL EISPACK ROUTINES TRED3, TINVIT, AND TRBAK3
1259 C WERE MODIFIED BY THE INTRODUCTION OF TWO ROUTINES FROM THE
1260 C BLAS LIBRARY - DDOT AND DAXPY.
1262 C IF TINVIT FAILS TO CONVERGE, TQL2 IS CALLED
1264 C SEE EISPACK USERS GUIDE, B. T. SMITH ET AL, SPRINGER-VERLAG
1265 C LECTURE NOTES IN COMPUTER SCIENCE, VOL. 6, 2-ND EDITION, 1976 .
1266 C NOTE THAT IMTQLV AND TINVTB HAVE INTERNAL MACHINE
1267 C DEPENDENT CONSTANTS.
1269 DATA ONE, ZERO /1.0D+00, 0.0D+00/
1270 CALL TRED3B(N,(N*N+N)/2,A,B(1,1),B(1,2),B(1,3))
1271 CALL IMTQLV(N,B(1,1),B(1,2),B(1,3),ROOT,INDB,IERR,B(1,4))
1272 IF (IERR .NE. 0) RETURN
1274 C TO REORDER ROOTS...
1284 IF (NVECT .LE. 0) RETURN
1285 CALL TINVTB(NV,N,B(1,1),B(1,2),B(1,3),NVECT,ROOT,INDB,VECT,IERR,
1286 + B(1,4),B(1,5),B(1,6),B(1,7),B(1,8))
1287 IF (IERR .EQ. 0) GO TO 160
1289 C IF INVERSE ITERATION GIVES AN ERROR IN DETERMINING THE
1290 C EIGENVECTORS, TRY THE QL ALGORITHM IF ALL THE EIGENVECTORS
1293 IF (NVECT .NE. N) RETURN
1300 CALL TQL2 (NV,N,B(1,1),B(1,2),VECT,IERR)
1304 IF (IERR .NE. 0) RETURN
1305 160 CALL TRBK3B(NV,N,(N*N+N)/2,A,NVECT,VECT)
1308 C*MODULE EIGEN *DECK GLDIAG
1309 SUBROUTINE GLDIAG(LDVECT,NVECT,N,H,WRK,EIG,VECTOR,IERR,IWRK)
1311 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1313 LOGICAL GOPARR,DSKWRK,MASWRK
1315 DIMENSION H(*),WRK(N,8),EIG(N),VECTOR(LDVECT,NVECT),IWRK(N)
1317 COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(400)
1318 COMMON /MACHSW/ KDIAG,ICORFL,IXDR
1319 COMMON /PAR / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
1321 C ----- GENERAL ROUTINE TO DIAGONALIZE A SYMMETRIC MATRIX -----
1322 C IF KDIAG = 0, USE A ROUTINE FROM THE VECTOR LIBRARY,
1323 C IF AVAILABLE (SEE THE SUBROUTINE 'GLDIAG'
1324 C IN VECTOR.SRC), OR EVVRSP OTHERWISE
1329 C N = DIMENSION (ORDER) OF MATRIX TO BE SOLVED
1330 C LDVECT = LEADING DIMENSION OF VECTOR
1331 C NVECT = NUMBER OF VECTORS DESIRED
1332 C H = MATRIX TO BE DIAGONALIZED
1333 C WRK = N*8 W.P. REAL WORDS OF SCRATCH SPACE
1334 C EIG = EIGENVECTORS (OUTPUT)
1335 C VECTOR = EIGENVECTORS (OUTPUT)
1336 C IERR = ERROR FLAG (OUTPUT)
1337 C IWRK = N INTEGER WORDS OF SCRATCH SPACE
1341 C ----- USE STEVE ELBERT'S ROUTINE -----
1343 IF(KDIAG.LE.1 .OR. KDIAG.GT.3) THEN
1346 CALL EVVRSP(IW,N,NVECT,LENH,LDVECT,H,WRK,IWRK,EIG,VECTOR
1350 C ----- USE MODIFIED EISPAK ROUTINE -----
1353 * CALL GIVEIS(N,NVECT,LDVECT,H,WRK,IWRK,EIG,VECTOR,IERR)
1355 C ----- USE JACOBI ROTATION ROUTINE -----
1359 CALL JACDG(H,VECTOR,EIG,IWRK,WRK,LDVECT,N)
1361 IF (MASWRK) WRITE(IW,9000) N,NVECT,LDVECT
1367 9000 FORMAT(1X,'IN -GLDIAG-, N,NVECT,LDVECT=',3I8/
1368 * 1X,'THE JACOBI CODE CANNOT COPE WITH N.NE.NVECT!'/
1369 * 1X,'SO THIS RUN DOES NOT PERMIT KDIAG=3.')
1371 C*MODULE EIGEN *DECK IMTQLV
1372 SUBROUTINE IMTQLV(N,D,E,E2,W,IND,IERR,RV1)
1373 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1375 DOUBLE PRECISION MACHEP
1376 DIMENSION D(N),E(N),E2(N),W(N),RV1(N),IND(N)
1378 C THIS ROUTINE IS A VARIANT OF IMTQL1 WHICH IS A TRANSLATION OF
1379 C ALGOL PROCEDURE IMTQL1, NUM. MATH. 12, 377-383(1968) BY MARTIN AND
1380 C WILKINSON, AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE.
1381 C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971).
1383 C THIS ROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC TRIDIAGONAL
1384 C MATRIX BY THE IMPLICIT QL METHOD AND ASSOCIATES WITH THEM
1385 C THEIR CORRESPONDING SUBMATRIX INDICES.
1389 C N IS THE ORDER OF THE MATRIX,
1391 C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
1393 C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
1394 C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY,
1396 C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
1397 C E2(1) IS ARBITRARY.
1401 C D AND E ARE UNALTERED,
1403 C ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED
1404 C AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE
1405 C MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES.
1406 C E2(1) IS ALSO SET TO ZERO,
1408 C W CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN
1409 C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND
1410 C ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE
1411 C THE SMALLEST EIGENVALUES,
1413 C IND CONTAINS THE SUBMATRIX INDICES ASSOCIATED WITH THE
1414 C CORRESPONDING EIGENVALUES IN W -- 1 FOR EIGENVALUES
1415 C BELONGING TO THE FIRST SUBMATRIX FROM THE TOP,
1416 C 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.,
1419 C ZERO FOR NORMAL RETURN,
1420 C J IF THE J-TH EIGENVALUE HAS NOT BEEN
1421 C DETERMINED AFTER 30 ITERATIONS,
1423 C RV1 IS A TEMPORARY STORAGE ARRAY.
1425 C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
1426 C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
1428 C ------------------------------------------------------------------
1430 C ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
1431 C THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
1434 MACHEP = 2.0D+00**(-50)
1442 IF (I .NE. 1) RV1(I-1) = E(I)
1450 C ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT **********
1452 IF (M .EQ. N) GO TO 160
1453 IF (ABS(RV1(M)) .LE. MACHEP * (ABS(W(M)) + ABS(W(M+1)))) GO TO
1455 C ********** GUARD AGAINST UNDERFLOWED ELEMENT OF E2 **********
1456 IF (E2(M+1) .EQ. 0.0D+00) GO TO 180
1459 160 IF (M .LE. K) GO TO 200
1460 IF (M .NE. N) E2(M+1) = 0.0D+00
1464 IF (M .EQ. L) GO TO 280
1465 IF (J .EQ. 30) GO TO 380
1467 C ********** FORM SHIFT **********
1468 G = (W(L+1) - P) / (2.0D+00 * RV1(L))
1469 R = SQRT(G*G+1.0D+00)
1470 G = W(M) - P + RV1(L) / (G + SIGN(R,G))
1475 C ********** FOR I=M-1 STEP -1 UNTIL L DO -- **********
1480 IF (ABS(F) .LT. ABS(G)) GO TO 220
1482 R = SQRT(C*C+1.0D+00)
1488 R = SQRT(S*S+1.0D+00)
1493 R = (W(I) - G) * S + 2.0D+00 * C * B
1503 C ********** ORDER EIGENVALUES **********
1504 280 IF (L .EQ. 1) GO TO 320
1505 C ********** FOR I=L STEP -1 UNTIL 2 DO -- **********
1508 IF (P .GE. W(I-1)) GO TO 340
1519 C ********** SET ERROR -- NO CONVERGENCE TO AN
1520 C EIGENVALUE AFTER 30 ITERATIONS **********
1523 C ********** LAST CARD OF IMTQLV **********
1525 C*MODULE EIGEN *DECK JACDG
1526 SUBROUTINE JACDG(A,VEC,EIG,JBIG,BIG,LDVEC,N)
1528 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1530 DIMENSION A(*),VEC(LDVEC,N),EIG(N),JBIG(N),BIG(N)
1532 PARAMETER (ONE=1.0D+00)
1534 C ----- JACOBI DIAGONALIZATION OF SYMMETRIC MATRIX -----
1535 C SYMMETRIC MATRIX -A- OF DIMENSION -N- IS DESTROYED ON EXIT.
1536 C ALL EIGENVECTORS ARE FOUND, SO -VEC- MUST BE SQUARE,
1537 C UNLESS SOMEONE TAKES THE TROUBLE TO LOOK AT -NMAX- BELOW.
1538 C -BIG- AND -JBIG- ARE SCRATCH WORK ARRAYS.
1540 CALL VCLR(VEC,1,LDVEC*N)
1546 NB2 = (NB1*NB1+NB1)/2
1550 CALL JACDIA(A,VEC,NB1,NB2,LDVEC,NMIN,NMAX,BIG,JBIG)
1553 EIG(I) = A((I*I+I)/2)
1556 CALL JACORD(VEC,EIG,NB1,LDVEC)
1559 C*MODULE EIGEN *DECK JACDIA
1560 SUBROUTINE JACDIA(F,VEC,NB1,NB2,LDVEC,NMIN,NMAX,BIG,JBIG)
1561 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1562 LOGICAL GOPARR,DSKWRK,MASWRK
1563 DIMENSION F(NB2),VEC(LDVEC,NB1),BIG(NB1),JBIG(NB1)
1565 COMMON /PAR / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
1567 PARAMETER (ROOT2=0.707106781186548D+00 )
1568 PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, D1050=1.05D+00,
1569 * D1500=1.5D+00, D3875=3.875D+00,
1570 * D0500=0.5D+00, D1375=1.375D+00, D0250=0.25D+00 )
1571 PARAMETER (C2=1.0D-12, C3=4.0D-16,
1572 * C4=2.0D-16, C5=8.0D-09, C6=3.0D-06 )
1574 C F IS THE MATRIX TO BE DIAGONALIZED, F IS STORED TRIANGULAR
1575 C VEC IS THE ARRAY OF EIGENVECTORS, DIMENSION NB1*NB1
1576 C BIG AND JBIG ARE TEMPORARY SCRATCH AREAS OF DIMENSION NB1
1577 C THE ROTATIONS AMONG THE FIRST NMIN BASIS FUNCTIONS ARE NOT
1579 C THE ROTATIONS AMONG THE LAST NB1-NMAX BASIS FUNCTIONS ARE NOT
1585 EPS = 64.0D+00*EPSLON(ONE)
1587 C LOOP OVER COLUMNS (K) OF TRIANGULAR MATRIX TO DETERMINE
1588 C LARGEST OFF-DIAGONAL ELEMENTS IN ROW(I).
1593 IF(I.LT.NMIN .OR. I.EQ.1) GO TO 20
1597 IF(ABS(BIG(I)).GE.ABS(F(II+K))) GO TO 10
1603 C ----- 2X2 JACOBI ITERATIONS BEGIN HERE -----
1605 MAXIT=MAX(NB2*20,500)
1610 C FIND SMALLEST DIAGONAL ELEMENT
1616 SD= MIN(SD,ABS(F(JJ)))
1618 TEST = MAX(EPS, C2*MAX(SD,C6))
1620 C FIND LARGEST OFF-DIAGONAL ELEMENT
1626 IF(T.GE.ABS(BIG(I))) GO TO 50
1631 C TEST FOR CONVERGENCE, THEN DETERMINE ROTATION.
1633 IF(T.LT.TEST) RETURN
1636 IF(ITER.GT.MAXIT) THEN
1638 WRITE(6,*) 'JACOBI DIAGONALIZATION FAILS, DIMENSION=',NB1
1639 WRITE(6,9020) ITER,T,TEST,SD
1648 DIF=F(IAA+IA)-F(IBB+IB)
1649 IF(ABS(DIF).GT.C3*T) GO TO 70
1655 IF(T2X25 . GT . C4) GO TO 80
1659 80 IF(T2X25 . GT . C5) GO TO 90
1660 SX=T2X2*(ONE-D1500*T2X25)
1663 90 IF(T2X25 . GT . C6) GO TO 100
1664 CX=ONE+T2X25*(T2X25*D1375 - D0500)
1665 SX= T2X2*(ONE + T2X25*(T2X25*D3875 - D1500))
1667 100 T=D0250 / SQRT(D0250 + T2X25)
1669 SX= SIGN( SQRT(D0500 - T),T2X2)
1675 F(IEAR)=F(IEAR)*CX+F(IEBR)*SX
1676 F(IEBR)=T-F(IEBR)*CX
1677 IF(IR-IA) 220,120,130
1683 IF(JBIG(IR)) 200,220,200
1687 IF(IR-IB) 180,150,160
1688 150 F(IEAA)=F(IEAA)*CX+F(IEAB)*SX
1689 F(IEAB)=TT*CX+F(IEBR)*SX
1690 F(IEBR)=TT*SX-F(IEBR)*CX
1693 160 IF( ABS(T) . GE . ABS(F(IEBR))) GO TO 170
1694 IF(IB.GT.NMAX) GO TO 170
1698 180 IF( ABS(T) . LT . ABS(BIG(IR))) GO TO 190
1702 190 IF(IA . NE . JBIG(IR) . AND . IB . NE . JBIG(IR)) GO TO 220
1708 IF(ABS(BIG(IR)) . GE . ABS(F(K))) GO TO 210
1716 T1=VEC(I,IA)*CX + VEC(I,IB)*SX
1717 T2=VEC(I,IA)*SX - VEC(I,IB)*CX
1723 9020 FORMAT(1X,'ITER=',I6,' T,TEST,SD=',1P,3E20.10)
1725 C*MODULE EIGEN *DECK JACORD
1726 SUBROUTINE JACORD(VEC,EIG,N,LDVEC)
1727 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1728 DIMENSION VEC(LDVEC,N),EIG(N)
1730 C ---- SORT EIGENDATA INTO ASCENDING ORDER -----
1735 IF (EIG(J) .LT. EIG(JJ)) JJ = J
1737 IF (JJ .EQ. I) GO TO 290
1743 VEC(J,JJ) = VEC(J,I)
1749 C*MODULE EIGEN *DECK TINVTB
1750 SUBROUTINE TINVTB(NM,N,D,E,E2,M,W,IND,Z,
1751 * IERR,RV1,RV2,RV3,RV4,RV6)
1752 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1753 DIMENSION D(N),E(N),E2(N),W(M),Z(NM,M),
1754 * RV1(N),RV2(N),RV3(N),RV4(N),RV6(N),IND(M)
1755 DOUBLE PRECISION MACHEP,NORM
1756 INTEGER P,Q,R,S,TAG,GROUP
1757 C ------------------------------------------------------------------
1759 C THIS ROUTINE IS A TRANSLATION OF THE INVERSE ITERATION TECH-
1760 C NIQUE IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON.
1761 C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971).
1763 C THIS ROUTINE FINDS THOSE EIGENVECTORS OF A TRIDIAGONAL
1764 C SYMMETRIC MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES,
1765 C USING INVERSE ITERATION.
1769 C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
1770 C ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
1771 C DIMENSION STATEMENT,
1773 C N IS THE ORDER OF THE MATRIX,
1775 C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
1777 C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
1778 C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY,
1780 C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E,
1781 C WITH ZEROS CORRESPONDING TO NEGLIGIBLE ELEMENTS OF E.
1782 C E(I) IS CONSIDERED NEGLIGIBLE IF IT IS NOT LARGER THAN
1783 C THE PRODUCT OF THE RELATIVE MACHINE PRECISION AND THE SUM
1784 C OF THE MAGNITUDES OF D(I) AND D(I-1). E2(1) MUST CONTAIN
1785 C 0.0 IF THE EIGENVALUES ARE IN ASCENDING ORDER, OR 2.0
1786 C IF THE EIGENVALUES ARE IN DESCENDING ORDER. IF BISECT,
1787 C TRIDIB, OR IMTQLV HAS BEEN USED TO FIND THE EIGENVALUES,
1788 C THEIR OUTPUT E2 ARRAY IS EXACTLY WHAT IS EXPECTED HERE,
1790 C M IS THE NUMBER OF SPECIFIED EIGENVALUES,
1792 C W CONTAINS THE M EIGENVALUES IN ASCENDING OR DESCENDING ORDER,
1794 C IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES
1795 C ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W --
1796 C 1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM
1797 C THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.
1801 C ALL INPUT ARRAYS ARE UNALTERED,
1803 C Z CONTAINS THE ASSOCIATED SET OF ORTHONORMAL EIGENVECTORS.
1804 C ANY VECTOR WHICH FAILS TO CONVERGE IS SET TO ZERO,
1807 C ZERO FOR NORMAL RETURN,
1808 C -R IF THE EIGENVECTOR CORRESPONDING TO THE R-TH
1809 C EIGENVALUE FAILS TO CONVERGE IN 5 ITERATIONS,
1811 C RV1, RV2, RV3, RV4, AND RV6 ARE TEMPORARY STORAGE ARRAYS.
1813 C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
1814 C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
1816 C ------------------------------------------------------------------
1818 C ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
1819 C THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
1822 MACHEP = 2.0D+00**(-50)
1825 IF (M .EQ. 0) GO TO 680
1827 ORDER = 1.0D+00 - E2(1)
1837 C ********** ESTABLISH AND PROCESS NEXT SUBMATRIX **********
1842 IF (Q .EQ. N) GO TO 140
1843 IF (E2(Q+1) .EQ. 0.0D+00) GO TO 140
1845 C ********** FIND VECTORS BY INVERSE ITERATION **********
1851 IF (IND(R) .NE. TAG) GO TO 660
1854 IF (S .NE. 0) GO TO 220
1855 C ********** CHECK FOR ISOLATED ROOT **********
1857 IF (P .NE. Q) GO TO 160
1860 160 NORM = ABS(D(P))
1863 180 NORM = NORM + ABS(D(I)) + ABS(E(I))
1864 C ********** EPS2 IS THE CRITERION FOR GROUPING,
1865 C EPS3 REPLACES ZERO PIVOTS AND EQUAL
1866 C ROOTS ARE MODIFIED BY EPS3,
1867 C EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW **********
1868 EPS2 = 1.0D-03 * NORM
1869 EPS3 = MACHEP * NORM
1872 UK = EPS4 / SQRT(UK)
1876 C ********** LOOK FOR CLOSE OR COINCIDENT ROOTS **********
1877 220 IF (ABS(X1-X0) .GE. EPS2) GO TO 200
1879 IF (ORDER * (X1 - X0) .LE. 0.0D+00) X1 = X0 + ORDER * EPS3
1880 C ********** ELIMINATION WITH INTERCHANGES AND
1881 C INITIALIZATION OF VECTOR **********
1886 IF (I .EQ. P) GO TO 280
1887 IF (ABS(E(I)) .LT. ABS(U)) GO TO 260
1888 C ********** WARNING -- A DIVIDE CHECK MAY OCCUR HERE IF
1889 C E2 ARRAY HAS NOT BEEN SPECIFIED CORRECTLY **********
1893 RV2(I-1) = D(I) - X1
1895 IF (I .NE. Q) RV3(I-1) = E(I+1)
1896 U = V - XU * RV2(I-1)
1904 280 U = D(I) - X1 - XU * V
1905 IF (I .NE. Q) V = E(I+1)
1908 IF (U .EQ. 0.0D+00) U = EPS3
1912 C ********** BACK SUBSTITUTION
1913 C FOR I=Q STEP -1 UNTIL P DO -- **********
1914 320 DO 340 II = P, Q
1916 RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I)
1920 C ********** ORTHOGONALIZE WITH RESPECT TO PREVIOUS
1921 C MEMBERS OF GROUP **********
1922 IF (GROUP .EQ. 0) GO TO 400
1925 DO 380 JJ = 1, GROUP
1927 IF (IND(J) .NE. TAG) GO TO 360
1928 XU = DDOT(IQMP,RV6(P),1,Z(P,J),1)
1930 CALL DAXPY(IQMP,-XU,Z(P,J),1,RV6(P),1)
1937 420 NORM = NORM + ABS(RV6(I))
1939 IF (NORM .GE. 1.0D+00) GO TO 560
1940 C ********** FORWARD SUBSTITUTION **********
1941 IF (ITS .EQ. 5) GO TO 540
1942 IF (NORM .NE. 0.0D+00) GO TO 440
1947 440 XU = EPS4 / NORM
1950 460 RV6(I) = RV6(I) * XU
1951 C ********** ELIMINATION OPERATIONS ON NEXT VECTOR
1952 C ITERATE **********
1953 480 DO 520 I = IP, Q
1955 C ********** IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE
1956 C WAS PERFORMED EARLIER IN THE
1957 C TRIANGULARIZATION PROCESS **********
1958 IF (RV1(I-1) .NE. E(I)) GO TO 500
1961 500 RV6(I) = U - RV4(I) * RV6(I-1)
1966 C ********** SET ERROR -- NON-CONVERGED EIGENVECTOR **********
1970 C ********** NORMALIZE SO THAT SUM OF SQUARES IS
1971 C 1 AND EXPAND TO FULL ORDER **********
1975 RV6(I) = RV6(I) / NORM
1976 580 U = U + RV6(I)**2
1978 XU = 1.0D+00 / SQRT(U)
1981 620 Z(I,R) = 0.0D+00
1984 640 Z(I,R) = RV6(I) * XU
1989 IF (Q .LT. N) GO TO 100
1991 C ********** LAST CARD OF TINVIT **********
1993 C*MODULE EIGEN *DECK TQL2
1995 C ------------------------------------------------------------------
1997 SUBROUTINE TQL2(NM,N,D,E,Z,IERR)
1998 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1999 DOUBLE PRECISION MACHEP
2000 DIMENSION D(N),E(N),Z(NM,N)
2002 C THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL2,
2003 C NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND
2005 C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971).
2007 C THIS ROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
2008 C OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE QL METHOD.
2009 C THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO
2010 C BE FOUND IF TRED2 HAS BEEN USED TO REDUCE THIS
2011 C FULL MATRIX TO TRIDIAGONAL FORM.
2015 C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
2016 C ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
2017 C DIMENSION STATEMENT,
2019 C N IS THE ORDER OF THE MATRIX,
2021 C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
2023 C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
2024 C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY,
2026 C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE
2027 C REDUCTION BY TRED2, IF PERFORMED. IF THE EIGENVECTORS
2028 C OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN
2029 C THE IDENTITY MATRIX.
2033 C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN
2034 C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT
2035 C UNORDERED FOR INDICES 1,2,...,IERR-1,
2037 C E HAS BEEN DESTROYED,
2039 C Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC
2040 C TRIDIAGONAL (OR FULL) MATRIX. IF AN ERROR EXIT IS MADE,
2041 C Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED
2045 C ZERO FOR NORMAL RETURN,
2046 C J IF THE J-TH EIGENVALUE HAS NOT BEEN
2047 C DETERMINED AFTER 30 ITERATIONS.
2049 C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
2050 C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
2052 C ------------------------------------------------------------------
2054 C ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
2055 C THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
2058 MACHEP = 2.0D+00**(-50)
2061 IF (N .EQ. 1) GO TO 400
2072 H = MACHEP * (ABS(D(L)) + ABS(E(L)))
2074 C ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT **********
2076 IF (ABS(E(M)) .LE. B) GO TO 140
2077 C ********** E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
2078 C THROUGH THE BOTTOM OF THE LOOP **********
2081 140 IF (M .EQ. L) GO TO 280
2082 160 IF (J .EQ. 30) GO TO 380
2084 C ********** FORM SHIFT **********
2087 P = (D(L1) - G) / (2.0D+00 * E(L))
2088 R = SQRT(P*P+1.0D+00)
2089 D(L) = E(L) / (P + SIGN(R,P))
2096 C ********** QL TRANSFORMATION **********
2101 C ********** FOR I=M-1 STEP -1 UNTIL L DO -- **********
2106 IF (ABS(P) .LT. ABS(E(I))) GO TO 200
2108 R = SQRT(C*C+1.0D+00)
2114 R = SQRT(C*C+1.0D+00)
2115 E(I+1) = S * E(I) * R
2118 220 P = C * D(I) - S * G
2119 D(I+1) = H + S * (C * G + S * D(I))
2120 C ********** FORM VECTOR **********
2121 CALL DROT(N,Z(1,I+1),1,Z(1,I),1,C,S)
2127 IF (ABS(E(L)) .GT. B) GO TO 160
2130 C ********** ORDER EIGENVALUES AND EIGENVECTORS **********
2137 IF (D(J) .GE. P) GO TO 320
2142 IF (K .EQ. I) GO TO 360
2146 CALL DSWAP(N,Z(1,I),1,Z(1,K),1)
2151 C ********** SET ERROR -- NO CONVERGENCE TO AN
2152 C EIGENVALUE AFTER 30 ITERATIONS **********
2155 C ********** LAST CARD OF TQL2 **********
2157 C*MODULE EIGEN *DECK TRBK3B
2159 C ------------------------------------------------------------------
2161 SUBROUTINE TRBK3B(NM,N,NV,A,M,Z)
2162 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2163 DIMENSION A(NV),Z(NM,M)
2165 C THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRBAK3,
2166 C NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
2167 C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
2169 C THIS ROUTINE FORMS THE EIGENVECTORS OF A REAL SYMMETRIC
2170 C MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
2171 C SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY TRED3B.
2175 C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
2176 C ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
2177 C DIMENSION STATEMENT,
2179 C N IS THE ORDER OF THE MATRIX,
2181 C NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A
2182 C AS DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT,
2184 C A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANSFORMATIONS
2185 C USED IN THE REDUCTION BY TRED3B IN ITS FIRST
2186 C N*(N+1)/2 POSITIONS,
2188 C M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED,
2190 C Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
2191 C IN ITS FIRST M COLUMNS.
2195 C Z CONTAINS THE TRANSFORMED EIGENVECTORS
2196 C IN ITS FIRST M COLUMNS.
2198 C NOTE THAT TRBAK3 PRESERVES VECTOR EUCLIDEAN NORMS.
2200 C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
2201 C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
2203 C ------------------------------------------------------------------
2205 IF (M .EQ. 0) GO TO 140
2206 IF (N .EQ. 1) GO TO 140
2213 IF (H .EQ. 0.0D+00) GO TO 120
2216 S = -DDOT(L,A(IZ+1),1,Z(1,J),1)
2218 C ********** DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW **********
2221 CALL DAXPY(L,S,A(IZ+1),1,Z(1,J),1)
2228 C ********** LAST CARD OF TRBAK3 **********
2230 C*MODULE EIGEN *DECK TRED3B
2232 C ------------------------------------------------------------------
2234 SUBROUTINE TRED3B(N,NV,A,D,E,E2)
2235 IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2236 DIMENSION A(NV),D(N),E(N),E2(N)
2238 C THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED3,
2239 C NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
2240 C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
2242 C THIS ROUTINE REDUCES A REAL SYMMETRIC MATRIX, STORED AS
2243 C A ONE-DIMENSIONAL ARRAY, TO A SYMMETRIC TRIDIAGONAL MATRIX
2244 C USING ORTHOGONAL SIMILARITY TRANSFORMATIONS.
2248 C N IS THE ORDER OF THE MATRIX,
2250 C NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A
2251 C AS DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT,
2253 C A CONTAINS THE LOWER TRIANGLE OF THE REAL SYMMETRIC
2254 C INPUT MATRIX, STORED ROW-WISE AS A ONE-DIMENSIONAL
2255 C ARRAY, IN ITS FIRST N*(N+1)/2 POSITIONS.
2259 C A CONTAINS INFORMATION ABOUT THE ORTHOGONAL
2260 C TRANSFORMATIONS USED IN THE REDUCTION,
2262 C D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX,
2264 C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
2265 C MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO,
2267 C E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
2268 C E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
2270 C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
2271 C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
2273 C ------------------------------------------------------------------
2275 C ********** FOR I=N STEP -1 UNTIL 1 DO -- **********
2282 IF (L .LT. 1) GO TO 120
2283 C ********** SCALE ROW (ALGOL TOL THEN NOT NEEDED) **********
2287 SCALE = SCALE + ABS(D(K))
2290 IF (SCALE .NE. 0.0D+00) GO TO 140
2300 E2(I) = SCALE * SCALE * H
2302 G = -SIGN(SQRT(H),F)
2306 A(IZ) = SCALE * D(L)
2307 IF (L .EQ. 1) GO TO 280
2315 C ********** FORM ELEMENT OF A*U **********
2316 IF (JM1 .EQ. 0) GO TO 200
2318 E(K) = E(K) + DT * A(JK)
2319 G = G + D(K) * A(JK)
2322 200 E(J) = G + A(JK) * DT
2324 C ********** FORM ELEMENT OF P **********
2334 C ********** FORM REDUCED A **********
2342 A(JK) = A(JK) - F * E(K) - G * D(K)
2346 A(IZ+1) = SCALE * SQRT(H)
2350 C ********** LAST CARD OF TRED3 **********