1      SUBROUTINE COMQR2(NM,N,LOW,IGH,ORTR,ORTI,HR,HI,WR,WI,ZR,ZI,IERR)
2C***BEGIN PROLOGUE  COMQR2
3C***DATE WRITTEN   760101   (YYMMDD)
4C***REVISION DATE  830518   (YYMMDD)
5C***CATEGORY NO.  D4C2B
6C***KEYWORDS  EIGENVALUES,EIGENVECTORS,EISPACK
7C***AUTHOR  SMITH, B. T., ET AL.
8C***PURPOSE  Computes eigenvalues and eigenvectors of complex upper
9C            Hessenberg matrix.
10C***DESCRIPTION
11C
12C     This subroutine is a translation of a unitary analogue of the
13C     ALGOL procedure  COMLR2, NUM. MATH. 16, 181-204(1970) by Peters
14C     and Wilkinson.
15C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
16C     The unitary analogue substitutes the QR algorithm of Francis
17C     (COMP. JOUR. 4, 332-345(1962)) for the LR algorithm.
18C
19C     This subroutine finds the eigenvalues and eigenvectors
20C     of a COMPLEX UPPER Hessenberg matrix by the QR
21C     method.  The eigenvectors of a COMPLEX GENERAL matrix
22C     can also be found if  CORTH  has been used to reduce
23C     this general matrix to Hessenberg form.
24C
25C     On INPUT
26C
27C        NM must be set to the row dimension of two-dimensional
28C          array parameters as declared in the calling program
29C          dimension statement.
30C
31C        N is the order of the matrix.
32C
33C        LOW and IGH are integers determined by the balancing
34C          subroutine  CBAL.  If  CBAL  has not been used,
35C          set LOW=1, IGH=N.
36C
37C        ORTR and ORTI contain information about the unitary trans-
38C          formations used in the reduction by  CORTH, if performed.
39C          Only elements LOW through IGH are used.  If the eigenvectors
40C          of the Hessenberg matrix are desired, set ORTR(J) and
41C          ORTI(J) to 0.0E0 for these elements.
42C
43C        HR and HI contain the real and imaginary parts,
44C          respectively, of the complex upper Hessenberg matrix.
45C          Their lower triangles below the subdiagonal contain further
46C          information about the transformations which were used in the
47C          reduction by  CORTH, if performed.  If the eigenvectors of
48C          the Hessenberg matrix are desired, these elements may be
49C          arbitrary.
50C
51C     On OUTPUT
52C
53C        ORTR, ORTI, and the upper Hessenberg portions of HR and HI
54C          have been destroyed.
55C
56C        WR and WI contain the real and imaginary parts,
57C          respectively, of the eigenvalues.  If an error
58C          exit is made, the eigenvalues should be correct
59C          for indices IERR+1,...,N.
60C
61C        ZR and ZI contain the real and imaginary parts,
62C          respectively, of the eigenvectors.  The eigenvectors
63C          are unnormalized.  If an error exit is made, none of
64C          the eigenvectors has been found.
65C
66C        IERR is set to
67C          Zero       for normal return,
68C          J          if the J-th eigenvalue has not been
69C                     determined after a total of 30*N iterations.
70C
71C     Calls CSROOT for complex square root.
72C     Calls PYTHAG(A,B) for sqrt(A**2 + B**2).
73C     Calls CDIV for complex division.
74C
75C     Questions and comments should be directed to B. S. Garbow,
76C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
77C     ------------------------------------------------------------------
78C***REFERENCES  B. T. SMITH, J. M. BOYLE, J. J. DONGARRA, B. S. GARBOW,
79C                 Y. IKEBE, V. C. KLEMA, C. B. MOLER, *MATRIX EIGEN-
80C                 SYSTEM ROUTINES - EISPACK GUIDE*, SPRINGER-VERLAG,
81C                 1976.
82C***ROUTINES CALLED  CDIV,CSROOT,PYTHAG
83C***END PROLOGUE  COMQR2
84C
85      INTEGER I,J,K,L,M,N,EN,II,JJ,LL,NM,NN,IGH,IP1
86      INTEGER ITN,ITS,LOW,LP1,ENM1,IEND,IERR
87      REAL HR(NM,N),HI(NM,N),WR(N),WI(N),ZR(NM,N),ZI(NM,N)
88      REAL ORTR(IGH),ORTI(IGH)
89      REAL SI,SR,TI,TR,XI,XR,YI,YR,ZZI,ZZR,NORM,S1,S2
90      REAL PYTHAG
91C
92C***FIRST EXECUTABLE STATEMENT  COMQR2
93      IERR = 0
94C     .......... INITIALIZE EIGENVECTOR MATRIX ..........
95      DO 100 I = 1, N
96C
97         DO 100 J = 1, N
98            ZR(I,J) = 0.0E0
99            ZI(I,J) = 0.0E0
100            IF (I .EQ. J) ZR(I,J) = 1.0E0
101  100 CONTINUE
102C     .......... FORM THE MATRIX OF ACCUMULATED TRANSFORMATIONS
103C                FROM THE INFORMATION LEFT BY CORTH ..........
104      IEND = IGH - LOW - 1
105      IF (IEND) 180, 150, 105
106C     .......... FOR I=IGH-1 STEP -1 UNTIL LOW+1 DO -- ..........
107  105 DO 140 II = 1, IEND
108         I = IGH - II
109         IF (ORTR(I) .EQ. 0.0E0 .AND. ORTI(I) .EQ. 0.0E0) GO TO 140
110         IF (HR(I,I-1) .EQ. 0.0E0 .AND. HI(I,I-1) .EQ. 0.0E0) GO TO 140
111C     .......... NORM BELOW IS NEGATIVE OF H FORMED IN CORTH ..........
112         NORM = HR(I,I-1) * ORTR(I) + HI(I,I-1) * ORTI(I)
113         IP1 = I + 1
114C
115         DO 110 K = IP1, IGH
116            ORTR(K) = HR(K,I-1)
117            ORTI(K) = HI(K,I-1)
118  110    CONTINUE
119C
120         DO 130 J = I, IGH
121            SR = 0.0E0
122            SI = 0.0E0
123C
124            DO 115 K = I, IGH
125               SR = SR + ORTR(K) * ZR(K,J) + ORTI(K) * ZI(K,J)
126               SI = SI + ORTR(K) * ZI(K,J) - ORTI(K) * ZR(K,J)
127  115       CONTINUE
128C
129            SR = SR / NORM
130            SI = SI / NORM
131C
132            DO 120 K = I, IGH
133               ZR(K,J) = ZR(K,J) + SR * ORTR(K) - SI * ORTI(K)
134               ZI(K,J) = ZI(K,J) + SR * ORTI(K) + SI * ORTR(K)
135  120       CONTINUE
136C
137  130    CONTINUE
138C
139  140 CONTINUE
140C     .......... CREATE REAL SUBDIAGONAL ELEMENTS ..........
141  150 L = LOW + 1
142C
143      DO 170 I = L, IGH
144         LL = MIN0(I+1,IGH)
145         IF (HI(I,I-1) .EQ. 0.0E0) GO TO 170
146         NORM = PYTHAG(HR(I,I-1),HI(I,I-1))
147         YR = HR(I,I-1) / NORM
148         YI = HI(I,I-1) / NORM
149         HR(I,I-1) = NORM
150         HI(I,I-1) = 0.0E0
151C
152         DO 155 J = I, N
153            SI = YR * HI(I,J) - YI * HR(I,J)
154            HR(I,J) = YR * HR(I,J) + YI * HI(I,J)
155            HI(I,J) = SI
156  155    CONTINUE
157C
158         DO 160 J = 1, LL
159            SI = YR * HI(J,I) + YI * HR(J,I)
160            HR(J,I) = YR * HR(J,I) - YI * HI(J,I)
161            HI(J,I) = SI
162  160    CONTINUE
163C
164         DO 165 J = LOW, IGH
165            SI = YR * ZI(J,I) + YI * ZR(J,I)
166            ZR(J,I) = YR * ZR(J,I) - YI * ZI(J,I)
167            ZI(J,I) = SI
168  165    CONTINUE
169C
170  170 CONTINUE
171C     .......... STORE ROOTS ISOLATED BY CBAL ..........
172  180 DO 200 I = 1, N
173         IF (I .GE. LOW .AND. I .LE. IGH) GO TO 200
174         WR(I) = HR(I,I)
175         WI(I) = HI(I,I)
176  200 CONTINUE
177C
178      EN = IGH
179      TR = 0.0E0
180      TI = 0.0E0
181      ITN = 30*N
182C     .......... SEARCH FOR NEXT EIGENVALUE ..........
183  220 IF (EN .LT. LOW) GO TO 680
184      ITS = 0
185      ENM1 = EN - 1
186C     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
187C                FOR L=EN STEP -1 UNTIL LOW DO -- ..........
188  240 DO 260 LL = LOW, EN
189         L = EN + LOW - LL
190         IF (L .EQ. LOW) GO TO 300
191         S1 = ABS(HR(L-1,L-1)) + ABS(HI(L-1,L-1))
192     1             + ABS(HR(L,L)) +ABS(HI(L,L))
193         S2 = S1 + ABS(HR(L,L-1))
194         IF (S2 .EQ. S1) GO TO 300
195  260 CONTINUE
196C     .......... FORM SHIFT ..........
197  300 IF (L .EQ. EN) GO TO 660
198      IF (ITN .EQ. 0) GO TO 1000
199      IF (ITS .EQ. 10 .OR. ITS .EQ. 20) GO TO 320
200      SR = HR(EN,EN)
201      SI = HI(EN,EN)
202      XR = HR(ENM1,EN) * HR(EN,ENM1)
203      XI = HI(ENM1,EN) * HR(EN,ENM1)
204      IF (XR .EQ. 0.0E0 .AND. XI .EQ. 0.0E0) GO TO 340
205      YR = (HR(ENM1,ENM1) - SR) / 2.0E0
206      YI = (HI(ENM1,ENM1) - SI) / 2.0E0
207      CALL CSROOT(YR**2-YI**2+XR,2.0E0*YR*YI+XI,ZZR,ZZI)
208      IF (YR * ZZR + YI * ZZI .GE. 0.0E0) GO TO 310
209      ZZR = -ZZR
210      ZZI = -ZZI
211  310 CALL CDIV(XR,XI,YR+ZZR,YI+ZZI,XR,XI)
212      SR = SR - XR
213      SI = SI - XI
214      GO TO 340
215C     .......... FORM EXCEPTIONAL SHIFT ..........
216  320 SR = ABS(HR(EN,ENM1)) + ABS(HR(ENM1,EN-2))
217      SI = 0.0E0
218C
219  340 DO 360 I = LOW, EN
220         HR(I,I) = HR(I,I) - SR
221         HI(I,I) = HI(I,I) - SI
222  360 CONTINUE
223C
224      TR = TR + SR
225      TI = TI + SI
226      ITS = ITS + 1
227      ITN = ITN - 1
228C     .......... REDUCE TO TRIANGLE (ROWS) ..........
229      LP1 = L + 1
230C
231      DO 500 I = LP1, EN
232         SR = HR(I,I-1)
233         HR(I,I-1) = 0.0E0
234         NORM = PYTHAG(PYTHAG(HR(I-1,I-1),HI(I-1,I-1)),SR)
235         XR = HR(I-1,I-1) / NORM
236         WR(I-1) = XR
237         XI = HI(I-1,I-1) / NORM
238         WI(I-1) = XI
239         HR(I-1,I-1) = NORM
240         HI(I-1,I-1) = 0.0E0
241         HI(I,I-1) = SR / NORM
242C
243         DO 490 J = I, N
244            YR = HR(I-1,J)
245            YI = HI(I-1,J)
246            ZZR = HR(I,J)
247            ZZI = HI(I,J)
248            HR(I-1,J) = XR * YR + XI * YI + HI(I,I-1) * ZZR
249            HI(I-1,J) = XR * YI - XI * YR + HI(I,I-1) * ZZI
250            HR(I,J) = XR * ZZR - XI * ZZI - HI(I,I-1) * YR
251            HI(I,J) = XR * ZZI + XI * ZZR - HI(I,I-1) * YI
252  490    CONTINUE
253C
254  500 CONTINUE
255C
256      SI = HI(EN,EN)
257      IF (SI .EQ. 0.0E0) GO TO 540
258      NORM = PYTHAG(HR(EN,EN),SI)
259      SR = HR(EN,EN) / NORM
260      SI = SI / NORM
261      HR(EN,EN) = NORM
262      HI(EN,EN) = 0.0E0
263      IF (EN .EQ. N) GO TO 540
264      IP1 = EN + 1
265C
266      DO 520 J = IP1, N
267         YR = HR(EN,J)
268         YI = HI(EN,J)
269         HR(EN,J) = SR * YR + SI * YI
270         HI(EN,J) = SR * YI - SI * YR
271  520 CONTINUE
272C     .......... INVERSE OPERATION (COLUMNS) ..........
273  540 DO 600 J = LP1, EN
274         XR = WR(J-1)
275         XI = WI(J-1)
276C
277         DO 580 I = 1, J
278            YR = HR(I,J-1)
279            YI = 0.0E0
280            ZZR = HR(I,J)
281            ZZI = HI(I,J)
282            IF (I .EQ. J) GO TO 560
283            YI = HI(I,J-1)
284            HI(I,J-1) = XR * YI + XI * YR + HI(J,J-1) * ZZI
285  560       HR(I,J-1) = XR * YR - XI * YI + HI(J,J-1) * ZZR
286            HR(I,J) = XR * ZZR + XI * ZZI - HI(J,J-1) * YR
287            HI(I,J) = XR * ZZI - XI * ZZR - HI(J,J-1) * YI
288  580    CONTINUE
289C
290         DO 590 I = LOW, IGH
291            YR = ZR(I,J-1)
292            YI = ZI(I,J-1)
293            ZZR = ZR(I,J)
294            ZZI = ZI(I,J)
295            ZR(I,J-1) = XR * YR - XI * YI + HI(J,J-1) * ZZR
296            ZI(I,J-1) = XR * YI + XI * YR + HI(J,J-1) * ZZI
297            ZR(I,J) = XR * ZZR + XI * ZZI - HI(J,J-1) * YR
298            ZI(I,J) = XR * ZZI - XI * ZZR - HI(J,J-1) * YI
299  590    CONTINUE
300C
301  600 CONTINUE
302C
303      IF (SI .EQ. 0.0E0) GO TO 240
304C
305      DO 630 I = 1, EN
306         YR = HR(I,EN)
307         YI = HI(I,EN)
308         HR(I,EN) = SR * YR - SI * YI
309         HI(I,EN) = SR * YI + SI * YR
310  630 CONTINUE
311C
312      DO 640 I = LOW, IGH
313         YR = ZR(I,EN)
314         YI = ZI(I,EN)
315         ZR(I,EN) = SR * YR - SI * YI
316         ZI(I,EN) = SR * YI + SI * YR
317  640 CONTINUE
318C
319      GO TO 240
320C     .......... A ROOT FOUND ..........
321  660 HR(EN,EN) = HR(EN,EN) + TR
322      WR(EN) = HR(EN,EN)
323      HI(EN,EN) = HI(EN,EN) + TI
324      WI(EN) = HI(EN,EN)
325      EN = ENM1
326      GO TO 220
327C     .......... ALL ROOTS FOUND.  BACKSUBSTITUTE TO FIND
328C                VECTORS OF UPPER TRIANGULAR FORM ..........
329  680 NORM = 0.0E0
330C
331      DO 720 I = 1, N
332C
333         DO 720 J = I, N
334            NORM = NORM + ABS(HR(I,J)) + ABS(HI(I,J))
335  720 CONTINUE
336C
337      IF (N .EQ. 1 .OR. NORM .EQ. 0.0E0) GO TO 1001
338C     .......... FOR EN=N STEP -1 UNTIL 2 DO -- ..........
339      DO 800 NN = 2, N
340         EN = N + 2 - NN
341         XR = WR(EN)
342         XI = WI(EN)
343         ENM1 = EN - 1
344C     .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- ..........
345         DO 780 II = 1, ENM1
346            I = EN - II
347            ZZR = HR(I,EN)
348            ZZI = HI(I,EN)
349            IF (I .EQ. ENM1) GO TO 760
350            IP1 = I + 1
351C
352            DO 740 J = IP1, ENM1
353               ZZR = ZZR + HR(I,J) * HR(J,EN) - HI(I,J) * HI(J,EN)
354               ZZI = ZZI + HR(I,J) * HI(J,EN) + HI(I,J) * HR(J,EN)
355  740       CONTINUE
356C
357  760       YR = XR - WR(I)
358            YI = XI - WI(I)
359            IF (YR .NE. 0.0E0 .OR. YI .NE. 0.0E0) GO TO 775
360            YR = NORM
361  770       YR = 0.5E0*YR
362            IF (NORM + YR .GT. NORM) GO TO 770
363            YR = 2.0E0*YR
364  775       CALL CDIV(ZZR,ZZI,YR,YI,HR(I,EN),HI(I,EN))
365  780    CONTINUE
366C
367  800 CONTINUE
368C     .......... END BACKSUBSTITUTION ..........
369      ENM1 = N - 1
370C     .......... VECTORS OF ISOLATED ROOTS ..........
371      DO  840 I = 1, ENM1
372         IF (I .GE. LOW .AND. I .LE. IGH) GO TO 840
373         IP1 = I + 1
374C
375         DO 820 J = IP1, N
376            ZR(I,J) = HR(I,J)
377            ZI(I,J) = HI(I,J)
378  820    CONTINUE
379C
380  840 CONTINUE
381C     .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE
382C                VECTORS OF ORIGINAL FULL MATRIX.
383C                FOR J=N STEP -1 UNTIL LOW+1 DO -- ..........
384      DO 880 JJ = LOW, ENM1
385         J = N + LOW - JJ
386         M = MIN0(J-1,IGH)
387C
388         DO 880 I = LOW, IGH
389            ZZR = ZR(I,J)
390            ZZI = ZI(I,J)
391C
392            DO 860 K = LOW, M
393               ZZR = ZZR + ZR(I,K) * HR(K,J) - ZI(I,K) * HI(K,J)
394               ZZI = ZZI + ZR(I,K) * HI(K,J) + ZI(I,K) * HR(K,J)
395  860       CONTINUE
396C
397            ZR(I,J) = ZZR
398            ZI(I,J) = ZZI
399  880 CONTINUE
400C
401      GO TO 1001
402C     .......... SET ERROR -- NO CONVERGENCE TO AN
403C                EIGENVALUE AFTER 30*N ITERATIONS ..........
404 1000 IERR = EN
405 1001 RETURN
406      END
407