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