LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
zlaqr4.f
Go to the documentation of this file.
1*> \brief \b ZLAQR4 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/zlaqr4.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaqr4.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaqr4.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZLAQR4( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
22* 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* COMPLEX*16 H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> ZLAQR4 implements one level of recursion for ZLAQR0.
39*> It is a complete implementation of the small bulge multi-shift
40*> QR algorithm. It may be called by ZLAQR0 and, for large enough
41*> deflation window size, it may be called by ZLAQR3. This
42*> subroutine is identical to ZLAQR0 except that it calls ZLAQR2
44*>
45*> ZLAQR4 computes the eigenvalues of a Hessenberg matrix H
46*> and, optionally, the matrices T and Z from the Schur decomposition
47*> H = Z T Z**H, where T is an upper triangular matrix (the
48*> Schur form), and Z is the unitary matrix of Schur vectors.
49*>
50*> Optionally Z may be postmultiplied into an input unitary
51*> matrix Q so that this routine can give the Schur factorization
52*> of a matrix A which has been reduced to the Hessenberg form H
53*> by the unitary matrix Q: A = Q*H*Q**H = (QZ)*H*(QZ)**H.
54*> \endverbatim
55*
56* Arguments:
57* ==========
58*
59*> \param[in] WANTT
60*> \verbatim
61*> WANTT is LOGICAL
62*> = .TRUE. : the full Schur form T is required;
63*> = .FALSE.: only eigenvalues are required.
64*> \endverbatim
65*>
66*> \param[in] WANTZ
67*> \verbatim
68*> WANTZ is LOGICAL
69*> = .TRUE. : the matrix of Schur vectors Z is required;
70*> = .FALSE.: Schur vectors are not required.
71*> \endverbatim
72*>
73*> \param[in] N
74*> \verbatim
75*> N is INTEGER
76*> The order of the matrix H. N >= 0.
77*> \endverbatim
78*>
79*> \param[in] ILO
80*> \verbatim
81*> ILO is INTEGER
82*> \endverbatim
83*>
84*> \param[in] IHI
85*> \verbatim
86*> IHI is INTEGER
87*> It is assumed that H is already upper triangular in rows
88*> and columns 1:ILO-1 and IHI+1:N and, if ILO > 1,
89*> H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
90*> previous call to ZGEBAL, and then passed to ZGEHRD when the
91*> matrix output by ZGEBAL is reduced to Hessenberg form.
92*> Otherwise, ILO and IHI should be set to 1 and N,
93*> respectively. If N > 0, then 1 <= ILO <= IHI <= N.
94*> If N = 0, then ILO = 1 and IHI = 0.
95*> \endverbatim
96*>
97*> \param[in,out] H
98*> \verbatim
99*> H is COMPLEX*16 array, dimension (LDH,N)
100*> On entry, the upper Hessenberg matrix H.
101*> On exit, if INFO = 0 and WANTT is .TRUE., then H
102*> contains the upper triangular matrix T from the Schur
103*> decomposition (the Schur form). If INFO = 0 and WANT is
104*> .FALSE., then the contents of H are unspecified on exit.
105*> (The output value of H when INFO > 0 is given under the
106*> description of INFO below.)
107*>
108*> This subroutine may explicitly set H(i,j) = 0 for i > j and
109*> j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N.
110*> \endverbatim
111*>
112*> \param[in] LDH
113*> \verbatim
114*> LDH is INTEGER
115*> The leading dimension of the array H. LDH >= max(1,N).
116*> \endverbatim
117*>
118*> \param[out] W
119*> \verbatim
120*> W is COMPLEX*16 array, dimension (N)
121*> The computed eigenvalues of H(ILO:IHI,ILO:IHI) are stored
122*> in W(ILO:IHI). If WANTT is .TRUE., then the eigenvalues are
123*> stored in the same order as on the diagonal of the Schur
124*> form returned in H, with W(i) = H(i,i).
125*> \endverbatim
126*>
127*> \param[in] ILOZ
128*> \verbatim
129*> ILOZ is INTEGER
130*> \endverbatim
131*>
132*> \param[in] IHIZ
133*> \verbatim
134*> IHIZ is INTEGER
135*> Specify the rows of Z to which transformations must be
136*> applied if WANTZ is .TRUE..
137*> 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
138*> \endverbatim
139*>
140*> \param[in,out] Z
141*> \verbatim
142*> Z is COMPLEX*16 array, dimension (LDZ,IHI)
143*> If WANTZ is .FALSE., then Z is not referenced.
144*> If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is
145*> replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the
146*> orthogonal Schur factor of H(ILO:IHI,ILO:IHI).
147*> (The output value of Z when INFO > 0 is given under
148*> the description of INFO below.)
149*> \endverbatim
150*>
151*> \param[in] LDZ
152*> \verbatim
153*> LDZ is INTEGER
154*> The leading dimension of the array Z. if WANTZ is .TRUE.
155*> then LDZ >= MAX(1,IHIZ). Otherwise, LDZ >= 1.
156*> \endverbatim
157*>
158*> \param[out] WORK
159*> \verbatim
160*> WORK is COMPLEX*16 array, dimension LWORK
161*> On exit, if LWORK = -1, WORK(1) returns an estimate of
162*> the optimal value for LWORK.
163*> \endverbatim
164*>
165*> \param[in] LWORK
166*> \verbatim
167*> LWORK is INTEGER
168*> The dimension of the array WORK. LWORK >= max(1,N)
169*> is sufficient, but LWORK typically as large as 6*N may
170*> be required for optimal performance. A workspace query
171*> to determine the optimal workspace size is recommended.
172*>
173*> If LWORK = -1, then ZLAQR4 does a workspace query.
174*> In this case, ZLAQR4 checks the input parameters and
175*> estimates the optimal workspace size for the given
176*> values of N, ILO and IHI. The estimate is returned
177*> in WORK(1). No error message related to LWORK is
178*> issued by XERBLA. Neither H nor Z are accessed.
179*> \endverbatim
180*>
181*> \param[out] INFO
182*> \verbatim
183*> INFO is INTEGER
184*> = 0: successful exit
185*> > 0: if INFO = i, ZLAQR4 failed to compute all of
186*> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
187*> and WI contain those eigenvalues which have been
188*> successfully computed. (Failures are rare.)
189*>
190*> If INFO > 0 and WANT is .FALSE., then on exit,
191*> the remaining unconverged eigenvalues are the eigen-
192*> values of the upper Hessenberg matrix rows and
193*> columns ILO through INFO of the final, output
194*> value of H.
195*>
196*> If INFO > 0 and WANTT is .TRUE., then on exit
197*>
198*> (*) (initial value of H)*U = U*(final value of H)
199*>
200*> where U is a unitary matrix. The final
201*> value of H is upper Hessenberg and triangular in
202*> rows and columns INFO+1 through IHI.
203*>
204*> If INFO > 0 and WANTZ is .TRUE., then on exit
205*>
206*> (final value of Z(ILO:IHI,ILOZ:IHIZ)
207*> = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U
208*>
209*> where U is the unitary matrix in (*) (regard-
210*> less of the value of WANTT.)
211*>
212*> If INFO > 0 and WANTZ is .FALSE., then Z is not
213*> accessed.
214*> \endverbatim
215*
216* Authors:
217* ========
218*
219*> \author Univ. of Tennessee
220*> \author Univ. of California Berkeley
221*> \author Univ. of Colorado Denver
222*> \author NAG Ltd.
223*
224*> \ingroup complex16OTHERauxiliary
225*
226*> \par Contributors:
227* ==================
228*>
229*> Karen Braman and Ralph Byers, Department of Mathematics,
230*> University of Kansas, USA
231*
232*> \par References:
233* ================
234*>
235*> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
236*> Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
237*> Performance, SIAM Journal of Matrix Analysis, volume 23, pages
238*> 929--947, 2002.
239*> \n
240*> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
241*> Algorithm Part II: Aggressive Early Deflation, SIAM Journal
242*> of Matrix Analysis, volume 23, pages 948--973, 2002.
243*>
244* =====================================================================
245 SUBROUTINE zlaqr4( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
246 \$ IHIZ, Z, LDZ, WORK, LWORK, INFO )
247*
248* -- LAPACK auxiliary routine --
249* -- LAPACK is a software package provided by Univ. of Tennessee, --
250* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
251*
252* .. Scalar Arguments ..
253 INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N
254 LOGICAL WANTT, WANTZ
255* ..
256* .. Array Arguments ..
257 COMPLEX*16 H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * )
258* ..
259*
260* ================================================================
261*
262* .. Parameters ..
263*
264* ==== Matrices of order NTINY or smaller must be processed by
265* . ZLAHQR because of insufficient subdiagonal scratch space.
266* . (This is a hard limit.) ====
267 INTEGER NTINY
268 parameter( ntiny = 15 )
269*
270* ==== Exceptional deflation windows: try to cure rare
271* . slow convergence by varying the size of the
272* . deflation window after KEXNW iterations. ====
273 INTEGER KEXNW
274 parameter( kexnw = 5 )
275*
276* ==== Exceptional shifts: try to cure rare slow convergence
277* . with ad-hoc exceptional shifts every KEXSH iterations.
278* . ====
279 INTEGER KEXSH
280 parameter( kexsh = 6 )
281*
282* ==== The constant WILK1 is used to form the exceptional
283* . shifts. ====
284 DOUBLE PRECISION WILK1
285 parameter( wilk1 = 0.75d0 )
286 COMPLEX*16 ZERO, ONE
287 parameter( zero = ( 0.0d0, 0.0d0 ),
288 \$ one = ( 1.0d0, 0.0d0 ) )
289 DOUBLE PRECISION TWO
290 parameter( two = 2.0d0 )
291* ..
292* .. Local Scalars ..
293 COMPLEX*16 AA, BB, CC, CDUM, DD, DET, RTDISC, SWAP, TR2
294 DOUBLE PRECISION S
295 INTEGER I, INF, IT, ITMAX, K, KACC22, KBOT, KDU, KS,
296 \$ kt, ktop, ku, kv, kwh, kwtop, kwv, ld, ls,
297 \$ lwkopt, ndec, ndfl, nh, nho, nibble, nmin, ns,
298 \$ nsmax, nsr, nve, nw, nwmax, nwr, nwupbd
299 LOGICAL SORTED
300 CHARACTER JBCMPZ*2
301* ..
302* .. External Functions ..
303 INTEGER ILAENV
304 EXTERNAL ilaenv
305* ..
306* .. Local Arrays ..
307 COMPLEX*16 ZDUM( 1, 1 )
308* ..
309* .. External Subroutines ..
310 EXTERNAL zlacpy, zlahqr, zlaqr2, zlaqr5
311* ..
312* .. Intrinsic Functions ..
313 INTRINSIC abs, dble, dcmplx, dimag, int, max, min, mod,
314 \$ sqrt
315* ..
316* .. Statement Functions ..
317 DOUBLE PRECISION CABS1
318* ..
319* .. Statement Function definitions ..
320 cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
321* ..
322* .. Executable Statements ..
323 info = 0
324*
325* ==== Quick return for N = 0: nothing to do. ====
326*
327 IF( n.EQ.0 ) THEN
328 work( 1 ) = one
329 RETURN
330 END IF
331*
332 IF( n.LE.ntiny ) THEN
333*
334* ==== Tiny matrices must use ZLAHQR. ====
335*
336 lwkopt = 1
337 IF( lwork.NE.-1 )
338 \$ CALL zlahqr( wantt, wantz, n, ilo, ihi, h, ldh, w, iloz,
339 \$ ihiz, z, ldz, info )
340 ELSE
341*
342* ==== Use small bulge multi-shift QR with aggressive early
343* . deflation on larger-than-tiny matrices. ====
344*
345* ==== Hope for the best. ====
346*
347 info = 0
348*
349* ==== Set up job flags for ILAENV. ====
350*
351 IF( wantt ) THEN
352 jbcmpz( 1: 1 ) = 'S'
353 ELSE
354 jbcmpz( 1: 1 ) = 'E'
355 END IF
356 IF( wantz ) THEN
357 jbcmpz( 2: 2 ) = 'V'
358 ELSE
359 jbcmpz( 2: 2 ) = 'N'
360 END IF
361*
362* ==== NWR = recommended deflation window size. At this
363* . point, N .GT. NTINY = 15, so there is enough
364* . subdiagonal workspace for NWR.GE.2 as required.
365* . (In fact, there is enough subdiagonal space for
366* . NWR.GE.4.) ====
367*
368 nwr = ilaenv( 13, 'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
369 nwr = max( 2, nwr )
370 nwr = min( ihi-ilo+1, ( n-1 ) / 3, nwr )
371*
372* ==== NSR = recommended number of simultaneous shifts.
373* . At this point N .GT. NTINY = 15, so there is at
374* . enough subdiagonal workspace for NSR to be even
375* . and greater than or equal to two as required. ====
376*
377 nsr = ilaenv( 15, 'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
378 nsr = min( nsr, ( n-3 ) / 6, ihi-ilo )
379 nsr = max( 2, nsr-mod( nsr, 2 ) )
380*
381* ==== Estimate optimal workspace ====
382*
383* ==== Workspace query call to ZLAQR2 ====
384*
385 CALL zlaqr2( wantt, wantz, n, ilo, ihi, nwr+1, h, ldh, iloz,
386 \$ ihiz, z, ldz, ls, ld, w, h, ldh, n, h, ldh, n, h,
387 \$ ldh, work, -1 )
388*
389* ==== Optimal workspace = MAX(ZLAQR5, ZLAQR2) ====
390*
391 lwkopt = max( 3*nsr / 2, int( work( 1 ) ) )
392*
393* ==== Quick return in case of workspace query. ====
394*
395 IF( lwork.EQ.-1 ) THEN
396 work( 1 ) = dcmplx( lwkopt, 0 )
397 RETURN
398 END IF
399*
400* ==== ZLAHQR/ZLAQR0 crossover point ====
401*
402 nmin = ilaenv( 12, 'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
403 nmin = max( ntiny, nmin )
404*
405* ==== Nibble crossover point ====
406*
407 nibble = ilaenv( 14, 'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
408 nibble = max( 0, nibble )
409*
410* ==== Accumulate reflections during ttswp? Use block
411* . 2-by-2 structure during matrix-matrix multiply? ====
412*
413 kacc22 = ilaenv( 16, 'ZLAQR4', jbcmpz, n, ilo, ihi, lwork )
414 kacc22 = max( 0, kacc22 )
415 kacc22 = min( 2, kacc22 )
416*
417* ==== NWMAX = the largest possible deflation window for
418* . which there is sufficient workspace. ====
419*
420 nwmax = min( ( n-1 ) / 3, lwork / 2 )
421 nw = nwmax
422*
423* ==== NSMAX = the Largest number of simultaneous shifts
424* . for which there is sufficient workspace. ====
425*
426 nsmax = min( ( n-3 ) / 6, 2*lwork / 3 )
427 nsmax = nsmax - mod( nsmax, 2 )
428*
429* ==== NDFL: an iteration count restarted at deflation. ====
430*
431 ndfl = 1
432*
433* ==== ITMAX = iteration limit ====
434*
435 itmax = max( 30, 2*kexsh )*max( 10, ( ihi-ilo+1 ) )
436*
437* ==== Last row and column in the active block ====
438*
439 kbot = ihi
440*
441* ==== Main Loop ====
442*
443 DO 70 it = 1, itmax
444*
445* ==== Done when KBOT falls below ILO ====
446*
447 IF( kbot.LT.ilo )
448 \$ GO TO 80
449*
450* ==== Locate active block ====
451*
452 DO 10 k = kbot, ilo + 1, -1
453 IF( h( k, k-1 ).EQ.zero )
454 \$ GO TO 20
455 10 CONTINUE
456 k = ilo
457 20 CONTINUE
458 ktop = k
459*
460* ==== Select deflation window size:
461* . Typical Case:
462* . If possible and advisable, nibble the entire
463* . active block. If not, use size MIN(NWR,NWMAX)
464* . or MIN(NWR+1,NWMAX) depending upon which has
465* . the smaller corresponding subdiagonal entry
466* . (a heuristic).
467* .
468* . Exceptional Case:
469* . If there have been no deflations in KEXNW or
470* . more iterations, then vary the deflation window
471* . size. At first, because, larger windows are,
472* . in general, more powerful than smaller ones,
473* . rapidly increase the window to the maximum possible.
474* . Then, gradually reduce the window size. ====
475*
476 nh = kbot - ktop + 1
477 nwupbd = min( nh, nwmax )
478 IF( ndfl.LT.kexnw ) THEN
479 nw = min( nwupbd, nwr )
480 ELSE
481 nw = min( nwupbd, 2*nw )
482 END IF
483 IF( nw.LT.nwmax ) THEN
484 IF( nw.GE.nh-1 ) THEN
485 nw = nh
486 ELSE
487 kwtop = kbot - nw + 1
488 IF( cabs1( h( kwtop, kwtop-1 ) ).GT.
489 \$ cabs1( h( kwtop-1, kwtop-2 ) ) )nw = nw + 1
490 END IF
491 END IF
492 IF( ndfl.LT.kexnw ) THEN
493 ndec = -1
494 ELSE IF( ndec.GE.0 .OR. nw.GE.nwupbd ) THEN
495 ndec = ndec + 1
496 IF( nw-ndec.LT.2 )
497 \$ ndec = 0
498 nw = nw - ndec
499 END IF
500*
501* ==== Aggressive early deflation:
502* . split workspace under the subdiagonal into
503* . - an nw-by-nw work array V in the lower
504* . left-hand-corner,
505* . - an NW-by-at-least-NW-but-more-is-better
506* . (NW-by-NHO) horizontal work array along
507* . the bottom edge,
508* . - an at-least-NW-but-more-is-better (NHV-by-NW)
509* . vertical work array along the left-hand-edge.
510* . ====
511*
512 kv = n - nw + 1
513 kt = nw + 1
514 nho = ( n-nw-1 ) - kt + 1
515 kwv = nw + 2
516 nve = ( n-nw ) - kwv + 1
517*
518* ==== Aggressive early deflation ====
519*
520 CALL zlaqr2( wantt, wantz, n, ktop, kbot, nw, h, ldh, iloz,
521 \$ ihiz, z, ldz, ls, ld, w, h( kv, 1 ), ldh, nho,
522 \$ h( kv, kt ), ldh, nve, h( kwv, 1 ), ldh, work,
523 \$ lwork )
524*
525* ==== Adjust KBOT accounting for new deflations. ====
526*
527 kbot = kbot - ld
528*
529* ==== KS points to the shifts. ====
530*
531 ks = kbot - ls + 1
532*
533* ==== Skip an expensive QR sweep if there is a (partly
534* . heuristic) reason to expect that many eigenvalues
535* . will deflate without it. Here, the QR sweep is
536* . skipped if many eigenvalues have just been deflated
537* . or if the remaining active block is small.
538*
539 IF( ( ld.EQ.0 ) .OR. ( ( 100*ld.LE.nw*nibble ) .AND. ( kbot-
540 \$ ktop+1.GT.min( nmin, nwmax ) ) ) ) THEN
541*
542* ==== NS = nominal number of simultaneous shifts.
543* . This may be lowered (slightly) if ZLAQR2
544* . did not provide that many shifts. ====
545*
546 ns = min( nsmax, nsr, max( 2, kbot-ktop ) )
547 ns = ns - mod( ns, 2 )
548*
549* ==== If there have been no deflations
550* . in a multiple of KEXSH iterations,
551* . then try exceptional shifts.
552* . Otherwise use shifts provided by
553* . ZLAQR2 above or from the eigenvalues
554* . of a trailing principal submatrix. ====
555*
556 IF( mod( ndfl, kexsh ).EQ.0 ) THEN
557 ks = kbot - ns + 1
558 DO 30 i = kbot, ks + 1, -2
559 w( i ) = h( i, i ) + wilk1*cabs1( h( i, i-1 ) )
560 w( i-1 ) = w( i )
561 30 CONTINUE
562 ELSE
563*
564* ==== Got NS/2 or fewer shifts? Use ZLAHQR
565* . on a trailing principal submatrix to
566* . get more. (Since NS.LE.NSMAX.LE.(N-3)/6,
567* . there is enough space below the subdiagonal
568* . to fit an NS-by-NS scratch array.) ====
569*
570 IF( kbot-ks+1.LE.ns / 2 ) THEN
571 ks = kbot - ns + 1
572 kt = n - ns + 1
573 CALL zlacpy( 'A', ns, ns, h( ks, ks ), ldh,
574 \$ h( kt, 1 ), ldh )
575 CALL zlahqr( .false., .false., ns, 1, ns,
576 \$ h( kt, 1 ), ldh, w( ks ), 1, 1, zdum,
577 \$ 1, inf )
578 ks = ks + inf
579*
580* ==== In case of a rare QR failure use
581* . eigenvalues of the trailing 2-by-2
582* . principal submatrix. Scale to avoid
583* . overflows, underflows and subnormals.
584* . (The scale factor S can not be zero,
585* . because H(KBOT,KBOT-1) is nonzero.) ====
586*
587 IF( ks.GE.kbot ) THEN
588 s = cabs1( h( kbot-1, kbot-1 ) ) +
589 \$ cabs1( h( kbot, kbot-1 ) ) +
590 \$ cabs1( h( kbot-1, kbot ) ) +
591 \$ cabs1( h( kbot, kbot ) )
592 aa = h( kbot-1, kbot-1 ) / s
593 cc = h( kbot, kbot-1 ) / s
594 bb = h( kbot-1, kbot ) / s
595 dd = h( kbot, kbot ) / s
596 tr2 = ( aa+dd ) / two
597 det = ( aa-tr2 )*( dd-tr2 ) - bb*cc
598 rtdisc = sqrt( -det )
599 w( kbot-1 ) = ( tr2+rtdisc )*s
600 w( kbot ) = ( tr2-rtdisc )*s
601*
602 ks = kbot - 1
603 END IF
604 END IF
605*
606 IF( kbot-ks+1.GT.ns ) THEN
607*
608* ==== Sort the shifts (Helps a little) ====
609*
610 sorted = .false.
611 DO 50 k = kbot, ks + 1, -1
612 IF( sorted )
613 \$ GO TO 60
614 sorted = .true.
615 DO 40 i = ks, k - 1
616 IF( cabs1( w( i ) ).LT.cabs1( w( i+1 ) ) )
617 \$ THEN
618 sorted = .false.
619 swap = w( i )
620 w( i ) = w( i+1 )
621 w( i+1 ) = swap
622 END IF
623 40 CONTINUE
624 50 CONTINUE
625 60 CONTINUE
626 END IF
627 END IF
628*
629* ==== If there are only two shifts, then use
630* . only one. ====
631*
632 IF( kbot-ks+1.EQ.2 ) THEN
633 IF( cabs1( w( kbot )-h( kbot, kbot ) ).LT.
634 \$ cabs1( w( kbot-1 )-h( kbot, kbot ) ) ) THEN
635 w( kbot-1 ) = w( kbot )
636 ELSE
637 w( kbot ) = w( kbot-1 )
638 END IF
639 END IF
640*
641* ==== Use up to NS of the the smallest magnitude
642* . shifts. If there aren't NS shifts available,
643* . then use them all, possibly dropping one to
644* . make the number of shifts even. ====
645*
646 ns = min( ns, kbot-ks+1 )
647 ns = ns - mod( ns, 2 )
648 ks = kbot - ns + 1
649*
650* ==== Small-bulge multi-shift QR sweep:
651* . split workspace under the subdiagonal into
652* . - a KDU-by-KDU work array U in the lower
653* . left-hand-corner,
654* . - a KDU-by-at-least-KDU-but-more-is-better
655* . (KDU-by-NHo) horizontal work array WH along
656* . the bottom edge,
657* . - and an at-least-KDU-but-more-is-better-by-KDU
658* . (NVE-by-KDU) vertical work WV arrow along
659* . the left-hand-edge. ====
660*
661 kdu = 2*ns
662 ku = n - kdu + 1
663 kwh = kdu + 1
664 nho = ( n-kdu+1-4 ) - ( kdu+1 ) + 1
665 kwv = kdu + 4
666 nve = n - kdu - kwv + 1
667*
668* ==== Small-bulge multi-shift QR sweep ====
669*
670 CALL zlaqr5( wantt, wantz, kacc22, n, ktop, kbot, ns,
671 \$ w( ks ), h, ldh, iloz, ihiz, z, ldz, work,
672 \$ 3, h( ku, 1 ), ldh, nve, h( kwv, 1 ), ldh,
673 \$ nho, h( ku, kwh ), ldh )
674 END IF
675*
676* ==== Note progress (or the lack of it). ====
677*
678 IF( ld.GT.0 ) THEN
679 ndfl = 1
680 ELSE
681 ndfl = ndfl + 1
682 END IF
683*
684* ==== End of main loop ====
685 70 CONTINUE
686*
687* ==== Iteration limit exceeded. Set INFO to show where
688* . the problem occurred and exit. ====
689*
690 info = kbot
691 80 CONTINUE
692 END IF
693*
694* ==== Return the optimal value of LWORK. ====
695*
696 work( 1 ) = dcmplx( lwkopt, 0 )
697*
698* ==== End of ZLAQR4 ====
699*
700 END
subroutine zlaqr5(WANTT, WANTZ, KACC22, N, KTOP, KBOT, NSHFTS, S, H, LDH, ILOZ, IHIZ, Z, LDZ, V, LDV, U, LDU, NV, WV, LDWV, NH, WH, LDWH)
ZLAQR5 performs a single small-bulge multi-shift QR sweep.
Definition: zlaqr5.f:257
subroutine zlahqr(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, INFO)
ZLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix,...
Definition: zlahqr.f:195
subroutine zlaqr2(WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SH, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK)
ZLAQR2 performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fu...
Definition: zlaqr2.f:270
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
subroutine zlaqr4(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)
ZLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...
Definition: zlaqr4.f:247