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