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