LAPACK  3.10.0
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 >= 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 > 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 > 0, then 1 <= ILO <= IHI <= 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 > 0 is given under the
107 *> description of INFO below.)
108 *>
109 *> This subroutine may explicitly set H(i,j) = 0 for i > 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 >= 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 <= ILOZ <= ILO; IHI <= IHIZ <= 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 > 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 >= MAX(1,IHIZ). Otherwise, LDZ >= 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 >= 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 *> INFO is INTEGER
185 *> = 0: successful exit
186 *> > 0: if INFO = i, CLAQR4 failed to compute all of
187 *> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
188 *> and WI contain those eigenvalues which have been
189 *> successfully computed. (Failures are rare.)
190 *>
191 *> If INFO > 0 and WANT is .FALSE., then on exit,
192 *> the remaining unconverged eigenvalues are the eigen-
193 *> values of the upper Hessenberg matrix rows and
194 *> columns ILO through INFO of the final, output
195 *> value of H.
196 *>
197 *> If INFO > 0 and WANTT is .TRUE., then on exit
198 *>
199 *> (*) (initial value of H)*U = U*(final value of H)
200 *>
201 *> where U is a unitary matrix. The final
202 *> value of H is upper Hessenberg and triangular in
203 *> rows and columns INFO+1 through IHI.
204 *>
205 *> If INFO > 0 and WANTZ is .TRUE., then on exit
206 *>
207 *> (final value of Z(ILO:IHI,ILOZ:IHIZ)
208 *> = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U
209 *>
210 *> where U is the unitary matrix in (*) (regard-
211 *> less of the value of WANTT.)
212 *>
213 *> If INFO > 0 and WANTZ is .FALSE., then Z is not
214 *> accessed.
215 *> \endverbatim
216 *
217 * Authors:
218 * ========
219 *
220 *> \author Univ. of Tennessee
221 *> \author Univ. of California Berkeley
222 *> \author Univ. of Colorado Denver
223 *> \author NAG Ltd.
224 *
225 *> \ingroup complexOTHERauxiliary
226 *
227 *> \par Contributors:
228 * ==================
229 *>
230 *> Karen Braman and Ralph Byers, Department of Mathematics,
231 *> University of Kansas, USA
232 *
233 *> \par References:
234 * ================
235 *>
236 *> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
237 *> Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
238 *> Performance, SIAM Journal of Matrix Analysis, volume 23, pages
239 *> 929--947, 2002.
240 *> \n
241 *> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
242 *> Algorithm Part II: Aggressive Early Deflation, SIAM Journal
243 *> of Matrix Analysis, volume 23, pages 948--973, 2002.
244 *>
245 * =====================================================================
246  SUBROUTINE claqr4( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
247  $ IHIZ, Z, LDZ, WORK, LWORK, INFO )
248 *
249 * -- LAPACK auxiliary routine --
250 * -- LAPACK is a software package provided by Univ. of Tennessee, --
251 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
252 *
253 * .. Scalar Arguments ..
254  INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N
255  LOGICAL WANTT, WANTZ
256 * ..
257 * .. Array Arguments ..
258  COMPLEX H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * )
259 * ..
260 *
261 *
262 * ================================================================
263 *
264 * .. Parameters ..
265 *
266 * ==== Matrices of order NTINY or smaller must be processed by
267 * . CLAHQR because of insufficient subdiagonal scratch space.
268 * . (This is a hard limit.) ====
269  INTEGER NTINY
270  parameter( ntiny = 15 )
271 *
272 * ==== Exceptional deflation windows: try to cure rare
273 * . slow convergence by varying the size of the
274 * . deflation window after KEXNW iterations. ====
275  INTEGER KEXNW
276  parameter( kexnw = 5 )
277 *
278 * ==== Exceptional shifts: try to cure rare slow convergence
279 * . with ad-hoc exceptional shifts every KEXSH iterations.
280 * . ====
281  INTEGER KEXSH
282  parameter( kexsh = 6 )
283 *
284 * ==== The constant WILK1 is used to form the exceptional
285 * . shifts. ====
286  REAL WILK1
287  parameter( wilk1 = 0.75e0 )
288  COMPLEX ZERO, ONE
289  parameter( zero = ( 0.0e0, 0.0e0 ),
290  $ one = ( 1.0e0, 0.0e0 ) )
291  REAL TWO
292  parameter( two = 2.0e0 )
293 * ..
294 * .. Local Scalars ..
295  COMPLEX AA, BB, CC, CDUM, DD, DET, RTDISC, SWAP, TR2
296  REAL S
297  INTEGER I, INF, IT, ITMAX, K, KACC22, KBOT, KDU, KS,
298  $ kt, ktop, ku, kv, kwh, kwtop, kwv, ld, ls,
299  $ lwkopt, ndec, ndfl, nh, nho, nibble, nmin, ns,
300  $ nsmax, nsr, nve, nw, nwmax, nwr, nwupbd
301  LOGICAL SORTED
302  CHARACTER JBCMPZ*2
303 * ..
304 * .. External Functions ..
305  INTEGER ILAENV
306  EXTERNAL ilaenv
307 * ..
308 * .. Local Arrays ..
309  COMPLEX ZDUM( 1, 1 )
310 * ..
311 * .. External Subroutines ..
312  EXTERNAL clacpy, clahqr, claqr2, claqr5
313 * ..
314 * .. Intrinsic Functions ..
315  INTRINSIC abs, aimag, cmplx, int, max, min, mod, real,
316  $ sqrt
317 * ..
318 * .. Statement Functions ..
319  REAL CABS1
320 * ..
321 * .. Statement Function definitions ..
322  cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
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 CLAHQR. ====
337 *
338  lwkopt = 1
339  IF( lwork.NE.-1 )
340  $ CALL clahqr( wantt, wantz, n, ilo, ihi, h, ldh, w, iloz,
341  $ 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 = 15, 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.4.) ====
369 *
370  nwr = ilaenv( 13, 'CLAQR4', 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 = 15, 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, 'CLAQR4', jbcmpz, n, ilo, ihi, lwork )
380  nsr = min( nsr, ( n-3 ) / 6, ihi-ilo )
381  nsr = max( 2, nsr-mod( nsr, 2 ) )
382 *
383 * ==== Estimate optimal workspace ====
384 *
385 * ==== Workspace query call to CLAQR2 ====
386 *
387  CALL claqr2( wantt, wantz, n, ilo, ihi, nwr+1, h, ldh, iloz,
388  $ ihiz, z, ldz, ls, ld, w, h, ldh, n, h, ldh, n, h,
389  $ ldh, work, -1 )
390 *
391 * ==== Optimal workspace = MAX(CLAQR5, CLAQR2) ====
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 ) = cmplx( lwkopt, 0 )
399  RETURN
400  END IF
401 *
402 * ==== CLAHQR/CLAQR0 crossover point ====
403 *
404  nmin = ilaenv( 12, 'CLAQR4', jbcmpz, n, ilo, ihi, lwork )
405  nmin = max( ntiny, nmin )
406 *
407 * ==== Nibble crossover point ====
408 *
409  nibble = ilaenv( 14, 'CLAQR4', 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, 'CLAQR4', 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-3 ) / 6, 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 70 it = 1, itmax
446 *
447 * ==== Done when KBOT falls below ILO ====
448 *
449  IF( kbot.LT.ilo )
450  $ GO TO 80
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( cabs1( h( kwtop, kwtop-1 ) ).GT.
491  $ cabs1( 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 claqr2( wantt, wantz, n, ktop, kbot, nw, h, ldh, iloz,
523  $ ihiz, z, ldz, ls, ld, w, h( kv, 1 ), ldh, nho,
524  $ h( kv, kt ), ldh, nve, h( kwv, 1 ), ldh, work,
525  $ 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 CLAQR2
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 * . CLAQR2 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, ks + 1, -2
561  w( i ) = h( i, i ) + wilk1*cabs1( h( i, i-1 ) )
562  w( i-1 ) = w( i )
563  30 CONTINUE
564  ELSE
565 *
566 * ==== Got NS/2 or fewer shifts? Use CLAHQR
567 * . on a trailing principal submatrix to
568 * . get more. (Since NS.LE.NSMAX.LE.(N-3)/6,
569 * . there is enough space below the subdiagonal
570 * . to fit an NS-by-NS scratch array.) ====
571 *
572  IF( kbot-ks+1.LE.ns / 2 ) THEN
573  ks = kbot - ns + 1
574  kt = n - ns + 1
575  CALL clacpy( 'A', ns, ns, h( ks, ks ), ldh,
576  $ h( kt, 1 ), ldh )
577  CALL clahqr( .false., .false., ns, 1, ns,
578  $ h( kt, 1 ), ldh, w( ks ), 1, 1, zdum,
579  $ 1, inf )
580  ks = ks + inf
581 *
582 * ==== In case of a rare QR failure use
583 * . eigenvalues of the trailing 2-by-2
584 * . principal submatrix. Scale to avoid
585 * . overflows, underflows and subnormals.
586 * . (The scale factor S can not be zero,
587 * . because H(KBOT,KBOT-1) is nonzero.) ====
588 *
589  IF( ks.GE.kbot ) THEN
590  s = cabs1( h( kbot-1, kbot-1 ) ) +
591  $ cabs1( h( kbot, kbot-1 ) ) +
592  $ cabs1( h( kbot-1, kbot ) ) +
593  $ cabs1( h( kbot, kbot ) )
594  aa = h( kbot-1, kbot-1 ) / s
595  cc = h( kbot, kbot-1 ) / s
596  bb = h( kbot-1, kbot ) / s
597  dd = h( kbot, kbot ) / s
598  tr2 = ( aa+dd ) / two
599  det = ( aa-tr2 )*( dd-tr2 ) - bb*cc
600  rtdisc = sqrt( -det )
601  w( kbot-1 ) = ( tr2+rtdisc )*s
602  w( kbot ) = ( tr2-rtdisc )*s
603 *
604  ks = kbot - 1
605  END IF
606  END IF
607 *
608  IF( kbot-ks+1.GT.ns ) THEN
609 *
610 * ==== Sort the shifts (Helps a little) ====
611 *
612  sorted = .false.
613  DO 50 k = kbot, ks + 1, -1
614  IF( sorted )
615  $ GO TO 60
616  sorted = .true.
617  DO 40 i = ks, k - 1
618  IF( cabs1( w( i ) ).LT.cabs1( w( i+1 ) ) )
619  $ THEN
620  sorted = .false.
621  swap = w( i )
622  w( i ) = w( i+1 )
623  w( i+1 ) = swap
624  END IF
625  40 CONTINUE
626  50 CONTINUE
627  60 CONTINUE
628  END IF
629  END IF
630 *
631 * ==== If there are only two shifts, then use
632 * . only one. ====
633 *
634  IF( kbot-ks+1.EQ.2 ) THEN
635  IF( cabs1( w( kbot )-h( kbot, kbot ) ).LT.
636  $ cabs1( w( kbot-1 )-h( kbot, kbot ) ) ) THEN
637  w( kbot-1 ) = w( kbot )
638  ELSE
639  w( kbot ) = w( kbot-1 )
640  END IF
641  END IF
642 *
643 * ==== Use up to NS of the the smallest magnitude
644 * . shifts. If there aren't NS shifts available,
645 * . then use them all, possibly dropping one to
646 * . make the number of shifts even. ====
647 *
648  ns = min( ns, kbot-ks+1 )
649  ns = ns - mod( ns, 2 )
650  ks = kbot - ns + 1
651 *
652 * ==== Small-bulge multi-shift QR sweep:
653 * . split workspace under the subdiagonal into
654 * . - a KDU-by-KDU work array U in the lower
655 * . left-hand-corner,
656 * . - a KDU-by-at-least-KDU-but-more-is-better
657 * . (KDU-by-NHo) horizontal work array WH along
658 * . the bottom edge,
659 * . - and an at-least-KDU-but-more-is-better-by-KDU
660 * . (NVE-by-KDU) vertical work WV arrow along
661 * . the left-hand-edge. ====
662 *
663  kdu = 2*ns
664  ku = n - kdu + 1
665  kwh = kdu + 1
666  nho = ( n-kdu+1-4 ) - ( kdu+1 ) + 1
667  kwv = kdu + 4
668  nve = n - kdu - kwv + 1
669 *
670 * ==== Small-bulge multi-shift QR sweep ====
671 *
672  CALL claqr5( wantt, wantz, kacc22, n, ktop, kbot, ns,
673  $ w( ks ), h, ldh, iloz, ihiz, z, ldz, work,
674  $ 3, h( ku, 1 ), ldh, nve, h( kwv, 1 ), ldh,
675  $ nho, h( ku, kwh ), ldh )
676  END IF
677 *
678 * ==== Note progress (or the lack of it). ====
679 *
680  IF( ld.GT.0 ) THEN
681  ndfl = 1
682  ELSE
683  ndfl = ndfl + 1
684  END IF
685 *
686 * ==== End of main loop ====
687  70 CONTINUE
688 *
689 * ==== Iteration limit exceeded. Set INFO to show where
690 * . the problem occurred and exit. ====
691 *
692  info = kbot
693  80 CONTINUE
694  END IF
695 *
696 * ==== Return the optimal value of LWORK. ====
697 *
698  work( 1 ) = cmplx( lwkopt, 0 )
699 *
700 * ==== End of CLAQR4 ====
701 *
702  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 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:269
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 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