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