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