e4088eed6b3ab8cf9551eaa1f8556bb69c8300cb
[unres.git] / source / unres / src_MD / src / eigen.f
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
23 C
24 C*MODULE EIGEN   *DECK EINVIT
25       SUBROUTINE EINVIT(NM,N,D,E,E2,M,W,IND,Z,IERR,RV1,RV2,RV3,RV4,RV6)
26 C*
27 C*    AUTHORS-
28 C*       THIS IS A MODIFICATION OF TINVIT FROM EISPACK EDITION 3
29 C*       DATED AUGUST 1983.
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)
34 C*
35 C*    PURPOSE -
36 C*       THIS ROUTINE FINDS THOSE EIGENVECTORS OF A TRIDIAGONAL
37 C*       SYMMETRIC MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES.
38 C*
39 C*    METHOD -
40 C*       INVERSE ITERATION.
41 C*
42 C*    ON ENTRY -
43 C*       NM     - INTEGER
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.
47 C*       N      - INTEGER
48 C*       D      - W.P. REAL (N)
49 C*                CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
50 C*       E      - W.P. REAL (N)
51 C*                CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
52 C*                IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
53 C*       E2     - W.P. REAL (N)
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.
64 C*       M      - INTEGER
65 C*                THE NUMBER OF SPECIFIED EIGENVECTORS.
66 C*       W      - W.P. REAL (M)
67 C*                CONTAINS THE M EIGENVALUES IN ASCENDING
68 C*                OR DESCENDING ORDER.
69 C*       IND    - INTEGER (M)
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
74 C*                SUBMATRIX, ETC.
75 C*       IERR   - INTEGER (LOGICAL UNIT NUMBER)
76 C*                LOGICAL UNIT FOR ERROR MESSAGES
77 C*
78 C*    ON EXIT -
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.
84 C*       IERR   - INTEGER
85 C*                SET TO
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)
90 C*
91 C*       RV1, RV2, RV3, RV4, AND RV6 ARE TEMPORARY STORAGE ARRAYS.
92 C*
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
103 C*
104 C*    DIFFERENCES FROM EISPACK 3 -
105 C*       EPS3 IS SCALED BY  EPSCAL  (ENHANCES CONVERGENCE, BUT
106 C*          LOWERS ACCURACY)!
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
114 C*       USE LEVEL 1 BLAS
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
119 C*
120 C*    NOTE -
121 C*       QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
122 C*       B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
123 C*
124 C
125       LOGICAL CONVGD,GOPARR,DSKWRK,MASWRK
126 C
127       INTEGER GROUP,I,IERR,ITS,J,JJ,M,N,NM,P,Q,R,S,SUBMAT,TAG
128       INTEGER IND(M)
129 C
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
136 C
137       COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
138 C
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)
141 C
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)')
145 C
146 C-----------------------------------------------------------------------
147 C
148       LUEMSG = IERR
149       IERR = 0
150       X0 = ZERO
151       UK = ZERO
152       NORM = ZERO
153       EPS2 = ZERO
154       EPS3 = ZERO
155       EPS4 = ZERO
156       GROUP = 0
157       TAG = 0
158       ORDER = ONE - E2(1)
159       Q = 0
160       DO 930 SUBMAT = 1, N
161          P = Q + 1
162 C
163 C        .......... ESTABLISH AND PROCESS NEXT SUBMATRIX ..........
164 C
165          DO 120 Q = P, N-1
166             IF (E2(Q+1) .EQ. ZERO) GO TO 140
167   120    CONTINUE
168          Q = N
169 C
170 C        .......... FIND VECTORS BY INVERSE ITERATION ..........
171 C
172   140    CONTINUE
173          TAG = TAG + 1
174          ANORM = ZERO
175          S = 0
176 C
177          DO 920 R = 1, M
178             IF (IND(R) .NE. TAG) GO TO 920
179             ITS = 1
180             X1 = W(R)
181             IF (S .NE. 0) GO TO 510
182 C
183 C        .......... CHECK FOR ISOLATED ROOT ..........
184 C
185             XU = ONE
186             IF (P .EQ. Q) THEN
187                RV6(P) = ONE
188                CONVGD = .TRUE.
189                GO TO 860
190 C
191             END IF
192             NORM = ABS(D(P))
193             DO 500 I = P+1, Q
194                NORM = MAX( NORM, ABS(D(I)) + ABS(E(I)) )
195   500       CONTINUE
196 C
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 .........
201 C
202             EPS2 = GRPTOL * NORM
203             EPS3 = EPSCAL * EPSLON(NORM)
204             UK = Q - P + 1
205             EPS4 = UK * EPS3
206             UK = EPS4 / SQRT(UK)
207             S = P
208             GROUP = 0
209             GO TO 520
210 C
211 C        .......... LOOK FOR CLOSE OR COINCIDENT ROOTS ..........
212 C
213   510       IF (ABS(X1-X0) .GE. EPS2) THEN
214 C
215 C                 ROOTS ARE SEPERATE
216 C
217                GROUP = 0
218             ELSE
219 C
220 C                 ROOTS ARE CLOSE
221 C
222                GROUP = GROUP + 1
223                IF (ORDER * (X1 - X0) .LE. EPS3) X1 = X0 + ORDER * EPS3
224             END IF
225 C
226 C        .......... ELIMINATION WITH INTERCHANGES AND
227 C                   INITIALIZATION OF VECTOR ..........
228 C
229   520       CONTINUE
230 C
231             U = D(P) - X1
232             V = E(P+1)
233             RV6(P) = UK
234             DO 550 I = P+1, Q
235                RV6(I) = UK
236                IF (ABS(E(I)) .GT. ABS(U)) THEN
237 C
238 C                 EXCHANGE ROWS BEFORE ELIMINATION
239 C
240 C                  *** WARNING -- A DIVIDE CHECK MAY OCCUR HERE IF
241 C                      E2 ARRAY HAS NOT BEEN SPECIFIED CORRECTLY .......
242 C
243                   XU = U / E(I)
244                   RV4(I) = XU
245                   RV1(I-1) = E(I)
246                   RV2(I-1) = D(I) - X1
247                   RV3(I-1) = E(I+1)
248                   U = V - XU * RV2(I-1)
249                   V = -XU * RV3(I-1)
250 C
251                ELSE
252 C
253 C                    STRAIGHT ELIMINATION
254 C
255                   XU = E(I) / U
256                   RV4(I) = XU
257                   RV1(I-1) = U
258                   RV2(I-1) = V
259                   RV3(I-1) = ZERO
260                   U = D(I) - X1 - XU * V
261                   V = E(I+1)
262                END IF
263   550       CONTINUE
264 C
265             IF (ABS(U) .LE. EPS3) U = EPS3
266             RV1(Q) = U
267             RV2(Q) = ZERO
268             RV3(Q) = ZERO
269 C
270 C              DO INVERSE ITERATIONS
271 C
272             CONVGD = .FALSE.
273             DO 800 ITS = 1, 5
274                IF (ITS .EQ. 1) GO TO 600
275 C
276 C                    .......... FORWARD SUBSTITUTION ..........
277 C
278                   IF (NORM .EQ. ZERO) THEN
279                      RV6(S) = EPS4
280                      S = S + 1
281                      IF (S .GT. Q) S = P
282                   ELSE
283                      XU = EPS4 / NORM
284                      CALL DSCAL (Q-P+1, XU, RV6(P), 1)
285                   END IF
286 C
287 C                     ... ELIMINATION OPERATIONS ON NEXT VECTOR
288 C
289                   DO 590 I = P+1, Q
290                      U = RV6(I)
291 C
292 C                         IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE
293 C                         WAS PERFORMED EARLIER IN THE
294 C                         TRIANGULARIZATION PROCESS ..........
295 C
296                      IF (RV1(I-1) .EQ. E(I)) THEN
297                         U = RV6(I-1)
298                         RV6(I-1) = RV6(I)
299                      ELSE
300                         U = RV6(I)
301                      END IF
302                      RV6(I) = U - RV4(I) * RV6(I-1)
303   590             CONTINUE
304   600          CONTINUE
305 C
306 C           .......... BACK SUBSTITUTION
307 C
308                RV6(Q) = RV6(Q) / RV1(Q)
309                V = U
310                U = RV6(Q)
311                NORM = ABS(U)
312                DO 620 I = Q-1, P, -1
313                   RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I)
314                   V = U
315                   U = RV6(I)
316                   NORM = NORM + ABS(U)
317   620          CONTINUE
318                IF (GROUP .EQ. 0) GO TO 700
319 C
320 C                 ....... ORTHOGONALIZE WITH RESPECT TO PREVIOUS
321 C                         MEMBERS OF GROUP ..........
322 C
323                   J = R
324                   DO 680 JJ = 1, GROUP
325   630                J = J - 1
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),
328      *                          Z(P,J),1,RV6(P),1)
329   680             CONTINUE
330                   NORM = DASUM(Q-P+1, RV6(P), 1)
331   700          CONTINUE
332 C
333                IF (CONVGD) GO TO 840
334                IF (NORM .GE. ONE) CONVGD = .TRUE.
335   800       CONTINUE
336 C
337 C        .......... NORMALIZE SO THAT SUM OF SQUARES IS
338 C                   1 AND EXPAND TO FULL ORDER ..........
339 C
340   840       CONTINUE
341 C
342             XU = ONE / DNRM2(Q-P+1,RV6(P),1)
343 C
344   860       CONTINUE
345             DO 870 I = 1, P-1
346                Z(I,R) = ZERO
347   870       CONTINUE
348             DO 890 I = P,Q
349                Z(I,R) = RV6(I) * XU
350   890       CONTINUE
351             DO 900 I = Q+1, N
352                Z(I,R) = ZERO
353   900       CONTINUE
354 C
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
359 C
360 C               *** SET ERROR -- NON-CONVERGED EIGENVECTOR ..........
361 C
362                IF (RHO .GT. HUNDRD) IERR = -R
363             END IF
364 C
365             X0 = X1
366   920    CONTINUE
367 C
368          IF (Q .EQ. N) GO TO 940
369   930 CONTINUE
370   940 CONTINUE
371       RETURN
372       END
373 C*MODULE EIGEN   *DECK ELAUM
374       SUBROUTINE ELAU(HINV,L,D,A,E)
375 C
376       DOUBLE PRECISION A(*)
377       DOUBLE PRECISION D(L)
378       DOUBLE PRECISION E(L)
379       DOUBLE PRECISION F
380       DOUBLE PRECISION G
381       DOUBLE PRECISION HALF
382       DOUBLE PRECISION HH
383       DOUBLE PRECISION HINV
384       DOUBLE PRECISION ZERO
385 C
386       PARAMETER (ZERO = 0.0D+00, HALF = 0.5D+00)
387 C
388       JL = L
389       E(1) = A(1) * D(1)
390       JK = 2
391       DO 210 J = 2, JL
392          F = D(J)
393          G = ZERO
394          JM1 = J - 1
395 C
396          DO 200 K = 1, JM1
397             G = G + A(JK) * D(K)
398             E(K) = E(K) + A(JK) * F
399             JK = JK + 1
400   200    CONTINUE
401 C
402          E(J) = G + A(JK) * F
403          JK = JK + 1
404   210 CONTINUE
405 C
406 C        .......... FORM P ..........
407 C
408       F = ZERO
409       DO 245 J = 1, L
410          E(J) = E(J) * HINV
411          F = F + E(J) * D(J)
412   245 CONTINUE
413 C
414 C     .......... FORM Q ..........
415 C
416       HH = F * HALF * HINV
417       DO 250 J = 1, L
418   250 E(J) = E(J) - HH * D(J)
419 C
420       RETURN
421       END
422 C*MODULE EIGEN   *DECK EPSLON
423       DOUBLE PRECISION FUNCTION EPSLON (X)
424 C*
425 C*    AUTHORS -
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
428 C*
429 C*    PURPOSE -
430 C*       ESTIMATE UNIT ROUNDOFF IN QUANTITIES OF SIZE X.
431 C*
432 C*    ON ENTRY -
433 C*       X      - WORKING PRECISION REAL
434 C*                VALUES TO FIND EPSLON FOR
435 C*
436 C*    ON EXIT -
437 C*       EPSLON - WORKING PRECISION REAL
438 C*                SMALLEST POSITIVE VALUE SUCH THAT X+EPSLON .NE. ZERO
439 C*
440 C*    QUALIFICATIONS -
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
450 C*       ASSUMPTION 2.
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.
459 C*
460 C*    DIFFERENCES FROM EISPACK 3 -
461 C*       USE IS MADE OF PARAMETER STATEMENTS AND INTRINSIC FUNCTIONS
462 C*       --NO EXECUTEABLE CODE CHANGES--
463 C*
464 C*    NOTE -
465 C*       QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
466 C*       B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
467 C
468       DOUBLE PRECISION A,B,C,EPS,X
469       DOUBLE PRECISION ZERO, ONE, THREE, FOUR
470 C
471       PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, THREE=3.0D+00, FOUR=4.0D+00)
472 C
473 C-----------------------------------------------------------------------
474 C
475       A = FOUR/THREE
476    10 B = A - ONE
477       C = B + B + B
478       EPS = ABS(C - ONE)
479       IF (EPS .EQ. ZERO) GO TO 10
480       EPSLON = EPS*ABS(X)
481       RETURN
482       END
483 C*MODULE EIGEN   *DECK EQLRAT
484       SUBROUTINE EQLRAT(N,DIAG,E,E2IN,D,IND,IERR,E2)
485 C*
486 C*    AUTHORS -
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)
492 C*
493 C*    PURPOSE -
494 C*       THIS ROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC
495 C*       TRIDIAGONAL MATRIX
496 C*
497 C*    METHOD -
498 C*       RATIONAL QL
499 C*
500 C*    ON ENTRY -
501 C*       N      - INTEGER
502 C*                THE ORDER OF THE MATRIX.
503 C*       D      - W.P. REAL (N)
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.
509 C*
510 C*     ON EXIT -
511 C*       D      - W.P. REAL (N)
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)
517 C*                DESTROYED.
518 C*       IERR   - INTEGER
519 C*                SET TO
520 C*                ZERO       FOR NORMAL RETURN,
521 C*                J          IF THE J-TH EIGENVALUE HAS NOT BEEN
522 C*                           DETERMINED AFTER 30 ITERATIONS.
523 C*
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
529 C*
530 C*    NOTE -
531 C*       QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
532 C*       B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
533 C
534       INTEGER I,J,L,M,N,II,L1,IERR
535       INTEGER IND(N)
536 C
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
540 C
541       PARAMETER (ZERO = 0.0D+00, SCALE= 1.0D+00/64.0D+00, ONE = 1.0D+00)
542 C
543 C-----------------------------------------------------------------------
544       IERR = 0
545       D(1)=DIAG(1)
546       IND(1) = 1
547       K = 0
548       ITAG = 0
549       IF (N .EQ. 1) GO TO 1001
550 C
551       DO 100 I = 2, N
552          D(I)=DIAG(I)
553   100 E2(I-1) = E2IN(I)
554 C
555       F = ZERO
556       T = ZERO
557       B = EPSLON(ONE)
558       C = B *B
559       B = B * SCALE
560       E2(N) = ZERO
561 C
562       DO 290 L = 1, N
563          H = ABS(D(L)) + ABS(E(L))
564          IF (T .GE. H) GO TO 105
565             T = H
566             B = EPSLON(T)
567             C = B * B
568             B = B * SCALE
569   105    CONTINUE
570 C     .......... LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT ..........
571          M = L - 1
572   110    M = M + 1
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 ..........
576 C
577          IF (M .LE. K) GO TO 125
578             IF (M .NE. N) E2IN(M+1) = ZERO
579             K = M
580             ITAG = ITAG + 1
581   125    CONTINUE
582          IF (M .EQ. L) GO TO 210
583 C
584 C           ITERATE
585 C
586          DO 205 J = 1, 30
587 C              .......... FORM SHIFT ..........
588             L1 = L + 1
589             S = SQRT(E2(L))
590             G = D(L)
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))
594             H = G - D(L)
595 C
596             DO 140 I = L1, N
597   140       D(I) = D(I) - H
598 C
599             F = F + H
600 C              .......... RATIONAL QL TRANSFORMATION ..........
601             G = D(M) + B
602             H = G
603             S = ZERO
604             DO 200 I = M-1,L,-1
605                P = G * H
606                R = P + E2(I)
607                E2(I+1) = S * R
608                S = E2(I) / R
609                D(I+1) = H + S * (H + D(I))
610                G = D(I) - E2(I) / G   + B
611                H = G * P / R
612   200       CONTINUE
613 C
614             E2(L) = S * G
615             D(L) = H
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
619             E2(L) = H * E2(L)
620             IF (E2(L) .EQ. ZERO) GO TO 210
621   205    CONTINUE
622 C     .......... SET ERROR -- NO CONVERGENCE TO AN
623 C                EIGENVALUE AFTER 30 ITERATIONS ..........
624       IERR = L
625       GO TO 1001
626 C
627 C           CONVERGED
628 C
629   210    P = D(L) + F
630 C           .......... ORDER EIGENVALUES ..........
631          I = 1
632          IF (L .EQ. 1) GO TO 250
633             IF (P .LT. D(1)) GO TO 230
634                I = L
635 C           .......... LOOP TO FIND ORDERED POSITION
636   220          I = I - 1
637                IF (P .LT. D(I)) GO TO 220
638 C
639                I = I + 1
640                IF (I .EQ. L) GO TO 250
641   230       CONTINUE
642             DO 240 II = L, I+1, -1
643                D(II) = D(II-1)
644                IND(II) = IND(II-1)
645   240       CONTINUE
646 C
647   250    CONTINUE
648          D(I) = P
649          IND(I) = ITAG
650   290 CONTINUE
651 C
652  1001 RETURN
653       END
654 C*MODULE EIGEN   *DECK ESTPI1
655       DOUBLE PRECISION FUNCTION ESTPI1 (N,EVAL,D,E,X,ANORM)
656 C*
657 C*    AUTHOR -
658 C*       STEPHEN T. ELBERT (AMES LABORATORY-USDOE) DATE: 5 DEC 1986
659 C*
660 C*    PURPOSE -
661 C*       EVALUATE SYMMETRIC TRIDIAGONAL MATRIX PERFORMANCE INDEX
662 C*       *        *         *                  *           *
663 C*       FOR 1 EIGENVECTOR
664 C*           *
665 C*
666 C*    METHOD -
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.
673 C*
674 C*    ON ENTRY -
675 C*       N      - INTEGER
676 C*                THE ORDER OF THE MATRIX  A.
677 C*       EVAL   - W.P. REAL
678 C*                THE EIGENVALUE CORRESPONDING TO VECTOR  X.
679 C*       D      - W.P. REAL (N)
680 C*                THE DIAGONAL VECTOR OF  A.
681 C*       E      - W.P. REAL (N)
682 C*                THE SUB-DIAGONAL VECTOR OF  A.
683 C*       X      - W.P. REAL (N)
684 C*                AN EIGENVECTOR OF  A.
685 C*       ANORM  - W.P. REAL
686 C*                THE NORM OF  A  IF IT HAS BEEN PREVIOUSLY COMPUTED.
687 C*
688 C*    ON EXIT -
689 C*       ANORM  - W.P. REAL
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
695 C*
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)
700 C
701       DOUBLE PRECISION ANORM,EVAL,RNORM,SIZE,XNORM
702       DOUBLE PRECISION D(N), E(N), X(N)
703       DOUBLE PRECISION EPSLON, ONE, ZERO
704 C
705       PARAMETER (ZERO = 0.0D+00, ONE = 1.0D+00)
706 C
707 C-----------------------------------------------------------------------
708 C
709       ESTPI1 = ZERO
710       IF( N .LE. 1 ) RETURN
711       SIZE = 10 * N
712       IF (ANORM .EQ. ZERO) THEN
713 C
714 C              COMPUTE NORM OF  A
715 C
716          ANORM = MAX( ABS(D(1)) + ABS(E(2))
717      *               ,ABS(D(N)) + ABS(E(N)))
718          DO 110 I = 2, N-1
719             ANORM = MAX( ANORM, ABS(E(I))+ABS(D(I))+ABS(E(I+1)))
720   110    CONTINUE
721          IF(ANORM .EQ. ZERO) ANORM = ONE
722       END IF
723 C
724 C           COMPUTE NORMS OF RESIDUAL AND EIGENVECTOR
725 C
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))
729       DO 120 I = 2, N-1
730          XNORM = XNORM + ABS(X(I))
731          RNORM = RNORM + ABS(E(I)*X(I-1) + (D(I)-EVAL)*X(I)
732      *                       + E(I+1)*X(I+1))
733   120 CONTINUE
734 C
735       ESTPI1 = RNORM / (EPSLON(SIZE)*ANORM*XNORM)
736       RETURN
737       END
738 C*MODULE EIGEN   *DECK ETRBK3
739       SUBROUTINE ETRBK3(NM,N,NV,A,M,Z)
740 C*
741 C*    AUTHORS-
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)
748 C*
749 C*    PURPOSE -
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.
753 C*
754 C*    METHOD -
755 C*       THE CALCULATION IS CARRIED OUT BY FORMING THE MATRIX PRODUCT
756 C*          Q*Z
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.
765 C*
766 C*
767 C*    COMPLEXITY -
768 C*       M*N**2
769 C*
770 C*    ON ENTRY-
771 C*       NM     - INTEGER
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.
775 C*       N      - INTEGER
776 C*                THE ORDER OF THE MATRIX  A.
777 C*       NV     - INTEGER
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.
784 C*       M      - INTEGER
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.
789 C*
790 C*    ON EXIT-
791 C*       Z      - W.P. REAL (NM,M)
792 C*                CONTAINS THE TRANSFORMED EIGENVECTORS
793 C*                IN ITS FIRST M COLUMNS.
794 C*
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.
800 C*
801 C*    NOTE -
802 C*       QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
803 C*       B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
804 C
805       INTEGER I,II,IM1,IZ,J,M,N,NM,NV
806 C
807       DOUBLE PRECISION A(NV),Z(NM,M)
808       DOUBLE PRECISION H,S,DDOT,ZERO
809 C
810       PARAMETER (ZERO = 0.0D+00)
811 C
812 C-----------------------------------------------------------------------
813 C
814       IF (M .EQ. 0) RETURN
815       IF (N .LE. 2) RETURN
816 C
817       II=3
818       DO 140 I = 3, N
819          IZ=II+1
820          II=II+I
821          H = A(II)
822          IF (H .EQ. ZERO) GO TO 140
823             IM1 = I - 1
824             DO 130 J = 1, M
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)
827   130       CONTINUE
828   140 CONTINUE
829       RETURN
830       END
831 C*MODULE EIGEN   *DECK ETRED3
832       SUBROUTINE ETRED3(N,NV,A,D,E,E2)
833 C*
834 C*    AUTHORS -
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
841 C*
842 C*    PURPOSE -
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.
847 C*
848 C*    METHOD -
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
855 C*       A SCALAR
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.
863 C*
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.
871 C*
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.
875 C*
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.
879 C*
880 C*    COMPLEXITY -
881 C*       2/3 N**3
882 C*
883 C*    ON ENTRY-
884 C*       N      - INTEGER
885 C*                THE ORDER OF THE MATRIX.
886 C*       NV     - INTEGER
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.
893 C*
894 C*    ON EXIT-
895 C*       A      - W.P. REAL (NV)
896 C*                CONTAINS INFORMATION ABOUT THE ORTHOGONAL
897 C*                TRANSFORMATIONS USED IN THE REDUCTION.
898 C*       D      - W.P. REAL (N)
899 C*                CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL
900 C*                MATRIX.
901 C*       E      - W.P. REAL (N)
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.
907 C*
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
913 C*       USE BLAS(1)
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
918 C*
919 C*    NOTE -
920 C*       QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
921 C*       B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
922 C
923       INTEGER I,IIA,IZ0,L,N,NV
924 C
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
929 C
930       PARAMETER (ZERO = 0.0D+00, ONE = 1.0D+00)
931 C
932 C-----------------------------------------------------------------------
933 C
934       IF (N .LE. 2) GO TO 310
935       IZ0 = (N*N+N)/2
936       AIIMAX = ABS(A(IZ0))
937       DO 300 I = N, 3, -1
938          L = I - 1
939          IIA = IZ0
940          IZ0 = IZ0 - I
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
944 C
945 C           THIS ROW IS ALREADY IN TRI-DIAGONAL FORM
946 C
947             D(I) = A(IIA)
948             IF (AIIMAX+D(I) .EQ. AIIMAX) D(I) = ZERO
949             E(I) = A(IIA-1)
950             IF (AIIMAX+E(I) .EQ. AIIMAX) E(I) = ZERO
951             E2(I) = E(I)*E(I)
952             A(IIA) = ZERO
953             GO TO 300
954 C
955          END IF
956 C
957          SCALEI = ONE / SCALE
958          CALL DSCAL(L,SCALEI,A(IZ0+1),1)
959          HROOT = DNRM2(L,A(IZ0+1),1)
960 C
961          F = A(IZ0+L)
962          G = -SIGN(HROOT,F)
963          E(I) = SCALE * G
964          E2(I) = E(I)*E(I)
965          H = HROOT*HROOT - F * G
966          A(IZ0+L) = F - G
967          D(I) = A(IIA)
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)
973 C
974   300 CONTINUE
975   310 CONTINUE
976       E(1) = ZERO
977       E2(1)= ZERO
978       D(1) = A(1)
979       IF(N.EQ.1) RETURN
980 C
981       E(2) = A(2)
982       E2(2)= A(2)*A(2)
983       D(2) = A(3)
984       RETURN
985       END
986 C*MODULE EIGEN   *DECK EVVRSP
987       SUBROUTINE EVVRSP(MSGFL,N,NVECT,LENA,NV,A,B,IND,ROOT,
988      *                  VECT,IORDER,IERR)
989 C*
990 C*    AUTHOR:  S. T. ELBERT, AMES LABORATORY-USDOE, JUNE 1985
991 C*
992 C*    PURPOSE -
993 C*       FINDS   (ALL) EIGENVALUES    AND    (SOME OR ALL) EIGENVECTORS
994 C*                     *    *                                   *
995 C*       OF A REAL SYMMETRIC PACKED MATRIX.
996 C*            *    *         *
997 C*
998 C*    METHOD -
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
1007 C*       ORIGINAL ARRAY.
1008 C*
1009 C*       THESE ROUTINES ARE MODIFICATIONS OF THE EISPACK 3
1010 C*       ROUTINES TRED3, TQLRAT, TINVIT AND TRBAK3
1011 C*
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)
1017 C*
1018 C*    ON ENTRY -
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.
1023 C*       N      - INTEGER
1024 C*                ORDER OF MATRIX A.
1025 C*       NVECT  - INTEGER
1026 C*                NUMBER OF VECTORS DESIRED.  0 .LE. NVECT .LE. N.
1027 C*       LENA   - INTEGER
1028 C*                DIMENSION OF  A  IN CALLING ROUTINE.  MUST NOT BE LESS
1029 C*                THAN (N*N+N)/2.
1030 C*       NV     - INTEGER
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.
1040 C*       IORDER - INTEGER
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.
1044 C*
1045 C*    ON EXIT -
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).
1053 C*       IERR   - INTEGER
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.)
1058 C*
1059 C
1060       LOGICAL GOPARR,DSKWRK,MASWRK
1061 C
1062       DOUBLE PRECISION A(LENA)
1063       DOUBLE PRECISION B(N,8)
1064       DOUBLE PRECISION ROOT(N)
1065       DOUBLE PRECISION T
1066       DOUBLE PRECISION VECT(NV,*)
1067 C
1068       INTEGER IND(N)
1069 C
1070       COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
1071 C
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')
1086 C
1087 C-----------------------------------------------------------------------
1088 C
1089       LMSGFL=MSGFL
1090       IF (MSGFL .EQ. 0) LMSGFL=6
1091       IERR = N - 1
1092       IF (N .LE. 0) GO TO 800
1093       IERR = N + 1
1094       IF ( (N*N+N)/2 .GT. LENA) GO TO 810
1095 C
1096 C        REDUCE REAL SYMMETRIC MATRIX A TO TRIDIAGONAL FORM
1097 C
1098       CALL ETRED3(N,LENA,A,B(1,1),B(1,2),B(1,3))
1099 C
1100 C        FIND ALL EIGENVALUES OF TRIDIAGONAL MATRIX
1101 C
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
1104 C
1105 C         CHECK THE DESIRED ORDER OF THE EIGENVALUES
1106 C
1107       B(1,3) = IORDER
1108       IF (IORDER .EQ. 0) GO TO 300
1109          IF (IORDER .NE. 2) GO TO 850
1110 C
1111 C         ORDER ROOTS IN DESCENDING ORDER (LARGEST FIRST)...
1112 C        TURN ROOT AND IND ARRAYS END FOR END
1113 C
1114          DO 210 I = 1, N/2
1115             J = N+1-I
1116             T = ROOT(I)
1117             ROOT(I) = ROOT(J)
1118             ROOT(J) = T
1119             L = IND(I)
1120             IND(I) = IND(J)
1121             IND(J) = L
1122   210    CONTINUE
1123 C
1124 C           FIND I AND J MARKING THE START AND END OF A SEQUENCE
1125 C           OF DEGENERATE ROOTS
1126 C
1127          I=0
1128   220    CONTINUE
1129             I = I+1
1130             IF (I .GT. N) GO TO 300
1131             DO 230 J=I,N
1132                IF (ROOT(J) .NE. ROOT(I)) GO TO 240
1133   230       CONTINUE
1134             J = N+1
1135   240       CONTINUE
1136             J = J-1
1137             IF (J .EQ. I) GO TO 220
1138 C
1139 C                    TURN AROUND IND BETWEEN I AND J
1140 C
1141             JSV = J
1142             KLIM = (J-I+1)/2
1143             DO 250 K=1,KLIM
1144                L = IND(J)
1145                IND(J) = IND(I)
1146                IND(I) = L
1147                I = I+1
1148                J = J-1
1149   250       CONTINUE
1150             I = JSV
1151          GO TO 220
1152 C
1153   300 CONTINUE
1154 C
1155       IF (NVECT .LE. 0) RETURN
1156       IF (NV .LT. N) GO TO 830
1157 C
1158 C        FIND EIGENVECTORS OF TRI-DIAGONAL MATRIX VIA INVERSE ITERATION
1159 C
1160       IERR = LMSGFL
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
1164 C
1165 C        FIND EIGENVECTORS OF SYMMETRIC MATRIX VIA BACK TRANSFORMATION
1166 C
1167   400 CONTINUE
1168       CALL ETRBK3(NV,N,LENA,A,NVECT,VECT)
1169       RETURN
1170 C
1171 C        ERROR MESSAGE SECTION
1172 C
1173   800 IF (LMSGFL .LT. 0) RETURN
1174       IF (MASWRK) WRITE(LMSGFL,906)
1175       GO TO 890
1176 C
1177   810 IF (LMSGFL .LT. 0) RETURN
1178       IF (MASWRK) WRITE(LMSGFL,901)
1179       GO TO 890
1180 C
1181   820 IF (LMSGFL .LT. 0) RETURN
1182       IF (MASWRK) WRITE(LMSGFL,902) IERR
1183       GO TO 890
1184 C
1185   830 IF (LMSGFL .LT. 0) RETURN
1186       IF (MASWRK) WRITE(LMSGFL,903)
1187       GO TO 890
1188 C
1189   840 CONTINUE
1190       IF ((LMSGFL .GT. 0).AND.MASWRK) WRITE(LMSGFL,904) -IERR
1191       GO TO 400
1192 C
1193   850 IERR=-1
1194       IF (LMSGFL .LT. 0) RETURN
1195       IF (MASWRK) WRITE(LMSGFL,905)
1196       GO TO 890
1197 C
1198   890 CONTINUE
1199       IF (MASWRK) WRITE(LMSGFL,900) N,NVECT,LENA,NV,IORDER,IERR
1200       RETURN
1201       END
1202 C*MODULE EIGEN   *DECK FREDA
1203       SUBROUTINE FREDA(L,D,A,E)
1204 C
1205       DOUBLE PRECISION A(*)
1206       DOUBLE PRECISION D(L)
1207       DOUBLE PRECISION E(L)
1208       DOUBLE PRECISION F
1209       DOUBLE PRECISION G
1210 C
1211       JK = 1
1212 C
1213 C     .......... FORM REDUCED A ..........
1214 C
1215       DO 280 J = 1, L
1216          F = D(J)
1217          G = E(J)
1218 C
1219          DO 260 K = 1, J
1220             A(JK) = A(JK) - F * E(K) - G * D(K)
1221             JK = JK + 1
1222   260    CONTINUE
1223 C
1224   280 CONTINUE
1225       RETURN
1226       END
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)
1231 C
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.
1235 C
1236 C     INPUT..
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
1245 C
1246 C     OUTPUT..
1247 C     A       DESTROYED .
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.)
1255 C
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.
1261 C
1262 C         IF TINVIT FAILS TO CONVERGE, TQL2 IS CALLED
1263 C
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.
1268 C
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
1273 C
1274 C     TO REORDER ROOTS...
1275 C     K = N/2
1276 C     B(1,3) = 2.0D+00
1277 C     DO 50 I = 1, K
1278 C        J = N+1-I
1279 C        T = ROOT(I)
1280 C        ROOT(I) = ROOT(J)
1281 C        ROOT(J) = T
1282 C 50  CONTINUE
1283 C
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
1288 C
1289 C      IF INVERSE ITERATION GIVES AN ERROR IN DETERMINING THE
1290 C      EIGENVECTORS, TRY THE QL ALGORITHM IF ALL THE EIGENVECTORS
1291 C      ARE DESIRED.
1292 C
1293       IF (NVECT .NE. N) RETURN
1294       DO 120 I = 1, NVECT
1295       DO 100 J = 1, N
1296       VECT(I,J) = ZERO
1297   100 CONTINUE
1298       VECT(I,I) = ONE
1299   120 CONTINUE
1300       CALL TQL2 (NV,N,B(1,1),B(1,2),VECT,IERR)
1301       DO 140 I = 1, NVECT
1302       ROOT(I) = B(I,1)
1303   140 CONTINUE
1304       IF (IERR .NE. 0) RETURN
1305   160 CALL TRBK3B(NV,N,(N*N+N)/2,A,NVECT,VECT)
1306       RETURN
1307       END
1308 C*MODULE EIGEN   *DECK GLDIAG
1309       SUBROUTINE GLDIAG(LDVECT,NVECT,N,H,WRK,EIG,VECTOR,IERR,IWRK)
1310 C
1311       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1312 C
1313       LOGICAL GOPARR,DSKWRK,MASWRK
1314 C
1315       DIMENSION H(*),WRK(N,8),EIG(N),VECTOR(LDVECT,NVECT),IWRK(N)
1316 C
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
1320 C
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
1325 C              = 1, USE EVVRSP
1326 C              = 2, USE GIVEIS
1327 C              = 3, USE JACOBI
1328 C
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
1338 C
1339       IERR = 0
1340 C
1341 C         ----- USE STEVE ELBERT'S ROUTINE -----
1342 C
1343       IF(KDIAG.LE.1  .OR.  KDIAG.GT.3) THEN
1344          LENH = (N*N+N)/2
1345          KORDER =0
1346          CALL EVVRSP(IW,N,NVECT,LENH,LDVECT,H,WRK,IWRK,EIG,VECTOR
1347      *              ,KORDER,IERR)
1348       END IF
1349 C
1350 C         ----- USE MODIFIED EISPAK ROUTINE -----
1351 C
1352       IF(KDIAG.EQ.2)
1353      *   CALL GIVEIS(N,NVECT,LDVECT,H,WRK,IWRK,EIG,VECTOR,IERR)
1354 C
1355 C         ----- USE JACOBI ROTATION ROUTINE -----
1356 C
1357       IF(KDIAG.EQ.3) THEN
1358          IF(NVECT.EQ.N) THEN
1359             CALL JACDG(H,VECTOR,EIG,IWRK,WRK,LDVECT,N)
1360          ELSE
1361             IF (MASWRK) WRITE(IW,9000) N,NVECT,LDVECT
1362             CALL ABRT
1363          END IF
1364       END IF
1365       RETURN
1366 C
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.')
1370       END
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)
1374       INTEGER TAG
1375       DOUBLE PRECISION MACHEP
1376       DIMENSION D(N),E(N),E2(N),W(N),RV1(N),IND(N)
1377 C
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).
1382 C
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.
1386 C
1387 C     ON INPUT-
1388 C
1389 C        N IS THE ORDER OF THE MATRIX,
1390 C
1391 C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
1392 C
1393 C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
1394 C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY,
1395 C
1396 C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
1397 C          E2(1) IS ARBITRARY.
1398 C
1399 C     ON OUTPUT-
1400 C
1401 C        D AND E ARE UNALTERED,
1402 C
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,
1407 C
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,
1412 C
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.,
1417 C
1418 C        IERR IS SET TO
1419 C          ZERO       FOR NORMAL RETURN,
1420 C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
1421 C                     DETERMINED AFTER 30 ITERATIONS,
1422 C
1423 C        RV1 IS A TEMPORARY STORAGE ARRAY.
1424 C
1425 C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
1426 C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
1427 C
1428 C     ------------------------------------------------------------------
1429 C
1430 C     ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
1431 C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
1432 C
1433 C                **********
1434       MACHEP = 2.0D+00**(-50)
1435 C
1436       IERR = 0
1437       K = 0
1438       TAG = 0
1439 C
1440       DO 100 I = 1, N
1441       W(I) = D(I)
1442       IF (I .NE. 1) RV1(I-1) = E(I)
1443   100 CONTINUE
1444 C
1445       E2(1) = 0.0D+00
1446       RV1(N) = 0.0D+00
1447 C
1448       DO 360 L = 1, N
1449       J = 0
1450 C     ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT **********
1451   120 DO 140 M = L, N
1452       IF (M .EQ. N) GO TO 160
1453       IF (ABS(RV1(M)) .LE. MACHEP * (ABS(W(M)) + ABS(W(M+1)))) GO TO
1454      +     160
1455 C     ********** GUARD AGAINST UNDERFLOWED ELEMENT OF E2 **********
1456       IF (E2(M+1) .EQ. 0.0D+00) GO TO 180
1457   140 CONTINUE
1458 C
1459   160 IF (M .LE. K) GO TO 200
1460       IF (M .NE. N) E2(M+1) = 0.0D+00
1461   180 K = M
1462       TAG = TAG + 1
1463   200 P = W(L)
1464       IF (M .EQ. L) GO TO 280
1465       IF (J .EQ. 30) GO TO 380
1466       J = J + 1
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))
1471       S = 1.0D+00
1472       C = 1.0D+00
1473       P = 0.0D+00
1474       MML = M - L
1475 C     ********** FOR I=M-1 STEP -1 UNTIL L DO -- **********
1476       DO 260 II = 1, MML
1477       I = M - II
1478       F = S * RV1(I)
1479       B = C * RV1(I)
1480       IF (ABS(F) .LT. ABS(G)) GO TO 220
1481       C = G / F
1482       R = SQRT(C*C+1.0D+00)
1483       RV1(I+1) = F * R
1484       S = 1.0D+00 / R
1485       C = C * S
1486       GO TO 240
1487   220 S = F / G
1488       R = SQRT(S*S+1.0D+00)
1489       RV1(I+1) = G * R
1490       C = 1.0D+00 / R
1491       S = S * C
1492   240 G = W(I+1) - P
1493       R = (W(I) - G) * S + 2.0D+00 * C * B
1494       P = S * R
1495       W(I+1) = G + P
1496       G = C * R - B
1497   260 CONTINUE
1498 C
1499       W(L) = W(L) - P
1500       RV1(L) = G
1501       RV1(M) = 0.0D+00
1502       GO TO 120
1503 C     ********** ORDER EIGENVALUES **********
1504   280 IF (L .EQ. 1) GO TO 320
1505 C     ********** FOR I=L STEP -1 UNTIL 2 DO -- **********
1506       DO 300 II = 2, L
1507       I = L + 2 - II
1508       IF (P .GE. W(I-1)) GO TO 340
1509       W(I) = W(I-1)
1510       IND(I) = IND(I-1)
1511   300 CONTINUE
1512 C
1513   320 I = 1
1514   340 W(I) = P
1515       IND(I) = TAG
1516   360 CONTINUE
1517 C
1518       GO TO 400
1519 C     ********** SET ERROR -- NO CONVERGENCE TO AN
1520 C                EIGENVALUE AFTER 30 ITERATIONS **********
1521   380 IERR = L
1522   400 RETURN
1523 C     ********** LAST CARD OF IMTQLV **********
1524       END
1525 C*MODULE EIGEN   *DECK JACDG
1526       SUBROUTINE JACDG(A,VEC,EIG,JBIG,BIG,LDVEC,N)
1527 C
1528       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1529 C
1530       DIMENSION A(*),VEC(LDVEC,N),EIG(N),JBIG(N),BIG(N)
1531 C
1532       PARAMETER (ONE=1.0D+00)
1533 C
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.
1539 C
1540       CALL VCLR(VEC,1,LDVEC*N)
1541       DO 20 I = 1,N
1542         VEC(I,I) = ONE
1543    20 CONTINUE
1544 C
1545       NB1 = N
1546       NB2 = (NB1*NB1+NB1)/2
1547       NMIN = 1
1548       NMAX = NB1
1549 C
1550       CALL JACDIA(A,VEC,NB1,NB2,LDVEC,NMIN,NMAX,BIG,JBIG)
1551 C
1552       DO 30 I=1,N
1553         EIG(I) = A((I*I+I)/2)
1554    30 CONTINUE
1555 C
1556       CALL JACORD(VEC,EIG,NB1,LDVEC)
1557       RETURN
1558       END
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)
1564 C
1565       COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
1566 C
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 )
1573 C
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
1578 C      ACCOUNTED FOR.
1579 C      THE ROTATIONS AMONG THE LAST NB1-NMAX BASIS FUNCTIONS ARE NOT
1580 C      ACCOUNTED FOR.
1581 C
1582       IEAA=0
1583       IEAB=0
1584       TT=ZERO
1585       EPS = 64.0D+00*EPSLON(ONE)
1586 C
1587 C      LOOP OVER COLUMNS (K) OF TRIANGULAR MATRIX TO DETERMINE
1588 C      LARGEST OFF-DIAGONAL ELEMENTS IN ROW(I).
1589 C
1590       DO 20 I=1,NB1
1591          BIG(I)=ZERO
1592          JBIG(I)=0
1593          IF(I.LT.NMIN  .OR.  I.EQ.1) GO TO 20
1594          II = (I*I-I)/2
1595          J=MIN(I-1,NMAX)
1596          DO 10 K=1,J
1597             IF(ABS(BIG(I)).GE.ABS(F(II+K))) GO TO 10
1598             BIG(I)=F(II+K)
1599             JBIG(I)=K
1600    10    CONTINUE
1601    20 CONTINUE
1602 C
1603 C     ----- 2X2 JACOBI ITERATIONS BEGIN HERE -----
1604 C
1605       MAXIT=MAX(NB2*20,500)
1606       ITER=0
1607    30 CONTINUE
1608       ITER=ITER+1
1609 C
1610 C      FIND SMALLEST DIAGONAL ELEMENT
1611 C
1612       SD=D1050
1613       JJ=0
1614       DO 40 J=1,NB1
1615          JJ=JJ+J
1616          SD= MIN(SD,ABS(F(JJ)))
1617    40 CONTINUE
1618       TEST = MAX(EPS, C2*MAX(SD,C6))
1619 C
1620 C      FIND LARGEST OFF-DIAGONAL ELEMENT
1621 C
1622       T=ZERO
1623       I1=MAX(2,NMIN)
1624       IB = I1
1625       DO 50 I=I1,NB1
1626          IF(T.GE.ABS(BIG(I))) GO TO 50
1627          T= ABS(BIG(I))
1628          IB=I
1629    50 CONTINUE
1630 C
1631 C      TEST FOR CONVERGENCE, THEN DETERMINE ROTATION.
1632 C
1633       IF(T.LT.TEST) RETURN
1634 C                   ******
1635 C
1636       IF(ITER.GT.MAXIT) THEN
1637          IF (MASWRK) THEN
1638             WRITE(6,*) 'JACOBI DIAGONALIZATION FAILS, DIMENSION=',NB1
1639             WRITE(6,9020) ITER,T,TEST,SD
1640          ENDIF
1641          CALL ABRT
1642          STOP
1643       END IF
1644 C
1645       IA=JBIG(IB)
1646       IAA=IA*(IA-1)/2
1647       IBB=IB*(IB-1)/2
1648       DIF=F(IAA+IA)-F(IBB+IB)
1649       IF(ABS(DIF).GT.C3*T) GO TO 70
1650       SX=ROOT2
1651       CX=ROOT2
1652       GO TO 110
1653    70 T2X2=BIG(IB)/DIF
1654       T2X25=T2X2*T2X2
1655       IF(T2X25 . GT . C4) GO TO 80
1656       CX=ONE
1657       SX=T2X2
1658       GO TO 110
1659    80 IF(T2X25 . GT . C5) GO TO 90
1660       SX=T2X2*(ONE-D1500*T2X25)
1661       CX=ONE-D0500*T2X25
1662       GO TO 110
1663    90 IF(T2X25 . GT . C6) GO TO 100
1664       CX=ONE+T2X25*(T2X25*D1375 - D0500)
1665       SX= T2X2*(ONE + T2X25*(T2X25*D3875 - D1500))
1666       GO TO 110
1667   100 T=D0250  / SQRT(D0250   + T2X25)
1668       CX= SQRT(D0500   + T)
1669       SX= SIGN( SQRT(D0500   - T),T2X2)
1670   110 IEAR=IAA+1
1671       IEBR=IBB+1
1672 C
1673       DO 230 IR=1,NB1
1674          T=F(IEAR)*SX
1675          F(IEAR)=F(IEAR)*CX+F(IEBR)*SX
1676          F(IEBR)=T-F(IEBR)*CX
1677          IF(IR-IA) 220,120,130
1678   120    TT=F(IEBR)
1679          IEAA=IEAR
1680          IEAB=IEBR
1681          F(IEBR)=BIG(IB)
1682          IEAR=IEAR+IR-1
1683          IF(JBIG(IR)) 200,220,200
1684   130    T=F(IEAR)
1685          IT=IA
1686          IEAR=IEAR+IR-1
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
1691          IEBR=IEBR+IR-1
1692          GO TO 200
1693   160    IF(  ABS(T) . GE .  ABS(F(IEBR))) GO TO 170
1694          IF(IB.GT.NMAX) GO TO 170
1695          T=F(IEBR)
1696          IT=IB
1697   170    IEBR=IEBR+IR-1
1698   180    IF(  ABS(T) . LT .  ABS(BIG(IR))) GO TO 190
1699          BIG(IR) = T
1700          JBIG(IR) = IT
1701          GO TO 220
1702   190    IF(IA . NE . JBIG(IR) . AND . IB . NE . JBIG(IR))  GO TO 220
1703   200    KQ=IEAR-IR-IA+1
1704          BIG(IR)=ZERO
1705          IR1=MIN(IR-1,NMAX)
1706          DO 210 I=1,IR1
1707             K=KQ+I
1708             IF(ABS(BIG(IR)) . GE . ABS(F(K)))  GO TO 210
1709             BIG(IR) = F(K)
1710             JBIG(IR)=I
1711   210    CONTINUE
1712   220    IEAR=IEAR+1
1713   230    IEBR=IEBR+1
1714 C
1715       DO 240 I=1,NB1
1716          T1=VEC(I,IA)*CX + VEC(I,IB)*SX
1717          T2=VEC(I,IA)*SX - VEC(I,IB)*CX
1718          VEC(I,IA)=T1
1719          VEC(I,IB)=T2
1720   240 CONTINUE
1721       GO TO 30
1722 C
1723  9020 FORMAT(1X,'ITER=',I6,' T,TEST,SD=',1P,3E20.10)
1724       END
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)
1729 C
1730 C     ---- SORT EIGENDATA INTO ASCENDING ORDER -----
1731 C
1732       DO 290 I = 1, N
1733          JJ = I
1734          DO 270 J = I, N
1735             IF (EIG(J) .LT. EIG(JJ)) JJ = J
1736   270    CONTINUE
1737          IF (JJ .EQ. I) GO TO 290
1738          T = EIG(JJ)
1739          EIG(JJ) = EIG(I)
1740          EIG(I) = T
1741          DO 280 J = 1, N
1742             T = VEC(J,JJ)
1743             VEC(J,JJ) = VEC(J,I)
1744             VEC(J,I) = T
1745   280    CONTINUE
1746   290 CONTINUE
1747       RETURN
1748       END
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     ------------------------------------------------------------------
1758 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).
1762 C
1763 C     THIS ROUTINE FINDS THOSE EIGENVECTORS OF A TRIDIAGONAL
1764 C     SYMMETRIC MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES,
1765 C     USING INVERSE ITERATION.
1766 C
1767 C     ON INPUT-
1768 C
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,
1772 C
1773 C        N IS THE ORDER OF THE MATRIX,
1774 C
1775 C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
1776 C
1777 C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
1778 C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY,
1779 C
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,
1789 C
1790 C        M IS THE NUMBER OF SPECIFIED EIGENVALUES,
1791 C
1792 C        W CONTAINS THE M EIGENVALUES IN ASCENDING OR DESCENDING ORDER,
1793 C
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.
1798 C
1799 C     ON OUTPUT-
1800 C
1801 C        ALL INPUT ARRAYS ARE UNALTERED,
1802 C
1803 C        Z CONTAINS THE ASSOCIATED SET OF ORTHONORMAL EIGENVECTORS.
1804 C          ANY VECTOR WHICH FAILS TO CONVERGE IS SET TO ZERO,
1805 C
1806 C        IERR IS SET TO
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,
1810 C
1811 C        RV1, RV2, RV3, RV4, AND RV6 ARE TEMPORARY STORAGE ARRAYS.
1812 C
1813 C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
1814 C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
1815 C
1816 C     ------------------------------------------------------------------
1817 C
1818 C     ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
1819 C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
1820 C
1821 C                **********
1822       MACHEP = 2.0D+00**(-50)
1823 C
1824       IERR = 0
1825       IF (M .EQ. 0) GO TO 680
1826       TAG = 0
1827       ORDER = 1.0D+00 - E2(1)
1828       XU = 0.0D+00
1829       UK = 0.0D+00
1830       X0 = 0.0D+00
1831       U  = 0.0D+00
1832       EPS2 = 0.0D+00
1833       EPS3 = 0.0D+00
1834       EPS4 = 0.0D+00
1835       GROUP = 0
1836       Q = 0
1837 C     ********** ESTABLISH AND PROCESS NEXT SUBMATRIX **********
1838   100 P = Q + 1
1839       IP = P + 1
1840 C
1841       DO 120 Q = P, N
1842       IF (Q .EQ. N) GO TO 140
1843       IF (E2(Q+1) .EQ. 0.0D+00) GO TO 140
1844   120 CONTINUE
1845 C     ********** FIND VECTORS BY INVERSE ITERATION **********
1846   140 TAG = TAG + 1
1847       IQMP = Q - P + 1
1848       S = 0
1849 C
1850       DO 660 R = 1, M
1851       IF (IND(R) .NE. TAG) GO TO 660
1852       ITS = 1
1853       X1 = W(R)
1854       IF (S .NE. 0) GO TO 220
1855 C     ********** CHECK FOR ISOLATED ROOT **********
1856       XU = 1.0D+00
1857       IF (P .NE. Q) GO TO 160
1858       RV6(P) = 1.0D+00
1859       GO TO 600
1860   160 NORM = ABS(D(P))
1861 C
1862       DO 180 I = IP, Q
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
1870       UK = IQMP
1871       EPS4 = UK * EPS3
1872       UK = EPS4 / SQRT(UK)
1873       S = P
1874   200 GROUP = 0
1875       GO TO 240
1876 C     ********** LOOK FOR CLOSE OR COINCIDENT ROOTS **********
1877   220 IF (ABS(X1-X0) .GE. EPS2) GO TO 200
1878       GROUP = GROUP + 1
1879       IF (ORDER * (X1 - X0) .LE. 0.0D+00) X1 = X0 + ORDER * EPS3
1880 C     ********** ELIMINATION WITH INTERCHANGES AND
1881 C                INITIALIZATION OF VECTOR **********
1882   240 V = 0.0D+00
1883 C
1884       DO 300 I = P, Q
1885       RV6(I) = UK
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 **********
1890       XU = U / E(I)
1891       RV4(I) = XU
1892       RV1(I-1) = E(I)
1893       RV2(I-1) = D(I) - X1
1894       RV3(I-1) = 0.0D+00
1895       IF (I .NE. Q) RV3(I-1) = E(I+1)
1896       U = V - XU * RV2(I-1)
1897       V = -XU * RV3(I-1)
1898       GO TO 300
1899   260 XU = E(I) / U
1900       RV4(I) = XU
1901       RV1(I-1) = U
1902       RV2(I-1) = V
1903       RV3(I-1) = 0.0D+00
1904   280 U = D(I) - X1 - XU * V
1905       IF (I .NE. Q) V = E(I+1)
1906   300 CONTINUE
1907 C
1908       IF (U .EQ. 0.0D+00) U = EPS3
1909       RV1(Q) = U
1910       RV2(Q) = 0.0D+00
1911       RV3(Q) = 0.0D+00
1912 C     ********** BACK SUBSTITUTION
1913 C                FOR I=Q STEP -1 UNTIL P DO -- **********
1914   320 DO 340 II = P, Q
1915       I = P + Q - II
1916       RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I)
1917       V = U
1918       U = RV6(I)
1919   340 CONTINUE
1920 C     ********** ORTHOGONALIZE WITH RESPECT TO PREVIOUS
1921 C                MEMBERS OF GROUP **********
1922       IF (GROUP .EQ. 0) GO TO 400
1923       J = R
1924 C
1925       DO 380 JJ = 1, GROUP
1926   360 J = J - 1
1927       IF (IND(J) .NE. TAG) GO TO 360
1928       XU = DDOT(IQMP,RV6(P),1,Z(P,J),1)
1929 C
1930       CALL DAXPY(IQMP,-XU,Z(P,J),1,RV6(P),1)
1931 C
1932   380 CONTINUE
1933 C
1934   400 NORM = 0.0D+00
1935 C
1936       DO 420 I = P, Q
1937   420 NORM = NORM + ABS(RV6(I))
1938 C
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
1943       RV6(S) = EPS4
1944       S = S + 1
1945       IF (S .GT. Q) S = P
1946       GO TO 480
1947   440 XU = EPS4 / NORM
1948 C
1949       DO 460 I = P, Q
1950   460 RV6(I) = RV6(I) * XU
1951 C     ********** ELIMINATION OPERATIONS ON NEXT VECTOR
1952 C                ITERATE **********
1953   480 DO 520 I = IP, Q
1954       U = RV6(I)
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
1959       U = RV6(I-1)
1960       RV6(I-1) = RV6(I)
1961   500 RV6(I) = U - RV4(I) * RV6(I-1)
1962   520 CONTINUE
1963 C
1964       ITS = ITS + 1
1965       GO TO 320
1966 C     ********** SET ERROR -- NON-CONVERGED EIGENVECTOR **********
1967   540 IERR = -R
1968       XU = 0.0D+00
1969       GO TO 600
1970 C     ********** NORMALIZE SO THAT SUM OF SQUARES IS
1971 C                1 AND EXPAND TO FULL ORDER **********
1972   560 U = 0.0D+00
1973 C
1974       DO 580 I = P, Q
1975       RV6(I) = RV6(I) / NORM
1976   580 U = U + RV6(I)**2
1977 C
1978       XU = 1.0D+00 / SQRT(U)
1979 C
1980   600 DO 620 I = 1, N
1981   620 Z(I,R) = 0.0D+00
1982 C
1983       DO 640 I = P, Q
1984   640 Z(I,R) = RV6(I) * XU
1985 C
1986       X0 = X1
1987   660 CONTINUE
1988 C
1989       IF (Q .LT. N) GO TO 100
1990   680 RETURN
1991 C     ********** LAST CARD OF TINVIT **********
1992       END
1993 C*MODULE EIGEN   *DECK TQL2
1994 C
1995 C     ------------------------------------------------------------------
1996 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)
2001 C
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
2004 C     WILKINSON.
2005 C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971).
2006 C
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.
2012 C
2013 C     ON INPUT-
2014 C
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,
2018 C
2019 C        N IS THE ORDER OF THE MATRIX,
2020 C
2021 C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
2022 C
2023 C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
2024 C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY,
2025 C
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.
2030 C
2031 C      ON OUTPUT-
2032 C
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,
2036 C
2037 C        E HAS BEEN DESTROYED,
2038 C
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
2042 C          EIGENVALUES,
2043 C
2044 C        IERR IS SET TO
2045 C          ZERO       FOR NORMAL RETURN,
2046 C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
2047 C                     DETERMINED AFTER 30 ITERATIONS.
2048 C
2049 C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
2050 C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
2051 C
2052 C     ------------------------------------------------------------------
2053 C
2054 C     ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
2055 C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
2056 C
2057 C                **********
2058       MACHEP = 2.0D+00**(-50)
2059 C
2060       IERR = 0
2061       IF (N .EQ. 1) GO TO 400
2062 C
2063       DO 100 I = 2, N
2064   100 E(I-1) = E(I)
2065 C
2066       F = 0.0D+00
2067       B = 0.0D+00
2068       E(N) = 0.0D+00
2069 C
2070       DO 300 L = 1, N
2071       J = 0
2072       H = MACHEP * (ABS(D(L)) + ABS(E(L)))
2073       IF (B .LT. H) B = H
2074 C     ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT **********
2075       DO 120 M = L, N
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 **********
2079   120 CONTINUE
2080 C
2081   140 IF (M .EQ. L) GO TO 280
2082   160 IF (J .EQ. 30) GO TO 380
2083       J = J + 1
2084 C     ********** FORM SHIFT **********
2085       L1 = L + 1
2086       G = D(L)
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))
2090       H = G - D(L)
2091 C
2092       DO 180 I = L1, N
2093   180 D(I) = D(I) - H
2094 C
2095       F = F + H
2096 C     ********** QL TRANSFORMATION **********
2097       P = D(M)
2098       C = 1.0D+00
2099       S = 0.0D+00
2100       MML = M - L
2101 C     ********** FOR I=M-1 STEP -1 UNTIL L DO -- **********
2102       DO 260 II = 1, MML
2103       I = M - II
2104       G = C * E(I)
2105       H = C * P
2106       IF (ABS(P) .LT. ABS(E(I))) GO TO 200
2107       C = E(I) / P
2108       R = SQRT(C*C+1.0D+00)
2109       E(I+1) = S * P * R
2110       S = C / R
2111       C = 1.0D+00 / R
2112       GO TO 220
2113   200 C = P / E(I)
2114       R = SQRT(C*C+1.0D+00)
2115       E(I+1) = S * E(I) * R
2116       S = 1.0D+00 / R
2117       C = C * S
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)
2122 C
2123   260 CONTINUE
2124 C
2125       E(L) = S * P
2126       D(L) = C * P
2127       IF (ABS(E(L)) .GT. B) GO TO 160
2128   280 D(L) = D(L) + F
2129   300 CONTINUE
2130 C     ********** ORDER EIGENVALUES AND EIGENVECTORS **********
2131       DO 360 II = 2, N
2132       I = II - 1
2133       K = I
2134       P = D(I)
2135 C
2136       DO 320 J = II, N
2137       IF (D(J) .GE. P) GO TO 320
2138       K = J
2139       P = D(J)
2140   320 CONTINUE
2141 C
2142       IF (K .EQ. I) GO TO 360
2143       D(K) = D(I)
2144       D(I) = P
2145 C
2146       CALL DSWAP(N,Z(1,I),1,Z(1,K),1)
2147 C
2148   360 CONTINUE
2149 C
2150       GO TO 400
2151 C     ********** SET ERROR -- NO CONVERGENCE TO AN
2152 C                EIGENVALUE AFTER 30 ITERATIONS **********
2153   380 IERR = L
2154   400 RETURN
2155 C     ********** LAST CARD OF TQL2 **********
2156       END
2157 C*MODULE EIGEN   *DECK TRBK3B
2158 C
2159 C     ------------------------------------------------------------------
2160 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)
2164 C
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).
2168 C
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.
2172 C
2173 C     ON INPUT-
2174 C
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,
2178 C
2179 C        N IS THE ORDER OF THE MATRIX,
2180 C
2181 C        NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A
2182 C          AS DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT,
2183 C
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,
2187 C
2188 C        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED,
2189 C
2190 C        Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
2191 C          IN ITS FIRST M COLUMNS.
2192 C
2193 C     ON OUTPUT-
2194 C
2195 C        Z CONTAINS THE TRANSFORMED EIGENVECTORS
2196 C          IN ITS FIRST M COLUMNS.
2197 C
2198 C     NOTE THAT TRBAK3 PRESERVES VECTOR EUCLIDEAN NORMS.
2199 C
2200 C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
2201 C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
2202 C
2203 C     ------------------------------------------------------------------
2204 C
2205       IF (M .EQ. 0) GO TO 140
2206       IF (N .EQ. 1) GO TO 140
2207 C
2208       DO 120 I = 2, N
2209       L = I - 1
2210       IZ = (I * L) / 2
2211       IK = IZ + I
2212       H = A(IK)
2213       IF (H .EQ. 0.0D+00) GO TO 120
2214 C
2215       DO 100 J = 1, M
2216       S = -DDOT(L,A(IZ+1),1,Z(1,J),1)
2217 C
2218 C     ********** DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW **********
2219       S = (S / H) / H
2220 C
2221       CALL DAXPY(L,S,A(IZ+1),1,Z(1,J),1)
2222 C
2223   100 CONTINUE
2224 C
2225   120 CONTINUE
2226 C
2227   140 RETURN
2228 C     ********** LAST CARD OF TRBAK3 **********
2229       END
2230 C*MODULE EIGEN   *DECK TRED3B
2231 C
2232 C     ------------------------------------------------------------------
2233 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)
2237 C
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).
2241 C
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.
2245 C
2246 C     ON INPUT-
2247 C
2248 C        N IS THE ORDER OF THE MATRIX,
2249 C
2250 C        NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A
2251 C          AS DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT,
2252 C
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.
2256 C
2257 C     ON OUTPUT-
2258 C
2259 C        A CONTAINS INFORMATION ABOUT THE ORTHOGONAL
2260 C          TRANSFORMATIONS USED IN THE REDUCTION,
2261 C
2262 C        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX,
2263 C
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,
2266 C
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.
2269 C
2270 C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
2271 C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
2272 C
2273 C     ------------------------------------------------------------------
2274 C
2275 C     ********** FOR I=N STEP -1 UNTIL 1 DO -- **********
2276       DO 300 II = 1, N
2277       I = N + 1 - II
2278       L = I - 1
2279       IZ = (I * L) / 2
2280       H = 0.0D+00
2281       SCALE = 0.0D+00
2282       IF (L .LT. 1) GO TO 120
2283 C     ********** SCALE ROW (ALGOL TOL THEN NOT NEEDED) **********
2284       DO 100 K = 1, L
2285       IZ = IZ + 1
2286       D(K) = A(IZ)
2287       SCALE = SCALE + ABS(D(K))
2288   100 CONTINUE
2289 C
2290       IF (SCALE .NE. 0.0D+00) GO TO 140
2291   120 E(I) = 0.0D+00
2292       E2(I) = 0.0D+00
2293       GO TO 280
2294 C
2295   140 DO 160 K = 1, L
2296       D(K) = D(K) / SCALE
2297       H = H + D(K) * D(K)
2298   160 CONTINUE
2299 C
2300       E2(I) = SCALE * SCALE * H
2301       F = D(L)
2302       G = -SIGN(SQRT(H),F)
2303       E(I) = SCALE * G
2304       H = H - F * G
2305       D(L) = F - G
2306       A(IZ) = SCALE * D(L)
2307       IF (L .EQ. 1) GO TO 280
2308       F = 0.0D+00
2309 C
2310       JK = 1
2311       DO 220 J = 1, L
2312       JM1 = J - 1
2313       DT = D(J)
2314       G = 0.0D+00
2315 C     ********** FORM ELEMENT OF A*U **********
2316       IF (JM1 .EQ. 0) GO TO 200
2317       DO 180 K = 1, JM1
2318       E(K) = E(K) + DT * A(JK)
2319       G = G + D(K) * A(JK)
2320       JK = JK + 1
2321   180 CONTINUE
2322   200 E(J) = G + A(JK) * DT
2323       JK = JK + 1
2324 C     ********** FORM ELEMENT OF P **********
2325   220 CONTINUE
2326       F = 0.0D+00
2327       DO 240 J = 1, L
2328       E(J) = E(J) / H
2329       F = F + E(J) * D(J)
2330   240 CONTINUE
2331 C
2332       HH = F / (H + H)
2333       JK = 0
2334 C     ********** FORM REDUCED A **********
2335       DO 260 J = 1, L
2336       F = D(J)
2337       G = E(J) - HH * F
2338       E(J) = G
2339 C
2340       DO 260 K = 1, J
2341       JK = JK + 1
2342       A(JK) = A(JK) - F * E(K) - G * D(K)
2343   260 CONTINUE
2344 C
2345   280 D(I) = A(IZ+1)
2346       A(IZ+1) = SCALE * SQRT(H)
2347   300 CONTINUE
2348 C
2349       RETURN
2350 C     ********** LAST CARD OF TRED3 **********
2351       END