1      SUBROUTINE COMQR(NM,N,LOW,IGH,HR,HI,WR,WI,IERR)
2C***BEGIN PROLOGUE  COMQR
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 of complex upper Hessenberg matrix
9C            using the QR method.
10C***DESCRIPTION
11C
12C     This subroutine is a translation of a unitary analogue of the
13C     ALGOL procedure  COMLR, NUM. MATH. 12, 369-376(1968) by Martin
14C     and Wilkinson.
15C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 396-403(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 of a COMPLEX
20C     upper Hessenberg matrix by the QR method.
21C
22C     On INPUT
23C
24C        NM must be set to the row dimension of two-dimensional
25C          array parameters as declared in the calling program
26C          dimension statement.
27C
28C        N is the order of the matrix.
29C
30C        LOW and IGH are integers determined by the balancing
31C          subroutine  CBAL.  If  CBAL  has not been used,
32C          set LOW=1, IGH=N.
33C
34C        HR and HI contain the real and imaginary parts,
35C          respectively, of the complex upper Hessenberg matrix.
36C          Their lower triangles below the subdiagonal contain
37C          information about the unitary transformations used in
38C          the reduction by  CORTH, if performed.
39C
40C     On OUTPUT
41C
42C        The upper Hessenberg portions of HR and HI have been
43C          destroyed.  Therefore, they must be saved before
44C          calling  COMQR  if subsequent calculation of
45C          eigenvectors is to be performed.
46C
47C        WR and WI contain the real and imaginary parts,
48C          respectively, of the eigenvalues.  If an error
49C          exit is made, the eigenvalues should be correct
50C          for indices IERR+1,...,N.
51C
52C        IERR is set to
53C          ZERO       for normal return,
54C          J          if the J-th eigenvalue has not been
55C                     determined after a total of 30*N iterations.
56C
57C     Calls CSROOT for complex square root.
58C     Calls PYTHAG(A,B) for sqrt(A**2 + B**2).
59C     Calls CDIV for complex division.
60C
61C     Questions and comments should be directed to B. S. Garbow,
62C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
63C     ------------------------------------------------------------------
64C***REFERENCES  B. T. SMITH, J. M. BOYLE, J. J. DONGARRA, B. S. GARBOW,
65C                 Y. IKEBE, V. C. KLEMA, C. B. MOLER, *MATRIX EIGEN-
66C                 SYSTEM ROUTINES - EISPACK GUIDE*, SPRINGER-VERLAG,
67C                 1976.
68C***ROUTINES CALLED  CDIV,CSROOT,PYTHAG
69C***END PROLOGUE  COMQR
70C
71      INTEGER I,J,L,N,EN,LL,NM,IGH,ITN,ITS,LOW,LP1,ENM1,IERR
72      REAL HR(NM,N),HI(NM,N),WR(N),WI(N)
73      REAL SI,SR,TI,TR,XI,XR,YI,YR,ZZI,ZZR,NORM,S1,S2
74      REAL PYTHAG
75C
76C***FIRST EXECUTABLE STATEMENT  COMQR
77      IERR = 0
78      IF (LOW .EQ. IGH) GO TO 180
79C     .......... CREATE REAL SUBDIAGONAL ELEMENTS ..........
80      L = LOW + 1
81C
82      DO 170 I = L, IGH
83         LL = MIN0(I+1,IGH)
84         IF (HI(I,I-1) .EQ. 0.0E0) GO TO 170
85         NORM = PYTHAG(HR(I,I-1),HI(I,I-1))
86         YR = HR(I,I-1) / NORM
87         YI = HI(I,I-1) / NORM
88         HR(I,I-1) = NORM
89         HI(I,I-1) = 0.0E0
90C
91         DO 155 J = I, IGH
92            SI = YR * HI(I,J) - YI * HR(I,J)
93            HR(I,J) = YR * HR(I,J) + YI * HI(I,J)
94            HI(I,J) = SI
95  155    CONTINUE
96C
97         DO 160 J = LOW, LL
98            SI = YR * HI(J,I) + YI * HR(J,I)
99            HR(J,I) = YR * HR(J,I) - YI * HI(J,I)
100            HI(J,I) = SI
101  160    CONTINUE
102C
103  170 CONTINUE
104C     .......... STORE ROOTS ISOLATED BY CBAL ..........
105  180 DO 200 I = 1, N
106         IF (I .GE. LOW .AND. I .LE. IGH) GO TO 200
107         WR(I) = HR(I,I)
108         WI(I) = HI(I,I)
109  200 CONTINUE
110C
111      EN = IGH
112      TR = 0.0E0
113      TI = 0.0E0
114      ITN = 30*N
115C     .......... SEARCH FOR NEXT EIGENVALUE ..........
116  220 IF (EN .LT. LOW) GO TO 1001
117      ITS = 0
118      ENM1 = EN - 1
119C     .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
120C                FOR L=EN STEP -1 UNTIL LOW E0 -- ..........
121  240 DO 260 LL = LOW, EN
122         L = EN + LOW - LL
123         IF (L .EQ. LOW) GO TO 300
124         S1 = ABS(HR(L-1,L-1)) + ABS(HI(L-1,L-1))
125     1             + ABS(HR(L,L)) +ABS(HI(L,L))
126         S2 = S1 + ABS(HR(L,L-1))
127         IF (S2 .EQ. S1) GO TO 300
128  260 CONTINUE
129C     .......... FORM SHIFT ..........
130  300 IF (L .EQ. EN) GO TO 660
131      IF (ITN .EQ. 0) GO TO 1000
132      IF (ITS .EQ. 10 .OR. ITS .EQ. 20) GO TO 320
133      SR = HR(EN,EN)
134      SI = HI(EN,EN)
135      XR = HR(ENM1,EN) * HR(EN,ENM1)
136      XI = HI(ENM1,EN) * HR(EN,ENM1)
137      IF (XR .EQ. 0.0E0 .AND. XI .EQ. 0.0E0) GO TO 340
138      YR = (HR(ENM1,ENM1) - SR) / 2.0E0
139      YI = (HI(ENM1,ENM1) - SI) / 2.0E0
140      CALL CSROOT(YR**2-YI**2+XR,2.0E0*YR*YI+XI,ZZR,ZZI)
141      IF (YR * ZZR + YI * ZZI .GE. 0.0E0) GO TO 310
142      ZZR = -ZZR
143      ZZI = -ZZI
144  310 CALL CDIV(XR,XI,YR+ZZR,YI+ZZI,XR,XI)
145      SR = SR - XR
146      SI = SI - XI
147      GO TO 340
148C     .......... FORM EXCEPTIONAL SHIFT ..........
149  320 SR = ABS(HR(EN,ENM1)) + ABS(HR(ENM1,EN-2))
150      SI = 0.0E0
151C
152  340 DO 360 I = LOW, EN
153         HR(I,I) = HR(I,I) - SR
154         HI(I,I) = HI(I,I) - SI
155  360 CONTINUE
156C
157      TR = TR + SR
158      TI = TI + SI
159      ITS = ITS + 1
160      ITN = ITN - 1
161C     .......... REDUCE TO TRIANGLE (ROWS) ..........
162      LP1 = L + 1
163C
164      DO 500 I = LP1, EN
165         SR = HR(I,I-1)
166         HR(I,I-1) = 0.0E0
167         NORM = PYTHAG(PYTHAG(HR(I-1,I-1),HI(I-1,I-1)),SR)
168         XR = HR(I-1,I-1) / NORM
169         WR(I-1) = XR
170         XI = HI(I-1,I-1) / NORM
171         WI(I-1) = XI
172         HR(I-1,I-1) = NORM
173         HI(I-1,I-1) = 0.0E0
174         HI(I,I-1) = SR / NORM
175C
176         DO 490 J = I, EN
177            YR = HR(I-1,J)
178            YI = HI(I-1,J)
179            ZZR = HR(I,J)
180            ZZI = HI(I,J)
181            HR(I-1,J) = XR * YR + XI * YI + HI(I,I-1) * ZZR
182            HI(I-1,J) = XR * YI - XI * YR + HI(I,I-1) * ZZI
183            HR(I,J) = XR * ZZR - XI * ZZI - HI(I,I-1) * YR
184            HI(I,J) = XR * ZZI + XI * ZZR - HI(I,I-1) * YI
185  490    CONTINUE
186C
187  500 CONTINUE
188C
189      SI = HI(EN,EN)
190      IF (SI .EQ. 0.0E0) GO TO 540
191      NORM = PYTHAG(HR(EN,EN),SI)
192      SR = HR(EN,EN) / NORM
193      SI = SI / NORM
194      HR(EN,EN) = NORM
195      HI(EN,EN) = 0.0E0
196C     .......... INVERSE OPERATION (COLUMNS) ..........
197  540 DO 600 J = LP1, EN
198         XR = WR(J-1)
199         XI = WI(J-1)
200C
201         DO 580 I = L, J
202            YR = HR(I,J-1)
203            YI = 0.0E0
204            ZZR = HR(I,J)
205            ZZI = HI(I,J)
206            IF (I .EQ. J) GO TO 560
207            YI = HI(I,J-1)
208            HI(I,J-1) = XR * YI + XI * YR + HI(J,J-1) * ZZI
209  560       HR(I,J-1) = XR * YR - XI * YI + HI(J,J-1) * ZZR
210            HR(I,J) = XR * ZZR + XI * ZZI - HI(J,J-1) * YR
211            HI(I,J) = XR * ZZI - XI * ZZR - HI(J,J-1) * YI
212  580    CONTINUE
213C
214  600 CONTINUE
215C
216      IF (SI .EQ. 0.0E0) GO TO 240
217C
218      DO 630 I = L, EN
219         YR = HR(I,EN)
220         YI = HI(I,EN)
221         HR(I,EN) = SR * YR - SI * YI
222         HI(I,EN) = SR * YI + SI * YR
223  630 CONTINUE
224C
225      GO TO 240
226C     .......... A ROOT FOUND ..........
227  660 WR(EN) = HR(EN,EN) + TR
228      WI(EN) = HI(EN,EN) + TI
229      EN = ENM1
230      GO TO 220
231C     .......... SET ERROR -- NO CONVERGENCE TO AN
232C                EIGENVALUE AFTER 30*N ITERATIONS ..........
233 1000 IERR = EN
234 1001 RETURN
235      END
236