LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
slaqr0.f
Go to the documentation of this file.
1 *> \brief \b SLAQR0 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 SLAQR0 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaqr0.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaqr0.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaqr0.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLAQR0( 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 *> SLAQR0 computes the eigenvalues of a Hessenberg matrix H
40 *> and, optionally, the matrices T and Z from the Schur decomposition
41 *> H = Z T Z**T, where T is an upper quasi-triangular matrix (the
42 *> Schur form), and Z is the orthogonal matrix of Schur vectors.
43 *>
44 *> Optionally Z may be postmultiplied into an input orthogonal
45 *> matrix Q so that this routine can give the Schur factorization
46 *> of a matrix A which has been reduced to the Hessenberg form H
47 *> by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T.
48 *> \endverbatim
49 *
50 * Arguments:
51 * ==========
52 *
53 *> \param[in] WANTT
54 *> \verbatim
55 *> WANTT is LOGICAL
56 *> = .TRUE. : the full Schur form T is required;
57 *> = .FALSE.: only eigenvalues are required.
58 *> \endverbatim
59 *>
60 *> \param[in] WANTZ
61 *> \verbatim
62 *> WANTZ is LOGICAL
63 *> = .TRUE. : the matrix of Schur vectors Z is required;
64 *> = .FALSE.: Schur vectors are not required.
65 *> \endverbatim
66 *>
67 *> \param[in] N
68 *> \verbatim
69 *> N is INTEGER
70 *> The order of the matrix H. N >= 0.
71 *> \endverbatim
72 *>
73 *> \param[in] ILO
74 *> \verbatim
75 *> ILO is INTEGER
76 *> \endverbatim
77 *>
78 *> \param[in] IHI
79 *> \verbatim
80 *> IHI is INTEGER
81 *> It is assumed that H is already upper triangular in rows
82 *> and columns 1:ILO-1 and IHI+1:N and, if ILO > 1,
83 *> H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
84 *> previous call to SGEBAL, and then passed to SGEHRD when the
85 *> matrix output by SGEBAL is reduced to Hessenberg form.
86 *> Otherwise, ILO and IHI should be set to 1 and N,
87 *> respectively. If N > 0, then 1 <= ILO <= IHI <= N.
88 *> If N = 0, then ILO = 1 and IHI = 0.
89 *> \endverbatim
90 *>
91 *> \param[in,out] H
92 *> \verbatim
93 *> H is REAL array, dimension (LDH,N)
94 *> On entry, the upper Hessenberg matrix H.
95 *> On exit, if INFO = 0 and WANTT is .TRUE., then H contains
96 *> the upper quasi-triangular matrix T from the Schur
97 *> decomposition (the Schur form); 2-by-2 diagonal blocks
98 *> (corresponding to complex conjugate pairs of eigenvalues)
99 *> are returned in standard form, with H(i,i) = H(i+1,i+1)
100 *> and H(i+1,i)*H(i,i+1) < 0. If INFO = 0 and WANTT is
101 *> .FALSE., then the contents of H are unspecified on exit.
102 *> (The output value of H when INFO > 0 is given under the
103 *> description of INFO below.)
104 *>
105 *> This subroutine may explicitly set H(i,j) = 0 for i > j and
106 *> j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N.
107 *> \endverbatim
108 *>
109 *> \param[in] LDH
110 *> \verbatim
111 *> LDH is INTEGER
112 *> The leading dimension of the array H. LDH >= max(1,N).
113 *> \endverbatim
114 *>
115 *> \param[out] WR
116 *> \verbatim
117 *> WR is REAL array, dimension (IHI)
118 *> \endverbatim
119 *>
120 *> \param[out] WI
121 *> \verbatim
122 *> WI is REAL array, dimension (IHI)
123 *> The real and imaginary parts, respectively, of the computed
124 *> eigenvalues of H(ILO:IHI,ILO:IHI) are stored in WR(ILO:IHI)
125 *> and WI(ILO:IHI). If two eigenvalues are computed as a
126 *> complex conjugate pair, they are stored in consecutive
127 *> elements of WR and WI, say the i-th and (i+1)th, with
128 *> WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., then
129 *> the eigenvalues are stored in the same order as on the
130 *> diagonal of the Schur form returned in H, with
131 *> WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal
132 *> block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and
133 *> WI(i+1) = -WI(i).
134 *> \endverbatim
135 *>
136 *> \param[in] ILOZ
137 *> \verbatim
138 *> ILOZ is INTEGER
139 *> \endverbatim
140 *>
141 *> \param[in] IHIZ
142 *> \verbatim
143 *> IHIZ is INTEGER
144 *> Specify the rows of Z to which transformations must be
145 *> applied if WANTZ is .TRUE..
146 *> 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
147 *> \endverbatim
148 *>
149 *> \param[in,out] Z
150 *> \verbatim
151 *> Z is REAL array, dimension (LDZ,IHI)
152 *> If WANTZ is .FALSE., then Z is not referenced.
153 *> If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is
154 *> replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the
155 *> orthogonal Schur factor of H(ILO:IHI,ILO:IHI).
156 *> (The output value of Z when INFO > 0 is given under
157 *> the description of INFO below.)
158 *> \endverbatim
159 *>
160 *> \param[in] LDZ
161 *> \verbatim
162 *> LDZ is INTEGER
163 *> The leading dimension of the array Z. if WANTZ is .TRUE.
164 *> then LDZ >= MAX(1,IHIZ). Otherwise, LDZ >= 1.
165 *> \endverbatim
166 *>
167 *> \param[out] WORK
168 *> \verbatim
169 *> WORK is REAL array, dimension LWORK
170 *> On exit, if LWORK = -1, WORK(1) returns an estimate of
171 *> the optimal value for LWORK.
172 *> \endverbatim
173 *>
174 *> \param[in] LWORK
175 *> \verbatim
176 *> LWORK is INTEGER
177 *> The dimension of the array WORK. LWORK >= max(1,N)
178 *> is sufficient, but LWORK typically as large as 6*N may
179 *> be required for optimal performance. A workspace query
180 *> to determine the optimal workspace size is recommended.
181 *>
182 *> If LWORK = -1, then SLAQR0 does a workspace query.
183 *> In this case, SLAQR0 checks the input parameters and
184 *> estimates the optimal workspace size for the given
185 *> values of N, ILO and IHI. The estimate is returned
186 *> in WORK(1). No error message related to LWORK is
187 *> issued by XERBLA. Neither H nor Z are accessed.
188 *> \endverbatim
189 *>
190 *> \param[out] INFO
191 *> \verbatim
192 *> INFO is INTEGER
193 *> = 0: successful exit
194 *> > 0: if INFO = i, SLAQR0 failed to compute all of
195 *> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
196 *> and WI contain those eigenvalues which have been
197 *> successfully computed. (Failures are rare.)
198 *>
199 *> If INFO > 0 and WANT is .FALSE., then on exit,
200 *> the remaining unconverged eigenvalues are the eigen-
201 *> values of the upper Hessenberg matrix rows and
202 *> columns ILO through INFO of the final, output
203 *> value of H.
204 *>
205 *> If INFO > 0 and WANTT is .TRUE., then on exit
206 *>
207 *> (*) (initial value of H)*U = U*(final value of H)
208 *>
209 *> where U is an orthogonal matrix. The final
210 *> value of H is upper Hessenberg and quasi-triangular
211 *> in rows and columns INFO+1 through IHI.
212 *>
213 *> If INFO > 0 and WANTZ is .TRUE., then on exit
214 *>
215 *> (final value of Z(ILO:IHI,ILOZ:IHIZ)
216 *> = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U
217 *>
218 *> where U is the orthogonal matrix in (*) (regard-
219 *> less of the value of WANTT.)
220 *>
221 *> If INFO > 0 and WANTZ is .FALSE., then Z is not
222 *> accessed.
223 *> \endverbatim
224 *
225 * Authors:
226 * ========
227 *
228 *> \author Univ. of Tennessee
229 *> \author Univ. of California Berkeley
230 *> \author Univ. of Colorado Denver
231 *> \author NAG Ltd.
232 *
233 *> \ingroup realOTHERauxiliary
234 *
235 *> \par Contributors:
236 * ==================
237 *>
238 *> Karen Braman and Ralph Byers, Department of Mathematics,
239 *> University of Kansas, USA
240 *
241 *> \par References:
242 * ================
243 *>
244 *> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
245 *> Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
246 *> Performance, SIAM Journal of Matrix Analysis, volume 23, pages
247 *> 929--947, 2002.
248 *> \n
249 *> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
250 *> Algorithm Part II: Aggressive Early Deflation, SIAM Journal
251 *> of Matrix Analysis, volume 23, pages 948--973, 2002.
252 *>
253 * =====================================================================
254  SUBROUTINE slaqr0( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
255  $ ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO )
256 *
257 * -- LAPACK auxiliary routine --
258 * -- LAPACK is a software package provided by Univ. of Tennessee, --
259 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
260 *
261 * .. Scalar Arguments ..
262  INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N
263  LOGICAL WANTT, WANTZ
264 * ..
265 * .. Array Arguments ..
266  REAL H( LDH, * ), WI( * ), WORK( * ), WR( * ),
267  $ z( ldz, * )
268 * ..
269 *
270 * ================================================================
271 * .. Parameters ..
272 *
273 * ==== Matrices of order NTINY or smaller must be processed by
274 * . SLAHQR because of insufficient subdiagonal scratch space.
275 * . (This is a hard limit.) ====
276  INTEGER NTINY
277  parameter( ntiny = 15 )
278 *
279 * ==== Exceptional deflation windows: try to cure rare
280 * . slow convergence by varying the size of the
281 * . deflation window after KEXNW iterations. ====
282  INTEGER KEXNW
283  parameter( kexnw = 5 )
284 *
285 * ==== Exceptional shifts: try to cure rare slow convergence
286 * . with ad-hoc exceptional shifts every KEXSH iterations.
287 * . ====
288  INTEGER KEXSH
289  parameter( kexsh = 6 )
290 *
291 * ==== The constants WILK1 and WILK2 are used to form the
292 * . exceptional shifts. ====
293  REAL WILK1, WILK2
294  parameter( wilk1 = 0.75e0, wilk2 = -0.4375e0 )
295  REAL ZERO, ONE
296  parameter( zero = 0.0e0, one = 1.0e0 )
297 * ..
298 * .. Local Scalars ..
299  REAL AA, BB, CC, CS, DD, SN, SS, SWAP
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  REAL ZDUM( 1, 1 )
313 * ..
314 * .. External Subroutines ..
315  EXTERNAL slacpy, slahqr, slanv2, slaqr3, slaqr4, slaqr5
316 * ..
317 * .. Intrinsic Functions ..
318  INTRINSIC abs, int, max, min, mod, real
319 * ..
320 * .. Executable Statements ..
321  info = 0
322 *
323 * ==== Quick return for N = 0: nothing to do. ====
324 *
325  IF( n.EQ.0 ) THEN
326  work( 1 ) = one
327  RETURN
328  END IF
329 *
330  IF( n.LE.ntiny ) THEN
331 *
332 * ==== Tiny matrices must use SLAHQR. ====
333 *
334  lwkopt = 1
335  IF( lwork.NE.-1 )
336  $ CALL slahqr( wantt, wantz, n, ilo, ihi, h, ldh, wr, wi,
337  $ iloz, ihiz, z, ldz, info )
338  ELSE
339 *
340 * ==== Use small bulge multi-shift QR with aggressive early
341 * . deflation on larger-than-tiny matrices. ====
342 *
343 * ==== Hope for the best. ====
344 *
345  info = 0
346 *
347 * ==== Set up job flags for ILAENV. ====
348 *
349  IF( wantt ) THEN
350  jbcmpz( 1: 1 ) = 'S'
351  ELSE
352  jbcmpz( 1: 1 ) = 'E'
353  END IF
354  IF( wantz ) THEN
355  jbcmpz( 2: 2 ) = 'V'
356  ELSE
357  jbcmpz( 2: 2 ) = 'N'
358  END IF
359 *
360 * ==== NWR = recommended deflation window size. At this
361 * . point, N .GT. NTINY = 15, so there is enough
362 * . subdiagonal workspace for NWR.GE.2 as required.
363 * . (In fact, there is enough subdiagonal space for
364 * . NWR.GE.4.) ====
365 *
366  nwr = ilaenv( 13, 'SLAQR0', jbcmpz, n, ilo, ihi, lwork )
367  nwr = max( 2, nwr )
368  nwr = min( ihi-ilo+1, ( n-1 ) / 3, nwr )
369 *
370 * ==== NSR = recommended number of simultaneous shifts.
371 * . At this point N .GT. NTINY = 15, so there is at
372 * . enough subdiagonal workspace for NSR to be even
373 * . and greater than or equal to two as required. ====
374 *
375  nsr = ilaenv( 15, 'SLAQR0', jbcmpz, n, ilo, ihi, lwork )
376  nsr = min( nsr, ( n-3 ) / 6, ihi-ilo )
377  nsr = max( 2, nsr-mod( nsr, 2 ) )
378 *
379 * ==== Estimate optimal workspace ====
380 *
381 * ==== Workspace query call to SLAQR3 ====
382 *
383  CALL slaqr3( wantt, wantz, n, ilo, ihi, nwr+1, h, ldh, iloz,
384  $ ihiz, z, ldz, ls, ld, wr, wi, h, ldh, n, h, ldh,
385  $ n, h, ldh, work, -1 )
386 *
387 * ==== Optimal workspace = MAX(SLAQR5, SLAQR3) ====
388 *
389  lwkopt = max( 3*nsr / 2, int( work( 1 ) ) )
390 *
391 * ==== Quick return in case of workspace query. ====
392 *
393  IF( lwork.EQ.-1 ) THEN
394  work( 1 ) = real( lwkopt )
395  RETURN
396  END IF
397 *
398 * ==== SLAHQR/SLAQR0 crossover point ====
399 *
400  nmin = ilaenv( 12, 'SLAQR0', jbcmpz, n, ilo, ihi, lwork )
401  nmin = max( ntiny, nmin )
402 *
403 * ==== Nibble crossover point ====
404 *
405  nibble = ilaenv( 14, 'SLAQR0', jbcmpz, n, ilo, ihi, lwork )
406  nibble = max( 0, nibble )
407 *
408 * ==== Accumulate reflections during ttswp? Use block
409 * . 2-by-2 structure during matrix-matrix multiply? ====
410 *
411  kacc22 = ilaenv( 16, 'SLAQR0', jbcmpz, n, ilo, ihi, lwork )
412  kacc22 = max( 0, kacc22 )
413  kacc22 = min( 2, kacc22 )
414 *
415 * ==== NWMAX = the largest possible deflation window for
416 * . which there is sufficient workspace. ====
417 *
418  nwmax = min( ( n-1 ) / 3, lwork / 2 )
419  nw = nwmax
420 *
421 * ==== NSMAX = the Largest number of simultaneous shifts
422 * . for which there is sufficient workspace. ====
423 *
424  nsmax = min( ( n-3 ) / 6, 2*lwork / 3 )
425  nsmax = nsmax - mod( nsmax, 2 )
426 *
427 * ==== NDFL: an iteration count restarted at deflation. ====
428 *
429  ndfl = 1
430 *
431 * ==== ITMAX = iteration limit ====
432 *
433  itmax = max( 30, 2*kexsh )*max( 10, ( ihi-ilo+1 ) )
434 *
435 * ==== Last row and column in the active block ====
436 *
437  kbot = ihi
438 *
439 * ==== Main Loop ====
440 *
441  DO 80 it = 1, itmax
442 *
443 * ==== Done when KBOT falls below ILO ====
444 *
445  IF( kbot.LT.ilo )
446  $ GO TO 90
447 *
448 * ==== Locate active block ====
449 *
450  DO 10 k = kbot, ilo + 1, -1
451  IF( h( k, k-1 ).EQ.zero )
452  $ GO TO 20
453  10 CONTINUE
454  k = ilo
455  20 CONTINUE
456  ktop = k
457 *
458 * ==== Select deflation window size:
459 * . Typical Case:
460 * . If possible and advisable, nibble the entire
461 * . active block. If not, use size MIN(NWR,NWMAX)
462 * . or MIN(NWR+1,NWMAX) depending upon which has
463 * . the smaller corresponding subdiagonal entry
464 * . (a heuristic).
465 * .
466 * . Exceptional Case:
467 * . If there have been no deflations in KEXNW or
468 * . more iterations, then vary the deflation window
469 * . size. At first, because, larger windows are,
470 * . in general, more powerful than smaller ones,
471 * . rapidly increase the window to the maximum possible.
472 * . Then, gradually reduce the window size. ====
473 *
474  nh = kbot - ktop + 1
475  nwupbd = min( nh, nwmax )
476  IF( ndfl.LT.kexnw ) THEN
477  nw = min( nwupbd, nwr )
478  ELSE
479  nw = min( nwupbd, 2*nw )
480  END IF
481  IF( nw.LT.nwmax ) THEN
482  IF( nw.GE.nh-1 ) THEN
483  nw = nh
484  ELSE
485  kwtop = kbot - nw + 1
486  IF( abs( h( kwtop, kwtop-1 ) ).GT.
487  $ abs( h( kwtop-1, kwtop-2 ) ) )nw = nw + 1
488  END IF
489  END IF
490  IF( ndfl.LT.kexnw ) THEN
491  ndec = -1
492  ELSE IF( ndec.GE.0 .OR. nw.GE.nwupbd ) THEN
493  ndec = ndec + 1
494  IF( nw-ndec.LT.2 )
495  $ ndec = 0
496  nw = nw - ndec
497  END IF
498 *
499 * ==== Aggressive early deflation:
500 * . split workspace under the subdiagonal into
501 * . - an nw-by-nw work array V in the lower
502 * . left-hand-corner,
503 * . - an NW-by-at-least-NW-but-more-is-better
504 * . (NW-by-NHO) horizontal work array along
505 * . the bottom edge,
506 * . - an at-least-NW-but-more-is-better (NHV-by-NW)
507 * . vertical work array along the left-hand-edge.
508 * . ====
509 *
510  kv = n - nw + 1
511  kt = nw + 1
512  nho = ( n-nw-1 ) - kt + 1
513  kwv = nw + 2
514  nve = ( n-nw ) - kwv + 1
515 *
516 * ==== Aggressive early deflation ====
517 *
518  CALL slaqr3( wantt, wantz, n, ktop, kbot, nw, h, ldh, iloz,
519  $ ihiz, z, ldz, ls, ld, wr, wi, h( kv, 1 ), ldh,
520  $ nho, h( kv, kt ), ldh, nve, h( kwv, 1 ), ldh,
521  $ work, lwork )
522 *
523 * ==== Adjust KBOT accounting for new deflations. ====
524 *
525  kbot = kbot - ld
526 *
527 * ==== KS points to the shifts. ====
528 *
529  ks = kbot - ls + 1
530 *
531 * ==== Skip an expensive QR sweep if there is a (partly
532 * . heuristic) reason to expect that many eigenvalues
533 * . will deflate without it. Here, the QR sweep is
534 * . skipped if many eigenvalues have just been deflated
535 * . or if the remaining active block is small.
536 *
537  IF( ( ld.EQ.0 ) .OR. ( ( 100*ld.LE.nw*nibble ) .AND. ( kbot-
538  $ ktop+1.GT.min( nmin, nwmax ) ) ) ) THEN
539 *
540 * ==== NS = nominal number of simultaneous shifts.
541 * . This may be lowered (slightly) if SLAQR3
542 * . did not provide that many shifts. ====
543 *
544  ns = min( nsmax, nsr, max( 2, kbot-ktop ) )
545  ns = ns - mod( ns, 2 )
546 *
547 * ==== If there have been no deflations
548 * . in a multiple of KEXSH iterations,
549 * . then try exceptional shifts.
550 * . Otherwise use shifts provided by
551 * . SLAQR3 above or from the eigenvalues
552 * . of a trailing principal submatrix. ====
553 *
554  IF( mod( ndfl, kexsh ).EQ.0 ) THEN
555  ks = kbot - ns + 1
556  DO 30 i = kbot, max( ks+1, ktop+2 ), -2
557  ss = abs( h( i, i-1 ) ) + abs( h( i-1, i-2 ) )
558  aa = wilk1*ss + h( i, i )
559  bb = ss
560  cc = wilk2*ss
561  dd = aa
562  CALL slanv2( aa, bb, cc, dd, wr( i-1 ), wi( i-1 ),
563  $ wr( i ), wi( i ), cs, sn )
564  30 CONTINUE
565  IF( ks.EQ.ktop ) THEN
566  wr( ks+1 ) = h( ks+1, ks+1 )
567  wi( ks+1 ) = zero
568  wr( ks ) = wr( ks+1 )
569  wi( ks ) = wi( ks+1 )
570  END IF
571  ELSE
572 *
573 * ==== Got NS/2 or fewer shifts? Use SLAQR4 or
574 * . SLAHQR on a trailing principal submatrix to
575 * . get more. (Since NS.LE.NSMAX.LE.(N-3)/6,
576 * . there is enough space below the subdiagonal
577 * . to fit an NS-by-NS scratch array.) ====
578 *
579  IF( kbot-ks+1.LE.ns / 2 ) THEN
580  ks = kbot - ns + 1
581  kt = n - ns + 1
582  CALL slacpy( 'A', ns, ns, h( ks, ks ), ldh,
583  $ h( kt, 1 ), ldh )
584  IF( ns.GT.nmin ) THEN
585  CALL slaqr4( .false., .false., ns, 1, ns,
586  $ h( kt, 1 ), ldh, wr( ks ),
587  $ wi( ks ), 1, 1, zdum, 1, work,
588  $ lwork, inf )
589  ELSE
590  CALL slahqr( .false., .false., ns, 1, ns,
591  $ h( kt, 1 ), ldh, wr( ks ),
592  $ wi( ks ), 1, 1, zdum, 1, inf )
593  END IF
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 slanv2( 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 slaqr5( 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 ) = real( lwkopt )
733 *
734 * ==== End of SLAQR0 ====
735 *
736  END
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
subroutine slanv2(A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN)
SLANV2 computes the Schur factorization of a real 2-by-2 nonsymmetric matrix in standard form.
Definition: slanv2.f:127
subroutine slaqr0(WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)
SLAQR0 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...
Definition: slaqr0.f:256
subroutine slaqr4(WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)
SLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...
Definition: slaqr4.f:265
subroutine slaqr3(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)
SLAQR3 performs the orthogonal similarity transformation of a Hessenberg matrix to detect and deflate...
Definition: slaqr3.f:275
subroutine slaqr5(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)
SLAQR5 performs a single small-bulge multi-shift QR sweep.
Definition: slaqr5.f:265
subroutine slahqr(WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, ILOZ, IHIZ, Z, LDZ, INFO)
SLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix,...
Definition: slahqr.f:207