LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
dlaqr4.f
Go to the documentation of this file.
1 *> \brief \b DLAQR4 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 DLAQR4 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqr4.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqr4.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqr4.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLAQR4( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
22 * ILOZ, 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 * DOUBLE PRECISION H( LDH, * ), WI( * ), WORK( * ), WR( * ),
30 * $ Z( LDZ, * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> DLAQR4 implements one level of recursion for DLAQR0.
40 *> It is a complete implementation of the small bulge multi-shift
41 *> QR algorithm. It may be called by DLAQR0 and, for large enough
42 *> deflation window size, it may be called by DLAQR3. This
43 *> subroutine is identical to DLAQR0 except that it calls DLAQR2
44 *> instead of DLAQR3.
45 *>
46 *> DLAQR4 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**T, where T is an upper quasi-triangular matrix (the
49 *> Schur form), and Z is the orthogonal matrix of Schur vectors.
50 *>
51 *> Optionally Z may be postmultiplied into an input orthogonal
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 orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T.
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 DGEBAL, and then passed to DGEHRD when the
92 *> matrix output by DGEBAL 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 DOUBLE PRECISION array, dimension (LDH,N)
101 *> On entry, the upper Hessenberg matrix H.
102 *> On exit, if INFO = 0 and WANTT is .TRUE., then H contains
103 *> the upper quasi-triangular matrix T from the Schur
104 *> decomposition (the Schur form); 2-by-2 diagonal blocks
105 *> (corresponding to complex conjugate pairs of eigenvalues)
106 *> are returned in standard form, with H(i,i) = H(i+1,i+1)
107 *> and H(i+1,i)*H(i,i+1) < 0. If INFO = 0 and WANTT is
108 *> .FALSE., then the contents of H are unspecified on exit.
109 *> (The output value of H when INFO > 0 is given under the
110 *> description of INFO below.)
111 *>
112 *> This subroutine may explicitly set H(i,j) = 0 for i > j and
113 *> j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N.
114 *> \endverbatim
115 *>
116 *> \param[in] LDH
117 *> \verbatim
118 *> LDH is INTEGER
119 *> The leading dimension of the array H. LDH >= max(1,N).
120 *> \endverbatim
121 *>
122 *> \param[out] WR
123 *> \verbatim
124 *> WR is DOUBLE PRECISION array, dimension (IHI)
125 *> \endverbatim
126 *>
127 *> \param[out] WI
128 *> \verbatim
129 *> WI is DOUBLE PRECISION array, dimension (IHI)
130 *> The real and imaginary parts, respectively, of the computed
131 *> eigenvalues of H(ILO:IHI,ILO:IHI) are stored in WR(ILO:IHI)
132 *> and WI(ILO:IHI). If two eigenvalues are computed as a
133 *> complex conjugate pair, they are stored in consecutive
134 *> elements of WR and WI, say the i-th and (i+1)th, with
135 *> WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., then
136 *> the eigenvalues are stored in the same order as on the
137 *> diagonal of the Schur form returned in H, with
138 *> WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal
139 *> block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and
140 *> WI(i+1) = -WI(i).
141 *> \endverbatim
142 *>
143 *> \param[in] ILOZ
144 *> \verbatim
145 *> ILOZ is INTEGER
146 *> \endverbatim
147 *>
148 *> \param[in] IHIZ
149 *> \verbatim
150 *> IHIZ is INTEGER
151 *> Specify the rows of Z to which transformations must be
152 *> applied if WANTZ is .TRUE..
153 *> 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
154 *> \endverbatim
155 *>
156 *> \param[in,out] Z
157 *> \verbatim
158 *> Z is DOUBLE PRECISION array, dimension (LDZ,IHI)
159 *> If WANTZ is .FALSE., then Z is not referenced.
160 *> If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is
161 *> replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the
162 *> orthogonal Schur factor of H(ILO:IHI,ILO:IHI).
163 *> (The output value of Z when INFO > 0 is given under
164 *> the description of INFO below.)
165 *> \endverbatim
166 *>
167 *> \param[in] LDZ
168 *> \verbatim
169 *> LDZ is INTEGER
170 *> The leading dimension of the array Z. if WANTZ is .TRUE.
171 *> then LDZ >= MAX(1,IHIZ). Otherwise, LDZ >= 1.
172 *> \endverbatim
173 *>
174 *> \param[out] WORK
175 *> \verbatim
176 *> WORK is DOUBLE PRECISION array, dimension LWORK
177 *> On exit, if LWORK = -1, WORK(1) returns an estimate of
178 *> the optimal value for LWORK.
179 *> \endverbatim
180 *>
181 *> \param[in] LWORK
182 *> \verbatim
183 *> LWORK is INTEGER
184 *> The dimension of the array WORK. LWORK >= max(1,N)
185 *> is sufficient, but LWORK typically as large as 6*N may
186 *> be required for optimal performance. A workspace query
187 *> to determine the optimal workspace size is recommended.
188 *>
189 *> If LWORK = -1, then DLAQR4 does a workspace query.
190 *> In this case, DLAQR4 checks the input parameters and
191 *> estimates the optimal workspace size for the given
192 *> values of N, ILO and IHI. The estimate is returned
193 *> in WORK(1). No error message related to LWORK is
194 *> issued by XERBLA. Neither H nor Z are accessed.
195 *> \endverbatim
196 *>
197 *> \param[out] INFO
198 *> \verbatim
199 *> INFO is INTEGER
200 *> = 0: successful exit
201 *> > 0: if INFO = i, DLAQR4 failed to compute all of
202 *> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
203 *> and WI contain those eigenvalues which have been
204 *> successfully computed. (Failures are rare.)
205 *>
206 *> If INFO > 0 and WANT is .FALSE., then on exit,
207 *> the remaining unconverged eigenvalues are the eigen-
208 *> values of the upper Hessenberg matrix rows and
209 *> columns ILO through INFO of the final, output
210 *> value of H.
211 *>
212 *> If INFO > 0 and WANTT is .TRUE., then on exit
213 *>
214 *> (*) (initial value of H)*U = U*(final value of H)
215 *>
216 *> where U is a orthogonal matrix. The final
217 *> value of H is upper Hessenberg and triangular in
218 *> rows and columns INFO+1 through IHI.
219 *>
220 *> If INFO > 0 and WANTZ is .TRUE., then on exit
221 *>
222 *> (final value of Z(ILO:IHI,ILOZ:IHIZ)
223 *> = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U
224 *>
225 *> where U is the orthogonal matrix in (*) (regard-
226 *> less of the value of WANTT.)
227 *>
228 *> If INFO > 0 and WANTZ is .FALSE., then Z is not
229 *> accessed.
230 *> \endverbatim
231 *
232 * Authors:
233 * ========
234 *
235 *> \author Univ. of Tennessee
236 *> \author Univ. of California Berkeley
237 *> \author Univ. of Colorado Denver
238 *> \author NAG Ltd.
239 *
240 *> \ingroup doubleOTHERauxiliary
241 *
242 *> \par Contributors:
243 * ==================
244 *>
245 *> Karen Braman and Ralph Byers, Department of Mathematics,
246 *> University of Kansas, USA
247 *
248 *> \par References:
249 * ================
250 *>
251 *> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
252 *> Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
253 *> Performance, SIAM Journal of Matrix Analysis, volume 23, pages
254 *> 929--947, 2002.
255 *> \n
256 *> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
257 *> Algorithm Part II: Aggressive Early Deflation, SIAM Journal
258 *> of Matrix Analysis, volume 23, pages 948--973, 2002.
259 *>
260 * =====================================================================
261  SUBROUTINE dlaqr4( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
262  $ ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO )
263 *
264 * -- LAPACK auxiliary routine --
265 * -- LAPACK is a software package provided by Univ. of Tennessee, --
266 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
267 *
268 * .. Scalar Arguments ..
269  INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N
270  LOGICAL WANTT, WANTZ
271 * ..
272 * .. Array Arguments ..
273  DOUBLE PRECISION H( LDH, * ), WI( * ), WORK( * ), WR( * ),
274  $ z( ldz, * )
275 * ..
276 *
277 * ================================================================
278 * .. Parameters ..
279 *
280 * ==== Matrices of order NTINY or smaller must be processed by
281 * . DLAHQR because of insufficient subdiagonal scratch space.
282 * . (This is a hard limit.) ====
283  INTEGER NTINY
284  parameter( ntiny = 15 )
285 *
286 * ==== Exceptional deflation windows: try to cure rare
287 * . slow convergence by varying the size of the
288 * . deflation window after KEXNW iterations. ====
289  INTEGER KEXNW
290  parameter( kexnw = 5 )
291 *
292 * ==== Exceptional shifts: try to cure rare slow convergence
293 * . with ad-hoc exceptional shifts every KEXSH iterations.
294 * . ====
295  INTEGER KEXSH
296  parameter( kexsh = 6 )
297 *
298 * ==== The constants WILK1 and WILK2 are used to form the
299 * . exceptional shifts. ====
300  DOUBLE PRECISION WILK1, WILK2
301  parameter( wilk1 = 0.75d0, wilk2 = -0.4375d0 )
302  DOUBLE PRECISION ZERO, ONE
303  parameter( zero = 0.0d0, one = 1.0d0 )
304 * ..
305 * .. Local Scalars ..
306  DOUBLE PRECISION AA, BB, CC, CS, DD, SN, SS, SWAP
307  INTEGER I, INF, IT, ITMAX, K, KACC22, KBOT, KDU, KS,
308  $ kt, ktop, ku, kv, kwh, kwtop, kwv, ld, ls,
309  $ lwkopt, ndec, ndfl, nh, nho, nibble, nmin, ns,
310  $ nsmax, nsr, nve, nw, nwmax, nwr, nwupbd
311  LOGICAL SORTED
312  CHARACTER JBCMPZ*2
313 * ..
314 * .. External Functions ..
315  INTEGER ILAENV
316  EXTERNAL ilaenv
317 * ..
318 * .. Local Arrays ..
319  DOUBLE PRECISION ZDUM( 1, 1 )
320 * ..
321 * .. External Subroutines ..
322  EXTERNAL dlacpy, dlahqr, dlanv2, dlaqr2, dlaqr5
323 * ..
324 * .. Intrinsic Functions ..
325  INTRINSIC abs, dble, int, max, min, mod
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 DLAHQR. ====
340 *
341  lwkopt = 1
342  IF( lwork.NE.-1 )
343  $ CALL dlahqr( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi,
344  $ iloz, 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 = 15, 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.4.) ====
372 *
373  nwr = ilaenv( 13, 'DLAQR4', 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 = 15, 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, 'DLAQR4', jbcmpz, n, ilo, ihi, lwork )
383  nsr = min( nsr, ( n-3 ) / 6, ihi-ilo )
384  nsr = max( 2, nsr-mod( nsr, 2 ) )
385 *
386 * ==== Estimate optimal workspace ====
387 *
388 * ==== Workspace query call to DLAQR2 ====
389 *
390  CALL dlaqr2( wantt, wantz, n, ilo, ihi, nwr+1, h, ldh, iloz,
391  $ ihiz, z, ldz, ls, ld, wr, wi, h, ldh, n, h, ldh,
392  $ n, h, ldh, work, -1 )
393 *
394 * ==== Optimal workspace = MAX(DLAQR5, DLAQR2) ====
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 ) = dble( lwkopt )
402  RETURN
403  END IF
404 *
405 * ==== DLAHQR/DLAQR0 crossover point ====
406 *
407  nmin = ilaenv( 12, 'DLAQR4', jbcmpz, n, ilo, ihi, lwork )
408  nmin = max( ntiny, nmin )
409 *
410 * ==== Nibble crossover point ====
411 *
412  nibble = ilaenv( 14, 'DLAQR4', 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, 'DLAQR4', 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-3 ) / 6, 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 80 it = 1, itmax
449 *
450 * ==== Done when KBOT falls below ILO ====
451 *
452  IF( kbot.LT.ilo )
453  $ GO TO 90
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( abs( h( kwtop, kwtop-1 ) ).GT.
494  $ abs( 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 dlaqr2( wantt, wantz, n, ktop, kbot, nw, h, ldh, iloz,
526  $ ihiz, z, ldz, ls, ld, wr, wi, h( kv, 1 ), ldh,
527  $ nho, h( kv, kt ), ldh, nve, h( kwv, 1 ), ldh,
528  $ work, 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 DLAQR2
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 * . DLAQR2 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, max( ks+1, ktop+2 ), -2
564  ss = abs( h( i, i-1 ) ) + abs( h( i-1, i-2 ) )
565  aa = wilk1*ss + h( i, i )
566  bb = ss
567  cc = wilk2*ss
568  dd = aa
569  CALL dlanv2( aa, bb, cc, dd, wr( i-1 ), wi( i-1 ),
570  $ wr( i ), wi( i ), cs, sn )
571  30 CONTINUE
572  IF( ks.EQ.ktop ) THEN
573  wr( ks+1 ) = h( ks+1, ks+1 )
574  wi( ks+1 ) = zero
575  wr( ks ) = wr( ks+1 )
576  wi( ks ) = wi( ks+1 )
577  END IF
578  ELSE
579 *
580 * ==== Got NS/2 or fewer shifts? Use DLAHQR
581 * . on a trailing principal submatrix to
582 * . get more. (Since NS.LE.NSMAX.LE.(N-3)/6,
583 * . there is enough space below the subdiagonal
584 * . to fit an NS-by-NS scratch array.) ====
585 *
586  IF( kbot-ks+1.LE.ns / 2 ) THEN
587  ks = kbot - ns + 1
588  kt = n - ns + 1
589  CALL dlacpy( 'A', ns, ns, h( ks, ks ), ldh,
590  $ h( kt, 1 ), ldh )
591  CALL dlahqr( .false., .false., ns, 1, ns,
592  $ h( kt, 1 ), ldh, wr( ks ), wi( ks ),
593  $ 1, 1, zdum, 1, inf )
594  ks = ks + inf
595 *
596 * ==== In case of a rare QR failure use
597 * . eigenvalues of the trailing 2-by-2
598 * . principal submatrix. ====
599 *
600  IF( ks.GE.kbot ) THEN
601  aa = h( kbot-1, kbot-1 )
602  cc = h( kbot, kbot-1 )
603  bb = h( kbot-1, kbot )
604  dd = h( kbot, kbot )
605  CALL dlanv2( aa, bb, cc, dd, wr( kbot-1 ),
606  $ wi( kbot-1 ), wr( kbot ),
607  $ wi( kbot ), cs, sn )
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 * . Bubble sort keeps complex conjugate
616 * . pairs together. ====
617 *
618  sorted = .false.
619  DO 50 k = kbot, ks + 1, -1
620  IF( sorted )
621  $ GO TO 60
622  sorted = .true.
623  DO 40 i = ks, k - 1
624  IF( abs( wr( i ) )+abs( wi( i ) ).LT.
625  $ abs( wr( i+1 ) )+abs( wi( i+1 ) ) ) THEN
626  sorted = .false.
627 *
628  swap = wr( i )
629  wr( i ) = wr( i+1 )
630  wr( i+1 ) = swap
631 *
632  swap = wi( i )
633  wi( i ) = wi( i+1 )
634  wi( i+1 ) = swap
635  END IF
636  40 CONTINUE
637  50 CONTINUE
638  60 CONTINUE
639  END IF
640 *
641 * ==== Shuffle shifts into pairs of real shifts
642 * . and pairs of complex conjugate shifts
643 * . assuming complex conjugate shifts are
644 * . already adjacent to one another. (Yes,
645 * . they are.) ====
646 *
647  DO 70 i = kbot, ks + 2, -2
648  IF( wi( i ).NE.-wi( i-1 ) ) THEN
649 *
650  swap = wr( i )
651  wr( i ) = wr( i-1 )
652  wr( i-1 ) = wr( i-2 )
653  wr( i-2 ) = swap
654 *
655  swap = wi( i )
656  wi( i ) = wi( i-1 )
657  wi( i-1 ) = wi( i-2 )
658  wi( i-2 ) = swap
659  END IF
660  70 CONTINUE
661  END IF
662 *
663 * ==== If there are only two shifts and both are
664 * . real, then use only one. ====
665 *
666  IF( kbot-ks+1.EQ.2 ) THEN
667  IF( wi( kbot ).EQ.zero ) THEN
668  IF( abs( wr( kbot )-h( kbot, kbot ) ).LT.
669  $ abs( wr( kbot-1 )-h( kbot, kbot ) ) ) THEN
670  wr( kbot-1 ) = wr( kbot )
671  ELSE
672  wr( kbot ) = wr( kbot-1 )
673  END IF
674  END IF
675  END IF
676 *
677 * ==== Use up to NS of the the smallest magnitude
678 * . shifts. If there aren't NS shifts available,
679 * . then use them all, possibly dropping one to
680 * . make the number of shifts even. ====
681 *
682  ns = min( ns, kbot-ks+1 )
683  ns = ns - mod( ns, 2 )
684  ks = kbot - ns + 1
685 *
686 * ==== Small-bulge multi-shift QR sweep:
687 * . split workspace under the subdiagonal into
688 * . - a KDU-by-KDU work array U in the lower
689 * . left-hand-corner,
690 * . - a KDU-by-at-least-KDU-but-more-is-better
691 * . (KDU-by-NHo) horizontal work array WH along
692 * . the bottom edge,
693 * . - and an at-least-KDU-but-more-is-better-by-KDU
694 * . (NVE-by-KDU) vertical work WV arrow along
695 * . the left-hand-edge. ====
696 *
697  kdu = 2*ns
698  ku = n - kdu + 1
699  kwh = kdu + 1
700  nho = ( n-kdu+1-4 ) - ( kdu+1 ) + 1
701  kwv = kdu + 4
702  nve = n - kdu - kwv + 1
703 *
704 * ==== Small-bulge multi-shift QR sweep ====
705 *
706  CALL dlaqr5( wantt, wantz, kacc22, n, ktop, kbot, ns,
707  $ wr( ks ), wi( ks ), h, ldh, iloz, ihiz, z,
708  $ ldz, work, 3, h( ku, 1 ), ldh, nve,
709  $ h( kwv, 1 ), ldh, nho, h( ku, kwh ), ldh )
710  END IF
711 *
712 * ==== Note progress (or the lack of it). ====
713 *
714  IF( ld.GT.0 ) THEN
715  ndfl = 1
716  ELSE
717  ndfl = ndfl + 1
718  END IF
719 *
720 * ==== End of main loop ====
721  80 CONTINUE
722 *
723 * ==== Iteration limit exceeded. Set INFO to show where
724 * . the problem occurred and exit. ====
725 *
726  info = kbot
727  90 CONTINUE
728  END IF
729 *
730 * ==== Return the optimal value of LWORK. ====
731 *
732  work( 1 ) = dble( lwkopt )
733 *
734 * ==== End of DLAQR4 ====
735 *
736  END
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine dlanv2(A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN)
DLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form.
Definition: dlanv2.f:127
subroutine dlaqr4(WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)
DLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...
Definition: dlaqr4.f:263
subroutine dlaqr5(WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, SR, SI, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV, WV, LDWV, NH, WH, LDWH)
DLAQR5 performs a single small-bulge multi-shift QR sweep.
Definition: dlaqr5.f:265
subroutine dlaqr2(WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK)
DLAQR2 performs the orthogonal similarity transformation of a Hessenberg matrix to detect and deflate...
Definition: dlaqr2.f:278
subroutine dlahqr(WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, INFO)
DLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix,...
Definition: dlahqr.f:207