LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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
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 slaqr2(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)
SLAQR2 performs the orthogonal similarity transformation of a Hessenberg matrix to detect and deflate...
Definition: slaqr2.f:280
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