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