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