LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
slaqr4.f
Go to the documentation of this file.
1 *> \brief \b SLAQR4 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 SLAQR4 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaqr4.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaqr4.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaqr4.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLAQR4( 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 * REAL H( LDH, * ), WI( * ), WORK( * ), WR( * ),
30 * $ Z( LDZ, * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> SLAQR4 implements one level of recursion for SLAQR0.
40 *> It is a complete implementation of the small bulge multi-shift
41 *> QR algorithm. It may be called by SLAQR0 and, for large enough
42 *> deflation window size, it may be called by SLAQR3. This
43 *> subroutine is identical to SLAQR0 except that it calls SLAQR2
44 *> instead of SLAQR3.
45 *>
46 *> SLAQR4 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 .GE. 0.
78 *> \endverbatim
79 *>
80 *> \param[in] ILO
81 *> \verbatim
82 *> ILO is INTEGER
83 *> \endverbatim
84 *>
85 *> \param[in] IHI
86 *> \verbatim
87 *> IHI is INTEGER
88 *> It is assumed that H is already upper triangular in rows
89 *> and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1,
90 *> H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
91 *> previous call to SGEBAL, and then passed to SGEHRD when the
92 *> matrix output by SGEBAL is reduced to Hessenberg form.
93 *> Otherwise, ILO and IHI should be set to 1 and N,
94 *> respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
95 *> If N = 0, then ILO = 1 and IHI = 0.
96 *> \endverbatim
97 *>
98 *> \param[in,out] H
99 *> \verbatim
100 *> H is REAL 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).LT.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.GT.0 is given under the
110 *> description of INFO below.)
111 *>
112 *> This subroutine may explicitly set H(i,j) = 0 for i.GT.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 .GE. max(1,N).
120 *> \endverbatim
121 *>
122 *> \param[out] WR
123 *> \verbatim
124 *> WR is REAL array, dimension (IHI)
125 *> \endverbatim
126 *>
127 *> \param[out] WI
128 *> \verbatim
129 *> WI is REAL 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) .GT. 0 and WI(i+1) .LT. 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 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N.
154 *> \endverbatim
155 *>
156 *> \param[in,out] Z
157 *> \verbatim
158 *> Z is REAL 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.GT.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.GE.MAX(1,IHIZ). Otherwize, LDZ.GE.1.
172 *> \endverbatim
173 *>
174 *> \param[out] WORK
175 *> \verbatim
176 *> WORK is REAL 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 .GE. 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 SLAQR4 does a workspace query.
190 *> In this case, SLAQR4 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 *> \verbatim
201 *> INFO is INTEGER
202 *> = 0: successful exit
203 *> .GT. 0: if INFO = i, SLAQR4 failed to compute all of
204 *> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
205 *> and WI contain those eigenvalues which have been
206 *> successfully computed. (Failures are rare.)
207 *>
208 *> If INFO .GT. 0 and WANT is .FALSE., then on exit,
209 *> the remaining unconverged eigenvalues are the eigen-
210 *> values of the upper Hessenberg matrix rows and
211 *> columns ILO through INFO of the final, output
212 *> value of H.
213 *>
214 *> If INFO .GT. 0 and WANTT is .TRUE., then on exit
215 *>
216 *> (*) (initial value of H)*U = U*(final value of H)
217 *>
218 *> where U is a orthogonal matrix. The final
219 *> value of H is upper Hessenberg and triangular in
220 *> rows and columns INFO+1 through IHI.
221 *>
222 *> If INFO .GT. 0 and WANTZ is .TRUE., then on exit
223 *>
224 *> (final value of Z(ILO:IHI,ILOZ:IHIZ)
225 *> = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U
226 *>
227 *> where U is the orthogonal matrix in (*) (regard-
228 *> less of the value of WANTT.)
229 *>
230 *> If INFO .GT. 0 and WANTZ is .FALSE., then Z is not
231 *> accessed.
232 *> \endverbatim
233 *
234 * Authors:
235 * ========
236 *
237 *> \author Univ. of Tennessee
238 *> \author Univ. of California Berkeley
239 *> \author Univ. of Colorado Denver
240 *> \author NAG Ltd.
241 *
242 *> \date September 2012
243 *
244 *> \ingroup realOTHERauxiliary
245 *
246 *> \par Contributors:
247 * ==================
248 *>
249 *> Karen Braman and Ralph Byers, Department of Mathematics,
250 *> University of Kansas, USA
251 *
252 *> \par References:
253 * ================
254 *>
255 *> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
256 *> Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
257 *> Performance, SIAM Journal of Matrix Analysis, volume 23, pages
258 *> 929--947, 2002.
259 *> \n
260 *> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
261 *> Algorithm Part II: Aggressive Early Deflation, SIAM Journal
262 *> of Matrix Analysis, volume 23, pages 948--973, 2002.
263 *>
264 * =====================================================================
265  SUBROUTINE slaqr4( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
266  $ iloz, ihiz, z, ldz, work, lwork, info )
267 *
268 * -- LAPACK auxiliary routine (version 3.4.2) --
269 * -- LAPACK is a software package provided by Univ. of Tennessee, --
270 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
271 * September 2012
272 *
273 * .. Scalar Arguments ..
274  INTEGER ihi, ihiz, ilo, iloz, info, ldh, ldz, lwork, n
275  LOGICAL wantt, wantz
276 * ..
277 * .. Array Arguments ..
278  REAL h( ldh, * ), wi( * ), work( * ), wr( * ),
279  $ z( ldz, * )
280 * ..
281 *
282 * ================================================================
283 *
284 * .. Parameters ..
285 *
286 * ==== Matrices of order NTINY or smaller must be processed by
287 * . SLAHQR because of insufficient subdiagonal scratch space.
288 * . (This is a hard limit.) ====
289  INTEGER ntiny
290  parameter( ntiny = 11 )
291 *
292 * ==== Exceptional deflation windows: try to cure rare
293 * . slow convergence by varying the size of the
294 * . deflation window after KEXNW iterations. ====
295  INTEGER kexnw
296  parameter( kexnw = 5 )
297 *
298 * ==== Exceptional shifts: try to cure rare slow convergence
299 * . with ad-hoc exceptional shifts every KEXSH iterations.
300 * . ====
301  INTEGER kexsh
302  parameter( kexsh = 6 )
303 *
304 * ==== The constants WILK1 and WILK2 are used to form the
305 * . exceptional shifts. ====
306  REAL wilk1, wilk2
307  parameter( wilk1 = 0.75e0, wilk2 = -0.4375e0 )
308  REAL zero, one
309  parameter( zero = 0.0e0, one = 1.0e0 )
310 * ..
311 * .. Local Scalars ..
312  REAL aa, bb, cc, cs, dd, sn, ss, swap
313  INTEGER i, inf, it, itmax, k, kacc22, kbot, kdu, ks,
314  $ kt, ktop, ku, kv, kwh, kwtop, kwv, ld, ls,
315  $ lwkopt, ndec, ndfl, nh, nho, nibble, nmin, ns,
316  $ nsmax, nsr, nve, nw, nwmax, nwr, nwupbd
317  LOGICAL sorted
318  CHARACTER jbcmpz*2
319 * ..
320 * .. External Functions ..
321  INTEGER ilaenv
322  EXTERNAL ilaenv
323 * ..
324 * .. Local Arrays ..
325  REAL zdum( 1, 1 )
326 * ..
327 * .. External Subroutines ..
328  EXTERNAL slacpy, slahqr, slanv2, slaqr2, slaqr5
329 * ..
330 * .. Intrinsic Functions ..
331  INTRINSIC abs, int, max, min, mod, real
332 * ..
333 * .. Executable Statements ..
334  info = 0
335 *
336 * ==== Quick return for N = 0: nothing to do. ====
337 *
338  IF( n.EQ.0 ) THEN
339  work( 1 ) = one
340  return
341  END IF
342 *
343  IF( n.LE.ntiny ) THEN
344 *
345 * ==== Tiny matrices must use SLAHQR. ====
346 *
347  lwkopt = 1
348  IF( lwork.NE.-1 )
349  $ CALL slahqr( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi,
350  $ iloz, ihiz, z, ldz, info )
351  ELSE
352 *
353 * ==== Use small bulge multi-shift QR with aggressive early
354 * . deflation on larger-than-tiny matrices. ====
355 *
356 * ==== Hope for the best. ====
357 *
358  info = 0
359 *
360 * ==== Set up job flags for ILAENV. ====
361 *
362  IF( wantt ) THEN
363  jbcmpz( 1: 1 ) = 'S'
364  ELSE
365  jbcmpz( 1: 1 ) = 'E'
366  END IF
367  IF( wantz ) THEN
368  jbcmpz( 2: 2 ) = 'V'
369  ELSE
370  jbcmpz( 2: 2 ) = 'N'
371  END IF
372 *
373 * ==== NWR = recommended deflation window size. At this
374 * . point, N .GT. NTINY = 11, so there is enough
375 * . subdiagonal workspace for NWR.GE.2 as required.
376 * . (In fact, there is enough subdiagonal space for
377 * . NWR.GE.3.) ====
378 *
379  nwr = ilaenv( 13, 'SLAQR4', jbcmpz, n, ilo, ihi, lwork )
380  nwr = max( 2, nwr )
381  nwr = min( ihi-ilo+1, ( n-1 ) / 3, nwr )
382 *
383 * ==== NSR = recommended number of simultaneous shifts.
384 * . At this point N .GT. NTINY = 11, so there is at
385 * . enough subdiagonal workspace for NSR to be even
386 * . and greater than or equal to two as required. ====
387 *
388  nsr = ilaenv( 15, 'SLAQR4', jbcmpz, n, ilo, ihi, lwork )
389  nsr = min( nsr, ( n+6 ) / 9, ihi-ilo )
390  nsr = max( 2, nsr-mod( nsr, 2 ) )
391 *
392 * ==== Estimate optimal workspace ====
393 *
394 * ==== Workspace query call to SLAQR2 ====
395 *
396  CALL slaqr2( wantt, wantz, n, ilo, ihi, nwr+1, h, ldh, iloz,
397  $ ihiz, z, ldz, ls, ld, wr, wi, h, ldh, n, h, ldh,
398  $ n, h, ldh, work, -1 )
399 *
400 * ==== Optimal workspace = MAX(SLAQR5, SLAQR2) ====
401 *
402  lwkopt = max( 3*nsr / 2, int( work( 1 ) ) )
403 *
404 * ==== Quick return in case of workspace query. ====
405 *
406  IF( lwork.EQ.-1 ) THEN
407  work( 1 ) = REAL( lwkopt )
408  return
409  END IF
410 *
411 * ==== SLAHQR/SLAQR0 crossover point ====
412 *
413  nmin = ilaenv( 12, 'SLAQR4', jbcmpz, n, ilo, ihi, lwork )
414  nmin = max( ntiny, nmin )
415 *
416 * ==== Nibble crossover point ====
417 *
418  nibble = ilaenv( 14, 'SLAQR4', jbcmpz, n, ilo, ihi, lwork )
419  nibble = max( 0, nibble )
420 *
421 * ==== Accumulate reflections during ttswp? Use block
422 * . 2-by-2 structure during matrix-matrix multiply? ====
423 *
424  kacc22 = ilaenv( 16, 'SLAQR4', jbcmpz, n, ilo, ihi, lwork )
425  kacc22 = max( 0, kacc22 )
426  kacc22 = min( 2, kacc22 )
427 *
428 * ==== NWMAX = the largest possible deflation window for
429 * . which there is sufficient workspace. ====
430 *
431  nwmax = min( ( n-1 ) / 3, lwork / 2 )
432  nw = nwmax
433 *
434 * ==== NSMAX = the Largest number of simultaneous shifts
435 * . for which there is sufficient workspace. ====
436 *
437  nsmax = min( ( n+6 ) / 9, 2*lwork / 3 )
438  nsmax = nsmax - mod( nsmax, 2 )
439 *
440 * ==== NDFL: an iteration count restarted at deflation. ====
441 *
442  ndfl = 1
443 *
444 * ==== ITMAX = iteration limit ====
445 *
446  itmax = max( 30, 2*kexsh )*max( 10, ( ihi-ilo+1 ) )
447 *
448 * ==== Last row and column in the active block ====
449 *
450  kbot = ihi
451 *
452 * ==== Main Loop ====
453 *
454  DO 80 it = 1, itmax
455 *
456 * ==== Done when KBOT falls below ILO ====
457 *
458  IF( kbot.LT.ilo )
459  $ go to 90
460 *
461 * ==== Locate active block ====
462 *
463  DO 10 k = kbot, ilo + 1, -1
464  IF( h( k, k-1 ).EQ.zero )
465  $ go to 20
466  10 continue
467  k = ilo
468  20 continue
469  ktop = k
470 *
471 * ==== Select deflation window size:
472 * . Typical Case:
473 * . If possible and advisable, nibble the entire
474 * . active block. If not, use size MIN(NWR,NWMAX)
475 * . or MIN(NWR+1,NWMAX) depending upon which has
476 * . the smaller corresponding subdiagonal entry
477 * . (a heuristic).
478 * .
479 * . Exceptional Case:
480 * . If there have been no deflations in KEXNW or
481 * . more iterations, then vary the deflation window
482 * . size. At first, because, larger windows are,
483 * . in general, more powerful than smaller ones,
484 * . rapidly increase the window to the maximum possible.
485 * . Then, gradually reduce the window size. ====
486 *
487  nh = kbot - ktop + 1
488  nwupbd = min( nh, nwmax )
489  IF( ndfl.LT.kexnw ) THEN
490  nw = min( nwupbd, nwr )
491  ELSE
492  nw = min( nwupbd, 2*nw )
493  END IF
494  IF( nw.LT.nwmax ) THEN
495  IF( nw.GE.nh-1 ) THEN
496  nw = nh
497  ELSE
498  kwtop = kbot - nw + 1
499  IF( abs( h( kwtop, kwtop-1 ) ).GT.
500  $ abs( h( kwtop-1, kwtop-2 ) ) )nw = nw + 1
501  END IF
502  END IF
503  IF( ndfl.LT.kexnw ) THEN
504  ndec = -1
505  ELSE IF( ndec.GE.0 .OR. nw.GE.nwupbd ) THEN
506  ndec = ndec + 1
507  IF( nw-ndec.LT.2 )
508  $ ndec = 0
509  nw = nw - ndec
510  END IF
511 *
512 * ==== Aggressive early deflation:
513 * . split workspace under the subdiagonal into
514 * . - an nw-by-nw work array V in the lower
515 * . left-hand-corner,
516 * . - an NW-by-at-least-NW-but-more-is-better
517 * . (NW-by-NHO) horizontal work array along
518 * . the bottom edge,
519 * . - an at-least-NW-but-more-is-better (NHV-by-NW)
520 * . vertical work array along the left-hand-edge.
521 * . ====
522 *
523  kv = n - nw + 1
524  kt = nw + 1
525  nho = ( n-nw-1 ) - kt + 1
526  kwv = nw + 2
527  nve = ( n-nw ) - kwv + 1
528 *
529 * ==== Aggressive early deflation ====
530 *
531  CALL slaqr2( wantt, wantz, n, ktop, kbot, nw, h, ldh, iloz,
532  $ ihiz, z, ldz, ls, ld, wr, wi, h( kv, 1 ), ldh,
533  $ nho, h( kv, kt ), ldh, nve, h( kwv, 1 ), ldh,
534  $ work, lwork )
535 *
536 * ==== Adjust KBOT accounting for new deflations. ====
537 *
538  kbot = kbot - ld
539 *
540 * ==== KS points to the shifts. ====
541 *
542  ks = kbot - ls + 1
543 *
544 * ==== Skip an expensive QR sweep if there is a (partly
545 * . heuristic) reason to expect that many eigenvalues
546 * . will deflate without it. Here, the QR sweep is
547 * . skipped if many eigenvalues have just been deflated
548 * . or if the remaining active block is small.
549 *
550  IF( ( ld.EQ.0 ) .OR. ( ( 100*ld.LE.nw*nibble ) .AND. ( kbot-
551  $ ktop+1.GT.min( nmin, nwmax ) ) ) ) THEN
552 *
553 * ==== NS = nominal number of simultaneous shifts.
554 * . This may be lowered (slightly) if SLAQR2
555 * . did not provide that many shifts. ====
556 *
557  ns = min( nsmax, nsr, max( 2, kbot-ktop ) )
558  ns = ns - mod( ns, 2 )
559 *
560 * ==== If there have been no deflations
561 * . in a multiple of KEXSH iterations,
562 * . then try exceptional shifts.
563 * . Otherwise use shifts provided by
564 * . SLAQR2 above or from the eigenvalues
565 * . of a trailing principal submatrix. ====
566 *
567  IF( mod( ndfl, kexsh ).EQ.0 ) THEN
568  ks = kbot - ns + 1
569  DO 30 i = kbot, max( ks+1, ktop+2 ), -2
570  ss = abs( h( i, i-1 ) ) + abs( h( i-1, i-2 ) )
571  aa = wilk1*ss + h( i, i )
572  bb = ss
573  cc = wilk2*ss
574  dd = aa
575  CALL slanv2( aa, bb, cc, dd, wr( i-1 ), wi( i-1 ),
576  $ wr( i ), wi( i ), cs, sn )
577  30 continue
578  IF( ks.EQ.ktop ) THEN
579  wr( ks+1 ) = h( ks+1, ks+1 )
580  wi( ks+1 ) = zero
581  wr( ks ) = wr( ks+1 )
582  wi( ks ) = wi( ks+1 )
583  END IF
584  ELSE
585 *
586 * ==== Got NS/2 or fewer shifts? Use SLAHQR
587 * . on a trailing principal submatrix to
588 * . get more. (Since NS.LE.NSMAX.LE.(N+6)/9,
589 * . there is enough space below the subdiagonal
590 * . to fit an NS-by-NS scratch array.) ====
591 *
592  IF( kbot-ks+1.LE.ns / 2 ) THEN
593  ks = kbot - ns + 1
594  kt = n - ns + 1
595  CALL slacpy( 'A', ns, ns, h( ks, ks ), ldh,
596  $ h( kt, 1 ), ldh )
597  CALL slahqr( .false., .false., ns, 1, ns,
598  $ h( kt, 1 ), ldh, wr( ks ), wi( ks ),
599  $ 1, 1, zdum, 1, inf )
600  ks = ks + inf
601 *
602 * ==== In case of a rare QR failure use
603 * . eigenvalues of the trailing 2-by-2
604 * . principal submatrix. ====
605 *
606  IF( ks.GE.kbot ) THEN
607  aa = h( kbot-1, kbot-1 )
608  cc = h( kbot, kbot-1 )
609  bb = h( kbot-1, kbot )
610  dd = h( kbot, kbot )
611  CALL slanv2( aa, bb, cc, dd, wr( kbot-1 ),
612  $ wi( kbot-1 ), wr( kbot ),
613  $ wi( kbot ), cs, sn )
614  ks = kbot - 1
615  END IF
616  END IF
617 *
618  IF( kbot-ks+1.GT.ns ) THEN
619 *
620 * ==== Sort the shifts (Helps a little)
621 * . Bubble sort keeps complex conjugate
622 * . pairs together. ====
623 *
624  sorted = .false.
625  DO 50 k = kbot, ks + 1, -1
626  IF( sorted )
627  $ go to 60
628  sorted = .true.
629  DO 40 i = ks, k - 1
630  IF( abs( wr( i ) )+abs( wi( i ) ).LT.
631  $ abs( wr( i+1 ) )+abs( wi( i+1 ) ) ) THEN
632  sorted = .false.
633 *
634  swap = wr( i )
635  wr( i ) = wr( i+1 )
636  wr( i+1 ) = swap
637 *
638  swap = wi( i )
639  wi( i ) = wi( i+1 )
640  wi( i+1 ) = swap
641  END IF
642  40 continue
643  50 continue
644  60 continue
645  END IF
646 *
647 * ==== Shuffle shifts into pairs of real shifts
648 * . and pairs of complex conjugate shifts
649 * . assuming complex conjugate shifts are
650 * . already adjacent to one another. (Yes,
651 * . they are.) ====
652 *
653  DO 70 i = kbot, ks + 2, -2
654  IF( wi( i ).NE.-wi( i-1 ) ) THEN
655 *
656  swap = wr( i )
657  wr( i ) = wr( i-1 )
658  wr( i-1 ) = wr( i-2 )
659  wr( i-2 ) = swap
660 *
661  swap = wi( i )
662  wi( i ) = wi( i-1 )
663  wi( i-1 ) = wi( i-2 )
664  wi( i-2 ) = swap
665  END IF
666  70 continue
667  END IF
668 *
669 * ==== If there are only two shifts and both are
670 * . real, then use only one. ====
671 *
672  IF( kbot-ks+1.EQ.2 ) THEN
673  IF( wi( kbot ).EQ.zero ) THEN
674  IF( abs( wr( kbot )-h( kbot, kbot ) ).LT.
675  $ abs( wr( kbot-1 )-h( kbot, kbot ) ) ) THEN
676  wr( kbot-1 ) = wr( kbot )
677  ELSE
678  wr( kbot ) = wr( kbot-1 )
679  END IF
680  END IF
681  END IF
682 *
683 * ==== Use up to NS of the the smallest magnatiude
684 * . shifts. If there aren't NS shifts available,
685 * . then use them all, possibly dropping one to
686 * . make the number of shifts even. ====
687 *
688  ns = min( ns, kbot-ks+1 )
689  ns = ns - mod( ns, 2 )
690  ks = kbot - ns + 1
691 *
692 * ==== Small-bulge multi-shift QR sweep:
693 * . split workspace under the subdiagonal into
694 * . - a KDU-by-KDU work array U in the lower
695 * . left-hand-corner,
696 * . - a KDU-by-at-least-KDU-but-more-is-better
697 * . (KDU-by-NHo) horizontal work array WH along
698 * . the bottom edge,
699 * . - and an at-least-KDU-but-more-is-better-by-KDU
700 * . (NVE-by-KDU) vertical work WV arrow along
701 * . the left-hand-edge. ====
702 *
703  kdu = 3*ns - 3
704  ku = n - kdu + 1
705  kwh = kdu + 1
706  nho = ( n-kdu+1-4 ) - ( kdu+1 ) + 1
707  kwv = kdu + 4
708  nve = n - kdu - kwv + 1
709 *
710 * ==== Small-bulge multi-shift QR sweep ====
711 *
712  CALL slaqr5( wantt, wantz, kacc22, n, ktop, kbot, ns,
713  $ wr( ks ), wi( ks ), h, ldh, iloz, ihiz, z,
714  $ ldz, work, 3, h( ku, 1 ), ldh, nve,
715  $ h( kwv, 1 ), ldh, nho, h( ku, kwh ), ldh )
716  END IF
717 *
718 * ==== Note progress (or the lack of it). ====
719 *
720  IF( ld.GT.0 ) THEN
721  ndfl = 1
722  ELSE
723  ndfl = ndfl + 1
724  END IF
725 *
726 * ==== End of main loop ====
727  80 continue
728 *
729 * ==== Iteration limit exceeded. Set INFO to show where
730 * . the problem occurred and exit. ====
731 *
732  info = kbot
733  90 continue
734  END IF
735 *
736 * ==== Return the optimal value of LWORK. ====
737 *
738  work( 1 ) = REAL( lwkopt )
739 *
740 * ==== End of SLAQR4 ====
741 *
742  END