LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 >= 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 > 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 > 0, then 1 <= ILO <= IHI <= 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 > 0 is given under the
99 *> description of INFO below.)
100 *>
101 *> This subroutine may explicitly set H(i,j) = 0 for i > 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 >= 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 <= ILOZ <= ILO; IHI <= IHIZ <= 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 > 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 >= MAX(1,IHIZ). Otherwise, LDZ >= 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 >= 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 *> > 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 > 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 > 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 > 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 > 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 *> \ingroup complexOTHERauxiliary
218 *
219 *> \par Contributors:
220 * ==================
221 *>
222 *> Karen Braman and Ralph Byers, Department of Mathematics,
223 *> University of Kansas, USA
224 *
225 *> \par References:
226 * ================
227 *>
228 *> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
229 *> Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
230 *> Performance, SIAM Journal of Matrix Analysis, volume 23, pages
231 *> 929--947, 2002.
232 *> \n
233 *> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
234 *> Algorithm Part II: Aggressive Early Deflation, SIAM Journal
235 *> of Matrix Analysis, volume 23, pages 948--973, 2002.
236 *>
237 * =====================================================================
238  SUBROUTINE claqr0( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
239  $ IHIZ, Z, LDZ, WORK, LWORK, INFO )
240 *
241 * -- LAPACK auxiliary routine --
242 * -- LAPACK is a software package provided by Univ. of Tennessee, --
243 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
244 *
245 * .. Scalar Arguments ..
246  INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N
247  LOGICAL WANTT, WANTZ
248 * ..
249 * .. Array Arguments ..
250  COMPLEX H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * )
251 * ..
252 *
253 * ================================================================
254 * .. Parameters ..
255 *
256 * ==== Matrices of order NTINY or smaller must be processed by
257 * . CLAHQR because of insufficient subdiagonal scratch space.
258 * . (This is a hard limit.) ====
259  INTEGER NTINY
260  parameter( ntiny = 15 )
261 *
262 * ==== Exceptional deflation windows: try to cure rare
263 * . slow convergence by varying the size of the
264 * . deflation window after KEXNW iterations. ====
265  INTEGER KEXNW
266  parameter( kexnw = 5 )
267 *
268 * ==== Exceptional shifts: try to cure rare slow convergence
269 * . with ad-hoc exceptional shifts every KEXSH iterations.
270 * . ====
271  INTEGER KEXSH
272  parameter( kexsh = 6 )
273 *
274 * ==== The constant WILK1 is used to form the exceptional
275 * . shifts. ====
276  REAL WILK1
277  parameter( wilk1 = 0.75e0 )
278  COMPLEX ZERO, ONE
279  parameter( zero = ( 0.0e0, 0.0e0 ),
280  $ one = ( 1.0e0, 0.0e0 ) )
281  REAL TWO
282  parameter( two = 2.0e0 )
283 * ..
284 * .. Local Scalars ..
285  COMPLEX AA, BB, CC, CDUM, DD, DET, RTDISC, SWAP, TR2
286  REAL S
287  INTEGER I, INF, IT, ITMAX, K, KACC22, KBOT, KDU, KS,
288  $ kt, ktop, ku, kv, kwh, kwtop, kwv, ld, ls,
289  $ lwkopt, ndec, ndfl, nh, nho, nibble, nmin, ns,
290  $ nsmax, nsr, nve, nw, nwmax, nwr, nwupbd
291  LOGICAL SORTED
292  CHARACTER JBCMPZ*2
293 * ..
294 * .. External Functions ..
295  INTEGER ILAENV
296  EXTERNAL ilaenv
297 * ..
298 * .. Local Arrays ..
299  COMPLEX ZDUM( 1, 1 )
300 * ..
301 * .. External Subroutines ..
302  EXTERNAL clacpy, clahqr, claqr3, claqr4, claqr5
303 * ..
304 * .. Intrinsic Functions ..
305  INTRINSIC abs, aimag, cmplx, int, max, min, mod, real,
306  $ sqrt
307 * ..
308 * .. Statement Functions ..
309  REAL CABS1
310 * ..
311 * .. Statement Function definitions ..
312  cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
313 * ..
314 * .. Executable Statements ..
315  info = 0
316 *
317 * ==== Quick return for N = 0: nothing to do. ====
318 *
319  IF( n.EQ.0 ) THEN
320  work( 1 ) = one
321  RETURN
322  END IF
323 *
324  IF( n.LE.ntiny ) THEN
325 *
326 * ==== Tiny matrices must use CLAHQR. ====
327 *
328  lwkopt = 1
329  IF( lwork.NE.-1 )
330  $ CALL clahqr( wantt, wantz, n, ilo, ihi, h, ldh, w, iloz,
331  $ ihiz, z, ldz, info )
332  ELSE
333 *
334 * ==== Use small bulge multi-shift QR with aggressive early
335 * . deflation on larger-than-tiny matrices. ====
336 *
337 * ==== Hope for the best. ====
338 *
339  info = 0
340 *
341 * ==== Set up job flags for ILAENV. ====
342 *
343  IF( wantt ) THEN
344  jbcmpz( 1: 1 ) = 'S'
345  ELSE
346  jbcmpz( 1: 1 ) = 'E'
347  END IF
348  IF( wantz ) THEN
349  jbcmpz( 2: 2 ) = 'V'
350  ELSE
351  jbcmpz( 2: 2 ) = 'N'
352  END IF
353 *
354 * ==== NWR = recommended deflation window size. At this
355 * . point, N .GT. NTINY = 15, so there is enough
356 * . subdiagonal workspace for NWR.GE.2 as required.
357 * . (In fact, there is enough subdiagonal space for
358 * . NWR.GE.4.) ====
359 *
360  nwr = ilaenv( 13, 'CLAQR0', jbcmpz, n, ilo, ihi, lwork )
361  nwr = max( 2, nwr )
362  nwr = min( ihi-ilo+1, ( n-1 ) / 3, nwr )
363 *
364 * ==== NSR = recommended number of simultaneous shifts.
365 * . At this point N .GT. NTINY = 15, so there is at
366 * . enough subdiagonal workspace for NSR to be even
367 * . and greater than or equal to two as required. ====
368 *
369  nsr = ilaenv( 15, 'CLAQR0', jbcmpz, n, ilo, ihi, lwork )
370  nsr = min( nsr, ( n-3 ) / 6, ihi-ilo )
371  nsr = max( 2, nsr-mod( nsr, 2 ) )
372 *
373 * ==== Estimate optimal workspace ====
374 *
375 * ==== Workspace query call to CLAQR3 ====
376 *
377  CALL claqr3( wantt, wantz, n, ilo, ihi, nwr+1, h, ldh, iloz,
378  $ ihiz, z, ldz, ls, ld, w, h, ldh, n, h, ldh, n, h,
379  $ ldh, work, -1 )
380 *
381 * ==== Optimal workspace = MAX(CLAQR5, CLAQR3) ====
382 *
383  lwkopt = max( 3*nsr / 2, int( work( 1 ) ) )
384 *
385 * ==== Quick return in case of workspace query. ====
386 *
387  IF( lwork.EQ.-1 ) THEN
388  work( 1 ) = cmplx( lwkopt, 0 )
389  RETURN
390  END IF
391 *
392 * ==== CLAHQR/CLAQR0 crossover point ====
393 *
394  nmin = ilaenv( 12, 'CLAQR0', jbcmpz, n, ilo, ihi, lwork )
395  nmin = max( ntiny, nmin )
396 *
397 * ==== Nibble crossover point ====
398 *
399  nibble = ilaenv( 14, 'CLAQR0', jbcmpz, n, ilo, ihi, lwork )
400  nibble = max( 0, nibble )
401 *
402 * ==== Accumulate reflections during ttswp? Use block
403 * . 2-by-2 structure during matrix-matrix multiply? ====
404 *
405  kacc22 = ilaenv( 16, 'CLAQR0', jbcmpz, n, ilo, ihi, lwork )
406  kacc22 = max( 0, kacc22 )
407  kacc22 = min( 2, kacc22 )
408 *
409 * ==== NWMAX = the largest possible deflation window for
410 * . which there is sufficient workspace. ====
411 *
412  nwmax = min( ( n-1 ) / 3, lwork / 2 )
413  nw = nwmax
414 *
415 * ==== NSMAX = the Largest number of simultaneous shifts
416 * . for which there is sufficient workspace. ====
417 *
418  nsmax = min( ( n-3 ) / 6, 2*lwork / 3 )
419  nsmax = nsmax - mod( nsmax, 2 )
420 *
421 * ==== NDFL: an iteration count restarted at deflation. ====
422 *
423  ndfl = 1
424 *
425 * ==== ITMAX = iteration limit ====
426 *
427  itmax = max( 30, 2*kexsh )*max( 10, ( ihi-ilo+1 ) )
428 *
429 * ==== Last row and column in the active block ====
430 *
431  kbot = ihi
432 *
433 * ==== Main Loop ====
434 *
435  DO 70 it = 1, itmax
436 *
437 * ==== Done when KBOT falls below ILO ====
438 *
439  IF( kbot.LT.ilo )
440  $ GO TO 80
441 *
442 * ==== Locate active block ====
443 *
444  DO 10 k = kbot, ilo + 1, -1
445  IF( h( k, k-1 ).EQ.zero )
446  $ GO TO 20
447  10 CONTINUE
448  k = ilo
449  20 CONTINUE
450  ktop = k
451 *
452 * ==== Select deflation window size:
453 * . Typical Case:
454 * . If possible and advisable, nibble the entire
455 * . active block. If not, use size MIN(NWR,NWMAX)
456 * . or MIN(NWR+1,NWMAX) depending upon which has
457 * . the smaller corresponding subdiagonal entry
458 * . (a heuristic).
459 * .
460 * . Exceptional Case:
461 * . If there have been no deflations in KEXNW or
462 * . more iterations, then vary the deflation window
463 * . size. At first, because, larger windows are,
464 * . in general, more powerful than smaller ones,
465 * . rapidly increase the window to the maximum possible.
466 * . Then, gradually reduce the window size. ====
467 *
468  nh = kbot - ktop + 1
469  nwupbd = min( nh, nwmax )
470  IF( ndfl.LT.kexnw ) THEN
471  nw = min( nwupbd, nwr )
472  ELSE
473  nw = min( nwupbd, 2*nw )
474  END IF
475  IF( nw.LT.nwmax ) THEN
476  IF( nw.GE.nh-1 ) THEN
477  nw = nh
478  ELSE
479  kwtop = kbot - nw + 1
480  IF( cabs1( h( kwtop, kwtop-1 ) ).GT.
481  $ cabs1( h( kwtop-1, kwtop-2 ) ) )nw = nw + 1
482  END IF
483  END IF
484  IF( ndfl.LT.kexnw ) THEN
485  ndec = -1
486  ELSE IF( ndec.GE.0 .OR. nw.GE.nwupbd ) THEN
487  ndec = ndec + 1
488  IF( nw-ndec.LT.2 )
489  $ ndec = 0
490  nw = nw - ndec
491  END IF
492 *
493 * ==== Aggressive early deflation:
494 * . split workspace under the subdiagonal into
495 * . - an nw-by-nw work array V in the lower
496 * . left-hand-corner,
497 * . - an NW-by-at-least-NW-but-more-is-better
498 * . (NW-by-NHO) horizontal work array along
499 * . the bottom edge,
500 * . - an at-least-NW-but-more-is-better (NHV-by-NW)
501 * . vertical work array along the left-hand-edge.
502 * . ====
503 *
504  kv = n - nw + 1
505  kt = nw + 1
506  nho = ( n-nw-1 ) - kt + 1
507  kwv = nw + 2
508  nve = ( n-nw ) - kwv + 1
509 *
510 * ==== Aggressive early deflation ====
511 *
512  CALL claqr3( wantt, wantz, n, ktop, kbot, nw, h, ldh, iloz,
513  $ ihiz, z, ldz, ls, ld, w, h( kv, 1 ), ldh, nho,
514  $ h( kv, kt ), ldh, nve, h( kwv, 1 ), ldh, work,
515  $ lwork )
516 *
517 * ==== Adjust KBOT accounting for new deflations. ====
518 *
519  kbot = kbot - ld
520 *
521 * ==== KS points to the shifts. ====
522 *
523  ks = kbot - ls + 1
524 *
525 * ==== Skip an expensive QR sweep if there is a (partly
526 * . heuristic) reason to expect that many eigenvalues
527 * . will deflate without it. Here, the QR sweep is
528 * . skipped if many eigenvalues have just been deflated
529 * . or if the remaining active block is small.
530 *
531  IF( ( ld.EQ.0 ) .OR. ( ( 100*ld.LE.nw*nibble ) .AND. ( kbot-
532  $ ktop+1.GT.min( nmin, nwmax ) ) ) ) THEN
533 *
534 * ==== NS = nominal number of simultaneous shifts.
535 * . This may be lowered (slightly) if CLAQR3
536 * . did not provide that many shifts. ====
537 *
538  ns = min( nsmax, nsr, max( 2, kbot-ktop ) )
539  ns = ns - mod( ns, 2 )
540 *
541 * ==== If there have been no deflations
542 * . in a multiple of KEXSH iterations,
543 * . then try exceptional shifts.
544 * . Otherwise use shifts provided by
545 * . CLAQR3 above or from the eigenvalues
546 * . of a trailing principal submatrix. ====
547 *
548  IF( mod( ndfl, kexsh ).EQ.0 ) THEN
549  ks = kbot - ns + 1
550  DO 30 i = kbot, ks + 1, -2
551  w( i ) = h( i, i ) + wilk1*cabs1( h( i, i-1 ) )
552  w( i-1 ) = w( i )
553  30 CONTINUE
554  ELSE
555 *
556 * ==== Got NS/2 or fewer shifts? Use CLAQR4 or
557 * . CLAHQR on a trailing principal submatrix to
558 * . get more. (Since NS.LE.NSMAX.LE.(N-3)/6,
559 * . there is enough space below the subdiagonal
560 * . to fit an NS-by-NS scratch array.) ====
561 *
562  IF( kbot-ks+1.LE.ns / 2 ) THEN
563  ks = kbot - ns + 1
564  kt = n - ns + 1
565  CALL clacpy( 'A', ns, ns, h( ks, ks ), ldh,
566  $ h( kt, 1 ), ldh )
567  IF( ns.GT.nmin ) THEN
568  CALL claqr4( .false., .false., ns, 1, ns,
569  $ h( kt, 1 ), ldh, w( ks ), 1, 1,
570  $ zdum, 1, work, lwork, inf )
571  ELSE
572  CALL clahqr( .false., .false., ns, 1, ns,
573  $ h( kt, 1 ), ldh, w( ks ), 1, 1,
574  $ zdum, 1, inf )
575  END IF
576  ks = ks + inf
577 *
578 * ==== In case of a rare QR failure use
579 * . eigenvalues of the trailing 2-by-2
580 * . principal submatrix. Scale to avoid
581 * . overflows, underflows and subnormals.
582 * . (The scale factor S can not be zero,
583 * . because H(KBOT,KBOT-1) is nonzero.) ====
584 *
585  IF( ks.GE.kbot ) THEN
586  s = cabs1( h( kbot-1, kbot-1 ) ) +
587  $ cabs1( h( kbot, kbot-1 ) ) +
588  $ cabs1( h( kbot-1, kbot ) ) +
589  $ cabs1( h( kbot, kbot ) )
590  aa = h( kbot-1, kbot-1 ) / s
591  cc = h( kbot, kbot-1 ) / s
592  bb = h( kbot-1, kbot ) / s
593  dd = h( kbot, kbot ) / s
594  tr2 = ( aa+dd ) / two
595  det = ( aa-tr2 )*( dd-tr2 ) - bb*cc
596  rtdisc = sqrt( -det )
597  w( kbot-1 ) = ( tr2+rtdisc )*s
598  w( kbot ) = ( tr2-rtdisc )*s
599 *
600  ks = kbot - 1
601  END IF
602  END IF
603 *
604  IF( kbot-ks+1.GT.ns ) THEN
605 *
606 * ==== Sort the shifts (Helps a little) ====
607 *
608  sorted = .false.
609  DO 50 k = kbot, ks + 1, -1
610  IF( sorted )
611  $ GO TO 60
612  sorted = .true.
613  DO 40 i = ks, k - 1
614  IF( cabs1( w( i ) ).LT.cabs1( w( i+1 ) ) )
615  $ THEN
616  sorted = .false.
617  swap = w( i )
618  w( i ) = w( i+1 )
619  w( i+1 ) = swap
620  END IF
621  40 CONTINUE
622  50 CONTINUE
623  60 CONTINUE
624  END IF
625  END IF
626 *
627 * ==== If there are only two shifts, then use
628 * . only one. ====
629 *
630  IF( kbot-ks+1.EQ.2 ) THEN
631  IF( cabs1( w( kbot )-h( kbot, kbot ) ).LT.
632  $ cabs1( w( kbot-1 )-h( kbot, kbot ) ) ) THEN
633  w( kbot-1 ) = w( kbot )
634  ELSE
635  w( kbot ) = w( kbot-1 )
636  END IF
637  END IF
638 *
639 * ==== Use up to NS of the the smallest magnitude
640 * . shifts. If there aren't NS shifts available,
641 * . then use them all, possibly dropping one to
642 * . make the number of shifts even. ====
643 *
644  ns = min( ns, kbot-ks+1 )
645  ns = ns - mod( ns, 2 )
646  ks = kbot - ns + 1
647 *
648 * ==== Small-bulge multi-shift QR sweep:
649 * . split workspace under the subdiagonal into
650 * . - a KDU-by-KDU work array U in the lower
651 * . left-hand-corner,
652 * . - a KDU-by-at-least-KDU-but-more-is-better
653 * . (KDU-by-NHo) horizontal work array WH along
654 * . the bottom edge,
655 * . - and an at-least-KDU-but-more-is-better-by-KDU
656 * . (NVE-by-KDU) vertical work WV arrow along
657 * . the left-hand-edge. ====
658 *
659  kdu = 2*ns
660  ku = n - kdu + 1
661  kwh = kdu + 1
662  nho = ( n-kdu+1-4 ) - ( kdu+1 ) + 1
663  kwv = kdu + 4
664  nve = n - kdu - kwv + 1
665 *
666 * ==== Small-bulge multi-shift QR sweep ====
667 *
668  CALL claqr5( wantt, wantz, kacc22, n, ktop, kbot, ns,
669  $ w( ks ), h, ldh, iloz, ihiz, z, ldz, work,
670  $ 3, h( ku, 1 ), ldh, nve, h( kwv, 1 ), ldh,
671  $ nho, h( ku, kwh ), ldh )
672  END IF
673 *
674 * ==== Note progress (or the lack of it). ====
675 *
676  IF( ld.GT.0 ) THEN
677  ndfl = 1
678  ELSE
679  ndfl = ndfl + 1
680  END IF
681 *
682 * ==== End of main loop ====
683  70 CONTINUE
684 *
685 * ==== Iteration limit exceeded. Set INFO to show where
686 * . the problem occurred and exit. ====
687 *
688  info = kbot
689  80 CONTINUE
690  END IF
691 *
692 * ==== Return the optimal value of LWORK. ====
693 *
694  work( 1 ) = cmplx( lwkopt, 0 )
695 *
696 * ==== End of CLAQR0 ====
697 *
698  END
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:257
subroutine claqr0(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)
CLAQR0 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...
Definition: claqr0.f:240
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:248
subroutine claqr3(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)
CLAQR3 performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fu...
Definition: claqr3.f:266
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,...
Definition: clahqr.f:195
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:103