LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
zlaqr4.f
Go to the documentation of this file.
1 *> \brief \b ZLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur decomposition.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZLAQR4 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaqr4.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaqr4.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaqr4.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZLAQR4( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
22 * IHIZ, Z, LDZ, WORK, LWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N
26 * LOGICAL WANTT, WANTZ
27 * ..
28 * .. Array Arguments ..
29 * COMPLEX*16 H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> ZLAQR4 implements one level of recursion for ZLAQR0.
39 *> It is a complete implementation of the small bulge multi-shift
40 *> QR algorithm. It may be called by ZLAQR0 and, for large enough
41 *> deflation window size, it may be called by ZLAQR3. This
42 *> subroutine is identical to ZLAQR0 except that it calls ZLAQR2
43 *> instead of ZLAQR3.
44 *>
45 *> ZLAQR4 computes the eigenvalues of a Hessenberg matrix H
46 *> and, optionally, the matrices T and Z from the Schur decomposition
47 *> H = Z T Z**H, where T is an upper triangular matrix (the
48 *> Schur form), and Z is the unitary matrix of Schur vectors.
49 *>
50 *> Optionally Z may be postmultiplied into an input unitary
51 *> matrix Q so that this routine can give the Schur factorization
52 *> of a matrix A which has been reduced to the Hessenberg form H
53 *> by the unitary matrix Q: A = Q*H*Q**H = (QZ)*H*(QZ)**H.
54 *> \endverbatim
55 *
56 * Arguments:
57 * ==========
58 *
59 *> \param[in] WANTT
60 *> \verbatim
61 *> WANTT is LOGICAL
62 *> = .TRUE. : the full Schur form T is required;
63 *> = .FALSE.: only eigenvalues are required.
64 *> \endverbatim
65 *>
66 *> \param[in] WANTZ
67 *> \verbatim
68 *> WANTZ is LOGICAL
69 *> = .TRUE. : the matrix of Schur vectors Z is required;
70 *> = .FALSE.: Schur vectors are not required.
71 *> \endverbatim
72 *>
73 *> \param[in] N
74 *> \verbatim
75 *> N is INTEGER
76 *> The order of the matrix H. N .GE. 0.
77 *> \endverbatim
78 *>
79 *> \param[in] ILO
80 *> \verbatim
81 *> ILO is INTEGER
82 *> \endverbatim
83 *>
84 *> \param[in] IHI
85 *> \verbatim
86 *> IHI is INTEGER
87 *> It is assumed that H is already upper triangular in rows
88 *> and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1,
89 *> H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
90 *> previous call to ZGEBAL, and then passed to ZGEHRD when the
91 *> matrix output by ZGEBAL is reduced to Hessenberg form.
92 *> Otherwise, ILO and IHI should be set to 1 and N,
93 *> respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
94 *> If N = 0, then ILO = 1 and IHI = 0.
95 *> \endverbatim
96 *>
97 *> \param[in,out] H
98 *> \verbatim
99 *> H is COMPLEX*16 array, dimension (LDH,N)
100 *> On entry, the upper Hessenberg matrix H.
101 *> On exit, if INFO = 0 and WANTT is .TRUE., then H
102 *> contains the upper triangular matrix T from the Schur
103 *> decomposition (the Schur form). If INFO = 0 and WANT is
104 *> .FALSE., then the contents of H are unspecified on exit.
105 *> (The output value of H when INFO.GT.0 is given under the
106 *> description of INFO below.)
107 *>
108 *> This subroutine may explicitly set H(i,j) = 0 for i.GT.j and
109 *> j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N.
110 *> \endverbatim
111 *>
112 *> \param[in] LDH
113 *> \verbatim
114 *> LDH is INTEGER
115 *> The leading dimension of the array H. LDH .GE. max(1,N).
116 *> \endverbatim
117 *>
118 *> \param[out] W
119 *> \verbatim
120 *> W is COMPLEX*16 array, dimension (N)
121 *> The computed eigenvalues of H(ILO:IHI,ILO:IHI) are stored
122 *> in W(ILO:IHI). If WANTT is .TRUE., then the eigenvalues are
123 *> stored in the same order as on the diagonal of the Schur
124 *> form returned in H, with W(i) = H(i,i).
125 *> \endverbatim
126 *>
127 *> \param[in] ILOZ
128 *> \verbatim
129 *> ILOZ is INTEGER
130 *> \endverbatim
131 *>
132 *> \param[in] IHIZ
133 *> \verbatim
134 *> IHIZ is INTEGER
135 *> Specify the rows of Z to which transformations must be
136 *> applied if WANTZ is .TRUE..
137 *> 1 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N.
138 *> \endverbatim
139 *>
140 *> \param[in,out] Z
141 *> \verbatim
142 *> Z is COMPLEX*16 array, dimension (LDZ,IHI)
143 *> If WANTZ is .FALSE., then Z is not referenced.
144 *> If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is
145 *> replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the
146 *> orthogonal Schur factor of H(ILO:IHI,ILO:IHI).
147 *> (The output value of Z when INFO.GT.0 is given under
148 *> the description of INFO below.)
149 *> \endverbatim
150 *>
151 *> \param[in] LDZ
152 *> \verbatim
153 *> LDZ is INTEGER
154 *> The leading dimension of the array Z. if WANTZ is .TRUE.
155 *> then LDZ.GE.MAX(1,IHIZ). Otherwize, LDZ.GE.1.
156 *> \endverbatim
157 *>
158 *> \param[out] WORK
159 *> \verbatim
160 *> WORK is COMPLEX*16 array, dimension LWORK
161 *> On exit, if LWORK = -1, WORK(1) returns an estimate of
162 *> the optimal value for LWORK.
163 *> \endverbatim
164 *>
165 *> \param[in] LWORK
166 *> \verbatim
167 *> LWORK is INTEGER
168 *> The dimension of the array WORK. LWORK .GE. max(1,N)
169 *> is sufficient, but LWORK typically as large as 6*N may
170 *> be required for optimal performance. A workspace query
171 *> to determine the optimal workspace size is recommended.
172 *>
173 *> If LWORK = -1, then ZLAQR4 does a workspace query.
174 *> In this case, ZLAQR4 checks the input parameters and
175 *> estimates the optimal workspace size for the given
176 *> values of N, ILO and IHI. The estimate is returned
177 *> in WORK(1). No error message related to LWORK is
178 *> issued by XERBLA. Neither H nor Z are accessed.
179 *> \endverbatim
180 *>
181 *> \param[out] INFO
182 *> \verbatim
183 *> INFO is INTEGER
184 *> = 0: successful exit
185 *> .GT. 0: if INFO = i, ZLAQR4 failed to compute all of
186 *> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
187 *> and WI contain those eigenvalues which have been
188 *> successfully computed. (Failures are rare.)
189 *>
190 *> If INFO .GT. 0 and WANT is .FALSE., then on exit,
191 *> the remaining unconverged eigenvalues are the eigen-
192 *> values of the upper Hessenberg matrix rows and
193 *> columns ILO through INFO of the final, output
194 *> value of H.
195 *>
196 *> If INFO .GT. 0 and WANTT is .TRUE., then on exit
197 *>
198 *> (*) (initial value of H)*U = U*(final value of H)
199 *>
200 *> where U is a unitary matrix. The final
201 *> value of H is upper Hessenberg and triangular in
202 *> rows and columns INFO+1 through IHI.
203 *>
204 *> If INFO .GT. 0 and WANTZ is .TRUE., then on exit
205 *>
206 *> (final value of Z(ILO:IHI,ILOZ:IHIZ)
207 *> = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U
208 *>
209 *> where U is the unitary matrix in (*) (regard-
210 *> less of the value of WANTT.)
211 *>
212 *> If INFO .GT. 0 and WANTZ is .FALSE., then Z is not
213 *> accessed.
214 *> \endverbatim
215 *
216 * Authors:
217 * ========
218 *
219 *> \author Univ. of Tennessee
220 *> \author Univ. of California Berkeley
221 *> \author Univ. of Colorado Denver
222 *> \author NAG Ltd.
223 *
224 *> \date September 2012
225 *
226 *> \ingroup complex16OTHERauxiliary
227 *
228 *> \par Contributors:
229 * ==================
230 *>
231 *> Karen Braman and Ralph Byers, Department of Mathematics,
232 *> University of Kansas, USA
233 *
234 *> \par References:
235 * ================
236 *>
237 *> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
238 *> Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
239 *> Performance, SIAM Journal of Matrix Analysis, volume 23, pages
240 *> 929--947, 2002.
241 *> \n
242 *> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
243 *> Algorithm Part II: Aggressive Early Deflation, SIAM Journal
244 *> of Matrix Analysis, volume 23, pages 948--973, 2002.
245 *>
246 * =====================================================================
247  SUBROUTINE zlaqr4( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
248  $ ihiz, z, ldz, work, lwork, info )
249 *
250 * -- LAPACK auxiliary routine (version 3.4.2) --
251 * -- LAPACK is a software package provided by Univ. of Tennessee, --
252 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
253 * September 2012
254 *
255 * .. Scalar Arguments ..
256  INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N
257  LOGICAL WANTT, WANTZ
258 * ..
259 * .. Array Arguments ..
260  COMPLEX*16 H( ldh, * ), W( * ), WORK( * ), Z( ldz, * )
261 * ..
262 *
263 * ================================================================
264 *
265 * .. Parameters ..
266 *
267 * ==== Matrices of order NTINY or smaller must be processed by
268 * . ZLAHQR because of insufficient subdiagonal scratch space.
269 * . (This is a hard limit.) ====
270  INTEGER NTINY
271  parameter ( ntiny = 11 )
272 *
273 * ==== Exceptional deflation windows: try to cure rare
274 * . slow convergence by varying the size of the
275 * . deflation window after KEXNW iterations. ====
276  INTEGER KEXNW
277  parameter ( kexnw = 5 )
278 *
279 * ==== Exceptional shifts: try to cure rare slow convergence
280 * . with ad-hoc exceptional shifts every KEXSH iterations.
281 * . ====
282  INTEGER KEXSH
283  parameter ( kexsh = 6 )
284 *
285 * ==== The constant WILK1 is used to form the exceptional
286 * . shifts. ====
287  DOUBLE PRECISION WILK1
288  parameter ( wilk1 = 0.75d0 )
289  COMPLEX*16 ZERO, ONE
290  parameter ( zero = ( 0.0d0, 0.0d0 ),
291  $ one = ( 1.0d0, 0.0d0 ) )
292  DOUBLE PRECISION TWO
293  parameter ( two = 2.0d0 )
294 * ..
295 * .. Local Scalars ..
296  COMPLEX*16 AA, BB, CC, CDUM, DD, DET, RTDISC, SWAP, TR2
297  DOUBLE PRECISION S
298  INTEGER I, INF, IT, ITMAX, K, KACC22, KBOT, KDU, KS,
299  $ kt, ktop, ku, kv, kwh, kwtop, kwv, ld, ls,
300  $ lwkopt, ndec, ndfl, nh, nho, nibble, nmin, ns,
301  $ nsmax, nsr, nve, nw, nwmax, nwr, nwupbd
302  LOGICAL SORTED
303  CHARACTER JBCMPZ*2
304 * ..
305 * .. External Functions ..
306  INTEGER ILAENV
307  EXTERNAL ilaenv
308 * ..
309 * .. Local Arrays ..
310  COMPLEX*16 ZDUM( 1, 1 )
311 * ..
312 * .. External Subroutines ..
313  EXTERNAL zlacpy, zlahqr, zlaqr2, zlaqr5
314 * ..
315 * .. Intrinsic Functions ..
316  INTRINSIC abs, dble, dcmplx, dimag, int, max, min, mod,
317  $ sqrt
318 * ..
319 * .. Statement Functions ..
320  DOUBLE PRECISION CABS1
321 * ..
322 * .. Statement Function definitions ..
323  cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
324 * ..
325 * .. Executable Statements ..
326  info = 0
327 *
328 * ==== Quick return for N = 0: nothing to do. ====
329 *
330  IF( n.EQ.0 ) THEN
331  work( 1 ) = one
332  RETURN
333  END IF
334 *
335  IF( n.LE.ntiny ) THEN
336 *
337 * ==== Tiny matrices must use ZLAHQR. ====
338 *
339  lwkopt = 1
340  IF( lwork.NE.-1 )
341  $ CALL zlahqr( wantt, wantz, n, ilo, ihi, h, ldh, w, iloz,
342  $ ihiz, z, ldz, info )
343  ELSE
344 *
345 * ==== Use small bulge multi-shift QR with aggressive early
346 * . deflation on larger-than-tiny matrices. ====
347 *
348 * ==== Hope for the best. ====
349 *
350  info = 0
351 *
352 * ==== Set up job flags for ILAENV. ====
353 *
354  IF( wantt ) THEN
355  jbcmpz( 1: 1 ) = 'S'
356  ELSE
357  jbcmpz( 1: 1 ) = 'E'
358  END IF
359  IF( wantz ) THEN
360  jbcmpz( 2: 2 ) = 'V'
361  ELSE
362  jbcmpz( 2: 2 ) = 'N'
363  END IF
364 *
365 * ==== NWR = recommended deflation window size. At this
366 * . point, N .GT. NTINY = 11, so there is enough
367 * . subdiagonal workspace for NWR.GE.2 as required.
368 * . (In fact, there is enough subdiagonal space for
369 * . NWR.GE.3.) ====
370 *
371  nwr = ilaenv( 13, 'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
372  nwr = max( 2, nwr )
373  nwr = min( ihi-ilo+1, ( n-1 ) / 3, nwr )
374 *
375 * ==== NSR = recommended number of simultaneous shifts.
376 * . At this point N .GT. NTINY = 11, so there is at
377 * . enough subdiagonal workspace for NSR to be even
378 * . and greater than or equal to two as required. ====
379 *
380  nsr = ilaenv( 15, 'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
381  nsr = min( nsr, ( n+6 ) / 9, ihi-ilo )
382  nsr = max( 2, nsr-mod( nsr, 2 ) )
383 *
384 * ==== Estimate optimal workspace ====
385 *
386 * ==== Workspace query call to ZLAQR2 ====
387 *
388  CALL zlaqr2( wantt, wantz, n, ilo, ihi, nwr+1, h, ldh, iloz,
389  $ ihiz, z, ldz, ls, ld, w, h, ldh, n, h, ldh, n, h,
390  $ ldh, work, -1 )
391 *
392 * ==== Optimal workspace = MAX(ZLAQR5, ZLAQR2) ====
393 *
394  lwkopt = max( 3*nsr / 2, int( work( 1 ) ) )
395 *
396 * ==== Quick return in case of workspace query. ====
397 *
398  IF( lwork.EQ.-1 ) THEN
399  work( 1 ) = dcmplx( lwkopt, 0 )
400  RETURN
401  END IF
402 *
403 * ==== ZLAHQR/ZLAQR0 crossover point ====
404 *
405  nmin = ilaenv( 12, 'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
406  nmin = max( ntiny, nmin )
407 *
408 * ==== Nibble crossover point ====
409 *
410  nibble = ilaenv( 14, 'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
411  nibble = max( 0, nibble )
412 *
413 * ==== Accumulate reflections during ttswp? Use block
414 * . 2-by-2 structure during matrix-matrix multiply? ====
415 *
416  kacc22 = ilaenv( 16, 'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
417  kacc22 = max( 0, kacc22 )
418  kacc22 = min( 2, kacc22 )
419 *
420 * ==== NWMAX = the largest possible deflation window for
421 * . which there is sufficient workspace. ====
422 *
423  nwmax = min( ( n-1 ) / 3, lwork / 2 )
424  nw = nwmax
425 *
426 * ==== NSMAX = the Largest number of simultaneous shifts
427 * . for which there is sufficient workspace. ====
428 *
429  nsmax = min( ( n+6 ) / 9, 2*lwork / 3 )
430  nsmax = nsmax - mod( nsmax, 2 )
431 *
432 * ==== NDFL: an iteration count restarted at deflation. ====
433 *
434  ndfl = 1
435 *
436 * ==== ITMAX = iteration limit ====
437 *
438  itmax = max( 30, 2*kexsh )*max( 10, ( ihi-ilo+1 ) )
439 *
440 * ==== Last row and column in the active block ====
441 *
442  kbot = ihi
443 *
444 * ==== Main Loop ====
445 *
446  DO 70 it = 1, itmax
447 *
448 * ==== Done when KBOT falls below ILO ====
449 *
450  IF( kbot.LT.ilo )
451  $ GO TO 80
452 *
453 * ==== Locate active block ====
454 *
455  DO 10 k = kbot, ilo + 1, -1
456  IF( h( k, k-1 ).EQ.zero )
457  $ GO TO 20
458  10 CONTINUE
459  k = ilo
460  20 CONTINUE
461  ktop = k
462 *
463 * ==== Select deflation window size:
464 * . Typical Case:
465 * . If possible and advisable, nibble the entire
466 * . active block. If not, use size MIN(NWR,NWMAX)
467 * . or MIN(NWR+1,NWMAX) depending upon which has
468 * . the smaller corresponding subdiagonal entry
469 * . (a heuristic).
470 * .
471 * . Exceptional Case:
472 * . If there have been no deflations in KEXNW or
473 * . more iterations, then vary the deflation window
474 * . size. At first, because, larger windows are,
475 * . in general, more powerful than smaller ones,
476 * . rapidly increase the window to the maximum possible.
477 * . Then, gradually reduce the window size. ====
478 *
479  nh = kbot - ktop + 1
480  nwupbd = min( nh, nwmax )
481  IF( ndfl.LT.kexnw ) THEN
482  nw = min( nwupbd, nwr )
483  ELSE
484  nw = min( nwupbd, 2*nw )
485  END IF
486  IF( nw.LT.nwmax ) THEN
487  IF( nw.GE.nh-1 ) THEN
488  nw = nh
489  ELSE
490  kwtop = kbot - nw + 1
491  IF( cabs1( h( kwtop, kwtop-1 ) ).GT.
492  $ cabs1( h( kwtop-1, kwtop-2 ) ) )nw = nw + 1
493  END IF
494  END IF
495  IF( ndfl.LT.kexnw ) THEN
496  ndec = -1
497  ELSE IF( ndec.GE.0 .OR. nw.GE.nwupbd ) THEN
498  ndec = ndec + 1
499  IF( nw-ndec.LT.2 )
500  $ ndec = 0
501  nw = nw - ndec
502  END IF
503 *
504 * ==== Aggressive early deflation:
505 * . split workspace under the subdiagonal into
506 * . - an nw-by-nw work array V in the lower
507 * . left-hand-corner,
508 * . - an NW-by-at-least-NW-but-more-is-better
509 * . (NW-by-NHO) horizontal work array along
510 * . the bottom edge,
511 * . - an at-least-NW-but-more-is-better (NHV-by-NW)
512 * . vertical work array along the left-hand-edge.
513 * . ====
514 *
515  kv = n - nw + 1
516  kt = nw + 1
517  nho = ( n-nw-1 ) - kt + 1
518  kwv = nw + 2
519  nve = ( n-nw ) - kwv + 1
520 *
521 * ==== Aggressive early deflation ====
522 *
523  CALL zlaqr2( wantt, wantz, n, ktop, kbot, nw, h, ldh, iloz,
524  $ ihiz, z, ldz, ls, ld, w, h( kv, 1 ), ldh, nho,
525  $ h( kv, kt ), ldh, nve, h( kwv, 1 ), ldh, work,
526  $ lwork )
527 *
528 * ==== Adjust KBOT accounting for new deflations. ====
529 *
530  kbot = kbot - ld
531 *
532 * ==== KS points to the shifts. ====
533 *
534  ks = kbot - ls + 1
535 *
536 * ==== Skip an expensive QR sweep if there is a (partly
537 * . heuristic) reason to expect that many eigenvalues
538 * . will deflate without it. Here, the QR sweep is
539 * . skipped if many eigenvalues have just been deflated
540 * . or if the remaining active block is small.
541 *
542  IF( ( ld.EQ.0 ) .OR. ( ( 100*ld.LE.nw*nibble ) .AND. ( kbot-
543  $ ktop+1.GT.min( nmin, nwmax ) ) ) ) THEN
544 *
545 * ==== NS = nominal number of simultaneous shifts.
546 * . This may be lowered (slightly) if ZLAQR2
547 * . did not provide that many shifts. ====
548 *
549  ns = min( nsmax, nsr, max( 2, kbot-ktop ) )
550  ns = ns - mod( ns, 2 )
551 *
552 * ==== If there have been no deflations
553 * . in a multiple of KEXSH iterations,
554 * . then try exceptional shifts.
555 * . Otherwise use shifts provided by
556 * . ZLAQR2 above or from the eigenvalues
557 * . of a trailing principal submatrix. ====
558 *
559  IF( mod( ndfl, kexsh ).EQ.0 ) THEN
560  ks = kbot - ns + 1
561  DO 30 i = kbot, ks + 1, -2
562  w( i ) = h( i, i ) + wilk1*cabs1( h( i, i-1 ) )
563  w( i-1 ) = w( i )
564  30 CONTINUE
565  ELSE
566 *
567 * ==== Got NS/2 or fewer shifts? Use ZLAHQR
568 * . on a trailing principal submatrix to
569 * . get more. (Since NS.LE.NSMAX.LE.(N+6)/9,
570 * . there is enough space below the subdiagonal
571 * . to fit an NS-by-NS scratch array.) ====
572 *
573  IF( kbot-ks+1.LE.ns / 2 ) THEN
574  ks = kbot - ns + 1
575  kt = n - ns + 1
576  CALL zlacpy( 'A', ns, ns, h( ks, ks ), ldh,
577  $ h( kt, 1 ), ldh )
578  CALL zlahqr( .false., .false., ns, 1, ns,
579  $ h( kt, 1 ), ldh, w( ks ), 1, 1, zdum,
580  $ 1, inf )
581  ks = ks + inf
582 *
583 * ==== In case of a rare QR failure use
584 * . eigenvalues of the trailing 2-by-2
585 * . principal submatrix. Scale to avoid
586 * . overflows, underflows and subnormals.
587 * . (The scale factor S can not be zero,
588 * . because H(KBOT,KBOT-1) is nonzero.) ====
589 *
590  IF( ks.GE.kbot ) THEN
591  s = cabs1( h( kbot-1, kbot-1 ) ) +
592  $ cabs1( h( kbot, kbot-1 ) ) +
593  $ cabs1( h( kbot-1, kbot ) ) +
594  $ cabs1( h( kbot, kbot ) )
595  aa = h( kbot-1, kbot-1 ) / s
596  cc = h( kbot, kbot-1 ) / s
597  bb = h( kbot-1, kbot ) / s
598  dd = h( kbot, kbot ) / s
599  tr2 = ( aa+dd ) / two
600  det = ( aa-tr2 )*( dd-tr2 ) - bb*cc
601  rtdisc = sqrt( -det )
602  w( kbot-1 ) = ( tr2+rtdisc )*s
603  w( kbot ) = ( tr2-rtdisc )*s
604 *
605  ks = kbot - 1
606  END IF
607  END IF
608 *
609  IF( kbot-ks+1.GT.ns ) THEN
610 *
611 * ==== Sort the shifts (Helps a little) ====
612 *
613  sorted = .false.
614  DO 50 k = kbot, ks + 1, -1
615  IF( sorted )
616  $ GO TO 60
617  sorted = .true.
618  DO 40 i = ks, k - 1
619  IF( cabs1( w( i ) ).LT.cabs1( w( i+1 ) ) )
620  $ THEN
621  sorted = .false.
622  swap = w( i )
623  w( i ) = w( i+1 )
624  w( i+1 ) = swap
625  END IF
626  40 CONTINUE
627  50 CONTINUE
628  60 CONTINUE
629  END IF
630  END IF
631 *
632 * ==== If there are only two shifts, then use
633 * . only one. ====
634 *
635  IF( kbot-ks+1.EQ.2 ) THEN
636  IF( cabs1( w( kbot )-h( kbot, kbot ) ).LT.
637  $ cabs1( w( kbot-1 )-h( kbot, kbot ) ) ) THEN
638  w( kbot-1 ) = w( kbot )
639  ELSE
640  w( kbot ) = w( kbot-1 )
641  END IF
642  END IF
643 *
644 * ==== Use up to NS of the the smallest magnatiude
645 * . shifts. If there aren't NS shifts available,
646 * . then use them all, possibly dropping one to
647 * . make the number of shifts even. ====
648 *
649  ns = min( ns, kbot-ks+1 )
650  ns = ns - mod( ns, 2 )
651  ks = kbot - ns + 1
652 *
653 * ==== Small-bulge multi-shift QR sweep:
654 * . split workspace under the subdiagonal into
655 * . - a KDU-by-KDU work array U in the lower
656 * . left-hand-corner,
657 * . - a KDU-by-at-least-KDU-but-more-is-better
658 * . (KDU-by-NHo) horizontal work array WH along
659 * . the bottom edge,
660 * . - and an at-least-KDU-but-more-is-better-by-KDU
661 * . (NVE-by-KDU) vertical work WV arrow along
662 * . the left-hand-edge. ====
663 *
664  kdu = 3*ns - 3
665  ku = n - kdu + 1
666  kwh = kdu + 1
667  nho = ( n-kdu+1-4 ) - ( kdu+1 ) + 1
668  kwv = kdu + 4
669  nve = n - kdu - kwv + 1
670 *
671 * ==== Small-bulge multi-shift QR sweep ====
672 *
673  CALL zlaqr5( wantt, wantz, kacc22, n, ktop, kbot, ns,
674  $ w( ks ), h, ldh, iloz, ihiz, z, ldz, work,
675  $ 3, h( ku, 1 ), ldh, nve, h( kwv, 1 ), ldh,
676  $ nho, h( ku, kwh ), ldh )
677  END IF
678 *
679 * ==== Note progress (or the lack of it). ====
680 *
681  IF( ld.GT.0 ) THEN
682  ndfl = 1
683  ELSE
684  ndfl = ndfl + 1
685  END IF
686 *
687 * ==== End of main loop ====
688  70 CONTINUE
689 *
690 * ==== Iteration limit exceeded. Set INFO to show where
691 * . the problem occurred and exit. ====
692 *
693  info = kbot
694  80 CONTINUE
695  END IF
696 *
697 * ==== Return the optimal value of LWORK. ====
698 *
699  work( 1 ) = dcmplx( lwkopt, 0 )
700 *
701 * ==== End of ZLAQR4 ====
702 *
703  END
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine zlaqr4(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)
ZLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...
Definition: zlaqr4.f:249
subroutine zlaqr5(WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, S, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV, WV, LDWV, NH, WH, LDWH)
ZLAQR5 performs a single small-bulge multi-shift QR sweep.
Definition: zlaqr5.f:253
subroutine zlaqr2(WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SH, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK)
ZLAQR2 performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fu...
Definition: zlaqr2.f:272
subroutine zlahqr(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, INFO)
ZLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.
Definition: zlahqr.f:197