Annotation of rpl/lapack/lapack/dlahqr.f, revision 1.20
1.11 bertrand 1: *> \brief \b DLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.
1.8 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.16 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.8 bertrand 7: *
8: *> \htmlonly
1.16 bertrand 9: *> Download DLAHQR + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlahqr.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlahqr.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlahqr.f">
1.8 bertrand 15: *> [TXT]</a>
1.16 bertrand 16: *> \endhtmlonly
1.8 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
22: * ILOZ, IHIZ, Z, LDZ, INFO )
1.16 bertrand 23: *
1.8 bertrand 24: * .. Scalar Arguments ..
25: * INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
26: * LOGICAL WANTT, WANTZ
27: * ..
28: * .. Array Arguments ..
29: * DOUBLE PRECISION H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * )
30: * ..
1.16 bertrand 31: *
1.8 bertrand 32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> DLAHQR is an auxiliary routine called by DHSEQR to update the
39: *> eigenvalues and Schur decomposition already computed by DHSEQR, by
40: *> dealing with the Hessenberg submatrix in rows and columns ILO to
41: *> IHI.
42: *> \endverbatim
43: *
44: * Arguments:
45: * ==========
46: *
47: *> \param[in] WANTT
48: *> \verbatim
49: *> WANTT is LOGICAL
50: *> = .TRUE. : the full Schur form T is required;
51: *> = .FALSE.: only eigenvalues are required.
52: *> \endverbatim
53: *>
54: *> \param[in] WANTZ
55: *> \verbatim
56: *> WANTZ is LOGICAL
57: *> = .TRUE. : the matrix of Schur vectors Z is required;
58: *> = .FALSE.: Schur vectors are not required.
59: *> \endverbatim
60: *>
61: *> \param[in] N
62: *> \verbatim
63: *> N is INTEGER
64: *> The order of the matrix H. N >= 0.
65: *> \endverbatim
66: *>
67: *> \param[in] ILO
68: *> \verbatim
69: *> ILO is INTEGER
70: *> \endverbatim
71: *>
72: *> \param[in] IHI
73: *> \verbatim
74: *> IHI is INTEGER
75: *> It is assumed that H is already upper quasi-triangular in
76: *> rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless
77: *> ILO = 1). DLAHQR works primarily with the Hessenberg
78: *> submatrix in rows and columns ILO to IHI, but applies
79: *> transformations to all of H if WANTT is .TRUE..
80: *> 1 <= ILO <= max(1,IHI); IHI <= N.
81: *> \endverbatim
82: *>
83: *> \param[in,out] H
84: *> \verbatim
85: *> H is DOUBLE PRECISION array, dimension (LDH,N)
86: *> On entry, the upper Hessenberg matrix H.
87: *> On exit, if INFO is zero and if WANTT is .TRUE., H is upper
88: *> quasi-triangular in rows and columns ILO:IHI, with any
89: *> 2-by-2 diagonal blocks in standard form. If INFO is zero
90: *> and WANTT is .FALSE., the contents of H are unspecified on
91: *> exit. The output state of H if INFO is nonzero is given
92: *> below under the description of INFO.
93: *> \endverbatim
94: *>
95: *> \param[in] LDH
96: *> \verbatim
97: *> LDH is INTEGER
98: *> The leading dimension of the array H. LDH >= max(1,N).
99: *> \endverbatim
100: *>
101: *> \param[out] WR
102: *> \verbatim
103: *> WR is DOUBLE PRECISION array, dimension (N)
104: *> \endverbatim
105: *>
106: *> \param[out] WI
107: *> \verbatim
108: *> WI is DOUBLE PRECISION array, dimension (N)
109: *> The real and imaginary parts, respectively, of the computed
110: *> eigenvalues ILO to IHI are stored in the corresponding
111: *> elements of WR and WI. If two eigenvalues are computed as a
112: *> complex conjugate pair, they are stored in consecutive
113: *> elements of WR and WI, say the i-th and (i+1)th, with
114: *> WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the
115: *> eigenvalues are stored in the same order as on the diagonal
116: *> of the Schur form returned in H, with WR(i) = H(i,i), and, if
117: *> H(i:i+1,i:i+1) is a 2-by-2 diagonal block,
118: *> WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(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 DOUBLE PRECISION array, dimension (LDZ,N)
137: *> If WANTZ is .TRUE., on entry Z must contain the current
138: *> matrix Z of transformations accumulated by DHSEQR, and on
139: *> exit Z has been updated; transformations are applied only to
140: *> the submatrix Z(ILOZ:IHIZ,ILO:IHI).
141: *> If WANTZ is .FALSE., Z is not referenced.
142: *> \endverbatim
143: *>
144: *> \param[in] LDZ
145: *> \verbatim
146: *> LDZ is INTEGER
147: *> The leading dimension of the array Z. LDZ >= max(1,N).
148: *> \endverbatim
149: *>
150: *> \param[out] INFO
151: *> \verbatim
152: *> INFO is INTEGER
1.19 bertrand 153: *> = 0: successful exit
154: *> > 0: If INFO = i, DLAHQR failed to compute all the
1.8 bertrand 155: *> eigenvalues ILO to IHI in a total of 30 iterations
156: *> per eigenvalue; elements i+1:ihi of WR and WI
157: *> contain those eigenvalues which have been
158: *> successfully computed.
159: *>
1.19 bertrand 160: *> If INFO > 0 and WANTT is .FALSE., then on exit,
1.8 bertrand 161: *> the remaining unconverged eigenvalues are the
162: *> eigenvalues of the upper Hessenberg matrix rows
1.19 bertrand 163: *> and columns ILO through INFO of the final, output
1.8 bertrand 164: *> value of H.
165: *>
1.19 bertrand 166: *> If INFO > 0 and WANTT is .TRUE., then on exit
1.8 bertrand 167: *> (*) (initial value of H)*U = U*(final value of H)
1.19 bertrand 168: *> where U is an orthogonal matrix. The final
1.8 bertrand 169: *> value of H is upper Hessenberg and triangular in
170: *> rows and columns INFO+1 through IHI.
171: *>
1.19 bertrand 172: *> If INFO > 0 and WANTZ is .TRUE., then on exit
1.8 bertrand 173: *> (final value of Z) = (initial value of Z)*U
174: *> where U is the orthogonal matrix in (*)
175: *> (regardless of the value of WANTT.)
176: *> \endverbatim
177: *
178: * Authors:
179: * ========
180: *
1.16 bertrand 181: *> \author Univ. of Tennessee
182: *> \author Univ. of California Berkeley
183: *> \author Univ. of Colorado Denver
184: *> \author NAG Ltd.
1.8 bertrand 185: *
186: *> \ingroup doubleOTHERauxiliary
187: *
188: *> \par Further Details:
189: * =====================
190: *>
191: *> \verbatim
192: *>
193: *> 02-96 Based on modifications by
194: *> David Day, Sandia National Laboratory, USA
195: *>
196: *> 12-04 Further modifications by
197: *> Ralph Byers, University of Kansas, USA
198: *> This is a modified version of DLAHQR from LAPACK version 3.0.
199: *> It is (1) more robust against overflow and underflow and
200: *> (2) adopts the more conservative Ahues & Tisseur stopping
201: *> criterion (LAWN 122, 1997).
202: *> \endverbatim
203: *>
204: * =====================================================================
1.1 bertrand 205: SUBROUTINE DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI,
206: $ ILOZ, IHIZ, Z, LDZ, INFO )
1.20 ! bertrand 207: IMPLICIT NONE
1.1 bertrand 208: *
1.20 ! bertrand 209: * -- LAPACK auxiliary routine --
1.8 bertrand 210: * -- LAPACK is a software package provided by Univ. of Tennessee, --
211: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.1 bertrand 212: *
213: * .. Scalar Arguments ..
214: INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, N
215: LOGICAL WANTT, WANTZ
216: * ..
217: * .. Array Arguments ..
218: DOUBLE PRECISION H( LDH, * ), WI( * ), WR( * ), Z( LDZ, * )
219: * ..
220: *
1.8 bertrand 221: * =========================================================
1.1 bertrand 222: *
223: * .. Parameters ..
224: DOUBLE PRECISION ZERO, ONE, TWO
225: PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0, TWO = 2.0d0 )
226: DOUBLE PRECISION DAT1, DAT2
227: PARAMETER ( DAT1 = 3.0d0 / 4.0d0, DAT2 = -0.4375d0 )
1.20 ! bertrand 228: INTEGER KEXSH
! 229: PARAMETER ( KEXSH = 10 )
1.1 bertrand 230: * ..
231: * .. Local Scalars ..
232: DOUBLE PRECISION AA, AB, BA, BB, CS, DET, H11, H12, H21, H21S,
233: $ H22, RT1I, RT1R, RT2I, RT2R, RTDISC, S, SAFMAX,
234: $ SAFMIN, SMLNUM, SN, SUM, T1, T2, T3, TR, TST,
235: $ ULP, V2, V3
1.20 ! bertrand 236: INTEGER I, I1, I2, ITS, ITMAX, J, K, L, M, NH, NR, NZ,
! 237: $ KDEFL
1.1 bertrand 238: * ..
239: * .. Local Arrays ..
240: DOUBLE PRECISION V( 3 )
241: * ..
242: * .. External Functions ..
243: DOUBLE PRECISION DLAMCH
244: EXTERNAL DLAMCH
245: * ..
246: * .. External Subroutines ..
247: EXTERNAL DCOPY, DLABAD, DLANV2, DLARFG, DROT
248: * ..
249: * .. Intrinsic Functions ..
250: INTRINSIC ABS, DBLE, MAX, MIN, SQRT
251: * ..
252: * .. Executable Statements ..
253: *
254: INFO = 0
255: *
256: * Quick return if possible
257: *
258: IF( N.EQ.0 )
259: $ RETURN
260: IF( ILO.EQ.IHI ) THEN
261: WR( ILO ) = H( ILO, ILO )
262: WI( ILO ) = ZERO
263: RETURN
264: END IF
265: *
266: * ==== clear out the trash ====
267: DO 10 J = ILO, IHI - 3
268: H( J+2, J ) = ZERO
269: H( J+3, J ) = ZERO
270: 10 CONTINUE
271: IF( ILO.LE.IHI-2 )
272: $ H( IHI, IHI-2 ) = ZERO
273: *
274: NH = IHI - ILO + 1
275: NZ = IHIZ - ILOZ + 1
276: *
277: * Set machine-dependent constants for the stopping criterion.
278: *
279: SAFMIN = DLAMCH( 'SAFE MINIMUM' )
280: SAFMAX = ONE / SAFMIN
281: CALL DLABAD( SAFMIN, SAFMAX )
282: ULP = DLAMCH( 'PRECISION' )
283: SMLNUM = SAFMIN*( DBLE( NH ) / ULP )
284: *
285: * I1 and I2 are the indices of the first row and last column of H
286: * to which transformations must be applied. If eigenvalues only are
287: * being computed, I1 and I2 are set inside the main loop.
288: *
289: IF( WANTT ) THEN
290: I1 = 1
291: I2 = N
292: END IF
293: *
1.14 bertrand 294: * ITMAX is the total number of QR iterations allowed.
295: *
1.16 bertrand 296: ITMAX = 30 * MAX( 10, NH )
1.14 bertrand 297: *
1.20 ! bertrand 298: * KDEFL counts the number of iterations since a deflation
! 299: *
! 300: KDEFL = 0
! 301: *
1.1 bertrand 302: * The main loop begins here. I is the loop index and decreases from
303: * IHI to ILO in steps of 1 or 2. Each iteration of the loop works
304: * with the active submatrix in rows and columns L to I.
305: * Eigenvalues I+1 to IHI have already converged. Either L = ILO or
306: * H(L,L-1) is negligible so that the matrix splits.
307: *
308: I = IHI
309: 20 CONTINUE
310: L = ILO
311: IF( I.LT.ILO )
312: $ GO TO 160
313: *
314: * Perform QR iterations on rows and columns ILO to I until a
315: * submatrix of order 1 or 2 splits off at the bottom because a
316: * subdiagonal element has become negligible.
317: *
318: DO 140 ITS = 0, ITMAX
319: *
320: * Look for a single small subdiagonal element.
321: *
322: DO 30 K = I, L + 1, -1
323: IF( ABS( H( K, K-1 ) ).LE.SMLNUM )
324: $ GO TO 40
325: TST = ABS( H( K-1, K-1 ) ) + ABS( H( K, K ) )
326: IF( TST.EQ.ZERO ) THEN
327: IF( K-2.GE.ILO )
328: $ TST = TST + ABS( H( K-1, K-2 ) )
329: IF( K+1.LE.IHI )
330: $ TST = TST + ABS( H( K+1, K ) )
331: END IF
332: * ==== The following is a conservative small subdiagonal
333: * . deflation criterion due to Ahues & Tisseur (LAWN 122,
334: * . 1997). It has better mathematical foundation and
335: * . improves accuracy in some cases. ====
336: IF( ABS( H( K, K-1 ) ).LE.ULP*TST ) THEN
337: AB = MAX( ABS( H( K, K-1 ) ), ABS( H( K-1, K ) ) )
338: BA = MIN( ABS( H( K, K-1 ) ), ABS( H( K-1, K ) ) )
339: AA = MAX( ABS( H( K, K ) ),
340: $ ABS( H( K-1, K-1 )-H( K, K ) ) )
341: BB = MIN( ABS( H( K, K ) ),
342: $ ABS( H( K-1, K-1 )-H( K, K ) ) )
343: S = AA + AB
344: IF( BA*( AB / S ).LE.MAX( SMLNUM,
345: $ ULP*( BB*( AA / S ) ) ) )GO TO 40
346: END IF
347: 30 CONTINUE
348: 40 CONTINUE
349: L = K
350: IF( L.GT.ILO ) THEN
351: *
352: * H(L,L-1) is negligible
353: *
354: H( L, L-1 ) = ZERO
355: END IF
356: *
357: * Exit from loop if a submatrix of order 1 or 2 has split off.
358: *
359: IF( L.GE.I-1 )
360: $ GO TO 150
1.20 ! bertrand 361: KDEFL = KDEFL + 1
1.1 bertrand 362: *
363: * Now the active submatrix is in rows and columns L to I. If
364: * eigenvalues only are being computed, only the active submatrix
365: * need be transformed.
366: *
367: IF( .NOT.WANTT ) THEN
368: I1 = L
369: I2 = I
370: END IF
371: *
1.20 ! bertrand 372: IF( MOD(KDEFL,2*KEXSH).EQ.0 ) THEN
1.1 bertrand 373: *
374: * Exceptional shift.
375: *
1.20 ! bertrand 376: S = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) )
! 377: H11 = DAT1*S + H( I, I )
1.1 bertrand 378: H12 = DAT2*S
379: H21 = S
380: H22 = H11
1.20 ! bertrand 381: ELSE IF( MOD(KDEFL,KEXSH).EQ.0 ) THEN
1.1 bertrand 382: *
383: * Exceptional shift.
384: *
1.20 ! bertrand 385: S = ABS( H( L+1, L ) ) + ABS( H( L+2, L+1 ) )
! 386: H11 = DAT1*S + H( L, L )
1.1 bertrand 387: H12 = DAT2*S
388: H21 = S
389: H22 = H11
390: ELSE
391: *
392: * Prepare to use Francis' double shift
393: * (i.e. 2nd degree generalized Rayleigh quotient)
394: *
395: H11 = H( I-1, I-1 )
396: H21 = H( I, I-1 )
397: H12 = H( I-1, I )
398: H22 = H( I, I )
399: END IF
400: S = ABS( H11 ) + ABS( H12 ) + ABS( H21 ) + ABS( H22 )
401: IF( S.EQ.ZERO ) THEN
402: RT1R = ZERO
403: RT1I = ZERO
404: RT2R = ZERO
405: RT2I = ZERO
406: ELSE
407: H11 = H11 / S
408: H21 = H21 / S
409: H12 = H12 / S
410: H22 = H22 / S
411: TR = ( H11+H22 ) / TWO
412: DET = ( H11-TR )*( H22-TR ) - H12*H21
413: RTDISC = SQRT( ABS( DET ) )
414: IF( DET.GE.ZERO ) THEN
415: *
416: * ==== complex conjugate shifts ====
417: *
418: RT1R = TR*S
419: RT2R = RT1R
420: RT1I = RTDISC*S
421: RT2I = -RT1I
422: ELSE
423: *
424: * ==== real shifts (use only one of them) ====
425: *
426: RT1R = TR + RTDISC
427: RT2R = TR - RTDISC
428: IF( ABS( RT1R-H22 ).LE.ABS( RT2R-H22 ) ) THEN
429: RT1R = RT1R*S
430: RT2R = RT1R
431: ELSE
432: RT2R = RT2R*S
433: RT1R = RT2R
434: END IF
435: RT1I = ZERO
436: RT2I = ZERO
437: END IF
438: END IF
439: *
440: * Look for two consecutive small subdiagonal elements.
441: *
442: DO 50 M = I - 2, L, -1
443: * Determine the effect of starting the double-shift QR
444: * iteration at row M, and see if this would make H(M,M-1)
445: * negligible. (The following uses scaling to avoid
446: * overflows and most underflows.)
447: *
448: H21S = H( M+1, M )
449: S = ABS( H( M, M )-RT2R ) + ABS( RT2I ) + ABS( H21S )
450: H21S = H( M+1, M ) / S
451: V( 1 ) = H21S*H( M, M+1 ) + ( H( M, M )-RT1R )*
452: $ ( ( H( M, M )-RT2R ) / S ) - RT1I*( RT2I / S )
453: V( 2 ) = H21S*( H( M, M )+H( M+1, M+1 )-RT1R-RT2R )
454: V( 3 ) = H21S*H( M+2, M+1 )
455: S = ABS( V( 1 ) ) + ABS( V( 2 ) ) + ABS( V( 3 ) )
456: V( 1 ) = V( 1 ) / S
457: V( 2 ) = V( 2 ) / S
458: V( 3 ) = V( 3 ) / S
459: IF( M.EQ.L )
460: $ GO TO 60
461: IF( ABS( H( M, M-1 ) )*( ABS( V( 2 ) )+ABS( V( 3 ) ) ).LE.
462: $ ULP*ABS( V( 1 ) )*( ABS( H( M-1, M-1 ) )+ABS( H( M,
463: $ M ) )+ABS( H( M+1, M+1 ) ) ) )GO TO 60
464: 50 CONTINUE
465: 60 CONTINUE
466: *
467: * Double-shift QR step
468: *
469: DO 130 K = M, I - 1
470: *
471: * The first iteration of this loop determines a reflection G
472: * from the vector V and applies it from left and right to H,
473: * thus creating a nonzero bulge below the subdiagonal.
474: *
475: * Each subsequent iteration determines a reflection G to
476: * restore the Hessenberg form in the (K-1)th column, and thus
477: * chases the bulge one step toward the bottom of the active
478: * submatrix. NR is the order of G.
479: *
480: NR = MIN( 3, I-K+1 )
481: IF( K.GT.M )
482: $ CALL DCOPY( NR, H( K, K-1 ), 1, V, 1 )
483: CALL DLARFG( NR, V( 1 ), V( 2 ), 1, T1 )
484: IF( K.GT.M ) THEN
485: H( K, K-1 ) = V( 1 )
486: H( K+1, K-1 ) = ZERO
487: IF( K.LT.I-1 )
488: $ H( K+2, K-1 ) = ZERO
489: ELSE IF( M.GT.L ) THEN
490: * ==== Use the following instead of
491: * . H( K, K-1 ) = -H( K, K-1 ) to
492: * . avoid a bug when v(2) and v(3)
493: * . underflow. ====
494: H( K, K-1 ) = H( K, K-1 )*( ONE-T1 )
495: END IF
496: V2 = V( 2 )
497: T2 = T1*V2
498: IF( NR.EQ.3 ) THEN
499: V3 = V( 3 )
500: T3 = T1*V3
501: *
502: * Apply G from the left to transform the rows of the matrix
503: * in columns K to I2.
504: *
505: DO 70 J = K, I2
506: SUM = H( K, J ) + V2*H( K+1, J ) + V3*H( K+2, J )
507: H( K, J ) = H( K, J ) - SUM*T1
508: H( K+1, J ) = H( K+1, J ) - SUM*T2
509: H( K+2, J ) = H( K+2, J ) - SUM*T3
510: 70 CONTINUE
511: *
512: * Apply G from the right to transform the columns of the
513: * matrix in rows I1 to min(K+3,I).
514: *
515: DO 80 J = I1, MIN( K+3, I )
516: SUM = H( J, K ) + V2*H( J, K+1 ) + V3*H( J, K+2 )
517: H( J, K ) = H( J, K ) - SUM*T1
518: H( J, K+1 ) = H( J, K+1 ) - SUM*T2
519: H( J, K+2 ) = H( J, K+2 ) - SUM*T3
520: 80 CONTINUE
521: *
522: IF( WANTZ ) THEN
523: *
524: * Accumulate transformations in the matrix Z
525: *
526: DO 90 J = ILOZ, IHIZ
527: SUM = Z( J, K ) + V2*Z( J, K+1 ) + V3*Z( J, K+2 )
528: Z( J, K ) = Z( J, K ) - SUM*T1
529: Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2
530: Z( J, K+2 ) = Z( J, K+2 ) - SUM*T3
531: 90 CONTINUE
532: END IF
533: ELSE IF( NR.EQ.2 ) THEN
534: *
535: * Apply G from the left to transform the rows of the matrix
536: * in columns K to I2.
537: *
538: DO 100 J = K, I2
539: SUM = H( K, J ) + V2*H( K+1, J )
540: H( K, J ) = H( K, J ) - SUM*T1
541: H( K+1, J ) = H( K+1, J ) - SUM*T2
542: 100 CONTINUE
543: *
544: * Apply G from the right to transform the columns of the
545: * matrix in rows I1 to min(K+3,I).
546: *
547: DO 110 J = I1, I
548: SUM = H( J, K ) + V2*H( J, K+1 )
549: H( J, K ) = H( J, K ) - SUM*T1
550: H( J, K+1 ) = H( J, K+1 ) - SUM*T2
551: 110 CONTINUE
552: *
553: IF( WANTZ ) THEN
554: *
555: * Accumulate transformations in the matrix Z
556: *
557: DO 120 J = ILOZ, IHIZ
558: SUM = Z( J, K ) + V2*Z( J, K+1 )
559: Z( J, K ) = Z( J, K ) - SUM*T1
560: Z( J, K+1 ) = Z( J, K+1 ) - SUM*T2
561: 120 CONTINUE
562: END IF
563: END IF
564: 130 CONTINUE
565: *
566: 140 CONTINUE
567: *
568: * Failure to converge in remaining number of iterations
569: *
570: INFO = I
571: RETURN
572: *
573: 150 CONTINUE
574: *
575: IF( L.EQ.I ) THEN
576: *
577: * H(I,I-1) is negligible: one eigenvalue has converged.
578: *
579: WR( I ) = H( I, I )
580: WI( I ) = ZERO
581: ELSE IF( L.EQ.I-1 ) THEN
582: *
583: * H(I-1,I-2) is negligible: a pair of eigenvalues have converged.
584: *
585: * Transform the 2-by-2 submatrix to standard Schur form,
586: * and compute and store the eigenvalues.
587: *
588: CALL DLANV2( H( I-1, I-1 ), H( I-1, I ), H( I, I-1 ),
589: $ H( I, I ), WR( I-1 ), WI( I-1 ), WR( I ), WI( I ),
590: $ CS, SN )
591: *
592: IF( WANTT ) THEN
593: *
594: * Apply the transformation to the rest of H.
595: *
596: IF( I2.GT.I )
597: $ CALL DROT( I2-I, H( I-1, I+1 ), LDH, H( I, I+1 ), LDH,
598: $ CS, SN )
599: CALL DROT( I-I1-1, H( I1, I-1 ), 1, H( I1, I ), 1, CS, SN )
600: END IF
601: IF( WANTZ ) THEN
602: *
603: * Apply the transformation to Z.
604: *
605: CALL DROT( NZ, Z( ILOZ, I-1 ), 1, Z( ILOZ, I ), 1, CS, SN )
606: END IF
607: END IF
1.20 ! bertrand 608: * reset deflation counter
! 609: KDEFL = 0
1.1 bertrand 610: *
611: * return to start of the main loop with new value of I.
612: *
613: I = L - 1
614: GO TO 20
615: *
616: 160 CONTINUE
617: RETURN
618: *
619: * End of DLAHQR
620: *
621: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlahqr.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack