LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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
subroutine clahqr(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, INFO)
CLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.
Definition: clahqr.f:197
subroutine claqr2(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)
CLAQR2 performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fu...
Definition: claqr2.f:271
subroutine claqr4(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)
CLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...
Definition: claqr4.f:251
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
subroutine claqr5(WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, S, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV, WV, LDWV, NH, WH, LDWH)
CLAQR5 performs a single small-bulge multi-shift QR sweep.
Definition: claqr5.f:253