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