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