LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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