1: *> \brief \b ZLAQR0 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 ZLAQR0 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaqr0.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaqr0.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaqr0.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZLAQR0( 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: *> ZLAQR0 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: *>
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 ZGEBAL, and then passed to ZGEHRD when the
85: *> matrix output by ZGEBAL 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 COMPLEX*16 array, dimension (LDH,N)
94: *> On entry, the upper Hessenberg matrix H.
95: *> On exit, if INFO = 0 and WANTT is .TRUE., then H
96: *> contains the upper triangular matrix T from the Schur
97: *> decomposition (the Schur form). If INFO = 0 and WANT is
98: *> .FALSE., then the contents of H are unspecified on exit.
99: *> (The output value of H when INFO > 0 is given under the
100: *> description of INFO below.)
101: *>
102: *> This subroutine may explicitly set H(i,j) = 0 for i > j and
103: *> j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N.
104: *> \endverbatim
105: *>
106: *> \param[in] LDH
107: *> \verbatim
108: *> LDH is INTEGER
109: *> The leading dimension of the array H. LDH >= max(1,N).
110: *> \endverbatim
111: *>
112: *> \param[out] W
113: *> \verbatim
114: *> W is COMPLEX*16 array, dimension (N)
115: *> The computed eigenvalues of H(ILO:IHI,ILO:IHI) are stored
116: *> in W(ILO:IHI). If WANTT is .TRUE., then the eigenvalues are
117: *> stored in the same order as on the diagonal of the Schur
118: *> form returned in H, with W(i) = H(i,i).
119: *> \endverbatim
120: *>
121: *> \param[in] ILOZ
122: *> \verbatim
123: *> ILOZ is INTEGER
124: *> \endverbatim
125: *>
126: *> \param[in] IHIZ
127: *> \verbatim
128: *> IHIZ is INTEGER
129: *> Specify the rows of Z to which transformations must be
130: *> applied if WANTZ is .TRUE..
131: *> 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
132: *> \endverbatim
133: *>
134: *> \param[in,out] Z
135: *> \verbatim
136: *> Z is COMPLEX*16 array, dimension (LDZ,IHI)
137: *> If WANTZ is .FALSE., then Z is not referenced.
138: *> If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is
139: *> replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the
140: *> orthogonal Schur factor of H(ILO:IHI,ILO:IHI).
141: *> (The output value of Z when INFO > 0 is given under
142: *> the description of INFO below.)
143: *> \endverbatim
144: *>
145: *> \param[in] LDZ
146: *> \verbatim
147: *> LDZ is INTEGER
148: *> The leading dimension of the array Z. if WANTZ is .TRUE.
149: *> then LDZ >= MAX(1,IHIZ). Otherwise, LDZ >= 1.
150: *> \endverbatim
151: *>
152: *> \param[out] WORK
153: *> \verbatim
154: *> WORK is COMPLEX*16 array, dimension LWORK
155: *> On exit, if LWORK = -1, WORK(1) returns an estimate of
156: *> the optimal value for LWORK.
157: *> \endverbatim
158: *>
159: *> \param[in] LWORK
160: *> \verbatim
161: *> LWORK is INTEGER
162: *> The dimension of the array WORK. LWORK >= max(1,N)
163: *> is sufficient, but LWORK typically as large as 6*N may
164: *> be required for optimal performance. A workspace query
165: *> to determine the optimal workspace size is recommended.
166: *>
167: *> If LWORK = -1, then ZLAQR0 does a workspace query.
168: *> In this case, ZLAQR0 checks the input parameters and
169: *> estimates the optimal workspace size for the given
170: *> values of N, ILO and IHI. The estimate is returned
171: *> in WORK(1). No error message related to LWORK is
172: *> issued by XERBLA. Neither H nor Z are accessed.
173: *> \endverbatim
174: *>
175: *> \param[out] INFO
176: *> \verbatim
177: *> INFO is INTEGER
178: *> = 0: successful exit
179: *> > 0: if INFO = i, ZLAQR0 failed to compute all of
180: *> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
181: *> and WI contain those eigenvalues which have been
182: *> successfully computed. (Failures are rare.)
183: *>
184: *> If INFO > 0 and WANT is .FALSE., then on exit,
185: *> the remaining unconverged eigenvalues are the eigen-
186: *> values of the upper Hessenberg matrix rows and
187: *> columns ILO through INFO of the final, output
188: *> value of H.
189: *>
190: *> If INFO > 0 and WANTT is .TRUE., then on exit
191: *>
192: *> (*) (initial value of H)*U = U*(final value of H)
193: *>
194: *> where U is a unitary matrix. The final
195: *> value of H is upper Hessenberg and triangular in
196: *> rows and columns INFO+1 through IHI.
197: *>
198: *> If INFO > 0 and WANTZ is .TRUE., then on exit
199: *>
200: *> (final value of Z(ILO:IHI,ILOZ:IHIZ)
201: *> = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U
202: *>
203: *> where U is the unitary matrix in (*) (regard-
204: *> less of the value of WANTT.)
205: *>
206: *> If INFO > 0 and WANTZ is .FALSE., then Z is not
207: *> accessed.
208: *> \endverbatim
209: *
210: * Authors:
211: * ========
212: *
213: *> \author Univ. of Tennessee
214: *> \author Univ. of California Berkeley
215: *> \author Univ. of Colorado Denver
216: *> \author NAG Ltd.
217: *
218: *> \ingroup complex16OTHERauxiliary
219: *
220: *> \par Contributors:
221: * ==================
222: *>
223: *> Karen Braman and Ralph Byers, Department of Mathematics,
224: *> University of Kansas, USA
225: *
226: *> \par References:
227: * ================
228: *>
229: *> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
230: *> Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
231: *> Performance, SIAM Journal of Matrix Analysis, volume 23, pages
232: *> 929--947, 2002.
233: *> \n
234: *> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
235: *> Algorithm Part II: Aggressive Early Deflation, SIAM Journal
236: *> of Matrix Analysis, volume 23, pages 948--973, 2002.
237: *>
238: * =====================================================================
239: SUBROUTINE ZLAQR0( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
240: $ IHIZ, Z, LDZ, WORK, LWORK, INFO )
241: *
242: * -- LAPACK auxiliary routine --
243: * -- LAPACK is a software package provided by Univ. of Tennessee, --
244: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
245: *
246: * .. Scalar Arguments ..
247: INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N
248: LOGICAL WANTT, WANTZ
249: * ..
250: * .. Array Arguments ..
251: COMPLEX*16 H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * )
252: * ..
253: *
254: * ================================================================
255: *
256: * .. Parameters ..
257: *
258: * ==== Matrices of order NTINY or smaller must be processed by
259: * . ZLAHQR because of insufficient subdiagonal scratch space.
260: * . (This is a hard limit.) ====
261: INTEGER NTINY
262: PARAMETER ( NTINY = 15 )
263: *
264: * ==== Exceptional deflation windows: try to cure rare
265: * . slow convergence by varying the size of the
266: * . deflation window after KEXNW iterations. ====
267: INTEGER KEXNW
268: PARAMETER ( KEXNW = 5 )
269: *
270: * ==== Exceptional shifts: try to cure rare slow convergence
271: * . with ad-hoc exceptional shifts every KEXSH iterations.
272: * . ====
273: INTEGER KEXSH
274: PARAMETER ( KEXSH = 6 )
275: *
276: * ==== The constant WILK1 is used to form the exceptional
277: * . shifts. ====
278: DOUBLE PRECISION WILK1
279: PARAMETER ( WILK1 = 0.75d0 )
280: COMPLEX*16 ZERO, ONE
281: PARAMETER ( ZERO = ( 0.0d0, 0.0d0 ),
282: $ ONE = ( 1.0d0, 0.0d0 ) )
283: DOUBLE PRECISION TWO
284: PARAMETER ( TWO = 2.0d0 )
285: * ..
286: * .. Local Scalars ..
287: COMPLEX*16 AA, BB, CC, CDUM, DD, DET, RTDISC, SWAP, TR2
288: DOUBLE PRECISION S
289: INTEGER I, INF, IT, ITMAX, K, KACC22, KBOT, KDU, KS,
290: $ KT, KTOP, KU, KV, KWH, KWTOP, KWV, LD, LS,
291: $ LWKOPT, NDEC, NDFL, NH, NHO, NIBBLE, NMIN, NS,
292: $ NSMAX, NSR, NVE, NW, NWMAX, NWR, NWUPBD
293: LOGICAL SORTED
294: CHARACTER JBCMPZ*2
295: * ..
296: * .. External Functions ..
297: INTEGER ILAENV
298: EXTERNAL ILAENV
299: * ..
300: * .. Local Arrays ..
301: COMPLEX*16 ZDUM( 1, 1 )
302: * ..
303: * .. External Subroutines ..
304: EXTERNAL ZLACPY, ZLAHQR, ZLAQR3, ZLAQR4, ZLAQR5
305: * ..
306: * .. Intrinsic Functions ..
307: INTRINSIC ABS, DBLE, DCMPLX, DIMAG, INT, MAX, MIN, MOD,
308: $ SQRT
309: * ..
310: * .. Statement Functions ..
311: DOUBLE PRECISION CABS1
312: * ..
313: * .. Statement Function definitions ..
314: CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
315: * ..
316: * .. Executable Statements ..
317: INFO = 0
318: *
319: * ==== Quick return for N = 0: nothing to do. ====
320: *
321: IF( N.EQ.0 ) THEN
322: WORK( 1 ) = ONE
323: RETURN
324: END IF
325: *
326: IF( N.LE.NTINY ) THEN
327: *
328: * ==== Tiny matrices must use ZLAHQR. ====
329: *
330: LWKOPT = 1
331: IF( LWORK.NE.-1 )
332: $ CALL ZLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ,
333: $ IHIZ, Z, LDZ, INFO )
334: ELSE
335: *
336: * ==== Use small bulge multi-shift QR with aggressive early
337: * . deflation on larger-than-tiny matrices. ====
338: *
339: * ==== Hope for the best. ====
340: *
341: INFO = 0
342: *
343: * ==== Set up job flags for ILAENV. ====
344: *
345: IF( WANTT ) THEN
346: JBCMPZ( 1: 1 ) = 'S'
347: ELSE
348: JBCMPZ( 1: 1 ) = 'E'
349: END IF
350: IF( WANTZ ) THEN
351: JBCMPZ( 2: 2 ) = 'V'
352: ELSE
353: JBCMPZ( 2: 2 ) = 'N'
354: END IF
355: *
356: * ==== NWR = recommended deflation window size. At this
357: * . point, N .GT. NTINY = 15, so there is enough
358: * . subdiagonal workspace for NWR.GE.2 as required.
359: * . (In fact, there is enough subdiagonal space for
360: * . NWR.GE.4.) ====
361: *
362: NWR = ILAENV( 13, 'ZLAQR0', JBCMPZ, N, ILO, IHI, LWORK )
363: NWR = MAX( 2, NWR )
364: NWR = MIN( IHI-ILO+1, ( N-1 ) / 3, NWR )
365: *
366: * ==== NSR = recommended number of simultaneous shifts.
367: * . At this point N .GT. NTINY = 15, so there is at
368: * . enough subdiagonal workspace for NSR to be even
369: * . and greater than or equal to two as required. ====
370: *
371: NSR = ILAENV( 15, 'ZLAQR0', JBCMPZ, N, ILO, IHI, LWORK )
372: NSR = MIN( NSR, ( N-3 ) / 6, IHI-ILO )
373: NSR = MAX( 2, NSR-MOD( NSR, 2 ) )
374: *
375: * ==== Estimate optimal workspace ====
376: *
377: * ==== Workspace query call to ZLAQR3 ====
378: *
379: CALL ZLAQR3( WANTT, WANTZ, N, ILO, IHI, NWR+1, H, LDH, ILOZ,
380: $ IHIZ, Z, LDZ, LS, LD, W, H, LDH, N, H, LDH, N, H,
381: $ LDH, WORK, -1 )
382: *
383: * ==== Optimal workspace = MAX(ZLAQR5, ZLAQR3) ====
384: *
385: LWKOPT = MAX( 3*NSR / 2, INT( WORK( 1 ) ) )
386: *
387: * ==== Quick return in case of workspace query. ====
388: *
389: IF( LWORK.EQ.-1 ) THEN
390: WORK( 1 ) = DCMPLX( LWKOPT, 0 )
391: RETURN
392: END IF
393: *
394: * ==== ZLAHQR/ZLAQR0 crossover point ====
395: *
396: NMIN = ILAENV( 12, 'ZLAQR0', JBCMPZ, N, ILO, IHI, LWORK )
397: NMIN = MAX( NTINY, NMIN )
398: *
399: * ==== Nibble crossover point ====
400: *
401: NIBBLE = ILAENV( 14, 'ZLAQR0', JBCMPZ, N, ILO, IHI, LWORK )
402: NIBBLE = MAX( 0, NIBBLE )
403: *
404: * ==== Accumulate reflections during ttswp? Use block
405: * . 2-by-2 structure during matrix-matrix multiply? ====
406: *
407: KACC22 = ILAENV( 16, 'ZLAQR0', JBCMPZ, N, ILO, IHI, LWORK )
408: KACC22 = MAX( 0, KACC22 )
409: KACC22 = MIN( 2, KACC22 )
410: *
411: * ==== NWMAX = the largest possible deflation window for
412: * . which there is sufficient workspace. ====
413: *
414: NWMAX = MIN( ( N-1 ) / 3, LWORK / 2 )
415: NW = NWMAX
416: *
417: * ==== NSMAX = the Largest number of simultaneous shifts
418: * . for which there is sufficient workspace. ====
419: *
420: NSMAX = MIN( ( N-3 ) / 6, 2*LWORK / 3 )
421: NSMAX = NSMAX - MOD( NSMAX, 2 )
422: *
423: * ==== NDFL: an iteration count restarted at deflation. ====
424: *
425: NDFL = 1
426: *
427: * ==== ITMAX = iteration limit ====
428: *
429: ITMAX = MAX( 30, 2*KEXSH )*MAX( 10, ( IHI-ILO+1 ) )
430: *
431: * ==== Last row and column in the active block ====
432: *
433: KBOT = IHI
434: *
435: * ==== Main Loop ====
436: *
437: DO 70 IT = 1, ITMAX
438: *
439: * ==== Done when KBOT falls below ILO ====
440: *
441: IF( KBOT.LT.ILO )
442: $ GO TO 80
443: *
444: * ==== Locate active block ====
445: *
446: DO 10 K = KBOT, ILO + 1, -1
447: IF( H( K, K-1 ).EQ.ZERO )
448: $ GO TO 20
449: 10 CONTINUE
450: K = ILO
451: 20 CONTINUE
452: KTOP = K
453: *
454: * ==== Select deflation window size:
455: * . Typical Case:
456: * . If possible and advisable, nibble the entire
457: * . active block. If not, use size MIN(NWR,NWMAX)
458: * . or MIN(NWR+1,NWMAX) depending upon which has
459: * . the smaller corresponding subdiagonal entry
460: * . (a heuristic).
461: * .
462: * . Exceptional Case:
463: * . If there have been no deflations in KEXNW or
464: * . more iterations, then vary the deflation window
465: * . size. At first, because, larger windows are,
466: * . in general, more powerful than smaller ones,
467: * . rapidly increase the window to the maximum possible.
468: * . Then, gradually reduce the window size. ====
469: *
470: NH = KBOT - KTOP + 1
471: NWUPBD = MIN( NH, NWMAX )
472: IF( NDFL.LT.KEXNW ) THEN
473: NW = MIN( NWUPBD, NWR )
474: ELSE
475: NW = MIN( NWUPBD, 2*NW )
476: END IF
477: IF( NW.LT.NWMAX ) THEN
478: IF( NW.GE.NH-1 ) THEN
479: NW = NH
480: ELSE
481: KWTOP = KBOT - NW + 1
482: IF( CABS1( H( KWTOP, KWTOP-1 ) ).GT.
483: $ CABS1( H( KWTOP-1, KWTOP-2 ) ) )NW = NW + 1
484: END IF
485: END IF
486: IF( NDFL.LT.KEXNW ) THEN
487: NDEC = -1
488: ELSE IF( NDEC.GE.0 .OR. NW.GE.NWUPBD ) THEN
489: NDEC = NDEC + 1
490: IF( NW-NDEC.LT.2 )
491: $ NDEC = 0
492: NW = NW - NDEC
493: END IF
494: *
495: * ==== Aggressive early deflation:
496: * . split workspace under the subdiagonal into
497: * . - an nw-by-nw work array V in the lower
498: * . left-hand-corner,
499: * . - an NW-by-at-least-NW-but-more-is-better
500: * . (NW-by-NHO) horizontal work array along
501: * . the bottom edge,
502: * . - an at-least-NW-but-more-is-better (NHV-by-NW)
503: * . vertical work array along the left-hand-edge.
504: * . ====
505: *
506: KV = N - NW + 1
507: KT = NW + 1
508: NHO = ( N-NW-1 ) - KT + 1
509: KWV = NW + 2
510: NVE = ( N-NW ) - KWV + 1
511: *
512: * ==== Aggressive early deflation ====
513: *
514: CALL ZLAQR3( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
515: $ IHIZ, Z, LDZ, LS, LD, W, H( KV, 1 ), LDH, NHO,
516: $ H( KV, KT ), LDH, NVE, H( KWV, 1 ), LDH, WORK,
517: $ LWORK )
518: *
519: * ==== Adjust KBOT accounting for new deflations. ====
520: *
521: KBOT = KBOT - LD
522: *
523: * ==== KS points to the shifts. ====
524: *
525: KS = KBOT - LS + 1
526: *
527: * ==== Skip an expensive QR sweep if there is a (partly
528: * . heuristic) reason to expect that many eigenvalues
529: * . will deflate without it. Here, the QR sweep is
530: * . skipped if many eigenvalues have just been deflated
531: * . or if the remaining active block is small.
532: *
533: IF( ( LD.EQ.0 ) .OR. ( ( 100*LD.LE.NW*NIBBLE ) .AND. ( KBOT-
534: $ KTOP+1.GT.MIN( NMIN, NWMAX ) ) ) ) THEN
535: *
536: * ==== NS = nominal number of simultaneous shifts.
537: * . This may be lowered (slightly) if ZLAQR3
538: * . did not provide that many shifts. ====
539: *
540: NS = MIN( NSMAX, NSR, MAX( 2, KBOT-KTOP ) )
541: NS = NS - MOD( NS, 2 )
542: *
543: * ==== If there have been no deflations
544: * . in a multiple of KEXSH iterations,
545: * . then try exceptional shifts.
546: * . Otherwise use shifts provided by
547: * . ZLAQR3 above or from the eigenvalues
548: * . of a trailing principal submatrix. ====
549: *
550: IF( MOD( NDFL, KEXSH ).EQ.0 ) THEN
551: KS = KBOT - NS + 1
552: DO 30 I = KBOT, KS + 1, -2
553: W( I ) = H( I, I ) + WILK1*CABS1( H( I, I-1 ) )
554: W( I-1 ) = W( I )
555: 30 CONTINUE
556: ELSE
557: *
558: * ==== Got NS/2 or fewer shifts? Use ZLAQR4 or
559: * . ZLAHQR on a trailing principal submatrix to
560: * . get more. (Since NS.LE.NSMAX.LE.(N-3)/6,
561: * . there is enough space below the subdiagonal
562: * . to fit an NS-by-NS scratch array.) ====
563: *
564: IF( KBOT-KS+1.LE.NS / 2 ) THEN
565: KS = KBOT - NS + 1
566: KT = N - NS + 1
567: CALL ZLACPY( 'A', NS, NS, H( KS, KS ), LDH,
568: $ H( KT, 1 ), LDH )
569: IF( NS.GT.NMIN ) THEN
570: CALL ZLAQR4( .false., .false., NS, 1, NS,
571: $ H( KT, 1 ), LDH, W( KS ), 1, 1,
572: $ ZDUM, 1, WORK, LWORK, INF )
573: ELSE
574: CALL ZLAHQR( .false., .false., NS, 1, NS,
575: $ H( KT, 1 ), LDH, W( KS ), 1, 1,
576: $ ZDUM, 1, INF )
577: END IF
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 ZLAQR0 ====
699: *
700: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlaqr0.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack