source: branches/stable/EISPACK/eispack.f

Last change on this file was 2, checked in by oldcode, 24 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 348.0 KB
Line 
1      SUBROUTINE CDIV(AR,AI,BR,BI,CR,CI)
2      DOUBLE PRECISION AR,AI,BR,BI,CR,CI
3C
4C     COMPLEX DIVISION, (CR,CI) = (AR,AI)/(BR,BI)
5C
6      DOUBLE PRECISION S,ARS,AIS,BRS,BIS
7      S = DABS(BR) + DABS(BI)
8      ARS = AR/S
9      AIS = AI/S
10      BRS = BR/S
11      BIS = BI/S
12      S = BRS**2 + BIS**2
13      CR = (ARS*BRS + AIS*BIS)/S
14      CI = (AIS*BRS - ARS*BIS)/S
15      RETURN
16      END
17      SUBROUTINE CSROOT(XR,XI,YR,YI)
18      DOUBLE PRECISION XR,XI,YR,YI
19C
20C     (YR,YI) = COMPLEX DSQRT(XR,XI) 
21C     BRANCH CHOSEN SO THAT YR .GE. 0.0 AND SIGN(YI) .EQ. SIGN(XI)
22C
23      DOUBLE PRECISION S,TR,TI,PYTHAG
24      TR = XR
25      TI = XI
26      S = DSQRT(0.5D0*(PYTHAG(TR,TI) + DABS(TR)))
27      IF (TR .GE. 0.0D0) YR = S
28      IF (TI .LT. 0.0D0) S = -S
29      IF (TR .LE. 0.0D0) YI = S
30      IF (TR .LT. 0.0D0) YR = 0.5D0*(TI/YI)
31      IF (TR .GT. 0.0D0) YI = 0.5D0*(TI/YR)
32      RETURN
33      END
34      DOUBLE PRECISION FUNCTION EPSLON (X)
35      DOUBLE PRECISION X
36C
37C     ESTIMATE UNIT ROUNDOFF IN QUANTITIES OF SIZE X.
38C
39      DOUBLE PRECISION A,B,C,EPS
40C
41C     THIS PROGRAM SHOULD FUNCTION PROPERLY ON ALL SYSTEMS
42C     SATISFYING THE FOLLOWING TWO ASSUMPTIONS,
43C        1.  THE BASE USED IN REPRESENTING FLOATING POINT
44C            NUMBERS IS NOT A POWER OF THREE.
45C        2.  THE QUANTITY  A  IN STATEMENT 10 IS REPRESENTED TO
46C            THE ACCURACY USED IN FLOATING POINT VARIABLES
47C            THAT ARE STORED IN MEMORY.
48C     THE STATEMENT NUMBER 10 AND THE GO TO 10 ARE INTENDED TO
49C     FORCE OPTIMIZING COMPILERS TO GENERATE CODE SATISFYING
50C     ASSUMPTION 2.
51C     UNDER THESE ASSUMPTIONS, IT SHOULD BE TRUE THAT,
52C            A  IS NOT EXACTLY EQUAL TO FOUR-THIRDS,
53C            B  HAS A ZERO FOR ITS LAST BIT OR DIGIT,
54C            C  IS NOT EXACTLY EQUAL TO ONE,
55C            EPS  MEASURES THE SEPARATION OF 1.0 FROM
56C                 THE NEXT LARGER FLOATING POINT NUMBER.
57C     THE DEVELOPERS OF EISPACK WOULD APPRECIATE BEING INFORMED
58C     ABOUT ANY SYSTEMS WHERE THESE ASSUMPTIONS DO NOT HOLD.
59C
60C     THIS VERSION DATED 4/6/83.
61C
62      A = 4.0D0/3.0D0
63   10 B = A - 1.0D0
64      C = B + B + B
65      EPS = DABS(C-1.0D0)
66      IF (EPS .EQ. 0.0D0) GO TO 10
67      EPSLON = EPS*DABS(X)
68      RETURN
69      END
70      DOUBLE PRECISION FUNCTION PYTHAG(A,B)
71      DOUBLE PRECISION A,B
72C
73C     FINDS DSQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW
74C
75      DOUBLE PRECISION P,R,S,T,U
76      P = DMAX1(DABS(A),DABS(B))
77      IF (P .EQ. 0.0D0) GO TO 20
78      R = (DMIN1(DABS(A),DABS(B))/P)**2
79   10 CONTINUE
80         T = 4.0D0 + R
81         IF (T .EQ. 4.0D0) GO TO 20
82         S = R/T
83         U = 1.0D0 + 2.0D0*S
84         P = U*P
85         R = (S/U)**2 * R
86      GO TO 10
87   20 PYTHAG = P
88      RETURN
89      END
90      SUBROUTINE BAKVEC(NM,N,T,E,M,Z,IERR)
91C
92      INTEGER I,J,M,N,NM,IERR
93      DOUBLE PRECISION T(NM,3),E(N),Z(NM,M)
94C
95C     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A NONSYMMETRIC
96C     TRIDIAGONAL MATRIX BY BACK TRANSFORMING THOSE OF THE
97C     CORRESPONDING SYMMETRIC MATRIX DETERMINED BY  FIGI.
98C
99C     ON INPUT
100C
101C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
102C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
103C          DIMENSION STATEMENT.
104C
105C        N IS THE ORDER OF THE MATRIX.
106C
107C        T CONTAINS THE NONSYMMETRIC MATRIX.  ITS SUBDIAGONAL IS
108C          STORED IN THE LAST N-1 POSITIONS OF THE FIRST COLUMN,
109C          ITS DIAGONAL IN THE N POSITIONS OF THE SECOND COLUMN,
110C          AND ITS SUPERDIAGONAL IN THE FIRST N-1 POSITIONS OF
111C          THE THIRD COLUMN.  T(1,1) AND T(N,3) ARE ARBITRARY.
112C
113C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE SYMMETRIC
114C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
115C
116C        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.
117C
118C        Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
119C          IN ITS FIRST M COLUMNS.
120C
121C     ON OUTPUT
122C
123C        T IS UNALTERED.
124C
125C        E IS DESTROYED.
126C
127C        Z CONTAINS THE TRANSFORMED EIGENVECTORS
128C          IN ITS FIRST M COLUMNS.
129C
130C        IERR IS SET TO
131C          ZERO       FOR NORMAL RETURN,
132C          2*N+I      IF E(I) IS ZERO WITH T(I,1) OR T(I-1,3) NON-ZERO.
133C                     IN THIS CASE, THE SYMMETRIC MATRIX IS NOT SIMILAR
134C                     TO THE ORIGINAL MATRIX, AND THE EIGENVECTORS
135C                     CANNOT BE FOUND BY THIS PROGRAM.
136C
137C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
138C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
139C
140C     THIS VERSION DATED AUGUST 1983.
141C
142C     ------------------------------------------------------------------
143C
144      IERR = 0
145      IF (M .EQ. 0) GO TO 1001
146      E(1) = 1.0D0
147      IF (N .EQ. 1) GO TO 1001
148C
149      DO 100 I = 2, N
150         IF (E(I) .NE. 0.0D0) GO TO 80
151         IF (T(I,1) .NE. 0.0D0 .OR. T(I-1,3) .NE. 0.0D0) GO TO 1000
152         E(I) = 1.0D0
153         GO TO 100
154   80    E(I) = E(I-1) * E(I) / T(I-1,3)
155  100 CONTINUE
156C
157      DO 120 J = 1, M
158C
159         DO 120 I = 2, N
160         Z(I,J) = Z(I,J) * E(I)
161  120 CONTINUE
162C
163      GO TO 1001
164C     .......... SET ERROR -- EIGENVECTORS CANNOT BE
165C                FOUND BY THIS PROGRAM ..........
166 1000 IERR = 2 * N + I
167 1001 RETURN
168      END
169      SUBROUTINE BALANC(NM,N,A,LOW,IGH,SCALE)
170C
171      INTEGER I,J,K,L,M,N,JJ,NM,IGH,LOW,IEXC
172      DOUBLE PRECISION A(NM,N),SCALE(N)
173      DOUBLE PRECISION C,F,G,R,S,B2,RADIX
174      LOGICAL NOCONV
175C
176C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BALANCE,
177C     NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH.
178C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971).
179C
180C     THIS SUBROUTINE BALANCES A REAL MATRIX AND ISOLATES
181C     EIGENVALUES WHENEVER POSSIBLE.
182C
183C     ON INPUT
184C
185C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
186C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
187C          DIMENSION STATEMENT.
188C
189C        N IS THE ORDER OF THE MATRIX.
190C
191C        A CONTAINS THE INPUT MATRIX TO BE BALANCED.
192C
193C     ON OUTPUT
194C
195C        A CONTAINS THE BALANCED MATRIX.
196C
197C        LOW AND IGH ARE TWO INTEGERS SUCH THAT A(I,J)
198C          IS EQUAL TO ZERO IF
199C           (1) I IS GREATER THAN J AND
200C           (2) J=1,...,LOW-1 OR I=IGH+1,...,N.
201C
202C        SCALE CONTAINS INFORMATION DETERMINING THE
203C           PERMUTATIONS AND SCALING FACTORS USED.
204C
205C     SUPPOSE THAT THE PRINCIPAL SUBMATRIX IN ROWS LOW THROUGH IGH
206C     HAS BEEN BALANCED, THAT P(J) DENOTES THE INDEX INTERCHANGED
207C     WITH J DURING THE PERMUTATION STEP, AND THAT THE ELEMENTS
208C     OF THE DIAGONAL MATRIX USED ARE DENOTED BY D(I,J).  THEN
209C        SCALE(J) = P(J),    FOR J = 1,...,LOW-1
210C                 = D(J,J),      J = LOW,...,IGH
211C                 = P(J)         J = IGH+1,...,N.
212C     THE ORDER IN WHICH THE INTERCHANGES ARE MADE IS N TO IGH+1,
213C     THEN 1 TO LOW-1.
214C
215C     NOTE THAT 1 IS RETURNED FOR IGH IF IGH IS ZERO FORMALLY.
216C
217C     THE ALGOL PROCEDURE EXC CONTAINED IN BALANCE APPEARS IN
218C     BALANC  IN LINE.  (NOTE THAT THE ALGOL ROLES OF IDENTIFIERS
219C     K,L HAVE BEEN REVERSED.)
220C
221C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
222C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
223C
224C     THIS VERSION DATED AUGUST 1983.
225C
226C     ------------------------------------------------------------------
227C
228      RADIX = 16.0D0
229C
230      B2 = RADIX * RADIX
231      K = 1
232      L = N
233      GO TO 100
234C     .......... IN-LINE PROCEDURE FOR ROW AND
235C                COLUMN EXCHANGE ..........
236   20 SCALE(M) = J
237      IF (J .EQ. M) GO TO 50
238C
239      DO 30 I = 1, L
240         F = A(I,J)
241         A(I,J) = A(I,M)
242         A(I,M) = F
243   30 CONTINUE
244C
245      DO 40 I = K, N
246         F = A(J,I)
247         A(J,I) = A(M,I)
248         A(M,I) = F
249   40 CONTINUE
250C
251   50 GO TO (80,130), IEXC
252C     .......... SEARCH FOR ROWS ISOLATING AN EIGENVALUE
253C                AND PUSH THEM DOWN ..........
254   80 IF (L .EQ. 1) GO TO 280
255      L = L - 1
256C     .......... FOR J=L STEP -1 UNTIL 1 DO -- ..........
257  100 DO 120 JJ = 1, L
258         J = L + 1 - JJ
259C
260         DO 110 I = 1, L
261            IF (I .EQ. J) GO TO 110
262            IF (A(J,I) .NE. 0.0D0) GO TO 120
263  110    CONTINUE
264C
265         M = L
266         IEXC = 1
267         GO TO 20
268  120 CONTINUE
269C
270      GO TO 140
271C     .......... SEARCH FOR COLUMNS ISOLATING AN EIGENVALUE
272C                AND PUSH THEM LEFT ..........
273  130 K = K + 1
274C
275  140 DO 170 J = K, L
276C
277         DO 150 I = K, L
278            IF (I .EQ. J) GO TO 150
279            IF (A(I,J) .NE. 0.0D0) GO TO 170
280  150    CONTINUE
281C
282         M = K
283         IEXC = 2
284         GO TO 20
285  170 CONTINUE
286C     .......... NOW BALANCE THE SUBMATRIX IN ROWS K TO L ..........
287      DO 180 I = K, L
288  180 SCALE(I) = 1.0D0
289C     .......... ITERATIVE LOOP FOR NORM REDUCTION ..........
290  190 NOCONV = .FALSE.
291C
292      DO 270 I = K, L
293         C = 0.0D0
294         R = 0.0D0
295C
296         DO 200 J = K, L
297            IF (J .EQ. I) GO TO 200
298            C = C + DABS(A(J,I))
299            R = R + DABS(A(I,J))
300  200    CONTINUE
301C     .......... GUARD AGAINST ZERO C OR R DUE TO UNDERFLOW ..........
302         IF (C .EQ. 0.0D0 .OR. R .EQ. 0.0D0) GO TO 270
303         G = R / RADIX
304         F = 1.0D0
305         S = C + R
306  210    IF (C .GE. G) GO TO 220
307         F = F * RADIX
308         C = C * B2
309         GO TO 210
310  220    G = R * RADIX
311  230    IF (C .LT. G) GO TO 240
312         F = F / RADIX
313         C = C / B2
314         GO TO 230
315C     .......... NOW BALANCE ..........
316  240    IF ((C + R) / F .GE. 0.95D0 * S) GO TO 270
317         G = 1.0D0 / F
318         SCALE(I) = SCALE(I) * F
319         NOCONV = .TRUE.
320C
321         DO 250 J = K, N
322  250    A(I,J) = A(I,J) * G
323C
324         DO 260 J = 1, L
325  260    A(J,I) = A(J,I) * F
326C
327  270 CONTINUE
328C
329      IF (NOCONV) GO TO 190
330C
331  280 LOW = K
332      IGH = L
333      RETURN
334      END
335      SUBROUTINE BALBAK(NM,N,LOW,IGH,SCALE,M,Z)
336C
337      INTEGER I,J,K,M,N,II,NM,IGH,LOW
338      DOUBLE PRECISION SCALE(N),Z(NM,M)
339      DOUBLE PRECISION S
340C
341C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BALBAK,
342C     NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH.
343C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971).
344C
345C     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A REAL GENERAL
346C     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
347C     BALANCED MATRIX DETERMINED BY  BALANC.
348C
349C     ON INPUT
350C
351C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
352C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
353C          DIMENSION STATEMENT.
354C
355C        N IS THE ORDER OF THE MATRIX.
356C
357C        LOW AND IGH ARE INTEGERS DETERMINED BY  BALANC.
358C
359C        SCALE CONTAINS INFORMATION DETERMINING THE PERMUTATIONS
360C          AND SCALING FACTORS USED BY  BALANC.
361C
362C        M IS THE NUMBER OF COLUMNS OF Z TO BE BACK TRANSFORMED.
363C
364C        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGEN-
365C          VECTORS TO BE BACK TRANSFORMED IN ITS FIRST M COLUMNS.
366C
367C     ON OUTPUT
368C
369C        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE
370C          TRANSFORMED EIGENVECTORS IN ITS FIRST M COLUMNS.
371C
372C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
373C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
374C
375C     THIS VERSION DATED AUGUST 1983.
376C
377C     ------------------------------------------------------------------
378C
379      IF (M .EQ. 0) GO TO 200
380      IF (IGH .EQ. LOW) GO TO 120
381C
382      DO 110 I = LOW, IGH
383         S = SCALE(I)
384C     .......... LEFT HAND EIGENVECTORS ARE BACK TRANSFORMED
385C                IF THE FOREGOING STATEMENT IS REPLACED BY
386C                S=1.0D0/SCALE(I). ..........
387         DO 100 J = 1, M
388  100    Z(I,J) = Z(I,J) * S
389C
390  110 CONTINUE
391C     ......... FOR I=LOW-1 STEP -1 UNTIL 1,
392C               IGH+1 STEP 1 UNTIL N DO -- ..........
393  120 DO 140 II = 1, N
394         I = II
395         IF (I .GE. LOW .AND. I .LE. IGH) GO TO 140
396         IF (I .LT. LOW) I = LOW - II
397         K = SCALE(I)
398         IF (K .EQ. I) GO TO 140
399C
400         DO 130 J = 1, M
401            S = Z(I,J)
402            Z(I,J) = Z(K,J)
403            Z(K,J) = S
404  130    CONTINUE
405C
406  140 CONTINUE
407C
408  200 RETURN
409      END
410      SUBROUTINE BANDR(NM,N,MB,A,D,E,E2,MATZ,Z)
411C
412      INTEGER J,K,L,N,R,I1,I2,J1,J2,KR,MB,MR,M1,NM,N2,R1,UGL,MAXL,MAXR
413      DOUBLE PRECISION A(NM,MB),D(N),E(N),E2(N),Z(NM,N)
414      DOUBLE PRECISION G,U,B1,B2,C2,F1,F2,S2,DMIN,DMINRT
415      LOGICAL MATZ
416C
417C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BANDRD,
418C     NUM. MATH. 12, 231-241(1968) BY SCHWARZ.
419C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 273-283(1971).
420C
421C     THIS SUBROUTINE REDUCES A REAL SYMMETRIC BAND MATRIX
422C     TO A SYMMETRIC TRIDIAGONAL MATRIX USING AND OPTIONALLY
423C     ACCUMULATING ORTHOGONAL SIMILARITY TRANSFORMATIONS.
424C
425C     ON INPUT
426C
427C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
428C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
429C          DIMENSION STATEMENT.
430C
431C        N IS THE ORDER OF THE MATRIX.
432C
433C        MB IS THE (HALF) BAND WIDTH OF THE MATRIX, DEFINED AS THE
434C          NUMBER OF ADJACENT DIAGONALS, INCLUDING THE PRINCIPAL
435C          DIAGONAL, REQUIRED TO SPECIFY THE NON-ZERO PORTION OF THE
436C          LOWER TRIANGLE OF THE MATRIX.
437C
438C        A CONTAINS THE LOWER TRIANGLE OF THE SYMMETRIC BAND INPUT
439C          MATRIX STORED AS AN N BY MB ARRAY.  ITS LOWEST SUBDIAGONAL
440C          IS STORED IN THE LAST N+1-MB POSITIONS OF THE FIRST COLUMN,
441C          ITS NEXT SUBDIAGONAL IN THE LAST N+2-MB POSITIONS OF THE
442C          SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND FINALLY
443C          ITS PRINCIPAL DIAGONAL IN THE N POSITIONS OF THE LAST COLUMN.
444C          CONTENTS OF STORAGES NOT PART OF THE MATRIX ARE ARBITRARY.
445C
446C        MATZ SHOULD BE SET TO .TRUE. IF THE TRANSFORMATION MATRIX IS
447C          TO BE ACCUMULATED, AND TO .FALSE. OTHERWISE.
448C
449C     ON OUTPUT
450C
451C        A HAS BEEN DESTROYED, EXCEPT FOR ITS LAST TWO COLUMNS WHICH
452C          CONTAIN A COPY OF THE TRIDIAGONAL MATRIX.
453C
454C        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX.
455C
456C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
457C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO.
458C
459C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
460C          E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
461C
462C        Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX PRODUCED IN
463C          THE REDUCTION IF MATZ HAS BEEN SET TO .TRUE.  OTHERWISE, Z
464C          IS NOT REFERENCED.
465C
466C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
467C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
468C
469C     THIS VERSION DATED AUGUST 1983.
470C
471C     ------------------------------------------------------------------
472C
473      DMIN = 2.0D0**(-64)
474      DMINRT = 2.0D0**(-32)
475C     .......... INITIALIZE DIAGONAL SCALING MATRIX ..........
476      DO 30 J = 1, N
477   30 D(J) = 1.0D0
478C
479      IF (.NOT. MATZ) GO TO 60
480C
481      DO 50 J = 1, N
482C
483         DO 40 K = 1, N
484   40    Z(J,K) = 0.0D0
485C
486         Z(J,J) = 1.0D0
487   50 CONTINUE
488C
489   60 M1 = MB - 1
490      IF (M1 - 1) 900, 800, 70
491   70 N2 = N - 2
492C
493      DO 700 K = 1, N2
494         MAXR = MIN0(M1,N-K)
495C     .......... FOR R=MAXR STEP -1 UNTIL 2 DO -- ..........
496         DO 600 R1 = 2, MAXR
497            R = MAXR + 2 - R1
498            KR = K + R
499            MR = MB - R
500            G = A(KR,MR)
501            A(KR-1,1) = A(KR-1,MR+1)
502            UGL = K
503C
504            DO 500 J = KR, N, M1
505               J1 = J - 1
506               J2 = J1 - 1
507               IF (G .EQ. 0.0D0) GO TO 600
508               B1 = A(J1,1) / G
509               B2 = B1 * D(J1) / D(J)
510               S2 = 1.0D0 / (1.0D0 + B1 * B2)
511               IF (S2 .GE. 0.5D0 ) GO TO 450
512               B1 = G / A(J1,1)
513               B2 = B1 * D(J) / D(J1)
514               C2 = 1.0D0 - S2
515               D(J1) = C2 * D(J1)
516               D(J) = C2 * D(J)
517               F1 = 2.0D0 * A(J,M1)
518               F2 = B1 * A(J1,MB)
519               A(J,M1) = -B2 * (B1 * A(J,M1) - A(J,MB)) - F2 + A(J,M1)
520               A(J1,MB) = B2 * (B2 * A(J,MB) + F1) + A(J1,MB)
521               A(J,MB) = B1 * (F2 - F1) + A(J,MB)
522C
523               DO 200 L = UGL, J2
524                  I2 = MB - J + L
525                  U = A(J1,I2+1) + B2 * A(J,I2)
526                  A(J,I2) = -B1 * A(J1,I2+1) + A(J,I2)
527                  A(J1,I2+1) = U
528  200          CONTINUE
529C
530               UGL = J
531               A(J1,1) = A(J1,1) + B2 * G
532               IF (J .EQ. N) GO TO 350
533               MAXL = MIN0(M1,N-J1)
534C
535               DO 300 L = 2, MAXL
536                  I1 = J1 + L
537                  I2 = MB - L
538                  U = A(I1,I2) + B2 * A(I1,I2+1)
539                  A(I1,I2+1) = -B1 * A(I1,I2) + A(I1,I2+1)
540                  A(I1,I2) = U
541  300          CONTINUE
542C
543               I1 = J + M1
544               IF (I1 .GT. N) GO TO 350
545               G = B2 * A(I1,1)
546  350          IF (.NOT. MATZ) GO TO 500
547C
548               DO 400 L = 1, N
549                  U = Z(L,J1) + B2 * Z(L,J)
550                  Z(L,J) = -B1 * Z(L,J1) + Z(L,J)
551                  Z(L,J1) = U
552  400          CONTINUE
553C
554               GO TO 500
555C
556  450          U = D(J1)
557               D(J1) = S2 * D(J)
558               D(J) = S2 * U
559               F1 = 2.0D0 * A(J,M1)
560               F2 = B1 * A(J,MB)
561               U = B1 * (F2 - F1) + A(J1,MB)
562               A(J,M1) = B2 * (B1 * A(J,M1) - A(J1,MB)) + F2 - A(J,M1)
563               A(J1,MB) = B2 * (B2 * A(J1,MB) + F1) + A(J,MB)
564               A(J,MB) = U
565C
566               DO 460 L = UGL, J2
567                  I2 = MB - J + L
568                  U = B2 * A(J1,I2+1) + A(J,I2)
569                  A(J,I2) = -A(J1,I2+1) + B1 * A(J,I2)
570                  A(J1,I2+1) = U
571  460          CONTINUE
572C
573               UGL = J
574               A(J1,1) = B2 * A(J1,1) + G
575               IF (J .EQ. N) GO TO 480
576               MAXL = MIN0(M1,N-J1)
577C
578               DO 470 L = 2, MAXL
579                  I1 = J1 + L
580                  I2 = MB - L
581                  U = B2 * A(I1,I2) + A(I1,I2+1)
582                  A(I1,I2+1) = -A(I1,I2) + B1 * A(I1,I2+1)
583                  A(I1,I2) = U
584  470          CONTINUE
585C
586               I1 = J + M1
587               IF (I1 .GT. N) GO TO 480
588               G = A(I1,1)
589               A(I1,1) = B1 * A(I1,1)
590  480          IF (.NOT. MATZ) GO TO 500
591C
592               DO 490 L = 1, N
593                  U = B2 * Z(L,J1) + Z(L,J)
594                  Z(L,J) = -Z(L,J1) + B1 * Z(L,J)
595                  Z(L,J1) = U
596  490          CONTINUE
597C
598  500       CONTINUE
599C
600  600    CONTINUE
601C
602         IF (MOD(K,64) .NE. 0) GO TO 700
603C     .......... RESCALE TO AVOID UNDERFLOW OR OVERFLOW ..........
604         DO 650 J = K, N
605            IF (D(J) .GE. DMIN) GO TO 650
606            MAXL = MAX0(1,MB+1-J)
607C
608            DO 610 L = MAXL, M1
609  610       A(J,L) = DMINRT * A(J,L)
610C
611            IF (J .EQ. N) GO TO 630
612            MAXL = MIN0(M1,N-J)
613C
614            DO 620 L = 1, MAXL
615               I1 = J + L
616               I2 = MB - L
617               A(I1,I2) = DMINRT * A(I1,I2)
618  620       CONTINUE
619C
620  630       IF (.NOT. MATZ) GO TO 645
621C
622            DO 640 L = 1, N
623  640       Z(L,J) = DMINRT * Z(L,J)
624C
625  645       A(J,MB) = DMIN * A(J,MB)
626            D(J) = D(J) / DMIN
627  650    CONTINUE
628C
629  700 CONTINUE
630C     .......... FORM SQUARE ROOT OF SCALING MATRIX ..........
631  800 DO 810 J = 2, N
632  810 E(J) = DSQRT(D(J))
633C
634      IF (.NOT. MATZ) GO TO 840
635C
636      DO 830 J = 1, N
637C
638         DO 820 K = 2, N
639  820    Z(J,K) = E(K) * Z(J,K)
640C
641  830 CONTINUE
642C
643  840 U = 1.0D0
644C
645      DO 850 J = 2, N
646         A(J,M1) = U * E(J) * A(J,M1)
647         U = E(J)
648         E2(J) = A(J,M1) ** 2
649         A(J,MB) = D(J) * A(J,MB)
650         D(J) = A(J,MB)
651         E(J) = A(J,M1)
652  850 CONTINUE
653C
654      D(1) = A(1,MB)
655      E(1) = 0.0D0
656      E2(1) = 0.0D0
657      GO TO 1001
658C
659  900 DO 950 J = 1, N
660         D(J) = A(J,MB)
661         E(J) = 0.0D0
662         E2(J) = 0.0D0
663  950 CONTINUE
664C
665 1001 RETURN
666      END
667      SUBROUTINE BANDV(NM,N,MBW,A,E21,M,W,Z,IERR,NV,RV,RV6)
668C
669      INTEGER I,J,K,M,N,R,II,IJ,JJ,KJ,MB,M1,NM,NV,IJ1,ITS,KJ1,MBW,M21,
670     X        IERR,MAXJ,MAXK,GROUP
671      DOUBLE PRECISION A(NM,MBW),W(M),Z(NM,M),RV(NV),RV6(N)
672      DOUBLE PRECISION U,V,UK,XU,X0,X1,E21,EPS2,EPS3,EPS4,NORM,ORDER,
673     X       EPSLON,PYTHAG
674C
675C     THIS SUBROUTINE FINDS THOSE EIGENVECTORS OF A REAL SYMMETRIC
676C     BAND MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES, USING INVERSE
677C     ITERATION.  THE SUBROUTINE MAY ALSO BE USED TO SOLVE SYSTEMS
678C     OF LINEAR EQUATIONS WITH A SYMMETRIC OR NON-SYMMETRIC BAND
679C     COEFFICIENT MATRIX.
680C
681C     ON INPUT
682C
683C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
684C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
685C          DIMENSION STATEMENT.
686C
687C        N IS THE ORDER OF THE MATRIX.
688C
689C        MBW IS THE NUMBER OF COLUMNS OF THE ARRAY A USED TO STORE THE
690C          BAND MATRIX.  IF THE MATRIX IS SYMMETRIC, MBW IS ITS (HALF)
691C          BAND WIDTH, DENOTED MB AND DEFINED AS THE NUMBER OF ADJACENT
692C          DIAGONALS, INCLUDING THE PRINCIPAL DIAGONAL, REQUIRED TO
693C          SPECIFY THE NON-ZERO PORTION OF THE LOWER TRIANGLE OF THE
694C          MATRIX.  IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS
695C          OF LINEAR EQUATIONS AND THE COEFFICIENT MATRIX IS NOT
696C          SYMMETRIC, IT MUST HOWEVER HAVE THE SAME NUMBER OF ADJACENT
697C          DIAGONALS ABOVE THE MAIN DIAGONAL AS BELOW, AND IN THIS
698C          CASE, MBW=2*MB-1.
699C
700C        A CONTAINS THE LOWER TRIANGLE OF THE SYMMETRIC BAND INPUT
701C          MATRIX STORED AS AN N BY MB ARRAY.  ITS LOWEST SUBDIAGONAL
702C          IS STORED IN THE LAST N+1-MB POSITIONS OF THE FIRST COLUMN,
703C          ITS NEXT SUBDIAGONAL IN THE LAST N+2-MB POSITIONS OF THE
704C          SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND FINALLY
705C          ITS PRINCIPAL DIAGONAL IN THE N POSITIONS OF COLUMN MB.
706C          IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS OF LINEAR
707C          EQUATIONS AND THE COEFFICIENT MATRIX IS NOT SYMMETRIC, A IS
708C          N BY 2*MB-1 INSTEAD WITH LOWER TRIANGLE AS ABOVE AND WITH
709C          ITS FIRST SUPERDIAGONAL STORED IN THE FIRST N-1 POSITIONS OF
710C          COLUMN MB+1, ITS SECOND SUPERDIAGONAL IN THE FIRST N-2
711C          POSITIONS OF COLUMN MB+2, FURTHER SUPERDIAGONALS SIMILARLY,
712C          AND FINALLY ITS HIGHEST SUPERDIAGONAL IN THE FIRST N+1-MB
713C          POSITIONS OF THE LAST COLUMN.
714C          CONTENTS OF STORAGES NOT PART OF THE MATRIX ARE ARBITRARY.
715C
716C        E21 SPECIFIES THE ORDERING OF THE EIGENVALUES AND CONTAINS
717C            0.0D0 IF THE EIGENVALUES ARE IN ASCENDING ORDER, OR
718C            2.0D0 IF THE EIGENVALUES ARE IN DESCENDING ORDER.
719C          IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS OF LINEAR
720C          EQUATIONS, E21 SHOULD BE SET TO 1.0D0 IF THE COEFFICIENT
721C          MATRIX IS SYMMETRIC AND TO -1.0D0 IF NOT.
722C
723C        M IS THE NUMBER OF SPECIFIED EIGENVALUES OR THE NUMBER OF
724C          SYSTEMS OF LINEAR EQUATIONS.
725C
726C        W CONTAINS THE M EIGENVALUES IN ASCENDING OR DESCENDING ORDER.
727C          IF THE SUBROUTINE IS BEING USED TO SOLVE SYSTEMS OF LINEAR
728C          EQUATIONS (A-W(R)*I)*X(R)=B(R), WHERE I IS THE IDENTITY
729C          MATRIX, W(R) SHOULD BE SET ACCORDINGLY, FOR R=1,2,...,M.
730C
731C        Z CONTAINS THE CONSTANT MATRIX COLUMNS (B(R),R=1,2,...,M), IF
732C          THE SUBROUTINE IS USED TO SOLVE SYSTEMS OF LINEAR EQUATIONS.
733C
734C        NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER RV
735C          AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT.
736C
737C     ON OUTPUT
738C
739C        A AND W ARE UNALTERED.
740C
741C        Z CONTAINS THE ASSOCIATED SET OF ORTHOGONAL EIGENVECTORS.
742C          ANY VECTOR WHICH FAILS TO CONVERGE IS SET TO ZERO.  IF THE
743C          SUBROUTINE IS USED TO SOLVE SYSTEMS OF LINEAR EQUATIONS,
744C          Z CONTAINS THE SOLUTION MATRIX COLUMNS (X(R),R=1,2,...,M).
745C
746C        IERR IS SET TO
747C          ZERO       FOR NORMAL RETURN,
748C          -R         IF THE EIGENVECTOR CORRESPONDING TO THE R-TH
749C                     EIGENVALUE FAILS TO CONVERGE, OR IF THE R-TH
750C                     SYSTEM OF LINEAR EQUATIONS IS NEARLY SINGULAR.
751C
752C        RV AND RV6 ARE TEMPORARY STORAGE ARRAYS.  NOTE THAT RV IS
753C          OF DIMENSION AT LEAST N*(2*MB-1).  IF THE SUBROUTINE
754C          IS BEING USED TO SOLVE SYSTEMS OF LINEAR EQUATIONS, THE
755C          DETERMINANT (UP TO SIGN) OF A-W(M)*I IS AVAILABLE, UPON
756C          RETURN, AS THE PRODUCT OF THE FIRST N ELEMENTS OF RV.
757C
758C     CALLS PYTHAG FOR  DSQRT(A*A + B*B) .
759C
760C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
761C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
762C
763C     THIS VERSION DATED AUGUST 1983.
764C
765C     ------------------------------------------------------------------
766C
767      IERR = 0
768      IF (M .EQ. 0) GO TO 1001
769      MB = MBW
770      IF (E21 .LT. 0.0D0) MB = (MBW + 1) / 2
771      M1 = MB - 1
772      M21 = M1 + MB
773      ORDER = 1.0D0 - DABS(E21)
774C     .......... FIND VECTORS BY INVERSE ITERATION ..........
775      DO 920 R = 1, M
776         ITS = 1
777         X1 = W(R)
778         IF (R .NE. 1) GO TO 100
779C     .......... COMPUTE NORM OF MATRIX ..........
780         NORM = 0.0D0
781C
782         DO 60 J = 1, MB
783            JJ = MB + 1 - J
784            KJ = JJ + M1
785            IJ = 1
786            V = 0.0D0
787C
788            DO 40 I = JJ, N
789               V = V + DABS(A(I,J))
790               IF (E21 .GE. 0.0D0) GO TO 40
791               V = V + DABS(A(IJ,KJ))
792               IJ = IJ + 1
793   40       CONTINUE
794C
795            NORM = DMAX1(NORM,V)
796   60    CONTINUE
797C
798         IF (E21 .LT. 0.0D0) NORM = 0.5D0 * NORM
799C     .......... EPS2 IS THE CRITERION FOR GROUPING,
800C                EPS3 REPLACES ZERO PIVOTS AND EQUAL
801C                ROOTS ARE MODIFIED BY EPS3,
802C                EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW ..........
803         IF (NORM .EQ. 0.0D0) NORM = 1.0D0
804         EPS2 = 1.0D-3 * NORM * DABS(ORDER)
805         EPS3 = EPSLON(NORM)
806         UK = N
807         UK = DSQRT(UK)
808         EPS4 = UK * EPS3
809   80    GROUP = 0
810         GO TO 120
811C     .......... LOOK FOR CLOSE OR COINCIDENT ROOTS ..........
812  100    IF (DABS(X1-X0) .GE. EPS2) GO TO 80
813         GROUP = GROUP + 1
814         IF (ORDER * (X1 - X0) .LE. 0.0D0) X1 = X0 + ORDER * EPS3
815C     .......... EXPAND MATRIX, SUBTRACT EIGENVALUE,
816C                AND INITIALIZE VECTOR ..........
817  120    DO 200 I = 1, N
818            IJ = I + MIN0(0,I-M1) * N
819            KJ = IJ + MB * N
820            IJ1 = KJ + M1 * N
821            IF (M1 .EQ. 0) GO TO 180
822C
823            DO 150 J = 1, M1
824               IF (IJ .GT. M1) GO TO 125
825               IF (IJ .GT. 0) GO TO 130
826               RV(IJ1) = 0.0D0
827               IJ1 = IJ1 + N
828               GO TO 130
829  125          RV(IJ) = A(I,J)
830  130          IJ = IJ + N
831               II = I + J
832               IF (II .GT. N) GO TO 150
833               JJ = MB - J
834               IF (E21 .GE. 0.0D0) GO TO 140
835               II = I
836               JJ = MB + J
837  140          RV(KJ) = A(II,JJ)
838               KJ = KJ + N
839  150       CONTINUE
840C
841  180       RV(IJ) = A(I,MB) - X1
842            RV6(I) = EPS4
843            IF (ORDER .EQ. 0.0D0) RV6(I) = Z(I,R)
844  200    CONTINUE
845C
846         IF (M1 .EQ. 0) GO TO 600
847C     .......... ELIMINATION WITH INTERCHANGES ..........
848         DO 580 I = 1, N
849            II = I + 1
850            MAXK = MIN0(I+M1-1,N)
851            MAXJ = MIN0(N-I,M21-2) * N
852C
853            DO 360 K = I, MAXK
854               KJ1 = K
855               J = KJ1 + N
856               JJ = J + MAXJ
857C
858               DO 340 KJ = J, JJ, N
859                  RV(KJ1) = RV(KJ)
860                  KJ1 = KJ
861  340          CONTINUE
862C
863               RV(KJ1) = 0.0D0
864  360       CONTINUE
865C
866            IF (I .EQ. N) GO TO 580
867            U = 0.0D0
868            MAXK = MIN0(I+M1,N)
869            MAXJ = MIN0(N-II,M21-2) * N
870C
871            DO 450 J = I, MAXK
872               IF (DABS(RV(J)) .LT. DABS(U)) GO TO 450
873               U = RV(J)
874               K = J
875  450       CONTINUE
876C
877            J = I + N
878            JJ = J + MAXJ
879            IF (K .EQ. I) GO TO 520
880            KJ = K
881C
882            DO 500 IJ = I, JJ, N
883               V = RV(IJ)
884               RV(IJ) = RV(KJ)
885               RV(KJ) = V
886               KJ = KJ + N
887  500       CONTINUE
888C
889            IF (ORDER .NE. 0.0D0) GO TO 520
890            V = RV6(I)
891            RV6(I) = RV6(K)
892            RV6(K) = V
893  520       IF (U .EQ. 0.0D0) GO TO 580
894C
895            DO 560 K = II, MAXK
896               V = RV(K) / U
897               KJ = K
898C
899               DO 540 IJ = J, JJ, N
900                  KJ = KJ + N
901                  RV(KJ) = RV(KJ) - V * RV(IJ)
902  540          CONTINUE
903C
904               IF (ORDER .EQ. 0.0D0) RV6(K) = RV6(K) - V * RV6(I)
905  560       CONTINUE
906C
907  580    CONTINUE
908C     .......... BACK SUBSTITUTION
909C                FOR I=N STEP -1 UNTIL 1 DO -- ..........
910  600    DO 630 II = 1, N
911            I = N + 1 - II
912            MAXJ = MIN0(II,M21)
913            IF (MAXJ .EQ. 1) GO TO 620
914            IJ1 = I
915            J = IJ1 + N
916            JJ = J + (MAXJ - 2) * N
917C
918            DO 610 IJ = J, JJ, N
919               IJ1 = IJ1 + 1
920               RV6(I) = RV6(I) - RV(IJ) * RV6(IJ1)
921  610       CONTINUE
922C
923  620       V = RV(I)
924            IF (DABS(V) .GE. EPS3) GO TO 625
925C     .......... SET ERROR -- NEARLY SINGULAR LINEAR SYSTEM ..........
926            IF (ORDER .EQ. 0.0D0) IERR = -R
927            V = DSIGN(EPS3,V)
928  625       RV6(I) = RV6(I) / V
929  630    CONTINUE
930C
931         XU = 1.0D0
932         IF (ORDER .EQ. 0.0D0) GO TO 870
933C     .......... ORTHOGONALIZE WITH RESPECT TO PREVIOUS
934C                MEMBERS OF GROUP ..........
935         IF (GROUP .EQ. 0) GO TO 700
936C
937         DO 680 JJ = 1, GROUP
938            J = R - GROUP - 1 + JJ
939            XU = 0.0D0
940C
941            DO 640 I = 1, N
942  640       XU = XU + RV6(I) * Z(I,J)
943C
944            DO 660 I = 1, N
945  660       RV6(I) = RV6(I) - XU * Z(I,J)
946C
947  680    CONTINUE
948C
949  700    NORM = 0.0D0
950C
951         DO 720 I = 1, N
952  720    NORM = NORM + DABS(RV6(I))
953C
954         IF (NORM .GE. 0.1D0) GO TO 840
955C     .......... IN-LINE PROCEDURE FOR CHOOSING
956C                A NEW STARTING VECTOR ..........
957         IF (ITS .GE. N) GO TO 830
958         ITS = ITS + 1
959         XU = EPS4 / (UK + 1.0D0)
960         RV6(1) = EPS4
961C
962         DO 760 I = 2, N
963  760    RV6(I) = XU
964C
965         RV6(ITS) = RV6(ITS) - EPS4 * UK
966         GO TO 600
967C     .......... SET ERROR -- NON-CONVERGED EIGENVECTOR ..........
968  830    IERR = -R
969         XU = 0.0D0
970         GO TO 870
971C     .......... NORMALIZE SO THAT SUM OF SQUARES IS
972C                1 AND EXPAND TO FULL ORDER ..........
973  840    U = 0.0D0
974C
975         DO 860 I = 1, N
976  860    U = PYTHAG(U,RV6(I))
977C
978         XU = 1.0D0 / U
979C
980  870    DO 900 I = 1, N
981  900    Z(I,R) = RV6(I) * XU
982C
983         X0 = X1
984  920 CONTINUE
985C
986 1001 RETURN
987      END
988      SUBROUTINE BISECT(N,EPS1,D,E,E2,LB,UB,MM,M,W,IND,IERR,RV4,RV5)
989C
990      INTEGER I,J,K,L,M,N,P,Q,R,S,II,MM,M1,M2,TAG,IERR,ISTURM
991      DOUBLE PRECISION D(N),E(N),E2(N),W(MM),RV4(N),RV5(N)
992      DOUBLE PRECISION U,V,LB,T1,T2,UB,XU,X0,X1,EPS1,TST1,TST2,EPSLON
993      INTEGER IND(MM)
994C
995C     THIS SUBROUTINE IS A TRANSLATION OF THE BISECTION TECHNIQUE
996C     IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON.
997C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971).
998C
999C     THIS SUBROUTINE FINDS THOSE EIGENVALUES OF A TRIDIAGONAL
1000C     SYMMETRIC MATRIX WHICH LIE IN A SPECIFIED INTERVAL,
1001C     USING BISECTION.
1002C
1003C     ON INPUT
1004C
1005C        N IS THE ORDER OF THE MATRIX.
1006C
1007C        EPS1 IS AN ABSOLUTE ERROR TOLERANCE FOR THE COMPUTED
1008C          EIGENVALUES.  IF THE INPUT EPS1 IS NON-POSITIVE,
1009C          IT IS RESET FOR EACH SUBMATRIX TO A DEFAULT VALUE,
1010C          NAMELY, MINUS THE PRODUCT OF THE RELATIVE MACHINE
1011C          PRECISION AND THE 1-NORM OF THE SUBMATRIX.
1012C
1013C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
1014C
1015C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
1016C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
1017C
1018C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
1019C          E2(1) IS ARBITRARY.
1020C
1021C        LB AND UB DEFINE THE INTERVAL TO BE SEARCHED FOR EIGENVALUES.
1022C          IF LB IS NOT LESS THAN UB, NO EIGENVALUES WILL BE FOUND.
1023C
1024C        MM SHOULD BE SET TO AN UPPER BOUND FOR THE NUMBER OF
1025C          EIGENVALUES IN THE INTERVAL.  WARNING. IF MORE THAN
1026C          MM EIGENVALUES ARE DETERMINED TO LIE IN THE INTERVAL,
1027C          AN ERROR RETURN IS MADE WITH NO EIGENVALUES FOUND.
1028C
1029C     ON OUTPUT
1030C
1031C        EPS1 IS UNALTERED UNLESS IT HAS BEEN RESET TO ITS
1032C          (LAST) DEFAULT VALUE.
1033C
1034C        D AND E ARE UNALTERED.
1035C
1036C        ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED
1037C          AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE
1038C          MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES.
1039C          E2(1) IS ALSO SET TO ZERO.
1040C
1041C        M IS THE NUMBER OF EIGENVALUES DETERMINED TO LIE IN (LB,UB).
1042C
1043C        W CONTAINS THE M EIGENVALUES IN ASCENDING ORDER.
1044C
1045C        IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES
1046C          ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W --
1047C          1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM
1048C          THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC..
1049C
1050C        IERR IS SET TO
1051C          ZERO       FOR NORMAL RETURN,
1052C          3*N+1      IF M EXCEEDS MM.
1053C
1054C        RV4 AND RV5 ARE TEMPORARY STORAGE ARRAYS.
1055C
1056C     THE ALGOL PROCEDURE STURMCNT CONTAINED IN TRISTURM
1057C     APPEARS IN BISECT IN-LINE.
1058C
1059C     NOTE THAT SUBROUTINE TQL1 OR IMTQL1 IS GENERALLY FASTER THAN
1060C     BISECT, IF MORE THAN N/4 EIGENVALUES ARE TO BE FOUND.
1061C
1062C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
1063C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
1064C
1065C     THIS VERSION DATED AUGUST 1983.
1066C
1067C     ------------------------------------------------------------------
1068C
1069      IERR = 0
1070      TAG = 0
1071      T1 = LB
1072      T2 = UB
1073C     .......... LOOK FOR SMALL SUB-DIAGONAL ENTRIES ..........
1074      DO 40 I = 1, N
1075         IF (I .EQ. 1) GO TO 20
1076         TST1 = DABS(D(I)) + DABS(D(I-1))
1077         TST2 = TST1 + DABS(E(I))
1078         IF (TST2 .GT. TST1) GO TO 40
1079   20    E2(I) = 0.0D0
1080   40 CONTINUE
1081C     .......... DETERMINE THE NUMBER OF EIGENVALUES
1082C                IN THE INTERVAL ..........
1083      P = 1
1084      Q = N
1085      X1 = UB
1086      ISTURM = 1
1087      GO TO 320
1088   60 M = S
1089      X1 = LB
1090      ISTURM = 2
1091      GO TO 320
1092   80 M = M - S
1093      IF (M .GT. MM) GO TO 980
1094      Q = 0
1095      R = 0
1096C     .......... ESTABLISH AND PROCESS NEXT SUBMATRIX, REFINING
1097C                INTERVAL BY THE GERSCHGORIN BOUNDS ..........
1098  100 IF (R .EQ. M) GO TO 1001
1099      TAG = TAG + 1
1100      P = Q + 1
1101      XU = D(P)
1102      X0 = D(P)
1103      U = 0.0D0
1104C
1105      DO 120 Q = P, N
1106         X1 = U
1107         U = 0.0D0
1108         V = 0.0D0
1109         IF (Q .EQ. N) GO TO 110
1110         U = DABS(E(Q+1))
1111         V = E2(Q+1)
1112  110    XU = DMIN1(D(Q)-(X1+U),XU)
1113         X0 = DMAX1(D(Q)+(X1+U),X0)
1114         IF (V .EQ. 0.0D0) GO TO 140
1115  120 CONTINUE
1116C
1117  140 X1 = EPSLON(DMAX1(DABS(XU),DABS(X0)))
1118      IF (EPS1 .LE. 0.0D0) EPS1 = -X1
1119      IF (P .NE. Q) GO TO 180
1120C     .......... CHECK FOR ISOLATED ROOT WITHIN INTERVAL ..........
1121      IF (T1 .GT. D(P) .OR. D(P) .GE. T2) GO TO 940
1122      M1 = P
1123      M2 = P
1124      RV5(P) = D(P)
1125      GO TO 900
1126  180 X1 = X1 * (Q - P + 1)
1127      LB = DMAX1(T1,XU-X1)
1128      UB = DMIN1(T2,X0+X1)
1129      X1 = LB
1130      ISTURM = 3
1131      GO TO 320
1132  200 M1 = S + 1
1133      X1 = UB
1134      ISTURM = 4
1135      GO TO 320
1136  220 M2 = S
1137      IF (M1 .GT. M2) GO TO 940
1138C     .......... FIND ROOTS BY BISECTION ..........
1139      X0 = UB
1140      ISTURM = 5
1141C
1142      DO 240 I = M1, M2
1143         RV5(I) = UB
1144         RV4(I) = LB
1145  240 CONTINUE
1146C     .......... LOOP FOR K-TH EIGENVALUE
1147C                FOR K=M2 STEP -1 UNTIL M1 DO --
1148C                (-DO- NOT USED TO LEGALIZE -COMPUTED GO TO-) ..........
1149      K = M2
1150  250    XU = LB
1151C     .......... FOR I=K STEP -1 UNTIL M1 DO -- ..........
1152         DO 260 II = M1, K
1153            I = M1 + K - II
1154            IF (XU .GE. RV4(I)) GO TO 260
1155            XU = RV4(I)
1156            GO TO 280
1157  260    CONTINUE
1158C
1159  280    IF (X0 .GT. RV5(K)) X0 = RV5(K)
1160C     .......... NEXT BISECTION STEP ..........
1161  300    X1 = (XU + X0) * 0.5D0
1162         IF ((X0 - XU) .LE. DABS(EPS1)) GO TO 420
1163         TST1 = 2.0D0 * (DABS(XU) + DABS(X0))
1164         TST2 = TST1 + (X0 - XU)
1165         IF (TST2 .EQ. TST1) GO TO 420
1166C     .......... IN-LINE PROCEDURE FOR STURM SEQUENCE ..........
1167  320    S = P - 1
1168         U = 1.0D0
1169C
1170         DO 340 I = P, Q
1171            IF (U .NE. 0.0D0) GO TO 325
1172            V = DABS(E(I)) / EPSLON(1.0D0)
1173            IF (E2(I) .EQ. 0.0D0) V = 0.0D0
1174            GO TO 330
1175  325       V = E2(I) / U
1176  330       U = D(I) - X1 - V
1177            IF (U .LT. 0.0D0) S = S + 1
1178  340    CONTINUE
1179C
1180         GO TO (60,80,200,220,360), ISTURM
1181C     .......... REFINE INTERVALS ..........
1182  360    IF (S .GE. K) GO TO 400
1183         XU = X1
1184         IF (S .GE. M1) GO TO 380
1185         RV4(M1) = X1
1186         GO TO 300
1187  380    RV4(S+1) = X1
1188         IF (RV5(S) .GT. X1) RV5(S) = X1
1189         GO TO 300
1190  400    X0 = X1
1191         GO TO 300
1192C     .......... K-TH EIGENVALUE FOUND ..........
1193  420    RV5(K) = X1
1194      K = K - 1
1195      IF (K .GE. M1) GO TO 250
1196C     .......... ORDER EIGENVALUES TAGGED WITH THEIR
1197C                SUBMATRIX ASSOCIATIONS ..........
1198  900 S = R
1199      R = R + M2 - M1 + 1
1200      J = 1
1201      K = M1
1202C
1203      DO 920 L = 1, R
1204         IF (J .GT. S) GO TO 910
1205         IF (K .GT. M2) GO TO 940
1206         IF (RV5(K) .GE. W(L)) GO TO 915
1207C
1208         DO 905 II = J, S
1209            I = L + S - II
1210            W(I+1) = W(I)
1211            IND(I+1) = IND(I)
1212  905    CONTINUE
1213C
1214  910    W(L) = RV5(K)
1215         IND(L) = TAG
1216         K = K + 1
1217         GO TO 920
1218  915    J = J + 1
1219  920 CONTINUE
1220C
1221  940 IF (Q .LT. N) GO TO 100
1222      GO TO 1001
1223C     .......... SET ERROR -- UNDERESTIMATE OF NUMBER OF
1224C                EIGENVALUES IN INTERVAL ..........
1225  980 IERR = 3 * N + 1
1226 1001 LB = T1
1227      UB = T2
1228      RETURN
1229      END
1230      SUBROUTINE BQR(NM,N,MB,A,T,R,IERR,NV,RV)
1231C
1232      INTEGER I,J,K,L,M,N,II,IK,JK,JM,KJ,KK,KM,LL,MB,MK,MN,MZ,
1233     X        M1,M2,M3,M4,NI,NM,NV,ITS,KJ1,M21,M31,IERR,IMULT
1234      DOUBLE PRECISION A(NM,MB),RV(NV)
1235      DOUBLE PRECISION F,G,Q,R,S,T,TST1,TST2,SCALE,PYTHAG
1236C
1237C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BQR,
1238C     NUM. MATH. 16, 85-92(1970) BY MARTIN, REINSCH, AND WILKINSON.
1239C     HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 266-272(1971).
1240C
1241C     THIS SUBROUTINE FINDS THE EIGENVALUE OF SMALLEST (USUALLY)
1242C     MAGNITUDE OF A REAL SYMMETRIC BAND MATRIX USING THE
1243C     QR ALGORITHM WITH SHIFTS OF ORIGIN.  CONSECUTIVE CALLS
1244C     CAN BE MADE TO FIND FURTHER EIGENVALUES.
1245C
1246C     ON INPUT
1247C
1248C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
1249C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
1250C          DIMENSION STATEMENT.
1251C
1252C        N IS THE ORDER OF THE MATRIX.
1253C
1254C        MB IS THE (HALF) BAND WIDTH OF THE MATRIX, DEFINED AS THE
1255C          NUMBER OF ADJACENT DIAGONALS, INCLUDING THE PRINCIPAL
1256C          DIAGONAL, REQUIRED TO SPECIFY THE NON-ZERO PORTION OF THE
1257C          LOWER TRIANGLE OF THE MATRIX.
1258C
1259C        A CONTAINS THE LOWER TRIANGLE OF THE SYMMETRIC BAND INPUT
1260C          MATRIX STORED AS AN N BY MB ARRAY.  ITS LOWEST SUBDIAGONAL
1261C          IS STORED IN THE LAST N+1-MB POSITIONS OF THE FIRST COLUMN,
1262C          ITS NEXT SUBDIAGONAL IN THE LAST N+2-MB POSITIONS OF THE
1263C          SECOND COLUMN, FURTHER SUBDIAGONALS SIMILARLY, AND FINALLY
1264C          ITS PRINCIPAL DIAGONAL IN THE N POSITIONS OF THE LAST COLUMN.
1265C          CONTENTS OF STORAGES NOT PART OF THE MATRIX ARE ARBITRARY.
1266C          ON A SUBSEQUENT CALL, ITS OUTPUT CONTENTS FROM THE PREVIOUS
1267C          CALL SHOULD BE PASSED.
1268C
1269C        T SPECIFIES THE SHIFT (OF EIGENVALUES) APPLIED TO THE DIAGONAL
1270C          OF A IN FORMING THE INPUT MATRIX. WHAT IS ACTUALLY DETERMINED
1271C          IS THE EIGENVALUE OF A+TI (I IS THE IDENTITY MATRIX) NEAREST
1272C          TO T.  ON A SUBSEQUENT CALL, THE OUTPUT VALUE OF T FROM THE
1273C          PREVIOUS CALL SHOULD BE PASSED IF THE NEXT NEAREST EIGENVALUE
1274C          IS SOUGHT.
1275C
1276C        R SHOULD BE SPECIFIED AS ZERO ON THE FIRST CALL, AND AS ITS
1277C          OUTPUT VALUE FROM THE PREVIOUS CALL ON A SUBSEQUENT CALL.
1278C          IT IS USED TO DETERMINE WHEN THE LAST ROW AND COLUMN OF
1279C          THE TRANSFORMED BAND MATRIX CAN BE REGARDED AS NEGLIGIBLE.
1280C
1281C        NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER RV
1282C          AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT.
1283C
1284C     ON OUTPUT
1285C
1286C        A CONTAINS THE TRANSFORMED BAND MATRIX.  THE MATRIX A+TI
1287C          DERIVED FROM THE OUTPUT PARAMETERS IS SIMILAR TO THE
1288C          INPUT A+TI TO WITHIN ROUNDING ERRORS.  ITS LAST ROW AND
1289C          COLUMN ARE NULL (IF IERR IS ZERO).
1290C
1291C        T CONTAINS THE COMPUTED EIGENVALUE OF A+TI (IF IERR IS ZERO).
1292C
1293C        R CONTAINS THE MAXIMUM OF ITS INPUT VALUE AND THE NORM OF THE
1294C          LAST COLUMN OF THE INPUT MATRIX A.
1295C
1296C        IERR IS SET TO
1297C          ZERO       FOR NORMAL RETURN,
1298C          N          IF THE EIGENVALUE HAS NOT BEEN
1299C                     DETERMINED AFTER 30 ITERATIONS.
1300C
1301C        RV IS A TEMPORARY STORAGE ARRAY OF DIMENSION AT LEAST
1302C          (2*MB**2+4*MB-3).  THE FIRST (3*MB-2) LOCATIONS CORRESPOND
1303C          TO THE ALGOL ARRAY B, THE NEXT (2*MB-1) LOCATIONS CORRESPOND
1304C          TO THE ALGOL ARRAY H, AND THE FINAL (2*MB**2-MB) LOCATIONS
1305C          CORRESPOND TO THE MB BY (2*MB-1) ALGOL ARRAY U.
1306C
1307C     NOTE. FOR A SUBSEQUENT CALL, N SHOULD BE REPLACED BY N-1, BUT
1308C     MB SHOULD NOT BE ALTERED EVEN WHEN IT EXCEEDS THE CURRENT N.
1309C
1310C     CALLS PYTHAG FOR  DSQRT(A*A + B*B) .
1311C
1312C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
1313C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
1314C
1315C     THIS VERSION DATED AUGUST 1983.
1316C
1317C     ------------------------------------------------------------------
1318C
1319      IERR = 0
1320      M1 = MIN0(MB,N)
1321      M = M1 - 1
1322      M2 = M + M
1323      M21 = M2 + 1
1324      M3 = M21 + M
1325      M31 = M3 + 1
1326      M4 = M31 + M2
1327      MN = M + N
1328      MZ = MB - M1
1329      ITS = 0
1330C     .......... TEST FOR CONVERGENCE ..........
1331   40 G = A(N,MB)
1332      IF (M .EQ. 0) GO TO 360
1333      F = 0.0D0
1334C
1335      DO 50 K = 1, M
1336         MK = K + MZ
1337         F = F + DABS(A(N,MK))
1338   50 CONTINUE
1339C
1340      IF (ITS .EQ. 0 .AND. F .GT. R) R = F
1341      TST1 = R
1342      TST2 = TST1 + F
1343      IF (TST2 .LE. TST1) GO TO 360
1344      IF (ITS .EQ. 30) GO TO 1000
1345      ITS = ITS + 1
1346C     .......... FORM SHIFT FROM BOTTOM 2 BY 2 MINOR ..........
1347      IF (F .GT. 0.25D0 * R .AND. ITS .LT. 5) GO TO 90
1348      F = A(N,MB-1)
1349      IF (F .EQ. 0.0D0) GO TO 70
1350      Q = (A(N-1,MB) - G) / (2.0D0 * F)
1351      S = PYTHAG(Q,1.0D0)
1352      G = G - F / (Q + DSIGN(S,Q))
1353   70 T = T + G
1354C
1355      DO 80 I = 1, N
1356   80 A(I,MB) = A(I,MB) - G
1357C
1358   90 DO 100 K = M31, M4
1359  100 RV(K) = 0.0D0
1360C
1361      DO 350 II = 1, MN
1362         I = II - M
1363         NI = N - II
1364         IF (NI .LT. 0) GO TO 230
1365C     .......... FORM COLUMN OF SHIFTED MATRIX A-G*I ..........
1366         L = MAX0(1,2-I)
1367C
1368         DO 110 K = 1, M3
1369  110    RV(K) = 0.0D0
1370C
1371         DO 120 K = L, M1
1372            KM = K + M
1373            MK = K + MZ
1374            RV(KM) = A(II,MK)
1375  120    CONTINUE
1376C
1377         LL = MIN0(M,NI)
1378         IF (LL .EQ. 0) GO TO 135
1379C
1380         DO 130 K = 1, LL
1381            KM = K + M21
1382            IK = II + K
1383            MK = MB - K
1384            RV(KM) = A(IK,MK)
1385  130    CONTINUE
1386C     .......... PRE-MULTIPLY WITH HOUSEHOLDER REFLECTIONS ..........
1387  135    LL = M2
1388         IMULT = 0
1389C     .......... MULTIPLICATION PROCEDURE ..........
1390  140    KJ = M4 - M1
1391C
1392         DO 170 J = 1, LL
1393            KJ = KJ + M1
1394            JM = J + M3
1395            IF (RV(JM) .EQ. 0.0D0) GO TO 170
1396            F = 0.0D0
1397C
1398            DO 150 K = 1, M1
1399               KJ = KJ + 1
1400               JK = J + K - 1
1401               F = F + RV(KJ) * RV(JK)
1402  150       CONTINUE
1403C
1404            F = F / RV(JM)
1405            KJ = KJ - M1
1406C
1407            DO 160 K = 1, M1
1408               KJ = KJ + 1
1409               JK = J + K - 1
1410               RV(JK) = RV(JK) - RV(KJ) * F
1411  160       CONTINUE
1412C
1413            KJ = KJ - M1
1414  170    CONTINUE
1415C
1416         IF (IMULT .NE. 0) GO TO 280
1417C     .......... HOUSEHOLDER REFLECTION ..........
1418         F = RV(M21)
1419         S = 0.0D0
1420         RV(M4) = 0.0D0
1421         SCALE = 0.0D0
1422C
1423         DO 180 K = M21, M3
1424  180    SCALE = SCALE + DABS(RV(K))
1425C
1426         IF (SCALE .EQ. 0.0D0) GO TO 210
1427C
1428         DO 190 K = M21, M3
1429  190    S = S + (RV(K)/SCALE)**2
1430C
1431         S = SCALE * SCALE * S
1432         G = -DSIGN(DSQRT(S),F)
1433         RV(M21) = G
1434         RV(M4) = S - F * G
1435         KJ = M4 + M2 * M1 + 1
1436         RV(KJ) = F - G
1437C
1438         DO 200 K = 2, M1
1439            KJ = KJ + 1
1440            KM = K + M2
1441            RV(KJ) = RV(KM)
1442  200    CONTINUE
1443C     .......... SAVE COLUMN OF TRIANGULAR FACTOR R ..........
1444  210    DO 220 K = L, M1
1445            KM = K + M
1446            MK = K + MZ
1447            A(II,MK) = RV(KM)
1448  220    CONTINUE
1449C
1450  230    L = MAX0(1,M1+1-I)
1451         IF (I .LE. 0) GO TO 300
1452C     .......... PERFORM ADDITIONAL STEPS ..........
1453         DO 240 K = 1, M21
1454  240    RV(K) = 0.0D0
1455C
1456         LL = MIN0(M1,NI+M1)
1457C     .......... GET ROW OF TRIANGULAR FACTOR R ..........
1458         DO 250 KK = 1, LL
1459            K = KK - 1
1460            KM = K + M1
1461            IK = I + K
1462            MK = MB - K
1463            RV(KM) = A(IK,MK)
1464  250    CONTINUE
1465C     .......... POST-MULTIPLY WITH HOUSEHOLDER REFLECTIONS ..........
1466         LL = M1
1467         IMULT = 1
1468         GO TO 140
1469C     .......... STORE COLUMN OF NEW A MATRIX ..........
1470  280    DO 290 K = L, M1
1471            MK = K + MZ
1472            A(I,MK) = RV(K)
1473  290    CONTINUE
1474C     .......... UPDATE HOUSEHOLDER REFLECTIONS ..........
1475  300    IF (L .GT. 1) L = L - 1
1476         KJ1 = M4 + L * M1
1477C
1478         DO 320 J = L, M2
1479            JM = J + M3
1480            RV(JM) = RV(JM+1)
1481C
1482            DO 320 K = 1, M1
1483               KJ1 = KJ1 + 1
1484               KJ = KJ1 - M1
1485               RV(KJ) = RV(KJ1)
1486  320    CONTINUE
1487C
1488  350 CONTINUE
1489C
1490      GO TO 40
1491C     .......... CONVERGENCE ..........
1492  360 T = T + G
1493C
1494      DO 380 I = 1, N
1495  380 A(I,MB) = A(I,MB) - G
1496C
1497      DO 400 K = 1, M1
1498         MK = K + MZ
1499         A(N,MK) = 0.0D0
1500  400 CONTINUE
1501C
1502      GO TO 1001
1503C     .......... SET ERROR -- NO CONVERGENCE TO
1504C                EIGENVALUE AFTER 30 ITERATIONS ..........
1505 1000 IERR = N
1506 1001 RETURN
1507      END
1508      SUBROUTINE CBABK2(NM,N,LOW,IGH,SCALE,M,ZR,ZI)
1509C
1510      INTEGER I,J,K,M,N,II,NM,IGH,LOW
1511      DOUBLE PRECISION SCALE(N),ZR(NM,M),ZI(NM,M)
1512      DOUBLE PRECISION S
1513C
1514C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE
1515C     CBABK2, WHICH IS A COMPLEX VERSION OF BALBAK,
1516C     NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH.
1517C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971).
1518C
1519C     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX GENERAL
1520C     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
1521C     BALANCED MATRIX DETERMINED BY  CBAL.
1522C
1523C     ON INPUT
1524C
1525C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
1526C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
1527C          DIMENSION STATEMENT.
1528C
1529C        N IS THE ORDER OF THE MATRIX.
1530C
1531C        LOW AND IGH ARE INTEGERS DETERMINED BY  CBAL.
1532C
1533C        SCALE CONTAINS INFORMATION DETERMINING THE PERMUTATIONS
1534C          AND SCALING FACTORS USED BY  CBAL.
1535C
1536C        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.
1537C
1538C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
1539C          RESPECTIVELY, OF THE EIGENVECTORS TO BE
1540C          BACK TRANSFORMED IN THEIR FIRST M COLUMNS.
1541C
1542C     ON OUTPUT
1543C
1544C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
1545C          RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS
1546C          IN THEIR FIRST M COLUMNS.
1547C
1548C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
1549C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
1550C
1551C     THIS VERSION DATED AUGUST 1983.
1552C
1553C     ------------------------------------------------------------------
1554C
1555      IF (M .EQ. 0) GO TO 200
1556      IF (IGH .EQ. LOW) GO TO 120
1557C
1558      DO 110 I = LOW, IGH
1559         S = SCALE(I)
1560C     .......... LEFT HAND EIGENVECTORS ARE BACK TRANSFORMED
1561C                IF THE FOREGOING STATEMENT IS REPLACED BY
1562C                S=1.0D0/SCALE(I). ..........
1563         DO 100 J = 1, M
1564            ZR(I,J) = ZR(I,J) * S
1565            ZI(I,J) = ZI(I,J) * S
1566  100    CONTINUE
1567C
1568  110 CONTINUE
1569C     .......... FOR I=LOW-1 STEP -1 UNTIL 1,
1570C                IGH+1 STEP 1 UNTIL N DO -- ..........
1571  120 DO 140 II = 1, N
1572         I = II
1573         IF (I .GE. LOW .AND. I .LE. IGH) GO TO 140
1574         IF (I .LT. LOW) I = LOW - II
1575         K = SCALE(I)
1576         IF (K .EQ. I) GO TO 140
1577C
1578         DO 130 J = 1, M
1579            S = ZR(I,J)
1580            ZR(I,J) = ZR(K,J)
1581            ZR(K,J) = S
1582            S = ZI(I,J)
1583            ZI(I,J) = ZI(K,J)
1584            ZI(K,J) = S
1585  130    CONTINUE
1586C
1587  140 CONTINUE
1588C
1589  200 RETURN
1590      END
1591      SUBROUTINE CBAL(NM,N,AR,AI,LOW,IGH,SCALE)
1592C
1593      INTEGER I,J,K,L,M,N,JJ,NM,IGH,LOW,IEXC
1594      DOUBLE PRECISION AR(NM,N),AI(NM,N),SCALE(N)
1595      DOUBLE PRECISION C,F,G,R,S,B2,RADIX
1596      LOGICAL NOCONV
1597C
1598C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE
1599C     CBALANCE, WHICH IS A COMPLEX VERSION OF BALANCE,
1600C     NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH.
1601C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971).
1602C
1603C     THIS SUBROUTINE BALANCES A COMPLEX MATRIX AND ISOLATES
1604C     EIGENVALUES WHENEVER POSSIBLE.
1605C
1606C     ON INPUT
1607C
1608C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
1609C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
1610C          DIMENSION STATEMENT.
1611C
1612C        N IS THE ORDER OF THE MATRIX.
1613C
1614C        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
1615C          RESPECTIVELY, OF THE COMPLEX MATRIX TO BE BALANCED.
1616C
1617C     ON OUTPUT
1618C
1619C        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
1620C          RESPECTIVELY, OF THE BALANCED MATRIX.
1621C
1622C        LOW AND IGH ARE TWO INTEGERS SUCH THAT AR(I,J) AND AI(I,J)
1623C          ARE EQUAL TO ZERO IF
1624C           (1) I IS GREATER THAN J AND
1625C           (2) J=1,...,LOW-1 OR I=IGH+1,...,N.
1626C
1627C        SCALE CONTAINS INFORMATION DETERMINING THE
1628C           PERMUTATIONS AND SCALING FACTORS USED.
1629C
1630C     SUPPOSE THAT THE PRINCIPAL SUBMATRIX IN ROWS LOW THROUGH IGH
1631C     HAS BEEN BALANCED, THAT P(J) DENOTES THE INDEX INTERCHANGED
1632C     WITH J DURING THE PERMUTATION STEP, AND THAT THE ELEMENTS
1633C     OF THE DIAGONAL MATRIX USED ARE DENOTED BY D(I,J).  THEN
1634C        SCALE(J) = P(J),    FOR J = 1,...,LOW-1
1635C                 = D(J,J)       J = LOW,...,IGH
1636C                 = P(J)         J = IGH+1,...,N.
1637C     THE ORDER IN WHICH THE INTERCHANGES ARE MADE IS N TO IGH+1,
1638C     THEN 1 TO LOW-1.
1639C
1640C     NOTE THAT 1 IS RETURNED FOR IGH IF IGH IS ZERO FORMALLY.
1641C
1642C     THE ALGOL PROCEDURE EXC CONTAINED IN CBALANCE APPEARS IN
1643C     CBAL  IN LINE.  (NOTE THAT THE ALGOL ROLES OF IDENTIFIERS
1644C     K,L HAVE BEEN REVERSED.)
1645C
1646C     ARITHMETIC IS REAL THROUGHOUT.
1647C
1648C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
1649C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
1650C
1651C     THIS VERSION DATED AUGUST 1983.
1652C
1653C     ------------------------------------------------------------------
1654C
1655      RADIX = 16.0D0
1656C
1657      B2 = RADIX * RADIX
1658      K = 1
1659      L = N
1660      GO TO 100
1661C     .......... IN-LINE PROCEDURE FOR ROW AND
1662C                COLUMN EXCHANGE ..........
1663   20 SCALE(M) = J
1664      IF (J .EQ. M) GO TO 50
1665C
1666      DO 30 I = 1, L
1667         F = AR(I,J)
1668         AR(I,J) = AR(I,M)
1669         AR(I,M) = F
1670         F = AI(I,J)
1671         AI(I,J) = AI(I,M)
1672         AI(I,M) = F
1673   30 CONTINUE
1674C
1675      DO 40 I = K, N
1676         F = AR(J,I)
1677         AR(J,I) = AR(M,I)
1678         AR(M,I) = F
1679         F = AI(J,I)
1680         AI(J,I) = AI(M,I)
1681         AI(M,I) = F
1682   40 CONTINUE
1683C
1684   50 GO TO (80,130), IEXC
1685C     .......... SEARCH FOR ROWS ISOLATING AN EIGENVALUE
1686C                AND PUSH THEM DOWN ..........
1687   80 IF (L .EQ. 1) GO TO 280
1688      L = L - 1
1689C     .......... FOR J=L STEP -1 UNTIL 1 DO -- ..........
1690  100 DO 120 JJ = 1, L
1691         J = L + 1 - JJ
1692C
1693         DO 110 I = 1, L
1694            IF (I .EQ. J) GO TO 110
1695            IF (AR(J,I) .NE. 0.0D0 .OR. AI(J,I) .NE. 0.0D0) GO TO 120
1696  110    CONTINUE
1697C
1698         M = L
1699         IEXC = 1
1700         GO TO 20
1701  120 CONTINUE
1702C
1703      GO TO 140
1704C     .......... SEARCH FOR COLUMNS ISOLATING AN EIGENVALUE
1705C                AND PUSH THEM LEFT ..........
1706  130 K = K + 1
1707C
1708  140 DO 170 J = K, L
1709C
1710         DO 150 I = K, L
1711            IF (I .EQ. J) GO TO 150
1712            IF (AR(I,J) .NE. 0.0D0 .OR. AI(I,J) .NE. 0.0D0) GO TO 170
1713  150    CONTINUE
1714C
1715         M = K
1716         IEXC = 2
1717         GO TO 20
1718  170 CONTINUE
1719C     .......... NOW BALANCE THE SUBMATRIX IN ROWS K TO L ..........
1720      DO 180 I = K, L
1721  180 SCALE(I) = 1.0D0
1722C     .......... ITERATIVE LOOP FOR NORM REDUCTION ..........
1723  190 NOCONV = .FALSE.
1724C
1725      DO 270 I = K, L
1726         C = 0.0D0
1727         R = 0.0D0
1728C
1729         DO 200 J = K, L
1730            IF (J .EQ. I) GO TO 200
1731            C = C + DABS(AR(J,I)) + DABS(AI(J,I))
1732            R = R + DABS(AR(I,J)) + DABS(AI(I,J))
1733  200    CONTINUE
1734C     .......... GUARD AGAINST ZERO C OR R DUE TO UNDERFLOW ..........
1735         IF (C .EQ. 0.0D0 .OR. R .EQ. 0.0D0) GO TO 270
1736         G = R / RADIX
1737         F = 1.0D0
1738         S = C + R
1739  210    IF (C .GE. G) GO TO 220
1740         F = F * RADIX
1741         C = C * B2
1742         GO TO 210
1743  220    G = R * RADIX
1744  230    IF (C .LT. G) GO TO 240
1745         F = F / RADIX
1746         C = C / B2
1747         GO TO 230
1748C     .......... NOW BALANCE ..........
1749  240    IF ((C + R) / F .GE. 0.95D0 * S) GO TO 270
1750         G = 1.0D0 / F
1751         SCALE(I) = SCALE(I) * F
1752         NOCONV = .TRUE.
1753C
1754         DO 250 J = K, N
1755            AR(I,J) = AR(I,J) * G
1756            AI(I,J) = AI(I,J) * G
1757  250    CONTINUE
1758C
1759         DO 260 J = 1, L
1760            AR(J,I) = AR(J,I) * F
1761            AI(J,I) = AI(J,I) * F
1762  260    CONTINUE
1763C
1764  270 CONTINUE
1765C
1766      IF (NOCONV) GO TO 190
1767C
1768  280 LOW = K
1769      IGH = L
1770      RETURN
1771      END
1772      SUBROUTINE CG(NM,N,AR,AI,WR,WI,MATZ,ZR,ZI,FV1,FV2,FV3,IERR)
1773C
1774      INTEGER N,NM,IS1,IS2,IERR,MATZ
1775      DOUBLE PRECISION AR(NM,N),AI(NM,N),WR(N),WI(N),ZR(NM,N),ZI(NM,N),
1776     X       FV1(N),FV2(N),FV3(N)
1777C
1778C     THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF
1779C     SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK)
1780C     TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED)
1781C     OF A COMPLEX GENERAL MATRIX.
1782C
1783C     ON INPUT
1784C
1785C        NM  MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
1786C        ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
1787C        DIMENSION STATEMENT.
1788C
1789C        N  IS THE ORDER OF THE MATRIX  A=(AR,AI).
1790C
1791C        AR  AND  AI  CONTAIN THE REAL AND IMAGINARY PARTS,
1792C        RESPECTIVELY, OF THE COMPLEX GENERAL MATRIX.
1793C
1794C        MATZ  IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF
1795C        ONLY EIGENVALUES ARE DESIRED.  OTHERWISE IT IS SET TO
1796C        ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS.
1797C
1798C     ON OUTPUT
1799C
1800C        WR  AND  WI  CONTAIN THE REAL AND IMAGINARY PARTS,
1801C        RESPECTIVELY, OF THE EIGENVALUES.
1802C
1803C        ZR  AND  ZI  CONTAIN THE REAL AND IMAGINARY PARTS,
1804C        RESPECTIVELY, OF THE EIGENVECTORS IF MATZ IS NOT ZERO.
1805C
1806C        IERR  IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR
1807C           COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR COMQR
1808C           AND COMQR2.  THE NORMAL COMPLETION CODE IS ZERO.
1809C
1810C        FV1, FV2, AND  FV3  ARE TEMPORARY STORAGE ARRAYS.
1811C
1812C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
1813C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
1814C
1815C     THIS VERSION DATED AUGUST 1983.
1816C
1817C     ------------------------------------------------------------------
1818C
1819      IF (N .LE. NM) GO TO 10
1820      IERR = 10 * N
1821      GO TO 50
1822C
1823   10 CALL  CBAL(NM,N,AR,AI,IS1,IS2,FV1)
1824      CALL  CORTH(NM,N,IS1,IS2,AR,AI,FV2,FV3)
1825      IF (MATZ .NE. 0) GO TO 20
1826C     .......... FIND EIGENVALUES ONLY ..........
1827      CALL  COMQR(NM,N,IS1,IS2,AR,AI,WR,WI,IERR)
1828      GO TO 50
1829C     .......... FIND BOTH EIGENVALUES AND EIGENVECTORS ..........
1830   20 CALL  COMQR2(NM,N,IS1,IS2,FV2,FV3,AR,AI,WR,WI,ZR,ZI,IERR)
1831      IF (IERR .NE. 0) GO TO 50
1832      CALL  CBABK2(NM,N,IS1,IS2,FV1,N,ZR,ZI)
1833   50 RETURN
1834      END
1835      SUBROUTINE CH(NM,N,AR,AI,W,MATZ,ZR,ZI,FV1,FV2,FM1,IERR)
1836C
1837      INTEGER I,J,N,NM,IERR,MATZ
1838      DOUBLE PRECISION AR(NM,N),AI(NM,N),W(N),ZR(NM,N),ZI(NM,N),
1839     X       FV1(N),FV2(N),FM1(2,N)
1840C
1841C     THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF
1842C     SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK)
1843C     TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED)
1844C     OF A COMPLEX HERMITIAN MATRIX.
1845C
1846C     ON INPUT
1847C
1848C        NM  MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
1849C        ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
1850C        DIMENSION STATEMENT.
1851C
1852C        N  IS THE ORDER OF THE MATRIX  A=(AR,AI).
1853C
1854C        AR  AND  AI  CONTAIN THE REAL AND IMAGINARY PARTS,
1855C        RESPECTIVELY, OF THE COMPLEX HERMITIAN MATRIX.
1856C
1857C        MATZ  IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF
1858C        ONLY EIGENVALUES ARE DESIRED.  OTHERWISE IT IS SET TO
1859C        ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS.
1860C
1861C     ON OUTPUT
1862C
1863C        W  CONTAINS THE EIGENVALUES IN ASCENDING ORDER.
1864C
1865C        ZR  AND  ZI  CONTAIN THE REAL AND IMAGINARY PARTS,
1866C        RESPECTIVELY, OF THE EIGENVECTORS IF MATZ IS NOT ZERO.
1867C
1868C        IERR  IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR
1869C           COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR TQLRAT
1870C           AND TQL2.  THE NORMAL COMPLETION CODE IS ZERO.
1871C
1872C        FV1, FV2, AND  FM1  ARE TEMPORARY STORAGE ARRAYS.
1873C
1874C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
1875C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
1876C
1877C     THIS VERSION DATED AUGUST 1983.
1878C
1879C     ------------------------------------------------------------------
1880C
1881      IF (N .LE. NM) GO TO 10
1882      IERR = 10 * N
1883      GO TO 50
1884C
1885   10 CALL  HTRIDI(NM,N,AR,AI,W,FV1,FV2,FM1)
1886      IF (MATZ .NE. 0) GO TO 20
1887C     .......... FIND EIGENVALUES ONLY ..........
1888      CALL  TQLRAT(N,W,FV2,IERR)
1889      GO TO 50
1890C     .......... FIND BOTH EIGENVALUES AND EIGENVECTORS ..........
1891   20 DO 40 I = 1, N
1892C
1893         DO 30 J = 1, N
1894            ZR(J,I) = 0.0D0
1895   30    CONTINUE
1896C
1897         ZR(I,I) = 1.0D0
1898   40 CONTINUE
1899C
1900      CALL  TQL2(NM,N,W,FV1,ZR,IERR)
1901      IF (IERR .NE. 0) GO TO 50
1902      CALL  HTRIBK(NM,N,AR,AI,FM1,N,ZR,ZI)
1903   50 RETURN
1904      END
1905      SUBROUTINE CINVIT(NM,N,AR,AI,WR,WI,SELECT,MM,M,ZR,ZI,
1906     X                  IERR,RM1,RM2,RV1,RV2)
1907C
1908      INTEGER I,J,K,M,N,S,II,MM,MP,NM,UK,IP1,ITS,KM1,IERR
1909      DOUBLE PRECISION AR(NM,N),AI(NM,N),WR(N),WI(N),ZR(NM,MM),
1910     X       ZI(NM,MM),RM1(N,N),RM2(N,N),RV1(N),RV2(N)
1911      DOUBLE PRECISION X,Y,EPS3,NORM,NORMV,EPSLON,GROWTO,ILAMBD,PYTHAG,
1912     X       RLAMBD,UKROOT
1913      LOGICAL SELECT(N)
1914C
1915C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE CX INVIT
1916C     BY PETERS AND WILKINSON.
1917C     HANDBOOK FOR AUTO. COMP. VOL.II-LINEAR ALGEBRA, 418-439(1971).
1918C
1919C     THIS SUBROUTINE FINDS THOSE EIGENVECTORS OF A COMPLEX UPPER
1920C     HESSENBERG MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES,
1921C     USING INVERSE ITERATION.
1922C
1923C     ON INPUT
1924C
1925C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
1926C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
1927C          DIMENSION STATEMENT.
1928C
1929C        N IS THE ORDER OF THE MATRIX.
1930C
1931C        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
1932C          RESPECTIVELY, OF THE HESSENBERG MATRIX.
1933C
1934C        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS, RESPECTIVELY,
1935C          OF THE EIGENVALUES OF THE MATRIX.  THE EIGENVALUES MUST BE
1936C          STORED IN A MANNER IDENTICAL TO THAT OF SUBROUTINE  COMLR,
1937C          WHICH RECOGNIZES POSSIBLE SPLITTING OF THE MATRIX.
1938C
1939C        SELECT SPECIFIES THE EIGENVECTORS TO BE FOUND.  THE
1940C          EIGENVECTOR CORRESPONDING TO THE J-TH EIGENVALUE IS
1941C          SPECIFIED BY SETTING SELECT(J) TO .TRUE..
1942C
1943C        MM SHOULD BE SET TO AN UPPER BOUND FOR THE NUMBER OF
1944C          EIGENVECTORS TO BE FOUND.
1945C
1946C     ON OUTPUT
1947C
1948C        AR, AI, WI, AND SELECT ARE UNALTERED.
1949C
1950C        WR MAY HAVE BEEN ALTERED SINCE CLOSE EIGENVALUES ARE PERTURBED
1951C          SLIGHTLY IN SEARCHING FOR INDEPENDENT EIGENVECTORS.
1952C
1953C        M IS THE NUMBER OF EIGENVECTORS ACTUALLY FOUND.
1954C
1955C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS, RESPECTIVELY,
1956C          OF THE EIGENVECTORS.  THE EIGENVECTORS ARE NORMALIZED
1957C          SO THAT THE COMPONENT OF LARGEST MAGNITUDE IS 1.
1958C          ANY VECTOR WHICH FAILS THE ACCEPTANCE TEST IS SET TO ZERO.
1959C
1960C        IERR IS SET TO
1961C          ZERO       FOR NORMAL RETURN,
1962C          -(2*N+1)   IF MORE THAN MM EIGENVECTORS HAVE BEEN SPECIFIED,
1963C          -K         IF THE ITERATION CORRESPONDING TO THE K-TH
1964C                     VALUE FAILS,
1965C          -(N+K)     IF BOTH ERROR SITUATIONS OCCUR.
1966C
1967C        RM1, RM2, RV1, AND RV2 ARE TEMPORARY STORAGE ARRAYS.
1968C
1969C     THE ALGOL PROCEDURE GUESSVEC APPEARS IN CINVIT IN LINE.
1970C
1971C     CALLS CDIV FOR COMPLEX DIVISION.
1972C     CALLS PYTHAG FOR  DSQRT(A*A + B*B) .
1973C
1974C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
1975C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
1976C
1977C     THIS VERSION DATED AUGUST 1983.
1978C
1979C     ------------------------------------------------------------------
1980C
1981      IERR = 0
1982      UK = 0
1983      S = 1
1984C
1985      DO 980 K = 1, N
1986         IF (.NOT. SELECT(K)) GO TO 980
1987         IF (S .GT. MM) GO TO 1000
1988         IF (UK .GE. K) GO TO 200
1989C     .......... CHECK FOR POSSIBLE SPLITTING ..........
1990         DO 120 UK = K, N
1991            IF (UK .EQ. N) GO TO 140
1992            IF (AR(UK+1,UK) .EQ. 0.0D0 .AND. AI(UK+1,UK) .EQ. 0.0D0)
1993     X         GO TO 140
1994  120    CONTINUE
1995C     .......... COMPUTE INFINITY NORM OF LEADING UK BY UK
1996C                (HESSENBERG) MATRIX ..........
1997  140    NORM = 0.0D0
1998         MP = 1
1999C
2000         DO 180 I = 1, UK
2001            X = 0.0D0
2002C
2003            DO 160 J = MP, UK
2004  160       X = X + PYTHAG(AR(I,J),AI(I,J))
2005C
2006            IF (X .GT. NORM) NORM = X
2007            MP = I
2008  180    CONTINUE
2009C     .......... EPS3 REPLACES ZERO PIVOT IN DECOMPOSITION
2010C                AND CLOSE ROOTS ARE MODIFIED BY EPS3 ..........
2011         IF (NORM .EQ. 0.0D0) NORM = 1.0D0
2012         EPS3 = EPSLON(NORM)
2013C     .......... GROWTO IS THE CRITERION FOR GROWTH ..........
2014         UKROOT = UK
2015         UKROOT = DSQRT(UKROOT)
2016         GROWTO = 0.1D0 / UKROOT
2017  200    RLAMBD = WR(K)
2018         ILAMBD = WI(K)
2019         IF (K .EQ. 1) GO TO 280
2020         KM1 = K - 1
2021         GO TO 240
2022C     .......... PERTURB EIGENVALUE IF IT IS CLOSE
2023C                TO ANY PREVIOUS EIGENVALUE ..........
2024  220    RLAMBD = RLAMBD + EPS3
2025C     .......... FOR I=K-1 STEP -1 UNTIL 1 DO -- ..........
2026  240    DO 260 II = 1, KM1
2027            I = K - II
2028            IF (SELECT(I) .AND. DABS(WR(I)-RLAMBD) .LT. EPS3 .AND.
2029     X         DABS(WI(I)-ILAMBD) .LT. EPS3) GO TO 220
2030  260    CONTINUE
2031C
2032         WR(K) = RLAMBD
2033C     .......... FORM UPPER HESSENBERG (AR,AI)-(RLAMBD,ILAMBD)*I
2034C                AND INITIAL COMPLEX VECTOR ..........
2035  280    MP = 1
2036C
2037         DO 320 I = 1, UK
2038C
2039            DO 300 J = MP, UK
2040               RM1(I,J) = AR(I,J)
2041               RM2(I,J) = AI(I,J)
2042  300       CONTINUE
2043C
2044            RM1(I,I) = RM1(I,I) - RLAMBD
2045            RM2(I,I) = RM2(I,I) - ILAMBD
2046            MP = I
2047            RV1(I) = EPS3
2048  320    CONTINUE
2049C     .......... TRIANGULAR DECOMPOSITION WITH INTERCHANGES,
2050C                REPLACING ZERO PIVOTS BY EPS3 ..........
2051         IF (UK .EQ. 1) GO TO 420
2052C
2053         DO 400 I = 2, UK
2054            MP = I - 1
2055            IF (PYTHAG(RM1(I,MP),RM2(I,MP)) .LE.
2056     X          PYTHAG(RM1(MP,MP),RM2(MP,MP))) GO TO 360
2057C
2058            DO 340 J = MP, UK
2059               Y = RM1(I,J)
2060               RM1(I,J) = RM1(MP,J)
2061               RM1(MP,J) = Y
2062               Y = RM2(I,J)
2063               RM2(I,J) = RM2(MP,J)
2064               RM2(MP,J) = Y
2065  340       CONTINUE
2066C
2067  360       IF (RM1(MP,MP) .EQ. 0.0D0 .AND. RM2(MP,MP) .EQ. 0.0D0)
2068     X         RM1(MP,MP) = EPS3
2069            CALL CDIV(RM1(I,MP),RM2(I,MP),RM1(MP,MP),RM2(MP,MP),X,Y)
2070            IF (X .EQ. 0.0D0 .AND. Y .EQ. 0.0D0) GO TO 400
2071C
2072            DO 380 J = I, UK
2073               RM1(I,J) = RM1(I,J) - X * RM1(MP,J) + Y * RM2(MP,J)
2074               RM2(I,J) = RM2(I,J) - X * RM2(MP,J) - Y * RM1(MP,J)
2075  380       CONTINUE
2076C
2077  400    CONTINUE
2078C
2079  420    IF (RM1(UK,UK) .EQ. 0.0D0 .AND. RM2(UK,UK) .EQ. 0.0D0)
2080     X      RM1(UK,UK) = EPS3
2081         ITS = 0
2082C     .......... BACK SUBSTITUTION
2083C                FOR I=UK STEP -1 UNTIL 1 DO -- ..........
2084  660    DO 720 II = 1, UK
2085            I = UK + 1 - II
2086            X = RV1(I)
2087            Y = 0.0D0
2088            IF (I .EQ. UK) GO TO 700
2089            IP1 = I + 1
2090C
2091            DO 680 J = IP1, UK
2092               X = X - RM1(I,J) * RV1(J) + RM2(I,J) * RV2(J)
2093               Y = Y - RM1(I,J) * RV2(J) - RM2(I,J) * RV1(J)
2094  680       CONTINUE
2095C
2096  700       CALL CDIV(X,Y,RM1(I,I),RM2(I,I),RV1(I),RV2(I))
2097  720    CONTINUE
2098C     .......... ACCEPTANCE TEST FOR EIGENVECTOR
2099C                AND NORMALIZATION ..........
2100         ITS = ITS + 1
2101         NORM = 0.0D0
2102         NORMV = 0.0D0
2103C
2104         DO 780 I = 1, UK
2105            X = PYTHAG(RV1(I),RV2(I))
2106            IF (NORMV .GE. X) GO TO 760
2107            NORMV = X
2108            J = I
2109  760       NORM = NORM + X
2110  780    CONTINUE
2111C
2112         IF (NORM .LT. GROWTO) GO TO 840
2113C     .......... ACCEPT VECTOR ..........
2114         X = RV1(J)
2115         Y = RV2(J)
2116C
2117         DO 820 I = 1, UK
2118            CALL CDIV(RV1(I),RV2(I),X,Y,ZR(I,S),ZI(I,S))
2119  820    CONTINUE
2120C
2121         IF (UK .EQ. N) GO TO 940
2122         J = UK + 1
2123         GO TO 900
2124C     .......... IN-LINE PROCEDURE FOR CHOOSING
2125C                A NEW STARTING VECTOR ..........
2126  840    IF (ITS .GE. UK) GO TO 880
2127         X = UKROOT
2128         Y = EPS3 / (X + 1.0D0)
2129         RV1(1) = EPS3
2130C
2131         DO 860 I = 2, UK
2132  860    RV1(I) = Y
2133C
2134         J = UK - ITS + 1
2135         RV1(J) = RV1(J) - EPS3 * X
2136         GO TO 660
2137C     .......... SET ERROR -- UNACCEPTED EIGENVECTOR ..........
2138  880    J = 1
2139         IERR = -K
2140C     .......... SET REMAINING VECTOR COMPONENTS TO ZERO ..........
2141  900    DO 920 I = J, N
2142            ZR(I,S) = 0.0D0
2143            ZI(I,S) = 0.0D0
2144  920    CONTINUE
2145C
2146  940    S = S + 1
2147  980 CONTINUE
2148C
2149      GO TO 1001
2150C     .......... SET ERROR -- UNDERESTIMATE OF EIGENVECTOR
2151C                SPACE REQUIRED ..........
2152 1000 IF (IERR .NE. 0) IERR = IERR - N
2153      IF (IERR .EQ. 0) IERR = -(2 * N + 1)
2154 1001 M = S - 1
2155      RETURN
2156      END
2157      SUBROUTINE COMBAK(NM,LOW,IGH,AR,AI,INT,M,ZR,ZI)
2158C
2159      INTEGER I,J,M,LA,MM,MP,NM,IGH,KP1,LOW,MP1
2160      DOUBLE PRECISION AR(NM,IGH),AI(NM,IGH),ZR(NM,M),ZI(NM,M)
2161      DOUBLE PRECISION XR,XI
2162      INTEGER INT(IGH)
2163C
2164C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMBAK,
2165C     NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON.
2166C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
2167C
2168C     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX GENERAL
2169C     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
2170C     UPPER HESSENBERG MATRIX DETERMINED BY  COMHES.
2171C
2172C     ON INPUT
2173C
2174C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
2175C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
2176C          DIMENSION STATEMENT.
2177C
2178C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
2179C          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED,
2180C          SET LOW=1 AND IGH EQUAL TO THE ORDER OF THE MATRIX.
2181C
2182C        AR AND AI CONTAIN THE MULTIPLIERS WHICH WERE USED IN THE
2183C          REDUCTION BY  COMHES  IN THEIR LOWER TRIANGLES
2184C          BELOW THE SUBDIAGONAL.
2185C
2186C        INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS
2187C          INTERCHANGED IN THE REDUCTION BY  COMHES.
2188C          ONLY ELEMENTS LOW THROUGH IGH ARE USED.
2189C
2190C        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.
2191C
2192C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
2193C          RESPECTIVELY, OF THE EIGENVECTORS TO BE
2194C          BACK TRANSFORMED IN THEIR FIRST M COLUMNS.
2195C
2196C     ON OUTPUT
2197C
2198C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
2199C          RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS
2200C          IN THEIR FIRST M COLUMNS.
2201C
2202C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
2203C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
2204C
2205C     THIS VERSION DATED AUGUST 1983.
2206C
2207C     ------------------------------------------------------------------
2208C
2209      IF (M .EQ. 0) GO TO 200
2210      LA = IGH - 1
2211      KP1 = LOW + 1
2212      IF (LA .LT. KP1) GO TO 200
2213C     .......... FOR MP=IGH-1 STEP -1 UNTIL LOW+1 DO -- ..........
2214      DO 140 MM = KP1, LA
2215         MP = LOW + IGH - MM
2216         MP1 = MP + 1
2217C
2218         DO 110 I = MP1, IGH
2219            XR = AR(I,MP-1)
2220            XI = AI(I,MP-1)
2221            IF (XR .EQ. 0.0D0 .AND. XI .EQ. 0.0D0) GO TO 110
2222C
2223            DO 100 J = 1, M
2224               ZR(I,J) = ZR(I,J) + XR * ZR(MP,J) - XI * ZI(MP,J)
2225               ZI(I,J) = ZI(I,J) + XR * ZI(MP,J) + XI * ZR(MP,J)
2226  100       CONTINUE
2227C
2228  110    CONTINUE
2229C
2230         I = INT(MP)
2231         IF (I .EQ. MP) GO TO 140
2232C
2233         DO 130 J = 1, M
2234            XR = ZR(I,J)
2235            ZR(I,J) = ZR(MP,J)
2236            ZR(MP,J) = XR
2237            XI = ZI(I,J)
2238            ZI(I,J) = ZI(MP,J)
2239            ZI(MP,J) = XI
2240  130    CONTINUE
2241C
2242  140 CONTINUE
2243C
2244  200 RETURN
2245      END
2246      SUBROUTINE COMHES(NM,N,LOW,IGH,AR,AI,INT)
2247C
2248      INTEGER I,J,M,N,LA,NM,IGH,KP1,LOW,MM1,MP1
2249      DOUBLE PRECISION AR(NM,N),AI(NM,N)
2250      DOUBLE PRECISION XR,XI,YR,YI
2251      INTEGER INT(IGH)
2252C
2253C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMHES,
2254C     NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON.
2255C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
2256C
2257C     GIVEN A COMPLEX GENERAL MATRIX, THIS SUBROUTINE
2258C     REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS
2259C     LOW THROUGH IGH TO UPPER HESSENBERG FORM BY
2260C     STABILIZED ELEMENTARY SIMILARITY TRANSFORMATIONS.
2261C
2262C     ON INPUT
2263C
2264C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
2265C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
2266C          DIMENSION STATEMENT.
2267C
2268C        N IS THE ORDER OF THE MATRIX.
2269C
2270C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
2271C          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED,
2272C          SET LOW=1, IGH=N.
2273C
2274C        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
2275C          RESPECTIVELY, OF THE COMPLEX INPUT MATRIX.
2276C
2277C     ON OUTPUT
2278C
2279C        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
2280C          RESPECTIVELY, OF THE HESSENBERG MATRIX.  THE
2281C          MULTIPLIERS WHICH WERE USED IN THE REDUCTION
2282C          ARE STORED IN THE REMAINING TRIANGLES UNDER THE
2283C          HESSENBERG MATRIX.
2284C
2285C        INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS
2286C          INTERCHANGED IN THE REDUCTION.
2287C          ONLY ELEMENTS LOW THROUGH IGH ARE USED.
2288C
2289C     CALLS CDIV FOR COMPLEX DIVISION.
2290C
2291C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
2292C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
2293C
2294C     THIS VERSION DATED AUGUST 1983.
2295C
2296C     ------------------------------------------------------------------
2297C
2298      LA = IGH - 1
2299      KP1 = LOW + 1
2300      IF (LA .LT. KP1) GO TO 200
2301C
2302      DO 180 M = KP1, LA
2303         MM1 = M - 1
2304         XR = 0.0D0
2305         XI = 0.0D0
2306         I = M
2307C
2308         DO 100 J = M, IGH
2309            IF (DABS(AR(J,MM1)) + DABS(AI(J,MM1))
2310     X         .LE. DABS(XR) + DABS(XI)) GO TO 100
2311            XR = AR(J,MM1)
2312            XI = AI(J,MM1)
2313            I = J
2314  100    CONTINUE
2315C
2316         INT(M) = I
2317         IF (I .EQ. M) GO TO 130
2318C     .......... INTERCHANGE ROWS AND COLUMNS OF AR AND AI ..........
2319         DO 110 J = MM1, N
2320            YR = AR(I,J)
2321            AR(I,J) = AR(M,J)
2322            AR(M,J) = YR
2323            YI = AI(I,J)
2324            AI(I,J) = AI(M,J)
2325            AI(M,J) = YI
2326  110    CONTINUE
2327C
2328         DO 120 J = 1, IGH
2329            YR = AR(J,I)
2330            AR(J,I) = AR(J,M)
2331            AR(J,M) = YR
2332            YI = AI(J,I)
2333            AI(J,I) = AI(J,M)
2334            AI(J,M) = YI
2335  120    CONTINUE
2336C     .......... END INTERCHANGE ..........
2337  130    IF (XR .EQ. 0.0D0 .AND. XI .EQ. 0.0D0) GO TO 180
2338         MP1 = M + 1
2339C
2340         DO 160 I = MP1, IGH
2341            YR = AR(I,MM1)
2342            YI = AI(I,MM1)
2343            IF (YR .EQ. 0.0D0 .AND. YI .EQ. 0.0D0) GO TO 160
2344            CALL CDIV(YR,YI,XR,XI,YR,YI)
2345            AR(I,MM1) = YR
2346            AI(I,MM1) = YI
2347C
2348            DO 140 J = M, N
2349               AR(I,J) = AR(I,J) - YR * AR(M,J) + YI * AI(M,J)
2350               AI(I,J) = AI(I,J) - YR * AI(M,J) - YI * AR(M,J)
2351  140       CONTINUE
2352C
2353            DO 150 J = 1, IGH
2354               AR(J,M) = AR(J,M) + YR * AR(J,I) - YI * AI(J,I)
2355               AI(J,M) = AI(J,M) + YR * AI(J,I) + YI * AR(J,I)
2356  150       CONTINUE
2357C
2358  160    CONTINUE
2359C
2360  180 CONTINUE
2361C
2362  200 RETURN
2363      END
2364      SUBROUTINE COMLR(NM,N,LOW,IGH,HR,HI,WR,WI,IERR)
2365C
2366      INTEGER I,J,L,M,N,EN,LL,MM,NM,IGH,IM1,ITN,ITS,LOW,MP1,ENM1,IERR
2367      DOUBLE PRECISION HR(NM,N),HI(NM,N),WR(N),WI(N)
2368      DOUBLE PRECISION SI,SR,TI,TR,XI,XR,YI,YR,ZZI,ZZR,TST1,TST2
2369C
2370C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMLR,
2371C     NUM. MATH. 12, 369-376(1968) BY MARTIN AND WILKINSON.
2372C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 396-403(1971).
2373C
2374C     THIS SUBROUTINE FINDS THE EIGENVALUES OF A COMPLEX
2375C     UPPER HESSENBERG MATRIX BY THE MODIFIED LR METHOD.
2376C
2377C     ON INPUT
2378C
2379C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
2380C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
2381C          DIMENSION STATEMENT.
2382C
2383C        N IS THE ORDER OF THE MATRIX.
2384C
2385C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
2386C          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED,
2387C          SET LOW=1, IGH=N.
2388C
2389C        HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS,
2390C          RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX.
2391C          THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN THE
2392C          MULTIPLIERS WHICH WERE USED IN THE REDUCTION BY  COMHES,
2393C          IF PERFORMED.
2394C
2395C     ON OUTPUT
2396C
2397C        THE UPPER HESSENBERG PORTIONS OF HR AND HI HAVE BEEN
2398C          DESTROYED.  THEREFORE, THEY MUST BE SAVED BEFORE
2399C          CALLING  COMLR  IF SUBSEQUENT CALCULATION OF
2400C          EIGENVECTORS IS TO BE PERFORMED.
2401C
2402C        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
2403C          RESPECTIVELY, OF THE EIGENVALUES.  IF AN ERROR
2404C          EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT
2405C          FOR INDICES IERR+1,...,N.
2406C
2407C        IERR IS SET TO
2408C          ZERO       FOR NORMAL RETURN,
2409C          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED
2410C                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT.
2411C
2412C     CALLS CDIV FOR COMPLEX DIVISION.
2413C     CALLS CSROOT FOR COMPLEX SQUARE ROOT.
2414C
2415C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
2416C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
2417C
2418C     THIS VERSION DATED AUGUST 1983.
2419C
2420C     ------------------------------------------------------------------
2421C
2422      IERR = 0
2423C     .......... STORE ROOTS ISOLATED BY CBAL ..........
2424      DO 200 I = 1, N
2425         IF (I .GE. LOW .AND. I .LE. IGH) GO TO 200
2426         WR(I) = HR(I,I)
2427         WI(I) = HI(I,I)
2428  200 CONTINUE
2429C
2430      EN = IGH
2431      TR = 0.0D0
2432      TI = 0.0D0
2433      ITN = 30*N
2434C     .......... SEARCH FOR NEXT EIGENVALUE ..........
2435  220 IF (EN .LT. LOW) GO TO 1001
2436      ITS = 0
2437      ENM1 = EN - 1
2438C     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
2439C                FOR L=EN STEP -1 UNTIL LOW D0 -- ..........
2440  240 DO 260 LL = LOW, EN
2441         L = EN + LOW - LL
2442         IF (L .EQ. LOW) GO TO 300
2443         TST1 = DABS(HR(L-1,L-1)) + DABS(HI(L-1,L-1))
2444     X            + DABS(HR(L,L)) + DABS(HI(L,L))
2445         TST2 = TST1 + DABS(HR(L,L-1)) + DABS(HI(L,L-1))
2446         IF (TST2 .EQ. TST1) GO TO 300
2447  260 CONTINUE
2448C     .......... FORM SHIFT ..........
2449  300 IF (L .EQ. EN) GO TO 660
2450      IF (ITN .EQ. 0) GO TO 1000
2451      IF (ITS .EQ. 10 .OR. ITS .EQ. 20) GO TO 320
2452      SR = HR(EN,EN)
2453      SI = HI(EN,EN)
2454      XR = HR(ENM1,EN) * HR(EN,ENM1) - HI(ENM1,EN) * HI(EN,ENM1)
2455      XI = HR(ENM1,EN) * HI(EN,ENM1) + HI(ENM1,EN) * HR(EN,ENM1)
2456      IF (XR .EQ. 0.0D0 .AND. XI .EQ. 0.0D0) GO TO 340
2457      YR = (HR(ENM1,ENM1) - SR) / 2.0D0
2458      YI = (HI(ENM1,ENM1) - SI) / 2.0D0
2459      CALL CSROOT(YR**2-YI**2+XR,2.0D0*YR*YI+XI,ZZR,ZZI)
2460      IF (YR * ZZR + YI * ZZI .GE. 0.0D0) GO TO 310
2461      ZZR = -ZZR
2462      ZZI = -ZZI
2463  310 CALL CDIV(XR,XI,YR+ZZR,YI+ZZI,XR,XI)
2464      SR = SR - XR
2465      SI = SI - XI
2466      GO TO 340
2467C     .......... FORM EXCEPTIONAL SHIFT ..........
2468  320 SR = DABS(HR(EN,ENM1)) + DABS(HR(ENM1,EN-2))
2469      SI = DABS(HI(EN,ENM1)) + DABS(HI(ENM1,EN-2))
2470C
2471  340 DO 360 I = LOW, EN
2472         HR(I,I) = HR(I,I) - SR
2473         HI(I,I) = HI(I,I) - SI
2474  360 CONTINUE
2475C
2476      TR = TR + SR
2477      TI = TI + SI
2478      ITS = ITS + 1
2479      ITN = ITN - 1
2480C     .......... LOOK FOR TWO CONSECUTIVE SMALL
2481C                SUB-DIAGONAL ELEMENTS ..........
2482      XR = DABS(HR(ENM1,ENM1)) + DABS(HI(ENM1,ENM1))
2483      YR = DABS(HR(EN,ENM1)) + DABS(HI(EN,ENM1))
2484      ZZR = DABS(HR(EN,EN)) + DABS(HI(EN,EN))
2485C     .......... FOR M=EN-1 STEP -1 UNTIL L DO -- ..........
2486      DO 380 MM = L, ENM1
2487         M = ENM1 + L - MM
2488         IF (M .EQ. L) GO TO 420
2489         YI = YR
2490         YR = DABS(HR(M,M-1)) + DABS(HI(M,M-1))
2491         XI = ZZR
2492         ZZR = XR
2493         XR = DABS(HR(M-1,M-1)) + DABS(HI(M-1,M-1))
2494         TST1 = ZZR / YI * (ZZR + XR + XI)
2495         TST2 = TST1 + YR
2496         IF (TST2 .EQ. TST1) GO TO 420
2497  380 CONTINUE
2498C     .......... TRIANGULAR DECOMPOSITION H=L*R ..........
2499  420 MP1 = M + 1
2500C
2501      DO 520 I = MP1, EN
2502         IM1 = I - 1
2503         XR = HR(IM1,IM1)
2504         XI = HI(IM1,IM1)
2505         YR = HR(I,IM1)
2506         YI = HI(I,IM1)
2507         IF (DABS(XR) + DABS(XI) .GE. DABS(YR) + DABS(YI)) GO TO 460
2508C     .......... INTERCHANGE ROWS OF HR AND HI ..........
2509         DO 440 J = IM1, EN
2510            ZZR = HR(IM1,J)
2511            HR(IM1,J) = HR(I,J)
2512            HR(I,J) = ZZR
2513            ZZI = HI(IM1,J)
2514            HI(IM1,J) = HI(I,J)
2515            HI(I,J) = ZZI
2516  440    CONTINUE
2517C
2518         CALL CDIV(XR,XI,YR,YI,ZZR,ZZI)
2519         WR(I) = 1.0D0
2520         GO TO 480
2521  460    CALL CDIV(YR,YI,XR,XI,ZZR,ZZI)
2522         WR(I) = -1.0D0
2523  480    HR(I,IM1) = ZZR
2524         HI(I,IM1) = ZZI
2525C
2526         DO 500 J = I, EN
2527            HR(I,J) = HR(I,J) - ZZR * HR(IM1,J) + ZZI * HI(IM1,J)
2528            HI(I,J) = HI(I,J) - ZZR * HI(IM1,J) - ZZI * HR(IM1,J)
2529  500    CONTINUE
2530C
2531  520 CONTINUE
2532C     .......... COMPOSITION R*L=H ..........
2533      DO 640 J = MP1, EN
2534         XR = HR(J,J-1)
2535         XI = HI(J,J-1)
2536         HR(J,J-1) = 0.0D0
2537         HI(J,J-1) = 0.0D0
2538C     .......... INTERCHANGE COLUMNS OF HR AND HI,
2539C                IF NECESSARY ..........
2540         IF (WR(J) .LE. 0.0D0) GO TO 580
2541C
2542         DO 540 I = L, J
2543            ZZR = HR(I,J-1)
2544            HR(I,J-1) = HR(I,J)
2545            HR(I,J) = ZZR
2546            ZZI = HI(I,J-1)
2547            HI(I,J-1) = HI(I,J)
2548            HI(I,J) = ZZI
2549  540    CONTINUE
2550C
2551  580    DO 600 I = L, J
2552            HR(I,J-1) = HR(I,J-1) + XR * HR(I,J) - XI * HI(I,J)
2553            HI(I,J-1) = HI(I,J-1) + XR * HI(I,J) + XI * HR(I,J)
2554  600    CONTINUE
2555C
2556  640 CONTINUE
2557C
2558      GO TO 240
2559C     .......... A ROOT FOUND ..........
2560  660 WR(EN) = HR(EN,EN) + TR
2561      WI(EN) = HI(EN,EN) + TI
2562      EN = ENM1
2563      GO TO 220
2564C     .......... SET ERROR -- ALL EIGENVALUES HAVE NOT
2565C                CONVERGED AFTER 30*N ITERATIONS ..........
2566 1000 IERR = EN
2567 1001 RETURN
2568      END
2569      SUBROUTINE COMLR2(NM,N,LOW,IGH,INT,HR,HI,WR,WI,ZR,ZI,IERR)
2570C
2571      INTEGER I,J,K,L,M,N,EN,II,JJ,LL,MM,NM,NN,IGH,IM1,IP1,
2572     X        ITN,ITS,LOW,MP1,ENM1,IEND,IERR
2573      DOUBLE PRECISION HR(NM,N),HI(NM,N),WR(N),WI(N),ZR(NM,N),ZI(NM,N)
2574      DOUBLE PRECISION SI,SR,TI,TR,XI,XR,YI,YR,ZZI,ZZR,NORM,TST1,TST2
2575      INTEGER INT(IGH)
2576C
2577C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE COMLR2,
2578C     NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON.
2579C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
2580C
2581C     THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
2582C     OF A COMPLEX UPPER HESSENBERG MATRIX BY THE MODIFIED LR
2583C     METHOD.  THE EIGENVECTORS OF A COMPLEX GENERAL MATRIX
2584C     CAN ALSO BE FOUND IF  COMHES  HAS BEEN USED TO REDUCE
2585C     THIS GENERAL MATRIX TO HESSENBERG FORM.
2586C
2587C     ON INPUT
2588C
2589C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
2590C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
2591C          DIMENSION STATEMENT.
2592C
2593C        N IS THE ORDER OF THE MATRIX.
2594C
2595C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
2596C          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED,
2597C          SET LOW=1, IGH=N.
2598C
2599C        INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS INTERCHANGED
2600C          IN THE REDUCTION BY  COMHES, IF PERFORMED.  ONLY ELEMENTS
2601C          LOW THROUGH IGH ARE USED.  IF THE EIGENVECTORS OF THE HESSEN-
2602C          BERG MATRIX ARE DESIRED, SET INT(J)=J FOR THESE ELEMENTS.
2603C
2604C        HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS,
2605C          RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX.
2606C          THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN THE
2607C          MULTIPLIERS WHICH WERE USED IN THE REDUCTION BY  COMHES,
2608C          IF PERFORMED.  IF THE EIGENVECTORS OF THE HESSENBERG
2609C          MATRIX ARE DESIRED, THESE ELEMENTS MUST BE SET TO ZERO.
2610C
2611C     ON OUTPUT
2612C
2613C        THE UPPER HESSENBERG PORTIONS OF HR AND HI HAVE BEEN
2614C          DESTROYED, BUT THE LOCATION HR(1,1) CONTAINS THE NORM
2615C          OF THE TRIANGULARIZED MATRIX.
2616C
2617C        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
2618C          RESPECTIVELY, OF THE EIGENVALUES.  IF AN ERROR
2619C          EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT
2620C          FOR INDICES IERR+1,...,N.
2621C
2622C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
2623C          RESPECTIVELY, OF THE EIGENVECTORS.  THE EIGENVECTORS
2624C          ARE UNNORMALIZED.  IF AN ERROR EXIT IS MADE, NONE OF
2625C          THE EIGENVECTORS HAS BEEN FOUND.
2626C
2627C        IERR IS SET TO
2628C          ZERO       FOR NORMAL RETURN,
2629C          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED
2630C                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT.
2631C
2632C     
2633C     CALLS CDIV FOR COMPLEX DIVISION.
2634C     CALLS CSROOT FOR COMPLEX SQUARE ROOT.
2635C
2636C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
2637C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
2638C
2639C     THIS VERSION DATED AUGUST 1983.
2640C
2641C     ------------------------------------------------------------------
2642C
2643      IERR = 0
2644C     .......... INITIALIZE EIGENVECTOR MATRIX ..........
2645      DO 100 I = 1, N
2646C
2647         DO 100 J = 1, N
2648            ZR(I,J) = 0.0D0
2649            ZI(I,J) = 0.0D0
2650            IF (I .EQ. J) ZR(I,J) = 1.0D0
2651  100 CONTINUE
2652C     .......... FORM THE MATRIX OF ACCUMULATED TRANSFORMATIONS
2653C                FROM THE INFORMATION LEFT BY COMHES ..........
2654      IEND = IGH - LOW - 1
2655      IF (IEND .LE. 0) GO TO 180
2656C     .......... FOR I=IGH-1 STEP -1 UNTIL LOW+1 DO -- ..........
2657      DO 160 II = 1, IEND
2658         I = IGH - II
2659         IP1 = I + 1
2660C
2661         DO 120 K = IP1, IGH
2662            ZR(K,I) = HR(K,I-1)
2663            ZI(K,I) = HI(K,I-1)
2664  120    CONTINUE
2665C
2666         J = INT(I)
2667         IF (I .EQ. J) GO TO 160
2668C
2669         DO 140 K = I, IGH
2670            ZR(I,K) = ZR(J,K)
2671            ZI(I,K) = ZI(J,K)
2672            ZR(J,K) = 0.0D0
2673            ZI(J,K) = 0.0D0
2674  140    CONTINUE
2675C
2676         ZR(J,I) = 1.0D0
2677  160 CONTINUE
2678C     .......... STORE ROOTS ISOLATED BY CBAL ..........
2679  180 DO 200 I = 1, N
2680         IF (I .GE. LOW .AND. I .LE. IGH) GO TO 200
2681         WR(I) = HR(I,I)
2682         WI(I) = HI(I,I)
2683  200 CONTINUE
2684C
2685      EN = IGH
2686      TR = 0.0D0
2687      TI = 0.0D0
2688      ITN = 30*N
2689C     .......... SEARCH FOR NEXT EIGENVALUE ..........
2690  220 IF (EN .LT. LOW) GO TO 680
2691      ITS = 0
2692      ENM1 = EN - 1
2693C     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
2694C                FOR L=EN STEP -1 UNTIL LOW DO -- ..........
2695  240 DO 260 LL = LOW, EN
2696         L = EN + LOW - LL
2697         IF (L .EQ. LOW) GO TO 300
2698         TST1 = DABS(HR(L-1,L-1)) + DABS(HI(L-1,L-1))
2699     X            + DABS(HR(L,L)) + DABS(HI(L,L))
2700         TST2 = TST1 + DABS(HR(L,L-1)) + DABS(HI(L,L-1))
2701         IF (TST2 .EQ. TST1) GO TO 300
2702  260 CONTINUE
2703C     .......... FORM SHIFT ..........
2704  300 IF (L .EQ. EN) GO TO 660
2705      IF (ITN .EQ. 0) GO TO 1000
2706      IF (ITS .EQ. 10 .OR. ITS .EQ. 20) GO TO 320
2707      SR = HR(EN,EN)
2708      SI = HI(EN,EN)
2709      XR = HR(ENM1,EN) * HR(EN,ENM1) - HI(ENM1,EN) * HI(EN,ENM1)
2710      XI = HR(ENM1,EN) * HI(EN,ENM1) + HI(ENM1,EN) * HR(EN,ENM1)
2711      IF (XR .EQ. 0.0D0 .AND. XI .EQ. 0.0D0) GO TO 340
2712      YR = (HR(ENM1,ENM1) - SR) / 2.0D0
2713      YI = (HI(ENM1,ENM1) - SI) / 2.0D0
2714      CALL CSROOT(YR**2-YI**2+XR,2.0D0*YR*YI+XI,ZZR,ZZI)
2715      IF (YR * ZZR + YI * ZZI .GE. 0.0D0) GO TO 310
2716      ZZR = -ZZR
2717      ZZI = -ZZI
2718  310 CALL CDIV(XR,XI,YR+ZZR,YI+ZZI,XR,XI)
2719      SR = SR - XR
2720      SI = SI - XI
2721      GO TO 340
2722C     .......... FORM EXCEPTIONAL SHIFT ..........
2723  320 SR = DABS(HR(EN,ENM1)) + DABS(HR(ENM1,EN-2))
2724      SI = DABS(HI(EN,ENM1)) + DABS(HI(ENM1,EN-2))
2725C
2726  340 DO 360 I = LOW, EN
2727         HR(I,I) = HR(I,I) - SR
2728         HI(I,I) = HI(I,I) - SI
2729  360 CONTINUE
2730C
2731      TR = TR + SR
2732      TI = TI + SI
2733      ITS = ITS + 1
2734      ITN = ITN - 1
2735C     .......... LOOK FOR TWO CONSECUTIVE SMALL
2736C                SUB-DIAGONAL ELEMENTS ..........
2737      XR = DABS(HR(ENM1,ENM1)) + DABS(HI(ENM1,ENM1))
2738      YR = DABS(HR(EN,ENM1)) + DABS(HI(EN,ENM1))
2739      ZZR = DABS(HR(EN,EN)) + DABS(HI(EN,EN))
2740C     .......... FOR M=EN-1 STEP -1 UNTIL L DO -- ..........
2741      DO 380 MM = L, ENM1
2742         M = ENM1 + L - MM
2743         IF (M .EQ. L) GO TO 420
2744         YI = YR
2745         YR = DABS(HR(M,M-1)) + DABS(HI(M,M-1))
2746         XI = ZZR
2747         ZZR = XR
2748         XR = DABS(HR(M-1,M-1)) + DABS(HI(M-1,M-1))
2749         TST1 = ZZR / YI * (ZZR + XR + XI)
2750         TST2 = TST1 + YR
2751         IF (TST2 .EQ. TST1) GO TO 420
2752  380 CONTINUE
2753C     .......... TRIANGULAR DECOMPOSITION H=L*R ..........
2754  420 MP1 = M + 1
2755C
2756      DO 520 I = MP1, EN
2757         IM1 = I - 1
2758         XR = HR(IM1,IM1)
2759         XI = HI(IM1,IM1)
2760         YR = HR(I,IM1)
2761         YI = HI(I,IM1)
2762         IF (DABS(XR) + DABS(XI) .GE. DABS(YR) + DABS(YI)) GO TO 460
2763C     .......... INTERCHANGE ROWS OF HR AND HI ..........
2764         DO 440 J = IM1, N
2765            ZZR = HR(IM1,J)
2766            HR(IM1,J) = HR(I,J)
2767            HR(I,J) = ZZR
2768            ZZI = HI(IM1,J)
2769            HI(IM1,J) = HI(I,J)
2770            HI(I,J) = ZZI
2771  440    CONTINUE
2772C
2773         CALL CDIV(XR,XI,YR,YI,ZZR,ZZI)
2774         WR(I) = 1.0D0
2775         GO TO 480
2776  460    CALL CDIV(YR,YI,XR,XI,ZZR,ZZI)
2777         WR(I) = -1.0D0
2778  480    HR(I,IM1) = ZZR
2779         HI(I,IM1) = ZZI
2780C
2781         DO 500 J = I, N
2782            HR(I,J) = HR(I,J) - ZZR * HR(IM1,J) + ZZI * HI(IM1,J)
2783            HI(I,J) = HI(I,J) - ZZR * HI(IM1,J) - ZZI * HR(IM1,J)
2784  500    CONTINUE
2785C
2786  520 CONTINUE
2787C     .......... COMPOSITION R*L=H ..........
2788      DO 640 J = MP1, EN
2789         XR = HR(J,J-1)
2790         XI = HI(J,J-1)
2791         HR(J,J-1) = 0.0D0
2792         HI(J,J-1) = 0.0D0
2793C     .......... INTERCHANGE COLUMNS OF HR, HI, ZR, AND ZI,
2794C                IF NECESSARY ..........
2795         IF (WR(J) .LE. 0.0D0) GO TO 580
2796C
2797         DO 540 I = 1, J
2798            ZZR = HR(I,J-1)
2799            HR(I,J-1) = HR(I,J)
2800            HR(I,J) = ZZR
2801            ZZI = HI(I,J-1)
2802            HI(I,J-1) = HI(I,J)
2803            HI(I,J) = ZZI
2804  540    CONTINUE
2805C
2806         DO 560 I = LOW, IGH
2807            ZZR = ZR(I,J-1)
2808            ZR(I,J-1) = ZR(I,J)
2809            ZR(I,J) = ZZR
2810            ZZI = ZI(I,J-1)
2811            ZI(I,J-1) = ZI(I,J)
2812            ZI(I,J) = ZZI
2813  560    CONTINUE
2814C
2815  580    DO 600 I = 1, J
2816            HR(I,J-1) = HR(I,J-1) + XR * HR(I,J) - XI * HI(I,J)
2817            HI(I,J-1) = HI(I,J-1) + XR * HI(I,J) + XI * HR(I,J)
2818  600    CONTINUE
2819C     .......... ACCUMULATE TRANSFORMATIONS ..........
2820         DO 620 I = LOW, IGH
2821            ZR(I,J-1) = ZR(I,J-1) + XR * ZR(I,J) - XI * ZI(I,J)
2822            ZI(I,J-1) = ZI(I,J-1) + XR * ZI(I,J) + XI * ZR(I,J)
2823  620    CONTINUE
2824C
2825  640 CONTINUE
2826C
2827      GO TO 240
2828C     .......... A ROOT FOUND ..........
2829  660 HR(EN,EN) = HR(EN,EN) + TR
2830      WR(EN) = HR(EN,EN)
2831      HI(EN,EN) = HI(EN,EN) + TI
2832      WI(EN) = HI(EN,EN)
2833      EN = ENM1
2834      GO TO 220
2835C     .......... ALL ROOTS FOUND.  BACKSUBSTITUTE TO FIND
2836C                VECTORS OF UPPER TRIANGULAR FORM ..........
2837  680 NORM = 0.0D0
2838C
2839      DO 720 I = 1, N
2840C
2841         DO 720 J = I, N
2842            TR = DABS(HR(I,J)) + DABS(HI(I,J))
2843            IF (TR .GT. NORM) NORM = TR
2844  720 CONTINUE
2845C
2846      HR(1,1) = NORM
2847      IF (N .EQ. 1 .OR. NORM .EQ. 0.0D0) GO TO 1001
2848C     .......... FOR EN=N STEP -1 UNTIL 2 DO -- ..........
2849      DO 800 NN = 2, N
2850         EN = N + 2 - NN
2851         XR = WR(EN)
2852         XI = WI(EN)
2853         HR(EN,EN) = 1.0D0
2854         HI(EN,EN) = 0.0D0
2855         ENM1 = EN - 1
2856C     .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- ..........
2857         DO 780 II = 1, ENM1
2858            I = EN - II
2859            ZZR = 0.0D0
2860            ZZI = 0.0D0
2861            IP1 = I + 1
2862C
2863            DO 740 J = IP1, EN
2864               ZZR = ZZR + HR(I,J) * HR(J,EN) - HI(I,J) * HI(J,EN)
2865               ZZI = ZZI + HR(I,J) * HI(J,EN) + HI(I,J) * HR(J,EN)
2866  740       CONTINUE
2867C
2868            YR = XR - WR(I)
2869            YI = XI - WI(I)
2870            IF (YR .NE. 0.0D0 .OR. YI .NE. 0.0D0) GO TO 765
2871               TST1 = NORM
2872               YR = TST1
2873  760          YR = 0.01D0 * YR
2874               TST2 = NORM + YR
2875               IF (TST2 .GT. TST1) GO TO 760
2876  765       CONTINUE
2877            CALL CDIV(ZZR,ZZI,YR,YI,HR(I,EN),HI(I,EN))
2878C     .......... OVERFLOW CONTROL ..........
2879            TR = DABS(HR(I,EN)) + DABS(HI(I,EN))
2880            IF (TR .EQ. 0.0D0) GO TO 780
2881            TST1 = TR
2882            TST2 = TST1 + 1.0D0/TST1
2883            IF (TST2 .GT. TST1) GO TO 780
2884            DO 770 J = I, EN
2885               HR(J,EN) = HR(J,EN)/TR
2886               HI(J,EN) = HI(J,EN)/TR
2887  770       CONTINUE
2888C
2889  780    CONTINUE
2890C
2891  800 CONTINUE
2892C     .......... END BACKSUBSTITUTION ..........
2893      ENM1 = N - 1
2894C     .......... VECTORS OF ISOLATED ROOTS ..........
2895      DO  840 I = 1, ENM1
2896         IF (I .GE. LOW .AND. I .LE. IGH) GO TO 840
2897         IP1 = I + 1
2898C
2899         DO 820 J = IP1, N
2900            ZR(I,J) = HR(I,J)
2901            ZI(I,J) = HI(I,J)
2902  820    CONTINUE
2903C
2904  840 CONTINUE
2905C     .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE
2906C                VECTORS OF ORIGINAL FULL MATRIX.
2907C                FOR J=N STEP -1 UNTIL LOW+1 DO -- ..........
2908      DO 880 JJ = LOW, ENM1
2909         J = N + LOW - JJ
2910         M = MIN0(J,IGH)
2911C
2912         DO 880 I = LOW, IGH
2913            ZZR = 0.0D0
2914            ZZI = 0.0D0
2915C
2916            DO 860 K = LOW, M
2917               ZZR = ZZR + ZR(I,K) * HR(K,J) - ZI(I,K) * HI(K,J)
2918               ZZI = ZZI + ZR(I,K) * HI(K,J) + ZI(I,K) * HR(K,J)
2919  860       CONTINUE
2920C
2921            ZR(I,J) = ZZR
2922            ZI(I,J) = ZZI
2923  880 CONTINUE
2924C
2925      GO TO 1001
2926C     .......... SET ERROR -- ALL EIGENVALUES HAVE NOT
2927C                CONVERGED AFTER 30*N ITERATIONS ..........
2928 1000 IERR = EN
2929 1001 RETURN
2930      END
2931      SUBROUTINE COMQR(NM,N,LOW,IGH,HR,HI,WR,WI,IERR)
2932C
2933      INTEGER I,J,L,N,EN,LL,NM,IGH,ITN,ITS,LOW,LP1,ENM1,IERR
2934      DOUBLE PRECISION HR(NM,N),HI(NM,N),WR(N),WI(N)
2935      DOUBLE PRECISION SI,SR,TI,TR,XI,XR,YI,YR,ZZI,ZZR,NORM,TST1,TST2,
2936     X       PYTHAG
2937C
2938C     THIS SUBROUTINE IS A TRANSLATION OF A UNITARY ANALOGUE OF THE
2939C     ALGOL PROCEDURE  COMLR, NUM. MATH. 12, 369-376(1968) BY MARTIN
2940C     AND WILKINSON.
2941C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 396-403(1971).
2942C     THE UNITARY ANALOGUE SUBSTITUTES THE QR ALGORITHM OF FRANCIS
2943C     (COMP. JOUR. 4, 332-345(1962)) FOR THE LR ALGORITHM.
2944C
2945C     THIS SUBROUTINE FINDS THE EIGENVALUES OF A COMPLEX
2946C     UPPER HESSENBERG MATRIX BY THE QR METHOD.
2947C
2948C     ON INPUT
2949C
2950C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
2951C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
2952C          DIMENSION STATEMENT.
2953C
2954C        N IS THE ORDER OF THE MATRIX.
2955C
2956C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
2957C          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED,
2958C          SET LOW=1, IGH=N.
2959C
2960C        HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS,
2961C          RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX.
2962C          THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN
2963C          INFORMATION ABOUT THE UNITARY TRANSFORMATIONS USED IN
2964C          THE REDUCTION BY  CORTH, IF PERFORMED.
2965C
2966C     ON OUTPUT
2967C
2968C        THE UPPER HESSENBERG PORTIONS OF HR AND HI HAVE BEEN
2969C          DESTROYED.  THEREFORE, THEY MUST BE SAVED BEFORE
2970C          CALLING  COMQR  IF SUBSEQUENT CALCULATION OF
2971C          EIGENVECTORS IS TO BE PERFORMED.
2972C
2973C        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
2974C          RESPECTIVELY, OF THE EIGENVALUES.  IF AN ERROR
2975C          EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT
2976C          FOR INDICES IERR+1,...,N.
2977C
2978C        IERR IS SET TO
2979C          ZERO       FOR NORMAL RETURN,
2980C          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED
2981C                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT.
2982C
2983C     CALLS CDIV FOR COMPLEX DIVISION.
2984C     CALLS CSROOT FOR COMPLEX SQUARE ROOT.
2985C     CALLS PYTHAG FOR  DSQRT(A*A + B*B) .
2986C
2987C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
2988C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
2989C
2990C     THIS VERSION DATED AUGUST 1983.
2991C
2992C     ------------------------------------------------------------------
2993C
2994      IERR = 0
2995      IF (LOW .EQ. IGH) GO TO 180
2996C     .......... CREATE REAL SUBDIAGONAL ELEMENTS ..........
2997      L = LOW + 1
2998C
2999      DO 170 I = L, IGH
3000         LL = MIN0(I+1,IGH)
3001         IF (HI(I,I-1) .EQ. 0.0D0) GO TO 170
3002         NORM = PYTHAG(HR(I,I-1),HI(I,I-1))
3003         YR = HR(I,I-1) / NORM
3004         YI = HI(I,I-1) / NORM
3005         HR(I,I-1) = NORM
3006         HI(I,I-1) = 0.0D0
3007C
3008         DO 155 J = I, IGH
3009            SI = YR * HI(I,J) - YI * HR(I,J)
3010            HR(I,J) = YR * HR(I,J) + YI * HI(I,J)
3011            HI(I,J) = SI
3012  155    CONTINUE
3013C
3014         DO 160 J = LOW, LL
3015            SI = YR * HI(J,I) + YI * HR(J,I)
3016            HR(J,I) = YR * HR(J,I) - YI * HI(J,I)
3017            HI(J,I) = SI
3018  160    CONTINUE
3019C
3020  170 CONTINUE
3021C     .......... STORE ROOTS ISOLATED BY CBAL ..........
3022  180 DO 200 I = 1, N
3023         IF (I .GE. LOW .AND. I .LE. IGH) GO TO 200
3024         WR(I) = HR(I,I)
3025         WI(I) = HI(I,I)
3026  200 CONTINUE
3027C
3028      EN = IGH
3029      TR = 0.0D0
3030      TI = 0.0D0
3031      ITN = 30*N
3032C     .......... SEARCH FOR NEXT EIGENVALUE ..........
3033  220 IF (EN .LT. LOW) GO TO 1001
3034      ITS = 0
3035      ENM1 = EN - 1
3036C     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
3037C                FOR L=EN STEP -1 UNTIL LOW D0 -- ..........
3038  240 DO 260 LL = LOW, EN
3039         L = EN + LOW - LL
3040         IF (L .EQ. LOW) GO TO 300
3041         TST1 = DABS(HR(L-1,L-1)) + DABS(HI(L-1,L-1))
3042     X            + DABS(HR(L,L)) + DABS(HI(L,L))
3043         TST2 = TST1 + DABS(HR(L,L-1))
3044         IF (TST2 .EQ. TST1) GO TO 300
3045  260 CONTINUE
3046C     .......... FORM SHIFT ..........
3047  300 IF (L .EQ. EN) GO TO 660
3048      IF (ITN .EQ. 0) GO TO 1000
3049      IF (ITS .EQ. 10 .OR. ITS .EQ. 20) GO TO 320
3050      SR = HR(EN,EN)
3051      SI = HI(EN,EN)
3052      XR = HR(ENM1,EN) * HR(EN,ENM1)
3053      XI = HI(ENM1,EN) * HR(EN,ENM1)
3054      IF (XR .EQ. 0.0D0 .AND. XI .EQ. 0.0D0) GO TO 340
3055      YR = (HR(ENM1,ENM1) - SR) / 2.0D0
3056      YI = (HI(ENM1,ENM1) - SI) / 2.0D0
3057      CALL CSROOT(YR**2-YI**2+XR,2.0D0*YR*YI+XI,ZZR,ZZI)
3058      IF (YR * ZZR + YI * ZZI .GE. 0.0D0) GO TO 310
3059      ZZR = -ZZR
3060      ZZI = -ZZI
3061  310 CALL CDIV(XR,XI,YR+ZZR,YI+ZZI,XR,XI)
3062      SR = SR - XR
3063      SI = SI - XI
3064      GO TO 340
3065C     .......... FORM EXCEPTIONAL SHIFT ..........
3066  320 SR = DABS(HR(EN,ENM1)) + DABS(HR(ENM1,EN-2))
3067      SI = 0.0D0
3068C
3069  340 DO 360 I = LOW, EN
3070         HR(I,I) = HR(I,I) - SR
3071         HI(I,I) = HI(I,I) - SI
3072  360 CONTINUE
3073C
3074      TR = TR + SR
3075      TI = TI + SI
3076      ITS = ITS + 1
3077      ITN = ITN - 1
3078C     .......... REDUCE TO TRIANGLE (ROWS) ..........
3079      LP1 = L + 1
3080C
3081      DO 500 I = LP1, EN
3082         SR = HR(I,I-1)
3083         HR(I,I-1) = 0.0D0
3084         NORM = PYTHAG(PYTHAG(HR(I-1,I-1),HI(I-1,I-1)),SR)
3085         XR = HR(I-1,I-1) / NORM
3086         WR(I-1) = XR
3087         XI = HI(I-1,I-1) / NORM
3088         WI(I-1) = XI
3089         HR(I-1,I-1) = NORM
3090         HI(I-1,I-1) = 0.0D0
3091         HI(I,I-1) = SR / NORM
3092C
3093         DO 490 J = I, EN
3094            YR = HR(I-1,J)
3095            YI = HI(I-1,J)
3096            ZZR = HR(I,J)
3097            ZZI = HI(I,J)
3098            HR(I-1,J) = XR * YR + XI * YI + HI(I,I-1) * ZZR
3099            HI(I-1,J) = XR * YI - XI * YR + HI(I,I-1) * ZZI
3100            HR(I,J) = XR * ZZR - XI * ZZI - HI(I,I-1) * YR
3101            HI(I,J) = XR * ZZI + XI * ZZR - HI(I,I-1) * YI
3102  490    CONTINUE
3103C
3104  500 CONTINUE
3105C
3106      SI = HI(EN,EN)
3107      IF (SI .EQ. 0.0D0) GO TO 540
3108      NORM = PYTHAG(HR(EN,EN),SI)
3109      SR = HR(EN,EN) / NORM
3110      SI = SI / NORM
3111      HR(EN,EN) = NORM
3112      HI(EN,EN) = 0.0D0
3113C     .......... INVERSE OPERATION (COLUMNS) ..........
3114  540 DO 600 J = LP1, EN
3115         XR = WR(J-1)
3116         XI = WI(J-1)
3117C
3118         DO 580 I = L, J
3119            YR = HR(I,J-1)
3120            YI = 0.0D0
3121            ZZR = HR(I,J)
3122            ZZI = HI(I,J)
3123            IF (I .EQ. J) GO TO 560
3124            YI = HI(I,J-1)
3125            HI(I,J-1) = XR * YI + XI * YR + HI(J,J-1) * ZZI
3126  560       HR(I,J-1) = XR * YR - XI * YI + HI(J,J-1) * ZZR
3127            HR(I,J) = XR * ZZR + XI * ZZI - HI(J,J-1) * YR
3128            HI(I,J) = XR * ZZI - XI * ZZR - HI(J,J-1) * YI
3129  580    CONTINUE
3130C
3131  600 CONTINUE
3132C
3133      IF (SI .EQ. 0.0D0) GO TO 240
3134C
3135      DO 630 I = L, EN
3136         YR = HR(I,EN)
3137         YI = HI(I,EN)
3138         HR(I,EN) = SR * YR - SI * YI
3139         HI(I,EN) = SR * YI + SI * YR
3140  630 CONTINUE
3141C
3142      GO TO 240
3143C     .......... A ROOT FOUND ..........
3144  660 WR(EN) = HR(EN,EN) + TR
3145      WI(EN) = HI(EN,EN) + TI
3146      EN = ENM1
3147      GO TO 220
3148C     .......... SET ERROR -- ALL EIGENVALUES HAVE NOT
3149C                CONVERGED AFTER 30*N ITERATIONS ..........
3150 1000 IERR = EN
3151 1001 RETURN
3152      END
3153      SUBROUTINE COMQR2(NM,N,LOW,IGH,ORTR,ORTI,HR,HI,WR,WI,ZR,ZI,IERR)
3154C
3155      INTEGER I,J,K,L,M,N,EN,II,JJ,LL,NM,NN,IGH,IP1,
3156     X        ITN,ITS,LOW,LP1,ENM1,IEND,IERR
3157      DOUBLE PRECISION HR(NM,N),HI(NM,N),WR(N),WI(N),ZR(NM,N),ZI(NM,N),
3158     X       ORTR(IGH),ORTI(IGH)
3159      DOUBLE PRECISION SI,SR,TI,TR,XI,XR,YI,YR,ZZI,ZZR,NORM,TST1,TST2,
3160     X       PYTHAG
3161C
3162C     THIS SUBROUTINE IS A TRANSLATION OF A UNITARY ANALOGUE OF THE
3163C     ALGOL PROCEDURE  COMLR2, NUM. MATH. 16, 181-204(1970) BY PETERS
3164C     AND WILKINSON.
3165C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
3166C     THE UNITARY ANALOGUE SUBSTITUTES THE QR ALGORITHM OF FRANCIS
3167C     (COMP. JOUR. 4, 332-345(1962)) FOR THE LR ALGORITHM.
3168C
3169C     THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
3170C     OF A COMPLEX UPPER HESSENBERG MATRIX BY THE QR
3171C     METHOD.  THE EIGENVECTORS OF A COMPLEX GENERAL MATRIX
3172C     CAN ALSO BE FOUND IF  CORTH  HAS BEEN USED TO REDUCE
3173C     THIS GENERAL MATRIX TO HESSENBERG FORM.
3174C
3175C     ON INPUT
3176C
3177C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
3178C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
3179C          DIMENSION STATEMENT.
3180C
3181C        N IS THE ORDER OF THE MATRIX.
3182C
3183C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
3184C          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED,
3185C          SET LOW=1, IGH=N.
3186C
3187C        ORTR AND ORTI CONTAIN INFORMATION ABOUT THE UNITARY TRANS-
3188C          FORMATIONS USED IN THE REDUCTION BY  CORTH, IF PERFORMED.
3189C          ONLY ELEMENTS LOW THROUGH IGH ARE USED.  IF THE EIGENVECTORS
3190C          OF THE HESSENBERG MATRIX ARE DESIRED, SET ORTR(J) AND
3191C          ORTI(J) TO 0.0D0 FOR THESE ELEMENTS.
3192C
3193C        HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS,
3194C          RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX.
3195C          THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN FURTHER
3196C          INFORMATION ABOUT THE TRANSFORMATIONS WHICH WERE USED IN THE
3197C          REDUCTION BY  CORTH, IF PERFORMED.  IF THE EIGENVECTORS OF
3198C          THE HESSENBERG MATRIX ARE DESIRED, THESE ELEMENTS MAY BE
3199C          ARBITRARY.
3200C
3201C     ON OUTPUT
3202C
3203C        ORTR, ORTI, AND THE UPPER HESSENBERG PORTIONS OF HR AND HI
3204C          HAVE BEEN DESTROYED.
3205C
3206C        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
3207C          RESPECTIVELY, OF THE EIGENVALUES.  IF AN ERROR
3208C          EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT
3209C          FOR INDICES IERR+1,...,N.
3210C
3211C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
3212C          RESPECTIVELY, OF THE EIGENVECTORS.  THE EIGENVECTORS
3213C          ARE UNNORMALIZED.  IF AN ERROR EXIT IS MADE, NONE OF
3214C          THE EIGENVECTORS HAS BEEN FOUND.
3215C
3216C        IERR IS SET TO
3217C          ZERO       FOR NORMAL RETURN,
3218C          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED
3219C                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT.
3220C
3221C     CALLS CDIV FOR COMPLEX DIVISION.
3222C     CALLS CSROOT FOR COMPLEX SQUARE ROOT.
3223C     CALLS PYTHAG FOR  DSQRT(A*A + B*B) .
3224C
3225C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
3226C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
3227C
3228C     THIS VERSION DATED AUGUST 1983.
3229C
3230C     ------------------------------------------------------------------
3231C
3232      IERR = 0
3233C     .......... INITIALIZE EIGENVECTOR MATRIX ..........
3234      DO 101 J = 1, N
3235C
3236         DO 100 I = 1, N
3237            ZR(I,J) = 0.0D0
3238            ZI(I,J) = 0.0D0
3239  100    CONTINUE
3240         ZR(J,J) = 1.0D0
3241  101 CONTINUE
3242C     .......... FORM THE MATRIX OF ACCUMULATED TRANSFORMATIONS
3243C                FROM THE INFORMATION LEFT BY CORTH ..........
3244      IEND = IGH - LOW - 1
3245      IF (IEND) 180, 150, 105
3246C     .......... FOR I=IGH-1 STEP -1 UNTIL LOW+1 DO -- ..........
3247  105 DO 140 II = 1, IEND
3248         I = IGH - II
3249         IF (ORTR(I) .EQ. 0.0D0 .AND. ORTI(I) .EQ. 0.0D0) GO TO 140
3250         IF (HR(I,I-1) .EQ. 0.0D0 .AND. HI(I,I-1) .EQ. 0.0D0) GO TO 140
3251C     .......... NORM BELOW IS NEGATIVE OF H FORMED IN CORTH ..........
3252         NORM = HR(I,I-1) * ORTR(I) + HI(I,I-1) * ORTI(I)
3253         IP1 = I + 1
3254C
3255         DO 110 K = IP1, IGH
3256            ORTR(K) = HR(K,I-1)
3257            ORTI(K) = HI(K,I-1)
3258  110    CONTINUE
3259C
3260         DO 130 J = I, IGH
3261            SR = 0.0D0
3262            SI = 0.0D0
3263C
3264            DO 115 K = I, IGH
3265               SR = SR + ORTR(K) * ZR(K,J) + ORTI(K) * ZI(K,J)
3266               SI = SI + ORTR(K) * ZI(K,J) - ORTI(K) * ZR(K,J)
3267  115       CONTINUE
3268C
3269            SR = SR / NORM
3270            SI = SI / NORM
3271C
3272            DO 120 K = I, IGH
3273               ZR(K,J) = ZR(K,J) + SR * ORTR(K) - SI * ORTI(K)
3274               ZI(K,J) = ZI(K,J) + SR * ORTI(K) + SI * ORTR(K)
3275  120       CONTINUE
3276C
3277  130    CONTINUE
3278C
3279  140 CONTINUE
3280C     .......... CREATE REAL SUBDIAGONAL ELEMENTS ..........
3281  150 L = LOW + 1
3282C
3283      DO 170 I = L, IGH
3284         LL = MIN0(I+1,IGH)
3285         IF (HI(I,I-1) .EQ. 0.0D0) GO TO 170
3286         NORM = PYTHAG(HR(I,I-1),HI(I,I-1))
3287         YR = HR(I,I-1) / NORM
3288         YI = HI(I,I-1) / NORM
3289         HR(I,I-1) = NORM
3290         HI(I,I-1) = 0.0D0
3291C
3292         DO 155 J = I, N
3293            SI = YR * HI(I,J) - YI * HR(I,J)
3294            HR(I,J) = YR * HR(I,J) + YI * HI(I,J)
3295            HI(I,J) = SI
3296  155    CONTINUE
3297C
3298         DO 160 J = 1, LL
3299            SI = YR * HI(J,I) + YI * HR(J,I)
3300            HR(J,I) = YR * HR(J,I) - YI * HI(J,I)
3301            HI(J,I) = SI
3302  160    CONTINUE
3303C
3304         DO 165 J = LOW, IGH
3305            SI = YR * ZI(J,I) + YI * ZR(J,I)
3306            ZR(J,I) = YR * ZR(J,I) - YI * ZI(J,I)
3307            ZI(J,I) = SI
3308  165    CONTINUE
3309C
3310  170 CONTINUE
3311C     .......... STORE ROOTS ISOLATED BY CBAL ..........
3312  180 DO 200 I = 1, N
3313         IF (I .GE. LOW .AND. I .LE. IGH) GO TO 200
3314         WR(I) = HR(I,I)
3315         WI(I) = HI(I,I)
3316  200 CONTINUE
3317C
3318      EN = IGH
3319      TR = 0.0D0
3320      TI = 0.0D0
3321      ITN = 30*N
3322C     .......... SEARCH FOR NEXT EIGENVALUE ..........
3323  220 IF (EN .LT. LOW) GO TO 680
3324      ITS = 0
3325      ENM1 = EN - 1
3326C     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
3327C                FOR L=EN STEP -1 UNTIL LOW DO -- ..........
3328  240 DO 260 LL = LOW, EN
3329         L = EN + LOW - LL
3330         IF (L .EQ. LOW) GO TO 300
3331         TST1 = DABS(HR(L-1,L-1)) + DABS(HI(L-1,L-1))
3332     X            + DABS(HR(L,L)) + DABS(HI(L,L))
3333         TST2 = TST1 + DABS(HR(L,L-1))
3334         IF (TST2 .EQ. TST1) GO TO 300
3335  260 CONTINUE
3336C     .......... FORM SHIFT ..........
3337  300 IF (L .EQ. EN) GO TO 660
3338      IF (ITN .EQ. 0) GO TO 1000
3339      IF (ITS .EQ. 10 .OR. ITS .EQ. 20) GO TO 320
3340      SR = HR(EN,EN)
3341      SI = HI(EN,EN)
3342      XR = HR(ENM1,EN) * HR(EN,ENM1)
3343      XI = HI(ENM1,EN) * HR(EN,ENM1)
3344      IF (XR .EQ. 0.0D0 .AND. XI .EQ. 0.0D0) GO TO 340
3345      YR = (HR(ENM1,ENM1) - SR) / 2.0D0
3346      YI = (HI(ENM1,ENM1) - SI) / 2.0D0
3347      CALL CSROOT(YR**2-YI**2+XR,2.0D0*YR*YI+XI,ZZR,ZZI)
3348      IF (YR * ZZR + YI * ZZI .GE. 0.0D0) GO TO 310
3349      ZZR = -ZZR
3350      ZZI = -ZZI
3351  310 CALL CDIV(XR,XI,YR+ZZR,YI+ZZI,XR,XI)
3352      SR = SR - XR
3353      SI = SI - XI
3354      GO TO 340
3355C     .......... FORM EXCEPTIONAL SHIFT ..........
3356  320 SR = DABS(HR(EN,ENM1)) + DABS(HR(ENM1,EN-2))
3357      SI = 0.0D0
3358C
3359  340 DO 360 I = LOW, EN
3360         HR(I,I) = HR(I,I) - SR
3361         HI(I,I) = HI(I,I) - SI
3362  360 CONTINUE
3363C
3364      TR = TR + SR
3365      TI = TI + SI
3366      ITS = ITS + 1
3367      ITN = ITN - 1
3368C     .......... REDUCE TO TRIANGLE (ROWS) ..........
3369      LP1 = L + 1
3370C
3371      DO 500 I = LP1, EN
3372         SR = HR(I,I-1)
3373         HR(I,I-1) = 0.0D0
3374         NORM = PYTHAG(PYTHAG(HR(I-1,I-1),HI(I-1,I-1)),SR)
3375         XR = HR(I-1,I-1) / NORM
3376         WR(I-1) = XR
3377         XI = HI(I-1,I-1) / NORM
3378         WI(I-1) = XI
3379         HR(I-1,I-1) = NORM
3380         HI(I-1,I-1) = 0.0D0
3381         HI(I,I-1) = SR / NORM
3382C
3383         DO 490 J = I, N
3384            YR = HR(I-1,J)
3385            YI = HI(I-1,J)
3386            ZZR = HR(I,J)
3387            ZZI = HI(I,J)
3388            HR(I-1,J) = XR * YR + XI * YI + HI(I,I-1) * ZZR
3389            HI(I-1,J) = XR * YI - XI * YR + HI(I,I-1) * ZZI
3390            HR(I,J) = XR * ZZR - XI * ZZI - HI(I,I-1) * YR
3391            HI(I,J) = XR * ZZI + XI * ZZR - HI(I,I-1) * YI
3392  490    CONTINUE
3393C
3394  500 CONTINUE
3395C
3396      SI = HI(EN,EN)
3397      IF (SI .EQ. 0.0D0) GO TO 540
3398      NORM = PYTHAG(HR(EN,EN),SI)
3399      SR = HR(EN,EN) / NORM
3400      SI = SI / NORM
3401      HR(EN,EN) = NORM
3402      HI(EN,EN) = 0.0D0
3403      IF (EN .EQ. N) GO TO 540
3404      IP1 = EN + 1
3405C
3406      DO 520 J = IP1, N
3407         YR = HR(EN,J)
3408         YI = HI(EN,J)
3409         HR(EN,J) = SR * YR + SI * YI
3410         HI(EN,J) = SR * YI - SI * YR
3411  520 CONTINUE
3412C     .......... INVERSE OPERATION (COLUMNS) ..........
3413  540 DO 600 J = LP1, EN
3414         XR = WR(J-1)
3415         XI = WI(J-1)
3416C
3417         DO 580 I = 1, J
3418            YR = HR(I,J-1)
3419            YI = 0.0D0
3420            ZZR = HR(I,J)
3421            ZZI = HI(I,J)
3422            IF (I .EQ. J) GO TO 560
3423            YI = HI(I,J-1)
3424            HI(I,J-1) = XR * YI + XI * YR + HI(J,J-1) * ZZI
3425  560       HR(I,J-1) = XR * YR - XI * YI + HI(J,J-1) * ZZR
3426            HR(I,J) = XR * ZZR + XI * ZZI - HI(J,J-1) * YR
3427            HI(I,J) = XR * ZZI - XI * ZZR - HI(J,J-1) * YI
3428  580    CONTINUE
3429C
3430         DO 590 I = LOW, IGH
3431            YR = ZR(I,J-1)
3432            YI = ZI(I,J-1)
3433            ZZR = ZR(I,J)
3434            ZZI = ZI(I,J)
3435            ZR(I,J-1) = XR * YR - XI * YI + HI(J,J-1) * ZZR
3436            ZI(I,J-1) = XR * YI + XI * YR + HI(J,J-1) * ZZI
3437            ZR(I,J) = XR * ZZR + XI * ZZI - HI(J,J-1) * YR
3438            ZI(I,J) = XR * ZZI - XI * ZZR - HI(J,J-1) * YI
3439  590    CONTINUE
3440C
3441  600 CONTINUE
3442C
3443      IF (SI .EQ. 0.0D0) GO TO 240
3444C
3445      DO 630 I = 1, EN
3446         YR = HR(I,EN)
3447         YI = HI(I,EN)
3448         HR(I,EN) = SR * YR - SI * YI
3449         HI(I,EN) = SR * YI + SI * YR
3450  630 CONTINUE
3451C
3452      DO 640 I = LOW, IGH
3453         YR = ZR(I,EN)
3454         YI = ZI(I,EN)
3455         ZR(I,EN) = SR * YR - SI * YI
3456         ZI(I,EN) = SR * YI + SI * YR
3457  640 CONTINUE
3458C
3459      GO TO 240
3460C     .......... A ROOT FOUND ..........
3461  660 HR(EN,EN) = HR(EN,EN) + TR
3462      WR(EN) = HR(EN,EN)
3463      HI(EN,EN) = HI(EN,EN) + TI
3464      WI(EN) = HI(EN,EN)
3465      EN = ENM1
3466      GO TO 220
3467C     .......... ALL ROOTS FOUND.  BACKSUBSTITUTE TO FIND
3468C                VECTORS OF UPPER TRIANGULAR FORM ..........
3469  680 NORM = 0.0D0
3470C
3471      DO 720 I = 1, N
3472C
3473         DO 720 J = I, N
3474            TR = DABS(HR(I,J)) + DABS(HI(I,J))
3475            IF (TR .GT. NORM) NORM = TR
3476  720 CONTINUE
3477C
3478      IF (N .EQ. 1 .OR. NORM .EQ. 0.0D0) GO TO 1001
3479C     .......... FOR EN=N STEP -1 UNTIL 2 DO -- ..........
3480      DO 800 NN = 2, N
3481         EN = N + 2 - NN
3482         XR = WR(EN)
3483         XI = WI(EN)
3484         HR(EN,EN) = 1.0D0
3485         HI(EN,EN) = 0.0D0
3486         ENM1 = EN - 1
3487C     .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- ..........
3488         DO 780 II = 1, ENM1
3489            I = EN - II
3490            ZZR = 0.0D0
3491            ZZI = 0.0D0
3492            IP1 = I + 1
3493C
3494            DO 740 J = IP1, EN
3495               ZZR = ZZR + HR(I,J) * HR(J,EN) - HI(I,J) * HI(J,EN)
3496               ZZI = ZZI + HR(I,J) * HI(J,EN) + HI(I,J) * HR(J,EN)
3497  740       CONTINUE
3498C
3499            YR = XR - WR(I)
3500            YI = XI - WI(I)
3501            IF (YR .NE. 0.0D0 .OR. YI .NE. 0.0D0) GO TO 765
3502               TST1 = NORM
3503               YR = TST1
3504  760          YR = 0.01D0 * YR
3505               TST2 = NORM + YR
3506               IF (TST2 .GT. TST1) GO TO 760
3507  765       CONTINUE
3508            CALL CDIV(ZZR,ZZI,YR,YI,HR(I,EN),HI(I,EN))
3509C     .......... OVERFLOW CONTROL ..........
3510            TR = DABS(HR(I,EN)) + DABS(HI(I,EN))
3511            IF (TR .EQ. 0.0D0) GO TO 780
3512            TST1 = TR
3513            TST2 = TST1 + 1.0D0/TST1
3514            IF (TST2 .GT. TST1) GO TO 780
3515            DO 770 J = I, EN
3516               HR(J,EN) = HR(J,EN)/TR
3517               HI(J,EN) = HI(J,EN)/TR
3518  770       CONTINUE
3519C
3520  780    CONTINUE
3521C
3522  800 CONTINUE
3523C     .......... END BACKSUBSTITUTION ..........
3524      ENM1 = N - 1
3525C     .......... VECTORS OF ISOLATED ROOTS ..........
3526      DO  840 I = 1, ENM1
3527         IF (I .GE. LOW .AND. I .LE. IGH) GO TO 840
3528         IP1 = I + 1
3529C
3530         DO 820 J = IP1, N
3531            ZR(I,J) = HR(I,J)
3532            ZI(I,J) = HI(I,J)
3533  820    CONTINUE
3534C
3535  840 CONTINUE
3536C     .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE
3537C                VECTORS OF ORIGINAL FULL MATRIX.
3538C                FOR J=N STEP -1 UNTIL LOW+1 DO -- ..........
3539      DO 880 JJ = LOW, ENM1
3540         J = N + LOW - JJ
3541         M = MIN0(J,IGH)
3542C
3543         DO 880 I = LOW, IGH
3544            ZZR = 0.0D0
3545            ZZI = 0.0D0
3546C
3547            DO 860 K = LOW, M
3548               ZZR = ZZR + ZR(I,K) * HR(K,J) - ZI(I,K) * HI(K,J)
3549               ZZI = ZZI + ZR(I,K) * HI(K,J) + ZI(I,K) * HR(K,J)
3550  860       CONTINUE
3551C
3552            ZR(I,J) = ZZR
3553            ZI(I,J) = ZZI
3554  880 CONTINUE
3555C
3556      GO TO 1001
3557C     .......... SET ERROR -- ALL EIGENVALUES HAVE NOT
3558C                CONVERGED AFTER 30*N ITERATIONS ..........
3559 1000 IERR = EN
3560 1001 RETURN
3561      END
3562      SUBROUTINE CORTB(NM,LOW,IGH,AR,AI,ORTR,ORTI,M,ZR,ZI)
3563C
3564      INTEGER I,J,M,LA,MM,MP,NM,IGH,KP1,LOW,MP1
3565      DOUBLE PRECISION AR(NM,IGH),AI(NM,IGH),ORTR(IGH),ORTI(IGH),
3566     X       ZR(NM,M),ZI(NM,M)
3567      DOUBLE PRECISION H,GI,GR
3568C
3569C     THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF
3570C     THE ALGOL PROCEDURE ORTBAK, NUM. MATH. 12, 349-368(1968)
3571C     BY MARTIN AND WILKINSON.
3572C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
3573C
3574C     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX GENERAL
3575C     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
3576C     UPPER HESSENBERG MATRIX DETERMINED BY  CORTH.
3577C
3578C     ON INPUT
3579C
3580C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
3581C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
3582C          DIMENSION STATEMENT.
3583C
3584C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
3585C          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED,
3586C          SET LOW=1 AND IGH EQUAL TO THE ORDER OF THE MATRIX.
3587C
3588C        AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY
3589C          TRANSFORMATIONS USED IN THE REDUCTION BY  CORTH
3590C          IN THEIR STRICT LOWER TRIANGLES.
3591C
3592C        ORTR AND ORTI CONTAIN FURTHER INFORMATION ABOUT THE
3593C          TRANSFORMATIONS USED IN THE REDUCTION BY  CORTH.
3594C          ONLY ELEMENTS LOW THROUGH IGH ARE USED.
3595C
3596C        M IS THE NUMBER OF COLUMNS OF ZR AND ZI TO BE BACK TRANSFORMED.
3597C
3598C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
3599C          RESPECTIVELY, OF THE EIGENVECTORS TO BE
3600C          BACK TRANSFORMED IN THEIR FIRST M COLUMNS.
3601C
3602C     ON OUTPUT
3603C
3604C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
3605C          RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS
3606C          IN THEIR FIRST M COLUMNS.
3607C
3608C        ORTR AND ORTI HAVE BEEN ALTERED.
3609C
3610C     NOTE THAT CORTB PRESERVES VECTOR EUCLIDEAN NORMS.
3611C
3612C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
3613C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
3614C
3615C     THIS VERSION DATED AUGUST 1983.
3616C
3617C     ------------------------------------------------------------------
3618C
3619      IF (M .EQ. 0) GO TO 200
3620      LA = IGH - 1
3621      KP1 = LOW + 1
3622      IF (LA .LT. KP1) GO TO 200
3623C     .......... FOR MP=IGH-1 STEP -1 UNTIL LOW+1 DO -- ..........
3624      DO 140 MM = KP1, LA
3625         MP = LOW + IGH - MM
3626         IF (AR(MP,MP-1) .EQ. 0.0D0 .AND. AI(MP,MP-1) .EQ. 0.0D0)
3627     X      GO TO 140
3628C     .......... H BELOW IS NEGATIVE OF H FORMED IN CORTH ..........
3629         H = AR(MP,MP-1) * ORTR(MP) + AI(MP,MP-1) * ORTI(MP)
3630         MP1 = MP + 1
3631C
3632         DO 100 I = MP1, IGH
3633            ORTR(I) = AR(I,MP-1)
3634            ORTI(I) = AI(I,MP-1)
3635  100    CONTINUE
3636C
3637         DO 130 J = 1, M
3638            GR = 0.0D0
3639            GI = 0.0D0
3640C
3641            DO 110 I = MP, IGH
3642               GR = GR + ORTR(I) * ZR(I,J) + ORTI(I) * ZI(I,J)
3643               GI = GI + ORTR(I) * ZI(I,J) - ORTI(I) * ZR(I,J)
3644  110       CONTINUE
3645C
3646            GR = GR / H
3647            GI = GI / H
3648C
3649            DO 120 I = MP, IGH
3650               ZR(I,J) = ZR(I,J) + GR * ORTR(I) - GI * ORTI(I)
3651               ZI(I,J) = ZI(I,J) + GR * ORTI(I) + GI * ORTR(I)
3652  120       CONTINUE
3653C
3654  130    CONTINUE
3655C
3656  140 CONTINUE
3657C
3658  200 RETURN
3659      END
3660      SUBROUTINE CORTH(NM,N,LOW,IGH,AR,AI,ORTR,ORTI)
3661C
3662      INTEGER I,J,M,N,II,JJ,LA,MP,NM,IGH,KP1,LOW
3663      DOUBLE PRECISION AR(NM,N),AI(NM,N),ORTR(IGH),ORTI(IGH)
3664      DOUBLE PRECISION F,G,H,FI,FR,SCALE,PYTHAG
3665C
3666C     THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF
3667C     THE ALGOL PROCEDURE ORTHES, NUM. MATH. 12, 349-368(1968)
3668C     BY MARTIN AND WILKINSON.
3669C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
3670C
3671C     GIVEN A COMPLEX GENERAL MATRIX, THIS SUBROUTINE
3672C     REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS
3673C     LOW THROUGH IGH TO UPPER HESSENBERG FORM BY
3674C     UNITARY SIMILARITY TRANSFORMATIONS.
3675C
3676C     ON INPUT
3677C
3678C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
3679C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
3680C          DIMENSION STATEMENT.
3681C
3682C        N IS THE ORDER OF THE MATRIX.
3683C
3684C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
3685C          SUBROUTINE  CBAL.  IF  CBAL  HAS NOT BEEN USED,
3686C          SET LOW=1, IGH=N.
3687C
3688C        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
3689C          RESPECTIVELY, OF THE COMPLEX INPUT MATRIX.
3690C
3691C     ON OUTPUT
3692C
3693C        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
3694C          RESPECTIVELY, OF THE HESSENBERG MATRIX.  INFORMATION
3695C          ABOUT THE UNITARY TRANSFORMATIONS USED IN THE REDUCTION
3696C          IS STORED IN THE REMAINING TRIANGLES UNDER THE
3697C          HESSENBERG MATRIX.
3698C
3699C        ORTR AND ORTI CONTAIN FURTHER INFORMATION ABOUT THE
3700C          TRANSFORMATIONS.  ONLY ELEMENTS LOW THROUGH IGH ARE USED.
3701C
3702C     CALLS PYTHAG FOR  DSQRT(A*A + B*B) .
3703C
3704C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
3705C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
3706C
3707C     THIS VERSION DATED AUGUST 1983.
3708C
3709C     ------------------------------------------------------------------
3710C
3711      LA = IGH - 1
3712      KP1 = LOW + 1
3713      IF (LA .LT. KP1) GO TO 200
3714C
3715      DO 180 M = KP1, LA
3716         H = 0.0D0
3717         ORTR(M) = 0.0D0
3718         ORTI(M) = 0.0D0
3719         SCALE = 0.0D0
3720C     .......... SCALE COLUMN (ALGOL TOL THEN NOT NEEDED) ..........
3721         DO 90 I = M, IGH
3722   90    SCALE = SCALE + DABS(AR(I,M-1)) + DABS(AI(I,M-1))
3723C
3724         IF (SCALE .EQ. 0.0D0) GO TO 180
3725         MP = M + IGH
3726C     .......... FOR I=IGH STEP -1 UNTIL M DO -- ..........
3727         DO 100 II = M, IGH
3728            I = MP - II
3729            ORTR(I) = AR(I,M-1) / SCALE
3730            ORTI(I) = AI(I,M-1) / SCALE
3731            H = H + ORTR(I) * ORTR(I) + ORTI(I) * ORTI(I)
3732  100    CONTINUE
3733C
3734         G = DSQRT(H)
3735         F = PYTHAG(ORTR(M),ORTI(M))
3736         IF (F .EQ. 0.0D0) GO TO 103
3737         H = H + F * G
3738         G = G / F
3739         ORTR(M) = (1.0D0 + G) * ORTR(M)
3740         ORTI(M) = (1.0D0 + G) * ORTI(M)
3741         GO TO 105
3742C
3743  103    ORTR(M) = G
3744         AR(M,M-1) = SCALE
3745C     .......... FORM (I-(U*UT)/H) * A ..........
3746  105    DO 130 J = M, N
3747            FR = 0.0D0
3748            FI = 0.0D0
3749C     .......... FOR I=IGH STEP -1 UNTIL M DO -- ..........
3750            DO 110 II = M, IGH
3751               I = MP - II
3752               FR = FR + ORTR(I) * AR(I,J) + ORTI(I) * AI(I,J)
3753               FI = FI + ORTR(I) * AI(I,J) - ORTI(I) * AR(I,J)
3754  110       CONTINUE
3755C
3756            FR = FR / H
3757            FI = FI / H
3758C
3759            DO 120 I = M, IGH
3760               AR(I,J) = AR(I,J) - FR * ORTR(I) + FI * ORTI(I)
3761               AI(I,J) = AI(I,J) - FR * ORTI(I) - FI * ORTR(I)
3762  120       CONTINUE
3763C
3764  130    CONTINUE
3765C     .......... FORM (I-(U*UT)/H)*A*(I-(U*UT)/H) ..........
3766         DO 160 I = 1, IGH
3767            FR = 0.0D0
3768            FI = 0.0D0
3769C     .......... FOR J=IGH STEP -1 UNTIL M DO -- ..........
3770            DO 140 JJ = M, IGH
3771               J = MP - JJ
3772               FR = FR + ORTR(J) * AR(I,J) - ORTI(J) * AI(I,J)
3773               FI = FI + ORTR(J) * AI(I,J) + ORTI(J) * AR(I,J)
3774  140       CONTINUE
3775C
3776            FR = FR / H
3777            FI = FI / H
3778C
3779            DO 150 J = M, IGH
3780               AR(I,J) = AR(I,J) - FR * ORTR(J) - FI * ORTI(J)
3781               AI(I,J) = AI(I,J) + FR * ORTI(J) - FI * ORTR(J)
3782  150       CONTINUE
3783C
3784  160    CONTINUE
3785C
3786         ORTR(M) = SCALE * ORTR(M)
3787         ORTI(M) = SCALE * ORTI(M)
3788         AR(M,M-1) = -G * AR(M,M-1)
3789         AI(M,M-1) = -G * AI(M,M-1)
3790  180 CONTINUE
3791C
3792  200 RETURN
3793      END
3794      SUBROUTINE ELMBAK(NM,LOW,IGH,A,INT,M,Z)
3795C
3796      INTEGER I,J,M,LA,MM,MP,NM,IGH,KP1,LOW,MP1
3797      DOUBLE PRECISION A(NM,IGH),Z(NM,M)
3798      DOUBLE PRECISION X
3799      INTEGER INT(IGH)
3800C
3801C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMBAK,
3802C     NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON.
3803C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
3804C
3805C     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A REAL GENERAL
3806C     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
3807C     UPPER HESSENBERG MATRIX DETERMINED BY  ELMHES.
3808C
3809C     ON INPUT
3810C
3811C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
3812C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
3813C          DIMENSION STATEMENT.
3814C
3815C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
3816C          SUBROUTINE  BALANC.  IF  BALANC  HAS NOT BEEN USED,
3817C          SET LOW=1 AND IGH EQUAL TO THE ORDER OF THE MATRIX.
3818C
3819C        A CONTAINS THE MULTIPLIERS WHICH WERE USED IN THE
3820C          REDUCTION BY  ELMHES  IN ITS LOWER TRIANGLE
3821C          BELOW THE SUBDIAGONAL.
3822C
3823C        INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS
3824C          INTERCHANGED IN THE REDUCTION BY  ELMHES.
3825C          ONLY ELEMENTS LOW THROUGH IGH ARE USED.
3826C
3827C        M IS THE NUMBER OF COLUMNS OF Z TO BE BACK TRANSFORMED.
3828C
3829C        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGEN-
3830C          VECTORS TO BE BACK TRANSFORMED IN ITS FIRST M COLUMNS.
3831C
3832C     ON OUTPUT
3833C
3834C        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE
3835C          TRANSFORMED EIGENVECTORS IN ITS FIRST M COLUMNS.
3836C
3837C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
3838C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
3839C
3840C     THIS VERSION DATED AUGUST 1983.
3841C
3842C     ------------------------------------------------------------------
3843C
3844      IF (M .EQ. 0) GO TO 200
3845      LA = IGH - 1
3846      KP1 = LOW + 1
3847      IF (LA .LT. KP1) GO TO 200
3848C     .......... FOR MP=IGH-1 STEP -1 UNTIL LOW+1 DO -- ..........
3849      DO 140 MM = KP1, LA
3850         MP = LOW + IGH - MM
3851         MP1 = MP + 1
3852C
3853         DO 110 I = MP1, IGH
3854            X = A(I,MP-1)
3855            IF (X .EQ. 0.0D0) GO TO 110
3856C
3857            DO 100 J = 1, M
3858  100       Z(I,J) = Z(I,J) + X * Z(MP,J)
3859C
3860  110    CONTINUE
3861C
3862         I = INT(MP)
3863         IF (I .EQ. MP) GO TO 140
3864C
3865         DO 130 J = 1, M
3866            X = Z(I,J)
3867            Z(I,J) = Z(MP,J)
3868            Z(MP,J) = X
3869  130    CONTINUE
3870C
3871  140 CONTINUE
3872C
3873  200 RETURN
3874      END
3875      SUBROUTINE ELMHES(NM,N,LOW,IGH,A,INT)
3876C
3877      INTEGER I,J,M,N,LA,NM,IGH,KP1,LOW,MM1,MP1
3878      DOUBLE PRECISION A(NM,N)
3879      DOUBLE PRECISION X,Y
3880      INTEGER INT(IGH)
3881C
3882C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMHES,
3883C     NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON.
3884C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
3885C
3886C     GIVEN A REAL GENERAL MATRIX, THIS SUBROUTINE
3887C     REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS
3888C     LOW THROUGH IGH TO UPPER HESSENBERG FORM BY
3889C     STABILIZED ELEMENTARY SIMILARITY TRANSFORMATIONS.
3890C
3891C     ON INPUT
3892C
3893C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
3894C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
3895C          DIMENSION STATEMENT.
3896C
3897C        N IS THE ORDER OF THE MATRIX.
3898C
3899C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
3900C          SUBROUTINE  BALANC.  IF  BALANC  HAS NOT BEEN USED,
3901C          SET LOW=1, IGH=N.
3902C
3903C        A CONTAINS THE INPUT MATRIX.
3904C
3905C     ON OUTPUT
3906C
3907C        A CONTAINS THE HESSENBERG MATRIX.  THE MULTIPLIERS
3908C          WHICH WERE USED IN THE REDUCTION ARE STORED IN THE
3909C          REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX.
3910C
3911C        INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS
3912C          INTERCHANGED IN THE REDUCTION.
3913C          ONLY ELEMENTS LOW THROUGH IGH ARE USED.
3914C
3915C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
3916C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
3917C
3918C     THIS VERSION DATED AUGUST 1983.
3919C
3920C     ------------------------------------------------------------------
3921C
3922      LA = IGH - 1
3923      KP1 = LOW + 1
3924      IF (LA .LT. KP1) GO TO 200
3925C
3926      DO 180 M = KP1, LA
3927         MM1 = M - 1
3928         X = 0.0D0
3929         I = M
3930C
3931         DO 100 J = M, IGH
3932            IF (DABS(A(J,MM1)) .LE. DABS(X)) GO TO 100
3933            X = A(J,MM1)
3934            I = J
3935  100    CONTINUE
3936C
3937         INT(M) = I
3938         IF (I .EQ. M) GO TO 130
3939C     .......... INTERCHANGE ROWS AND COLUMNS OF A ..........
3940         DO 110 J = MM1, N
3941            Y = A(I,J)
3942            A(I,J) = A(M,J)
3943            A(M,J) = Y
3944  110    CONTINUE
3945C
3946         DO 120 J = 1, IGH
3947            Y = A(J,I)
3948            A(J,I) = A(J,M)
3949            A(J,M) = Y
3950  120    CONTINUE
3951C     .......... END INTERCHANGE ..........
3952  130    IF (X .EQ. 0.0D0) GO TO 180
3953         MP1 = M + 1
3954C
3955         DO 160 I = MP1, IGH
3956            Y = A(I,MM1)
3957            IF (Y .EQ. 0.0D0) GO TO 160
3958            Y = Y / X
3959            A(I,MM1) = Y
3960C
3961            DO 140 J = M, N
3962  140       A(I,J) = A(I,J) - Y * A(M,J)
3963C
3964            DO 150 J = 1, IGH
3965  150       A(J,M) = A(J,M) + Y * A(J,I)
3966C
3967  160    CONTINUE
3968C
3969  180 CONTINUE
3970C
3971  200 RETURN
3972      END
3973      SUBROUTINE ELTRAN(NM,N,LOW,IGH,A,INT,Z)
3974C
3975      INTEGER I,J,N,KL,MM,MP,NM,IGH,LOW,MP1
3976      DOUBLE PRECISION A(NM,IGH),Z(NM,N)
3977      INTEGER INT(IGH)
3978C
3979C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMTRANS,
3980C     NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON.
3981C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
3982C
3983C     THIS SUBROUTINE ACCUMULATES THE STABILIZED ELEMENTARY
3984C     SIMILARITY TRANSFORMATIONS USED IN THE REDUCTION OF A
3985C     REAL GENERAL MATRIX TO UPPER HESSENBERG FORM BY  ELMHES.
3986C
3987C     ON INPUT
3988C
3989C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
3990C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
3991C          DIMENSION STATEMENT.
3992C
3993C        N IS THE ORDER OF THE MATRIX.
3994C
3995C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
3996C          SUBROUTINE  BALANC.  IF  BALANC  HAS NOT BEEN USED,
3997C          SET LOW=1, IGH=N.
3998C
3999C        A CONTAINS THE MULTIPLIERS WHICH WERE USED IN THE
4000C          REDUCTION BY  ELMHES  IN ITS LOWER TRIANGLE
4001C          BELOW THE SUBDIAGONAL.
4002C
4003C        INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS
4004C          INTERCHANGED IN THE REDUCTION BY  ELMHES.
4005C          ONLY ELEMENTS LOW THROUGH IGH ARE USED.
4006C
4007C     ON OUTPUT
4008C
4009C        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE
4010C          REDUCTION BY  ELMHES.
4011C
4012C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
4013C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
4014C
4015C     THIS VERSION DATED AUGUST 1983.
4016C
4017C     ------------------------------------------------------------------
4018C
4019C     .......... INITIALIZE Z TO IDENTITY MATRIX ..........
4020      DO 80 J = 1, N
4021C
4022         DO 60 I = 1, N
4023   60    Z(I,J) = 0.0D0
4024C
4025         Z(J,J) = 1.0D0
4026   80 CONTINUE
4027C
4028      KL = IGH - LOW - 1
4029      IF (KL .LT. 1) GO TO 200
4030C     .......... FOR MP=IGH-1 STEP -1 UNTIL LOW+1 DO -- ..........
4031      DO 140 MM = 1, KL
4032         MP = IGH - MM
4033         MP1 = MP + 1
4034C
4035         DO 100 I = MP1, IGH
4036  100    Z(I,MP) = A(I,MP-1)
4037C
4038         I = INT(MP)
4039         IF (I .EQ. MP) GO TO 140
4040C
4041         DO 130 J = MP, IGH
4042            Z(MP,J) = Z(I,J)
4043            Z(I,J) = 0.0D0
4044  130    CONTINUE
4045C
4046         Z(I,MP) = 1.0D0
4047  140 CONTINUE
4048C
4049  200 RETURN
4050      END
4051      SUBROUTINE FIGI(NM,N,T,D,E,E2,IERR)
4052C
4053      INTEGER I,N,NM,IERR
4054      DOUBLE PRECISION T(NM,3),D(N),E(N),E2(N)
4055C
4056C     GIVEN A NONSYMMETRIC TRIDIAGONAL MATRIX SUCH THAT THE PRODUCTS
4057C     OF CORRESPONDING PAIRS OF OFF-DIAGONAL ELEMENTS ARE ALL
4058C     NON-NEGATIVE, THIS SUBROUTINE REDUCES IT TO A SYMMETRIC
4059C     TRIDIAGONAL MATRIX WITH THE SAME EIGENVALUES.  IF, FURTHER,
4060C     A ZERO PRODUCT ONLY OCCURS WHEN BOTH FACTORS ARE ZERO,
4061C     THE REDUCED MATRIX IS SIMILAR TO THE ORIGINAL MATRIX.
4062C
4063C     ON INPUT
4064C
4065C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
4066C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
4067C          DIMENSION STATEMENT.
4068C
4069C        N IS THE ORDER OF THE MATRIX.
4070C
4071C        T CONTAINS THE INPUT MATRIX.  ITS SUBDIAGONAL IS
4072C          STORED IN THE LAST N-1 POSITIONS OF THE FIRST COLUMN,
4073C          ITS DIAGONAL IN THE N POSITIONS OF THE SECOND COLUMN,
4074C          AND ITS SUPERDIAGONAL IN THE FIRST N-1 POSITIONS OF
4075C          THE THIRD COLUMN.  T(1,1) AND T(N,3) ARE ARBITRARY.
4076C
4077C     ON OUTPUT
4078C
4079C        T IS UNALTERED.
4080C
4081C        D CONTAINS THE DIAGONAL ELEMENTS OF THE SYMMETRIC MATRIX.
4082C
4083C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE SYMMETRIC
4084C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS NOT SET.
4085C
4086C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
4087C          E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
4088C
4089C        IERR IS SET TO
4090C          ZERO       FOR NORMAL RETURN,
4091C          N+I        IF T(I,1)*T(I-1,3) IS NEGATIVE,
4092C          -(3*N+I)   IF T(I,1)*T(I-1,3) IS ZERO WITH ONE FACTOR
4093C                     NON-ZERO.  IN THIS CASE, THE EIGENVECTORS OF
4094C                     THE SYMMETRIC MATRIX ARE NOT SIMPLY RELATED
4095C                     TO THOSE OF  T  AND SHOULD NOT BE SOUGHT.
4096C
4097C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
4098C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
4099C
4100C     THIS VERSION DATED AUGUST 1983.
4101C
4102C     ------------------------------------------------------------------
4103C
4104      IERR = 0
4105C
4106      DO 100 I = 1, N
4107         IF (I .EQ. 1) GO TO 90
4108         E2(I) = T(I,1) * T(I-1,3)
4109         IF (E2(I)) 1000, 60, 80
4110   60    IF (T(I,1) .EQ. 0.0D0 .AND. T(I-1,3) .EQ. 0.0D0) GO TO 80
4111C     .......... SET ERROR -- PRODUCT OF SOME PAIR OF OFF-DIAGONAL
4112C                ELEMENTS IS ZERO WITH ONE MEMBER NON-ZERO ..........
4113         IERR = -(3 * N + I)
4114   80    E(I) = DSQRT(E2(I))
4115   90    D(I) = T(I,2)
4116  100 CONTINUE
4117C
4118      GO TO 1001
4119C     .......... SET ERROR -- PRODUCT OF SOME PAIR OF OFF-DIAGONAL
4120C                ELEMENTS IS NEGATIVE ..........
4121 1000 IERR = N + I
4122 1001 RETURN
4123      END
4124      SUBROUTINE FIGI2(NM,N,T,D,E,Z,IERR)
4125C
4126      INTEGER I,J,N,NM,IERR
4127      DOUBLE PRECISION T(NM,3),D(N),E(N),Z(NM,N)
4128      DOUBLE PRECISION H
4129C
4130C     GIVEN A NONSYMMETRIC TRIDIAGONAL MATRIX SUCH THAT THE PRODUCTS
4131C     OF CORRESPONDING PAIRS OF OFF-DIAGONAL ELEMENTS ARE ALL
4132C     NON-NEGATIVE, AND ZERO ONLY WHEN BOTH FACTORS ARE ZERO, THIS
4133C     SUBROUTINE REDUCES IT TO A SYMMETRIC TRIDIAGONAL MATRIX
4134C     USING AND ACCUMULATING DIAGONAL SIMILARITY TRANSFORMATIONS.
4135C
4136C     ON INPUT
4137C
4138C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
4139C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
4140C          DIMENSION STATEMENT.
4141C
4142C        N IS THE ORDER OF THE MATRIX.
4143C
4144C        T CONTAINS THE INPUT MATRIX.  ITS SUBDIAGONAL IS
4145C          STORED IN THE LAST N-1 POSITIONS OF THE FIRST COLUMN,
4146C          ITS DIAGONAL IN THE N POSITIONS OF THE SECOND COLUMN,
4147C          AND ITS SUPERDIAGONAL IN THE FIRST N-1 POSITIONS OF
4148C          THE THIRD COLUMN.  T(1,1) AND T(N,3) ARE ARBITRARY.
4149C
4150C     ON OUTPUT
4151C
4152C        T IS UNALTERED.
4153C
4154C        D CONTAINS THE DIAGONAL ELEMENTS OF THE SYMMETRIC MATRIX.
4155C
4156C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE SYMMETRIC
4157C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS NOT SET.
4158C
4159C        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN
4160C          THE REDUCTION.
4161C
4162C        IERR IS SET TO
4163C          ZERO       FOR NORMAL RETURN,
4164C          N+I        IF T(I,1)*T(I-1,3) IS NEGATIVE,
4165C          2*N+I      IF T(I,1)*T(I-1,3) IS ZERO WITH
4166C                     ONE FACTOR NON-ZERO.
4167C
4168C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
4169C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
4170C
4171C     THIS VERSION DATED AUGUST 1983.
4172C
4173C     ------------------------------------------------------------------
4174C
4175      IERR = 0
4176C
4177      DO 100 I = 1, N
4178C
4179         DO 50 J = 1, N
4180   50    Z(I,J) = 0.0D0
4181C
4182         IF (I .EQ. 1) GO TO 70
4183         H = T(I,1) * T(I-1,3)
4184         IF (H) 900, 60, 80
4185   60    IF (T(I,1) .NE. 0.0D0 .OR. T(I-1,3) .NE. 0.0D0) GO TO 1000
4186         E(I) = 0.0D0
4187   70    Z(I,I) = 1.0D0
4188         GO TO 90
4189   80    E(I) = DSQRT(H)
4190         Z(I,I) = Z(I-1,I-1) * E(I) / T(I-1,3)
4191   90    D(I) = T(I,2)
4192  100 CONTINUE
4193C
4194      GO TO 1001
4195C     .......... SET ERROR -- PRODUCT OF SOME PAIR OF OFF-DIAGONAL
4196C                ELEMENTS IS NEGATIVE ..........
4197  900 IERR = N + I
4198      GO TO 1001
4199C     .......... SET ERROR -- PRODUCT OF SOME PAIR OF OFF-DIAGONAL
4200C                ELEMENTS IS ZERO WITH ONE MEMBER NON-ZERO ..........
4201 1000 IERR = 2 * N + I
4202 1001 RETURN
4203      END
4204      SUBROUTINE HQR(NM,N,LOW,IGH,H,WR,WI,IERR)
4205C
4206      INTEGER I,J,K,L,M,N,EN,LL,MM,NA,NM,IGH,ITN,ITS,LOW,MP2,ENM2,IERR
4207      DOUBLE PRECISION H(NM,N),WR(N),WI(N)
4208      DOUBLE PRECISION P,Q,R,S,T,W,X,Y,ZZ,NORM,TST1,TST2
4209      LOGICAL NOTLAS
4210C
4211C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR,
4212C     NUM. MATH. 14, 219-231(1970) BY MARTIN, PETERS, AND WILKINSON.
4213C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 359-371(1971).
4214C
4215C     THIS SUBROUTINE FINDS THE EIGENVALUES OF A REAL
4216C     UPPER HESSENBERG MATRIX BY THE QR METHOD.
4217C
4218C     ON INPUT
4219C
4220C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
4221C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
4222C          DIMENSION STATEMENT.
4223C
4224C        N IS THE ORDER OF THE MATRIX.
4225C
4226C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
4227C          SUBROUTINE  BALANC.  IF  BALANC  HAS NOT BEEN USED,
4228C          SET LOW=1, IGH=N.
4229C
4230C        H CONTAINS THE UPPER HESSENBERG MATRIX.  INFORMATION ABOUT
4231C          THE TRANSFORMATIONS USED IN THE REDUCTION TO HESSENBERG
4232C          FORM BY  ELMHES  OR  ORTHES, IF PERFORMED, IS STORED
4233C          IN THE REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX.
4234C
4235C     ON OUTPUT
4236C
4237C        H HAS BEEN DESTROYED.  THEREFORE, IT MUST BE SAVED
4238C          BEFORE CALLING  HQR  IF SUBSEQUENT CALCULATION AND
4239C          BACK TRANSFORMATION OF EIGENVECTORS IS TO BE PERFORMED.
4240C
4241C        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
4242C          RESPECTIVELY, OF THE EIGENVALUES.  THE EIGENVALUES
4243C          ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS
4244C          OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE
4245C          HAVING THE POSITIVE IMAGINARY PART FIRST.  IF AN
4246C          ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT
4247C          FOR INDICES IERR+1,...,N.
4248C
4249C        IERR IS SET TO
4250C          ZERO       FOR NORMAL RETURN,
4251C          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED
4252C                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT.
4253C
4254C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
4255C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
4256C
4257C     THIS VERSION DATED AUGUST 1983.
4258C
4259C     ------------------------------------------------------------------
4260C
4261      IERR = 0
4262      NORM = 0.0D0
4263      K = 1
4264C     .......... STORE ROOTS ISOLATED BY BALANC
4265C                AND COMPUTE MATRIX NORM ..........
4266      DO 50 I = 1, N
4267C
4268         DO 40 J = K, N
4269   40    NORM = NORM + DABS(H(I,J))
4270C
4271         K = I
4272         IF (I .GE. LOW .AND. I .LE. IGH) GO TO 50
4273         WR(I) = H(I,I)
4274         WI(I) = 0.0D0
4275   50 CONTINUE
4276C
4277      EN = IGH
4278      T = 0.0D0
4279      ITN = 30*N
4280C     .......... SEARCH FOR NEXT EIGENVALUES ..........
4281   60 IF (EN .LT. LOW) GO TO 1001
4282      ITS = 0
4283      NA = EN - 1
4284      ENM2 = NA - 1
4285C     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
4286C                FOR L=EN STEP -1 UNTIL LOW DO -- ..........
4287   70 DO 80 LL = LOW, EN
4288         L = EN + LOW - LL
4289         IF (L .EQ. LOW) GO TO 100
4290         S = DABS(H(L-1,L-1)) + DABS(H(L,L))
4291         IF (S .EQ. 0.0D0) S = NORM
4292         TST1 = S
4293         TST2 = TST1 + DABS(H(L,L-1))
4294         IF (TST2 .EQ. TST1) GO TO 100
4295   80 CONTINUE
4296C     .......... FORM SHIFT ..........
4297  100 X = H(EN,EN)
4298      IF (L .EQ. EN) GO TO 270
4299      Y = H(NA,NA)
4300      W = H(EN,NA) * H(NA,EN)
4301      IF (L .EQ. NA) GO TO 280
4302      IF (ITN .EQ. 0) GO TO 1000
4303      IF (ITS .NE. 10 .AND. ITS .NE. 20) GO TO 130
4304C     .......... FORM EXCEPTIONAL SHIFT ..........
4305      T = T + X
4306C
4307      DO 120 I = LOW, EN
4308  120 H(I,I) = H(I,I) - X
4309C
4310      S = DABS(H(EN,NA)) + DABS(H(NA,ENM2))
4311      X = 0.75D0 * S
4312      Y = X
4313      W = -0.4375D0 * S * S
4314  130 ITS = ITS + 1
4315      ITN = ITN - 1
4316C     .......... LOOK FOR TWO CONSECUTIVE SMALL
4317C                SUB-DIAGONAL ELEMENTS.
4318C                FOR M=EN-2 STEP -1 UNTIL L DO -- ..........
4319      DO 140 MM = L, ENM2
4320         M = ENM2 + L - MM
4321         ZZ = H(M,M)
4322         R = X - ZZ
4323         S = Y - ZZ
4324         P = (R * S - W) / H(M+1,M) + H(M,M+1)
4325         Q = H(M+1,M+1) - ZZ - R - S
4326         R = H(M+2,M+1)
4327         S = DABS(P) + DABS(Q) + DABS(R)
4328         P = P / S
4329         Q = Q / S
4330         R = R / S
4331         IF (M .EQ. L) GO TO 150
4332         TST1 = DABS(P)*(DABS(H(M-1,M-1)) + DABS(ZZ) + DABS(H(M+1,M+1)))
4333         TST2 = TST1 + DABS(H(M,M-1))*(DABS(Q) + DABS(R))
4334         IF (TST2 .EQ. TST1) GO TO 150
4335  140 CONTINUE
4336C
4337  150 MP2 = M + 2
4338C
4339      DO 160 I = MP2, EN
4340         H(I,I-2) = 0.0D0
4341         IF (I .EQ. MP2) GO TO 160
4342         H(I,I-3) = 0.0D0
4343  160 CONTINUE
4344C     .......... DOUBLE QR STEP INVOLVING ROWS L TO EN AND
4345C                COLUMNS M TO EN ..........
4346      DO 260 K = M, NA
4347         NOTLAS = K .NE. NA
4348         IF (K .EQ. M) GO TO 170
4349         P = H(K,K-1)
4350         Q = H(K+1,K-1)
4351         R = 0.0D0
4352         IF (NOTLAS) R = H(K+2,K-1)
4353         X = DABS(P) + DABS(Q) + DABS(R)
4354         IF (X .EQ. 0.0D0) GO TO 260
4355         P = P / X
4356         Q = Q / X
4357         R = R / X
4358  170    S = DSIGN(DSQRT(P*P+Q*Q+R*R),P)
4359         IF (K .EQ. M) GO TO 180
4360         H(K,K-1) = -S * X
4361         GO TO 190
4362  180    IF (L .NE. M) H(K,K-1) = -H(K,K-1)
4363  190    P = P + S
4364         X = P / S
4365         Y = Q / S
4366         ZZ = R / S
4367         Q = Q / P
4368         R = R / P
4369         IF (NOTLAS) GO TO 225
4370C     .......... ROW MODIFICATION ..........
4371         DO 200 J = K, N
4372            P = H(K,J) + Q * H(K+1,J)
4373            H(K,J) = H(K,J) - P * X
4374            H(K+1,J) = H(K+1,J) - P * Y
4375  200    CONTINUE
4376C
4377         J = MIN0(EN,K+3)
4378C     .......... COLUMN MODIFICATION ..........
4379         DO 210 I = 1, J
4380            P = X * H(I,K) + Y * H(I,K+1)
4381            H(I,K) = H(I,K) - P
4382            H(I,K+1) = H(I,K+1) - P * Q
4383  210    CONTINUE
4384         GO TO 255
4385  225    CONTINUE
4386C     .......... ROW MODIFICATION ..........
4387         DO 230 J = K, N
4388            P = H(K,J) + Q * H(K+1,J) + R * H(K+2,J)
4389            H(K,J) = H(K,J) - P * X
4390            H(K+1,J) = H(K+1,J) - P * Y
4391            H(K+2,J) = H(K+2,J) - P * ZZ
4392  230    CONTINUE
4393C
4394         J = MIN0(EN,K+3)
4395C     .......... COLUMN MODIFICATION ..........
4396         DO 240 I = 1, J
4397            P = X * H(I,K) + Y * H(I,K+1) + ZZ * H(I,K+2)
4398            H(I,K) = H(I,K) - P
4399            H(I,K+1) = H(I,K+1) - P * Q
4400            H(I,K+2) = H(I,K+2) - P * R
4401  240    CONTINUE
4402  255    CONTINUE
4403C
4404  260 CONTINUE
4405C
4406      GO TO 70
4407C     .......... ONE ROOT FOUND ..........
4408  270 WR(EN) = X + T
4409      WI(EN) = 0.0D0
4410      EN = NA
4411      GO TO 60
4412C     .......... TWO ROOTS FOUND ..........
4413  280 P = (Y - X) / 2.0D0
4414      Q = P * P + W
4415      ZZ = DSQRT(DABS(Q))
4416      X = X + T
4417      IF (Q .LT. 0.0D0) GO TO 320
4418C     .......... REAL PAIR ..........
4419      ZZ = P + DSIGN(ZZ,P)
4420      WR(NA) = X + ZZ
4421      WR(EN) = WR(NA)
4422      IF (ZZ .NE. 0.0D0) WR(EN) = X - W / ZZ
4423      WI(NA) = 0.0D0
4424      WI(EN) = 0.0D0
4425      GO TO 330
4426C     .......... COMPLEX PAIR ..........
4427  320 WR(NA) = X + P
4428      WR(EN) = X + P
4429      WI(NA) = ZZ
4430      WI(EN) = -ZZ
4431  330 EN = ENM2
4432      GO TO 60
4433C     .......... SET ERROR -- ALL EIGENVALUES HAVE NOT
4434C                CONVERGED AFTER 30*N ITERATIONS ..........
4435 1000 IERR = EN
4436 1001 RETURN
4437      END
4438      SUBROUTINE HQR2(NM,N,LOW,IGH,H,WR,WI,Z,IERR)
4439C
4440      INTEGER I,J,K,L,M,N,EN,II,JJ,LL,MM,NA,NM,NN,
4441     X        IGH,ITN,ITS,LOW,MP2,ENM2,IERR
4442      DOUBLE PRECISION H(NM,N),WR(N),WI(N),Z(NM,N)
4443      DOUBLE PRECISION P,Q,R,S,T,W,X,Y,RA,SA,VI,VR,ZZ,NORM,TST1,TST2
4444      LOGICAL NOTLAS
4445C
4446C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR2,
4447C     NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON.
4448C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
4449C
4450C     THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
4451C     OF A REAL UPPER HESSENBERG MATRIX BY THE QR METHOD.  THE
4452C     EIGENVECTORS OF A REAL GENERAL MATRIX CAN ALSO BE FOUND
4453C     IF  ELMHES  AND  ELTRAN  OR  ORTHES  AND  ORTRAN  HAVE
4454C     BEEN USED TO REDUCE THIS GENERAL MATRIX TO HESSENBERG FORM
4455C     AND TO ACCUMULATE THE SIMILARITY TRANSFORMATIONS.
4456C
4457C     ON INPUT
4458C
4459C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
4460C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
4461C          DIMENSION STATEMENT.
4462C
4463C        N IS THE ORDER OF THE MATRIX.
4464C
4465C        LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
4466C          SUBROUTINE  BALANC.  IF  BALANC  HAS NOT BEEN USED,
4467C          SET LOW=1, IGH=N.
4468C
4469C        H CONTAINS THE UPPER HESSENBERG MATRIX.
4470C
4471C        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED BY  ELTRAN
4472C          AFTER THE REDUCTION BY  ELMHES, OR BY  ORTRAN  AFTER THE
4473C          REDUCTION BY  ORTHES, IF PERFORMED.  IF THE EIGENVECTORS
4474C          OF THE HESSENBERG MATRIX ARE DESIRED, Z MUST CONTAIN THE
4475C          IDENTITY MATRIX.
4476C
4477C     ON OUTPUT
4478C
4479C        H HAS BEEN DESTROYED.
4480C
4481C        WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
4482C          RESPECTIVELY, OF THE EIGENVALUES.  THE EIGENVALUES
4483C          ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS
4484C          OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE
4485C          HAVING THE POSITIVE IMAGINARY PART FIRST.  IF AN
4486C          ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT
4487C          FOR INDICES IERR+1,...,N.
4488C
4489C        Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS.
4490C          IF THE I-TH EIGENVALUE IS REAL, THE I-TH COLUMN OF Z
4491C          CONTAINS ITS EIGENVECTOR.  IF THE I-TH EIGENVALUE IS COMPLEX
4492C          WITH POSITIVE IMAGINARY PART, THE I-TH AND (I+1)-TH
4493C          COLUMNS OF Z CONTAIN THE REAL AND IMAGINARY PARTS OF ITS
4494C          EIGENVECTOR.  THE EIGENVECTORS ARE UNNORMALIZED.  IF AN
4495C          ERROR EXIT IS MADE, NONE OF THE EIGENVECTORS HAS BEEN FOUND.
4496C
4497C        IERR IS SET TO
4498C          ZERO       FOR NORMAL RETURN,
4499C          J          IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED
4500C                     WHILE THE J-TH EIGENVALUE IS BEING SOUGHT.
4501C
4502C     CALLS CDIV FOR COMPLEX DIVISION.
4503C
4504C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
4505C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
4506C
4507C     THIS VERSION DATED AUGUST 1983.
4508C
4509C     ------------------------------------------------------------------
4510C
4511      IERR = 0
4512      NORM = 0.0D0
4513      K = 1
4514C     .......... STORE ROOTS ISOLATED BY BALANC
4515C                AND COMPUTE MATRIX NORM ..........
4516      DO 50 I = 1, N
4517C
4518         DO 40 J = K, N
4519   40    NORM = NORM + DABS(H(I,J))
4520C
4521         K = I
4522         IF (I .GE. LOW .AND. I .LE. IGH) GO TO 50
4523         WR(I) = H(I,I)
4524         WI(I) = 0.0D0
4525   50 CONTINUE
4526C
4527      EN = IGH
4528      T = 0.0D0
4529      ITN = 30*N
4530C     .......... SEARCH FOR NEXT EIGENVALUES ..........
4531   60 IF (EN .LT. LOW) GO TO 340
4532      ITS = 0
4533      NA = EN - 1
4534      ENM2 = NA - 1
4535C     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
4536C                FOR L=EN STEP -1 UNTIL LOW DO -- ..........
4537   70 DO 80 LL = LOW, EN
4538         L = EN + LOW - LL
4539         IF (L .EQ. LOW) GO TO 100
4540         S = DABS(H(L-1,L-1)) + DABS(H(L,L))
4541         IF (S .EQ. 0.0D0) S = NORM
4542         TST1 = S
4543         TST2 = TST1 + DABS(H(L,L-1))
4544         IF (TST2 .EQ. TST1) GO TO 100
4545   80 CONTINUE
4546C     .......... FORM SHIFT ..........
4547  100 X = H(EN,EN)
4548      IF (L .EQ. EN) GO TO 270
4549      Y = H(NA,NA)
4550      W = H(EN,NA) * H(NA,EN)
4551      IF (L .EQ. NA) GO TO 280
4552      IF (ITN .EQ. 0) GO TO 1000
4553      IF (ITS .NE. 10 .AND. ITS .NE. 20) GO TO 130
4554C     .......... FORM EXCEPTIONAL SHIFT ..........
4555      T = T + X
4556C
4557      DO 120 I = LOW, EN
4558  120 H(I,I) = H(I,I) - X
4559C
4560      S = DABS(H(EN,NA)) + DABS(H(NA,ENM2))
4561      X = 0.75D0 * S
4562      Y = X
4563      W = -0.4375D0 * S * S
4564  130 ITS = ITS + 1
4565      ITN = ITN - 1
4566C     .......... LOOK FOR TWO CONSECUTIVE SMALL
4567C                SUB-DIAGONAL ELEMENTS.
4568C                FOR M=EN-2 STEP -1 UNTIL L DO -- ..........
4569      DO 140 MM = L, ENM2
4570         M = ENM2 + L - MM
4571         ZZ = H(M,M)
4572         R = X - ZZ
4573         S = Y - ZZ
4574         P = (R * S - W) / H(M+1,M) + H(M,M+1)
4575         Q = H(M+1,M+1) - ZZ - R - S
4576         R = H(M+2,M+1)
4577         S = DABS(P) + DABS(Q) + DABS(R)
4578         P = P / S
4579         Q = Q / S
4580         R = R / S
4581         IF (M .EQ. L) GO TO 150
4582         TST1 = DABS(P)*(DABS(H(M-1,M-1)) + DABS(ZZ) + DABS(H(M+1,M+1)))
4583         TST2 = TST1 + DABS(H(M,M-1))*(DABS(Q) + DABS(R))
4584         IF (TST2 .EQ. TST1) GO TO 150
4585  140 CONTINUE
4586C
4587  150 MP2 = M + 2
4588C
4589      DO 160 I = MP2, EN
4590         H(I,I-2) = 0.0D0
4591         IF (I .EQ. MP2) GO TO 160
4592         H(I,I-3) = 0.0D0
4593  160 CONTINUE
4594C     .......... DOUBLE QR STEP INVOLVING ROWS L TO EN AND
4595C                COLUMNS M TO EN ..........
4596      DO 260 K = M, NA
4597         NOTLAS = K .NE. NA
4598         IF (K .EQ. M) GO TO 170
4599         P = H(K,K-1)
4600         Q = H(K+1,K-1)
4601         R = 0.0D0
4602         IF (NOTLAS) R = H(K+2,K-1)
4603         X = DABS(P) + DABS(Q) + DABS(R)
4604         IF (X .EQ. 0.0D0) GO TO 260
4605         P = P / X
4606         Q = Q / X
4607         R = R / X
4608  170    S = DSIGN(DSQRT(P*P+Q*Q+R*R),P)
4609         IF (K .EQ. M) GO TO 180
4610         H(K,K-1) = -S * X
4611         GO TO 190
4612  180    IF (L .NE. M) H(K,K-1) = -H(K,K-1)
4613  190    P = P + S
4614         X = P / S
4615         Y = Q / S
4616         ZZ = R / S
4617         Q = Q / P
4618         R = R / P
4619         IF (NOTLAS) GO TO 225
4620C     .......... ROW MODIFICATION ..........
4621         DO 200 J = K, N
4622            P = H(K,J) + Q * H(K+1,J)
4623            H(K,J) = H(K,J) - P * X
4624            H(K+1,J) = H(K+1,J) - P * Y
4625  200    CONTINUE
4626C
4627         J = MIN0(EN,K+3)
4628C     .......... COLUMN MODIFICATION ..........
4629         DO 210 I = 1, J
4630            P = X * H(I,K) + Y * H(I,K+1)
4631            H(I,K) = H(I,K) - P
4632            H(I,K+1) = H(I,K+1) - P * Q
4633  210    CONTINUE
4634C     .......... ACCUMULATE TRANSFORMATIONS ..........
4635         DO 220 I = LOW, IGH
4636            P = X * Z(I,K) + Y * Z(I,K+1)
4637            Z(I,K) = Z(I,K) - P
4638            Z(I,K+1) = Z(I,K+1) - P * Q
4639  220    CONTINUE
4640         GO TO 255
4641  225    CONTINUE
4642C     .......... ROW MODIFICATION ..........
4643         DO 230 J = K, N
4644            P = H(K,J) + Q * H(K+1,J) + R * H(K+2,J)
4645            H(K,J) = H(K,J) - P * X
4646            H(K+1,J) = H(K+1,J) - P * Y
4647            H(K+2,J) = H(K+2,J) - P * ZZ
4648  230    CONTINUE
4649C
4650         J = MIN0(EN,K+3)
4651C     .......... COLUMN MODIFICATION ..........
4652         DO 240 I = 1, J
4653            P = X * H(I,K) + Y * H(I,K+1) + ZZ * H(I,K+2)
4654            H(I,K) = H(I,K) - P
4655            H(I,K+1) = H(I,K+1) - P * Q
4656            H(I,K+2) = H(I,K+2) - P * R
4657  240    CONTINUE
4658C     .......... ACCUMULATE TRANSFORMATIONS ..........
4659         DO 250 I = LOW, IGH
4660            P = X * Z(I,K) + Y * Z(I,K+1) + ZZ * Z(I,K+2)
4661            Z(I,K) = Z(I,K) - P
4662            Z(I,K+1) = Z(I,K+1) - P * Q
4663            Z(I,K+2) = Z(I,K+2) - P * R
4664  250    CONTINUE
4665  255    CONTINUE
4666C
4667  260 CONTINUE
4668C
4669      GO TO 70
4670C     .......... ONE ROOT FOUND ..........
4671  270 H(EN,EN) = X + T
4672      WR(EN) = H(EN,EN)
4673      WI(EN) = 0.0D0
4674      EN = NA
4675      GO TO 60
4676C     .......... TWO ROOTS FOUND ..........
4677  280 P = (Y - X) / 2.0D0
4678      Q = P * P + W
4679      ZZ = DSQRT(DABS(Q))
4680      H(EN,EN) = X + T
4681      X = H(EN,EN)
4682      H(NA,NA) = Y + T
4683      IF (Q .LT. 0.0D0) GO TO 320
4684C     .......... REAL PAIR ..........
4685      ZZ = P + DSIGN(ZZ,P)
4686      WR(NA) = X + ZZ
4687      WR(EN) = WR(NA)
4688      IF (ZZ .NE. 0.0D0) WR(EN) = X - W / ZZ
4689      WI(NA) = 0.0D0
4690      WI(EN) = 0.0D0
4691      X = H(EN,NA)
4692      S = DABS(X) + DABS(ZZ)
4693      P = X / S
4694      Q = ZZ / S
4695      R = DSQRT(P*P+Q*Q)
4696      P = P / R
4697      Q = Q / R
4698C     .......... ROW MODIFICATION ..........
4699      DO 290 J = NA, N
4700         ZZ = H(NA,J)
4701         H(NA,J) = Q * ZZ + P * H(EN,J)
4702         H(EN,J) = Q * H(EN,J) - P * ZZ
4703  290 CONTINUE
4704C     .......... COLUMN MODIFICATION ..........
4705      DO 300 I = 1, EN
4706         ZZ = H(I,NA)
4707         H(I,NA) = Q * ZZ + P * H(I,EN)
4708         H(I,EN) = Q * H(I,EN) - P * ZZ
4709  300 CONTINUE
4710C     .......... ACCUMULATE TRANSFORMATIONS ..........
4711      DO 310 I = LOW, IGH
4712         ZZ = Z(I,NA)
4713         Z(I,NA) = Q * ZZ + P * Z(I,EN)
4714         Z(I,EN) = Q * Z(I,EN) - P * ZZ
4715  310 CONTINUE
4716C
4717      GO TO 330
4718C     .......... COMPLEX PAIR ..........
4719  320 WR(NA) = X + P
4720      WR(EN) = X + P
4721      WI(NA) = ZZ
4722      WI(EN) = -ZZ
4723  330 EN = ENM2
4724      GO TO 60
4725C     .......... ALL ROOTS FOUND.  BACKSUBSTITUTE TO FIND
4726C                VECTORS OF UPPER TRIANGULAR FORM ..........
4727  340 IF (NORM .EQ. 0.0D0) GO TO 1001
4728C     .......... FOR EN=N STEP -1 UNTIL 1 DO -- ..........
4729      DO 800 NN = 1, N
4730         EN = N + 1 - NN
4731         P = WR(EN)
4732         Q = WI(EN)
4733         NA = EN - 1
4734         IF (Q) 710, 600, 800
4735C     .......... REAL VECTOR ..........
4736  600    M = EN
4737         H(EN,EN) = 1.0D0
4738         IF (NA .EQ. 0) GO TO 800
4739C     .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- ..........
4740         DO 700 II = 1, NA
4741            I = EN - II
4742            W = H(I,I) - P
4743            R = 0.0D0
4744C
4745            DO 610 J = M, EN
4746  610       R = R + H(I,J) * H(J,EN)
4747C
4748            IF (WI(I) .GE. 0.0D0) GO TO 630
4749            ZZ = W
4750            S = R
4751            GO TO 700
4752  630       M = I
4753            IF (WI(I) .NE. 0.0D0) GO TO 640
4754            T = W
4755            IF (T .NE. 0.0D0) GO TO 635
4756               TST1 = NORM
4757               T = TST1
4758  632          T = 0.01D0 * T
4759               TST2 = NORM + T
4760               IF (TST2 .GT. TST1) GO TO 632
4761  635       H(I,EN) = -R / T
4762            GO TO 680
4763C     .......... SOLVE REAL EQUATIONS ..........
4764  640       X = H(I,I+1)
4765            Y = H(I+1,I)
4766            Q = (WR(I) - P) * (WR(I) - P) + WI(I) * WI(I)
4767            T = (X * S - ZZ * R) / Q
4768            H(I,EN) = T
4769            IF (DABS(X) .LE. DABS(ZZ)) GO TO 650
4770            H(I+1,EN) = (-R - W * T) / X
4771            GO TO 680
4772  650       H(I+1,EN) = (-S - Y * T) / ZZ
4773C
4774C     .......... OVERFLOW CONTROL ..........
4775  680       T = DABS(H(I,EN))
4776            IF (T .EQ. 0.0D0) GO TO 700
4777            TST1 = T
4778            TST2 = TST1 + 1.0D0/TST1
4779            IF (TST2 .GT. TST1) GO TO 700
4780            DO 690 J = I, EN
4781               H(J,EN) = H(J,EN)/T
4782  690       CONTINUE
4783C
4784  700    CONTINUE
4785C     .......... END REAL VECTOR ..........
4786         GO TO 800
4787C     .......... COMPLEX VECTOR ..........
4788  710    M = NA
4789C     .......... LAST VECTOR COMPONENT CHOSEN IMAGINARY SO THAT
4790C                EIGENVECTOR MATRIX IS TRIANGULAR ..........
4791         IF (DABS(H(EN,NA)) .LE. DABS(H(NA,EN))) GO TO 720
4792         H(NA,NA) = Q / H(EN,NA)
4793         H(NA,EN) = -(H(EN,EN) - P) / H(EN,NA)
4794         GO TO 730
4795  720    CALL CDIV(0.0D0,-H(NA,EN),H(NA,NA)-P,Q,H(NA,NA),H(NA,EN))
4796  730    H(EN,NA) = 0.0D0
4797         H(EN,EN) = 1.0D0
4798         ENM2 = NA - 1
4799         IF (ENM2 .EQ. 0) GO TO 800
4800C     .......... FOR I=EN-2 STEP -1 UNTIL 1 DO -- ..........
4801         DO 795 II = 1, ENM2
4802            I = NA - II
4803            W = H(I,I) - P
4804            RA = 0.0D0
4805            SA = 0.0D0
4806C
4807            DO 760 J = M, EN
4808               RA = RA + H(I,J) * H(J,NA)
4809               SA = SA + H(I,J) * H(J,EN)
4810  760       CONTINUE
4811C
4812            IF (WI(I) .GE. 0.0D0) GO TO 770
4813            ZZ = W
4814            R = RA
4815            S = SA
4816            GO TO 795
4817  770       M = I
4818            IF (WI(I) .NE. 0.0D0) GO TO 780
4819            CALL CDIV(-RA,-SA,W,Q,H(I,NA),H(I,EN))
4820            GO TO 790
4821C     .......... SOLVE COMPLEX EQUATIONS ..........
4822  780       X = H(I,I+1)
4823            Y = H(I+1,I)
4824            VR = (WR(I) - P) * (WR(I) - P) + WI(I) * WI(I) - Q * Q
4825            VI = (WR(I) - P) * 2.0D0 * Q
4826            IF (VR .NE. 0.0D0 .OR. VI .NE. 0.0D0) GO TO 784
4827               TST1 = NORM * (DABS(W) + DABS(Q) + DABS(X)
4828     X                      + DABS(Y) + DABS(ZZ))
4829               VR = TST1
4830  783          VR = 0.01D0 * VR
4831               TST2 = TST1 + VR
4832               IF (TST2 .GT. TST1) GO TO 783
4833  784       CALL CDIV(X*R-ZZ*RA+Q*SA,X*S-ZZ*SA-Q*RA,VR,VI,
4834     X                H(I,NA),H(I,EN))
4835            IF (DABS(X) .LE. DABS(ZZ) + DABS(Q)) GO TO 785
4836            H(I+1,NA) = (-RA - W * H(I,NA) + Q * H(I,EN)) / X
4837            H(I+1,EN) = (-SA - W * H(I,EN) - Q * H(I,NA)) / X
4838            GO TO 790
4839  785       CALL CDIV(-R-Y*H(I,NA),-S-Y*H(I,EN),ZZ,Q,
4840     X                H(I+1,NA),H(I+1,EN))
4841C
4842C     .......... OVERFLOW CONTROL ..........
4843  790       T = DMAX1(DABS(H(I,NA)), DABS(H(I,EN)))
4844            IF (T .EQ. 0.0D0) GO TO 795
4845            TST1 = T
4846            TST2 = TST1 + 1.0D0/TST1
4847            IF (TST2 .GT. TST1) GO TO 795
4848            DO 792 J = I, EN
4849               H(J,NA) = H(J,NA)/T
4850               H(J,EN) = H(J,EN)/T
4851  792       CONTINUE
4852C
4853  795    CONTINUE
4854C     .......... END COMPLEX VECTOR ..........
4855  800 CONTINUE
4856C     .......... END BACK SUBSTITUTION.
4857C                VECTORS OF ISOLATED ROOTS ..........
4858      DO 840 I = 1, N
4859         IF (I .GE. LOW .AND. I .LE. IGH) GO TO 840
4860C
4861         DO 820 J = I, N
4862  820    Z(I,J) = H(I,J)
4863C
4864  840 CONTINUE
4865C     .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE
4866C                VECTORS OF ORIGINAL FULL MATRIX.
4867C                FOR J=N STEP -1 UNTIL LOW DO -- ..........
4868      DO 880 JJ = LOW, N
4869         J = N + LOW - JJ
4870         M = MIN0(J,IGH)
4871C
4872         DO 880 I = LOW, IGH
4873            ZZ = 0.0D0
4874C
4875            DO 860 K = LOW, M
4876  860       ZZ = ZZ + Z(I,K) * H(K,J)
4877C
4878            Z(I,J) = ZZ
4879  880 CONTINUE
4880C
4881      GO TO 1001
4882C     .......... SET ERROR -- ALL EIGENVALUES HAVE NOT
4883C                CONVERGED AFTER 30*N ITERATIONS ..........
4884 1000 IERR = EN
4885 1001 RETURN
4886      END
4887      SUBROUTINE HTRIB3(NM,N,A,TAU,M,ZR,ZI)
4888C
4889      INTEGER I,J,K,L,M,N,NM
4890      DOUBLE PRECISION A(NM,N),TAU(2,N),ZR(NM,M),ZI(NM,M)
4891      DOUBLE PRECISION H,S,SI
4892C
4893C     THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF
4894C     THE ALGOL PROCEDURE TRBAK3, NUM. MATH. 11, 181-195(1968)
4895C     BY MARTIN, REINSCH, AND WILKINSON.
4896C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
4897C
4898C     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX HERMITIAN
4899C     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
4900C     REAL SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY  HTRID3.
4901C
4902C     ON INPUT
4903C
4904C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
4905C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
4906C          DIMENSION STATEMENT.
4907C
4908C        N IS THE ORDER OF THE MATRIX.
4909C
4910C        A CONTAINS INFORMATION ABOUT THE UNITARY TRANSFORMATIONS
4911C          USED IN THE REDUCTION BY  HTRID3.
4912C
4913C        TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS.
4914C
4915C        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.
4916C
4917C        ZR CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
4918C          IN ITS FIRST M COLUMNS.
4919C
4920C     ON OUTPUT
4921C
4922C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
4923C          RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS
4924C          IN THEIR FIRST M COLUMNS.
4925C
4926C     NOTE THAT THE LAST COMPONENT OF EACH RETURNED VECTOR
4927C     IS REAL AND THAT VECTOR EUCLIDEAN NORMS ARE PRESERVED.
4928C
4929C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
4930C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
4931C
4932C     THIS VERSION DATED AUGUST 1983.
4933C
4934C     ------------------------------------------------------------------
4935C
4936      IF (M .EQ. 0) GO TO 200
4937C     .......... TRANSFORM THE EIGENVECTORS OF THE REAL SYMMETRIC
4938C                TRIDIAGONAL MATRIX TO THOSE OF THE HERMITIAN
4939C                TRIDIAGONAL MATRIX. ..........
4940      DO 50 K = 1, N
4941C
4942         DO 50 J = 1, M
4943            ZI(K,J) = -ZR(K,J) * TAU(2,K)
4944            ZR(K,J) = ZR(K,J) * TAU(1,K)
4945   50 CONTINUE
4946C
4947      IF (N .EQ. 1) GO TO 200
4948C     .......... RECOVER AND APPLY THE HOUSEHOLDER MATRICES ..........
4949      DO 140 I = 2, N
4950         L = I - 1
4951         H = A(I,I)
4952         IF (H .EQ. 0.0D0) GO TO 140
4953C
4954         DO 130 J = 1, M
4955            S = 0.0D0
4956            SI = 0.0D0
4957C
4958            DO 110 K = 1, L
4959               S = S + A(I,K) * ZR(K,J) - A(K,I) * ZI(K,J)
4960               SI = SI + A(I,K) * ZI(K,J) + A(K,I) * ZR(K,J)
4961  110       CONTINUE
4962C     .......... DOUBLE DIVISIONS AVOID POSSIBLE UNDERFLOW ..........
4963            S = (S / H) / H
4964            SI = (SI / H) / H
4965C
4966            DO 120 K = 1, L
4967               ZR(K,J) = ZR(K,J) - S * A(I,K) - SI * A(K,I)
4968               ZI(K,J) = ZI(K,J) - SI * A(I,K) + S * A(K,I)
4969  120       CONTINUE
4970C
4971  130    CONTINUE
4972C
4973  140 CONTINUE
4974C
4975  200 RETURN
4976      END
4977      SUBROUTINE HTRIBK(NM,N,AR,AI,TAU,M,ZR,ZI)
4978C
4979      INTEGER I,J,K,L,M,N,NM
4980      DOUBLE PRECISION AR(NM,N),AI(NM,N),TAU(2,N),ZR(NM,M),ZI(NM,M)
4981      DOUBLE PRECISION H,S,SI
4982C
4983C     THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF
4984C     THE ALGOL PROCEDURE TRBAK1, NUM. MATH. 11, 181-195(1968)
4985C     BY MARTIN, REINSCH, AND WILKINSON.
4986C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
4987C
4988C     THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX HERMITIAN
4989C     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
4990C     REAL SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY  HTRIDI.
4991C
4992C     ON INPUT
4993C
4994C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
4995C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
4996C          DIMENSION STATEMENT.
4997C
4998C        N IS THE ORDER OF THE MATRIX.
4999C
5000C        AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY TRANS-
5001C          FORMATIONS USED IN THE REDUCTION BY  HTRIDI  IN THEIR
5002C          FULL LOWER TRIANGLES EXCEPT FOR THE DIAGONAL OF AR.
5003C
5004C        TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS.
5005C
5006C        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.
5007C
5008C        ZR CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
5009C          IN ITS FIRST M COLUMNS.
5010C
5011C     ON OUTPUT
5012C
5013C        ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
5014C          RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS
5015C          IN THEIR FIRST M COLUMNS.
5016C
5017C     NOTE THAT THE LAST COMPONENT OF EACH RETURNED VECTOR
5018C     IS REAL AND THAT VECTOR EUCLIDEAN NORMS ARE PRESERVED.
5019C
5020C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
5021C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
5022C
5023C     THIS VERSION DATED AUGUST 1983.
5024C
5025C     ------------------------------------------------------------------
5026C
5027      IF (M .EQ. 0) GO TO 200
5028C     .......... TRANSFORM THE EIGENVECTORS OF THE REAL SYMMETRIC
5029C                TRIDIAGONAL MATRIX TO THOSE OF THE HERMITIAN
5030C                TRIDIAGONAL MATRIX. ..........
5031      DO 50 K = 1, N
5032C
5033         DO 50 J = 1, M
5034            ZI(K,J) = -ZR(K,J) * TAU(2,K)
5035            ZR(K,J) = ZR(K,J) * TAU(1,K)
5036   50 CONTINUE
5037C
5038      IF (N .EQ. 1) GO TO 200
5039C     .......... RECOVER AND APPLY THE HOUSEHOLDER MATRICES ..........
5040      DO 140 I = 2, N
5041         L = I - 1
5042         H = AI(I,I)
5043         IF (H .EQ. 0.0D0) GO TO 140
5044C
5045         DO 130 J = 1, M
5046            S = 0.0D0
5047            SI = 0.0D0
5048C
5049            DO 110 K = 1, L
5050               S = S + AR(I,K) * ZR(K,J) - AI(I,K) * ZI(K,J)
5051               SI = SI + AR(I,K) * ZI(K,J) + AI(I,K) * ZR(K,J)
5052  110       CONTINUE
5053C     .......... DOUBLE DIVISIONS AVOID POSSIBLE UNDERFLOW ..........
5054            S = (S / H) / H
5055            SI = (SI / H) / H
5056C
5057            DO 120 K = 1, L
5058               ZR(K,J) = ZR(K,J) - S * AR(I,K) - SI * AI(I,K)
5059               ZI(K,J) = ZI(K,J) - SI * AR(I,K) + S * AI(I,K)
5060  120       CONTINUE
5061C
5062  130    CONTINUE
5063C
5064  140 CONTINUE
5065C
5066  200 RETURN
5067      END
5068      SUBROUTINE HTRID3(NM,N,A,D,E,E2,TAU)
5069C
5070      INTEGER I,J,K,L,N,II,NM,JM1,JP1
5071      DOUBLE PRECISION A(NM,N),D(N),E(N),E2(N),TAU(2,N)
5072      DOUBLE PRECISION F,G,H,FI,GI,HH,SI,SCALE,PYTHAG
5073C
5074C     THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF
5075C     THE ALGOL PROCEDURE TRED3, NUM. MATH. 11, 181-195(1968)
5076C     BY MARTIN, REINSCH, AND WILKINSON.
5077C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
5078C
5079C     THIS SUBROUTINE REDUCES A COMPLEX HERMITIAN MATRIX, STORED AS
5080C     A SINGLE SQUARE ARRAY, TO A REAL SYMMETRIC TRIDIAGONAL MATRIX
5081C     USING UNITARY SIMILARITY TRANSFORMATIONS.
5082C
5083C     ON INPUT
5084C
5085C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
5086C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
5087C          DIMENSION STATEMENT.
5088C
5089C        N IS THE ORDER OF THE MATRIX.
5090C
5091C        A CONTAINS THE LOWER TRIANGLE OF THE COMPLEX HERMITIAN INPUT
5092C          MATRIX.  THE REAL PARTS OF THE MATRIX ELEMENTS ARE STORED
5093C          IN THE FULL LOWER TRIANGLE OF A, AND THE IMAGINARY PARTS
5094C          ARE STORED IN THE TRANSPOSED POSITIONS OF THE STRICT UPPER
5095C          TRIANGLE OF A.  NO STORAGE IS REQUIRED FOR THE ZERO
5096C          IMAGINARY PARTS OF THE DIAGONAL ELEMENTS.
5097C
5098C     ON OUTPUT
5099C
5100C        A CONTAINS INFORMATION ABOUT THE UNITARY TRANSFORMATIONS
5101C          USED IN THE REDUCTION.
5102C
5103C        D CONTAINS THE DIAGONAL ELEMENTS OF THE THE TRIDIAGONAL MATRIX.
5104C
5105C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
5106C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO.
5107C
5108C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
5109C          E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
5110C
5111C        TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS.
5112C
5113C     CALLS PYTHAG FOR  DSQRT(A*A + B*B) .
5114C
5115C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
5116C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
5117C
5118C     THIS VERSION DATED AUGUST 1983.
5119C
5120C     ------------------------------------------------------------------
5121C
5122      TAU(1,N) = 1.0D0
5123      TAU(2,N) = 0.0D0
5124C     .......... FOR I=N STEP -1 UNTIL 1 DO -- ..........
5125      DO 300 II = 1, N
5126         I = N + 1 - II
5127         L = I - 1
5128         H = 0.0D0
5129         SCALE = 0.0D0
5130         IF (L .LT. 1) GO TO 130
5131C     .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) ..........
5132         DO 120 K = 1, L
5133  120    SCALE = SCALE + DABS(A(I,K)) + DABS(A(K,I))
5134C
5135         IF (SCALE .NE. 0.0D0) GO TO 140
5136         TAU(1,L) = 1.0D0
5137         TAU(2,L) = 0.0D0
5138  130    E(I) = 0.0D0
5139         E2(I) = 0.0D0
5140         GO TO 290
5141C
5142  140    DO 150 K = 1, L
5143            A(I,K) = A(I,K) / SCALE
5144            A(K,I) = A(K,I) / SCALE
5145            H = H + A(I,K) * A(I,K) + A(K,I) * A(K,I)
5146  150    CONTINUE
5147C
5148         E2(I) = SCALE * SCALE * H
5149         G = DSQRT(H)
5150         E(I) = SCALE * G
5151         F = PYTHAG(A(I,L),A(L,I))
5152C     .......... FORM NEXT DIAGONAL ELEMENT OF MATRIX T ..........
5153         IF (F .EQ. 0.0D0) GO TO 160
5154         TAU(1,L) = (A(L,I) * TAU(2,I) - A(I,L) * TAU(1,I)) / F
5155         SI = (A(I,L) * TAU(2,I) + A(L,I) * TAU(1,I)) / F
5156         H = H + F * G
5157         G = 1.0D0 + G / F
5158         A(I,L) = G * A(I,L)
5159         A(L,I) = G * A(L,I)
5160         IF (L .EQ. 1) GO TO 270
5161         GO TO 170
5162  160    TAU(1,L) = -TAU(1,I)
5163         SI = TAU(2,I)
5164         A(I,L) = G
5165  170    F = 0.0D0
5166C
5167         DO 240 J = 1, L
5168            G = 0.0D0
5169            GI = 0.0D0
5170            IF (J .EQ. 1) GO TO 190
5171            JM1 = J - 1
5172C     .......... FORM ELEMENT OF A*U ..........
5173            DO 180 K = 1, JM1
5174               G = G + A(J,K) * A(I,K) + A(K,J) * A(K,I)
5175               GI = GI - A(J,K) * A(K,I) + A(K,J) * A(I,K)
5176  180       CONTINUE
5177C
5178  190       G = G + A(J,J) * A(I,J)
5179            GI = GI - A(J,J) * A(J,I)
5180            JP1 = J + 1
5181            IF (L .LT. JP1) GO TO 220
5182C
5183            DO 200 K = JP1, L
5184               G = G + A(K,J) * A(I,K) - A(J,K) * A(K,I)
5185               GI = GI - A(K,J) * A(K,I) - A(J,K) * A(I,K)
5186  200       CONTINUE
5187C     .......... FORM ELEMENT OF P ..........
5188  220       E(J) = G / H
5189            TAU(2,J) = GI / H
5190            F = F + E(J) * A(I,J) - TAU(2,J) * A(J,I)
5191  240    CONTINUE
5192C
5193         HH = F / (H + H)
5194C     .......... FORM REDUCED A ..........
5195         DO 260 J = 1, L
5196            F = A(I,J)
5197            G = E(J) - HH * F
5198            E(J) = G
5199            FI = -A(J,I)
5200            GI = TAU(2,J) - HH * FI
5201            TAU(2,J) = -GI
5202            A(J,J) = A(J,J) - 2.0D0 * (F * G + FI * GI)
5203            IF (J .EQ. 1) GO TO 260
5204            JM1 = J - 1
5205C
5206            DO 250 K = 1, JM1
5207               A(J,K) = A(J,K) - F * E(K) - G * A(I,K)
5208     X                         + FI * TAU(2,K) + GI * A(K,I)
5209               A(K,J) = A(K,J) - F * TAU(2,K) - G * A(K,I)
5210     X                         - FI * E(K) - GI * A(I,K)
5211  250       CONTINUE
5212C
5213  260    CONTINUE
5214C
5215  270    DO 280 K = 1, L
5216            A(I,K) = SCALE * A(I,K)
5217            A(K,I) = SCALE * A(K,I)
5218  280    CONTINUE
5219C
5220         TAU(2,L) = -SI
5221  290    D(I) = A(I,I)
5222         A(I,I) = SCALE * DSQRT(H)
5223  300 CONTINUE
5224C
5225      RETURN
5226      END
5227      SUBROUTINE HTRIDI(NM,N,AR,AI,D,E,E2,TAU)
5228C
5229      INTEGER I,J,K,L,N,II,NM,JP1
5230      DOUBLE PRECISION AR(NM,N),AI(NM,N),D(N),E(N),E2(N),TAU(2,N)
5231      DOUBLE PRECISION F,G,H,FI,GI,HH,SI,SCALE,PYTHAG
5232C
5233C     THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF
5234C     THE ALGOL PROCEDURE TRED1, NUM. MATH. 11, 181-195(1968)
5235C     BY MARTIN, REINSCH, AND WILKINSON.
5236C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
5237C
5238C     THIS SUBROUTINE REDUCES A COMPLEX HERMITIAN MATRIX
5239C     TO A REAL SYMMETRIC TRIDIAGONAL MATRIX USING
5240C     UNITARY SIMILARITY TRANSFORMATIONS.
5241C
5242C     ON INPUT
5243C
5244C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
5245C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
5246C          DIMENSION STATEMENT.
5247C
5248C        N IS THE ORDER OF THE MATRIX.
5249C
5250C        AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
5251C          RESPECTIVELY, OF THE COMPLEX HERMITIAN INPUT MATRIX.
5252C          ONLY THE LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.
5253C
5254C     ON OUTPUT
5255C
5256C        AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY TRANS-
5257C          FORMATIONS USED IN THE REDUCTION IN THEIR FULL LOWER
5258C          TRIANGLES.  THEIR STRICT UPPER TRIANGLES AND THE
5259C          DIAGONAL OF AR ARE UNALTERED.
5260C
5261C        D CONTAINS THE DIAGONAL ELEMENTS OF THE THE TRIDIAGONAL MATRIX.
5262C
5263C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
5264C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO.
5265C
5266C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
5267C          E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
5268C
5269C        TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS.
5270C
5271C     CALLS PYTHAG FOR  DSQRT(A*A + B*B) .
5272C
5273C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
5274C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
5275C
5276C     THIS VERSION DATED AUGUST 1983.
5277C
5278C     ------------------------------------------------------------------
5279C
5280      TAU(1,N) = 1.0D0
5281      TAU(2,N) = 0.0D0
5282C
5283      DO 100 I = 1, N
5284  100 D(I) = AR(I,I)
5285C     .......... FOR I=N STEP -1 UNTIL 1 DO -- ..........
5286      DO 300 II = 1, N
5287         I = N + 1 - II
5288         L = I - 1
5289         H = 0.0D0
5290         SCALE = 0.0D0
5291         IF (L .LT. 1) GO TO 130
5292C     .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) ..........
5293         DO 120 K = 1, L
5294  120    SCALE = SCALE + DABS(AR(I,K)) + DABS(AI(I,K))
5295C
5296         IF (SCALE .NE. 0.0D0) GO TO 140
5297         TAU(1,L) = 1.0D0
5298         TAU(2,L) = 0.0D0
5299  130    E(I) = 0.0D0
5300         E2(I) = 0.0D0
5301         GO TO 290
5302C
5303  140    DO 150 K = 1, L
5304            AR(I,K) = AR(I,K) / SCALE
5305            AI(I,K) = AI(I,K) / SCALE
5306            H = H + AR(I,K) * AR(I,K) + AI(I,K) * AI(I,K)
5307  150    CONTINUE
5308C
5309         E2(I) = SCALE * SCALE * H
5310         G = DSQRT(H)
5311         E(I) = SCALE * G
5312         F = PYTHAG(AR(I,L),AI(I,L))
5313C     .......... FORM NEXT DIAGONAL ELEMENT OF MATRIX T ..........
5314         IF (F .EQ. 0.0D0) GO TO 160
5315         TAU(1,L) = (AI(I,L) * TAU(2,I) - AR(I,L) * TAU(1,I)) / F
5316         SI = (AR(I,L) * TAU(2,I) + AI(I,L) * TAU(1,I)) / F
5317         H = H + F * G
5318         G = 1.0D0 + G / F
5319         AR(I,L) = G * AR(I,L)
5320         AI(I,L) = G * AI(I,L)
5321         IF (L .EQ. 1) GO TO 270
5322         GO TO 170
5323  160    TAU(1,L) = -TAU(1,I)
5324         SI = TAU(2,I)
5325         AR(I,L) = G
5326  170    F = 0.0D0
5327C
5328         DO 240 J = 1, L
5329            G = 0.0D0
5330            GI = 0.0D0
5331C     .......... FORM ELEMENT OF A*U ..........
5332            DO 180 K = 1, J
5333               G = G + AR(J,K) * AR(I,K) + AI(J,K) * AI(I,K)
5334               GI = GI - AR(J,K) * AI(I,K) + AI(J,K) * AR(I,K)
5335  180       CONTINUE
5336C
5337            JP1 = J + 1
5338            IF (L .LT. JP1) GO TO 220
5339C
5340            DO 200 K = JP1, L
5341               G = G + AR(K,J) * AR(I,K) - AI(K,J) * AI(I,K)
5342               GI = GI - AR(K,J) * AI(I,K) - AI(K,J) * AR(I,K)
5343  200       CONTINUE
5344C     .......... FORM ELEMENT OF P ..........
5345  220       E(J) = G / H
5346            TAU(2,J) = GI / H
5347            F = F + E(J) * AR(I,J) - TAU(2,J) * AI(I,J)
5348  240    CONTINUE
5349C
5350         HH = F / (H + H)
5351C     .......... FORM REDUCED A ..........
5352         DO 260 J = 1, L
5353            F = AR(I,J)
5354            G = E(J) - HH * F
5355            E(J) = G
5356            FI = -AI(I,J)
5357            GI = TAU(2,J) - HH * FI
5358            TAU(2,J) = -GI
5359C
5360            DO 260 K = 1, J
5361               AR(J,K) = AR(J,K) - F * E(K) - G * AR(I,K)
5362     X                           + FI * TAU(2,K) + GI * AI(I,K)
5363               AI(J,K) = AI(J,K) - F * TAU(2,K) - G * AI(I,K)
5364     X                           - FI * E(K) - GI * AR(I,K)
5365  260    CONTINUE
5366C
5367  270    DO 280 K = 1, L
5368            AR(I,K) = SCALE * AR(I,K)
5369            AI(I,K) = SCALE * AI(I,K)
5370  280    CONTINUE
5371C
5372         TAU(2,L) = -SI
5373  290    HH = D(I)
5374         D(I) = AR(I,I)
5375         AR(I,I) = HH
5376         AI(I,I) = SCALE * DSQRT(H)
5377  300 CONTINUE
5378C
5379      RETURN
5380      END
5381      SUBROUTINE IMTQL1(N,D,E,IERR)
5382C
5383      INTEGER I,J,L,M,N,II,MML,IERR
5384      DOUBLE PRECISION D(N),E(N)
5385      DOUBLE PRECISION B,C,F,G,P,R,S,TST1,TST2,PYTHAG
5386C
5387C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE IMTQL1,
5388C     NUM. MATH. 12, 377-383(1968) BY MARTIN AND WILKINSON,
5389C     AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE.
5390C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971).
5391C
5392C     THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC
5393C     TRIDIAGONAL MATRIX BY THE IMPLICIT QL METHOD.
5394C
5395C     ON INPUT
5396C
5397C        N IS THE ORDER OF THE MATRIX.
5398C
5399C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
5400C
5401C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
5402C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
5403C
5404C      ON OUTPUT
5405C
5406C        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
5407C          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND
5408C          ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE
5409C          THE SMALLEST EIGENVALUES.
5410C
5411C        E HAS BEEN DESTROYED.
5412C
5413C        IERR IS SET TO
5414C          ZERO       FOR NORMAL RETURN,
5415C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
5416C                     DETERMINED AFTER 30 ITERATIONS.
5417C
5418C     CALLS PYTHAG FOR  DSQRT(A*A + B*B) .
5419C
5420C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW,
5421C     MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY
5422C
5423C     THIS VERSION DATED AUGUST 1983.
5424C
5425C     ------------------------------------------------------------------
5426C
5427      IERR = 0
5428      IF (N .EQ. 1) GO TO 1001
5429C
5430      DO 100 I = 2, N
5431  100 E(I-1) = E(I)
5432C
5433      E(N) = 0.0D0
5434C
5435      DO 290 L = 1, N
5436         J = 0
5437C     .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT ..........
5438  105    DO 110 M = L, N
5439            IF (M .EQ. N) GO TO 120
5440            TST1 = DABS(D(M)) + DABS(D(M+1))
5441            TST2 = TST1 + DABS(E(M))
5442            IF (TST2 .EQ. TST1) GO TO 120
5443  110    CONTINUE
5444C
5445  120    P = D(L)
5446         IF (M .EQ. L) GO TO 215
5447         IF (J .EQ. 30) GO TO 1000
5448         J = J + 1
5449C     .......... FORM SHIFT ..........
5450         G = (D(L+1) - P) / (2.0D0 * E(L))
5451         R = PYTHAG(G,1.0D0)
5452         G = D(M) - P + E(L) / (G + DSIGN(R,G))
5453         S = 1.0D0
5454         C = 1.0D0
5455         P = 0.0D0
5456         MML = M - L
5457C     .......... FOR I=M-1 STEP -1 UNTIL L DO -- ..........
5458         DO 200 II = 1, MML
5459            I = M - II
5460            F = S * E(I)
5461            B = C * E(I)
5462            R = PYTHAG(F,G)
5463            E(I+1) = R
5464            IF (R .EQ. 0.0D0) GO TO 210
5465            S = F / R
5466            C = G / R
5467            G = D(I+1) - P
5468            R = (D(I) - G) * S + 2.0D0 * C * B
5469            P = S * R
5470            D(I+1) = G + P
5471            G = C * R - B
5472  200    CONTINUE
5473C
5474         D(L) = D(L) - P
5475         E(L) = G
5476         E(M) = 0.0D0
5477         GO TO 105
5478C     .......... RECOVER FROM UNDERFLOW ..........
5479  210    D(I+1) = D(I+1) - P
5480         E(M) = 0.0D0
5481         GO TO 105