Annotation of rpl/lapack/lapack/dlaqr3.f, revision 1.22
1.12 bertrand 1: *> \brief \b DLAQR3 performs the orthogonal similarity transformation of a Hessenberg matrix to detect and deflate fully converged eigenvalues from a trailing principal submatrix (aggressive early deflation).
1.9 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.17 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.9 bertrand 7: *
8: *> \htmlonly
1.17 bertrand 9: *> Download DLAQR3 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqr3.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqr3.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqr3.f">
1.9 bertrand 15: *> [TXT]</a>
1.17 bertrand 16: *> \endhtmlonly
1.9 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLAQR3( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
22: * IHIZ, Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T,
23: * LDT, NV, WV, LDWV, WORK, LWORK )
1.17 bertrand 24: *
1.9 bertrand 25: * .. Scalar Arguments ..
26: * INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
27: * $ LDZ, LWORK, N, ND, NH, NS, NV, NW
28: * LOGICAL WANTT, WANTZ
29: * ..
30: * .. Array Arguments ..
31: * DOUBLE PRECISION H( LDH, * ), SI( * ), SR( * ), T( LDT, * ),
32: * $ V( LDV, * ), WORK( * ), WV( LDWV, * ),
33: * $ Z( LDZ, * )
34: * ..
1.17 bertrand 35: *
1.9 bertrand 36: *
37: *> \par Purpose:
38: * =============
39: *>
40: *> \verbatim
41: *>
42: *> Aggressive early deflation:
43: *>
44: *> DLAQR3 accepts as input an upper Hessenberg matrix
45: *> H and performs an orthogonal similarity transformation
46: *> designed to detect and deflate fully converged eigenvalues from
47: *> a trailing principal submatrix. On output H has been over-
48: *> written by a new Hessenberg matrix that is a perturbation of
49: *> an orthogonal similarity transformation of H. It is to be
50: *> hoped that the final version of H has many zero subdiagonal
51: *> entries.
52: *> \endverbatim
53: *
54: * Arguments:
55: * ==========
56: *
57: *> \param[in] WANTT
58: *> \verbatim
59: *> WANTT is LOGICAL
60: *> If .TRUE., then the Hessenberg matrix H is fully updated
61: *> so that the quasi-triangular Schur factor may be
62: *> computed (in cooperation with the calling subroutine).
63: *> If .FALSE., then only enough of H is updated to preserve
64: *> the eigenvalues.
65: *> \endverbatim
66: *>
67: *> \param[in] WANTZ
68: *> \verbatim
69: *> WANTZ is LOGICAL
70: *> If .TRUE., then the orthogonal matrix Z is updated so
71: *> so that the orthogonal Schur factor may be computed
72: *> (in cooperation with the calling subroutine).
73: *> If .FALSE., then Z is not referenced.
74: *> \endverbatim
75: *>
76: *> \param[in] N
77: *> \verbatim
78: *> N is INTEGER
79: *> The order of the matrix H and (if WANTZ is .TRUE.) the
80: *> order of the orthogonal matrix Z.
81: *> \endverbatim
82: *>
83: *> \param[in] KTOP
84: *> \verbatim
85: *> KTOP is INTEGER
86: *> It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0.
87: *> KBOT and KTOP together determine an isolated block
88: *> along the diagonal of the Hessenberg matrix.
89: *> \endverbatim
90: *>
91: *> \param[in] KBOT
92: *> \verbatim
93: *> KBOT is INTEGER
94: *> It is assumed without a check that either
95: *> KBOT = N or H(KBOT+1,KBOT)=0. KBOT and KTOP together
96: *> determine an isolated block along the diagonal of the
97: *> Hessenberg matrix.
98: *> \endverbatim
99: *>
100: *> \param[in] NW
101: *> \verbatim
102: *> NW is INTEGER
1.21 bertrand 103: *> Deflation window size. 1 <= NW <= (KBOT-KTOP+1).
1.9 bertrand 104: *> \endverbatim
105: *>
106: *> \param[in,out] H
107: *> \verbatim
108: *> H is DOUBLE PRECISION array, dimension (LDH,N)
109: *> On input the initial N-by-N section of H stores the
110: *> Hessenberg matrix undergoing aggressive early deflation.
111: *> On output H has been transformed by an orthogonal
112: *> similarity transformation, perturbed, and the returned
113: *> to Hessenberg form that (it is to be hoped) has some
114: *> zero subdiagonal entries.
115: *> \endverbatim
116: *>
117: *> \param[in] LDH
118: *> \verbatim
1.19 bertrand 119: *> LDH is INTEGER
1.9 bertrand 120: *> Leading dimension of H just as declared in the calling
1.21 bertrand 121: *> subroutine. N <= LDH
1.9 bertrand 122: *> \endverbatim
123: *>
124: *> \param[in] ILOZ
125: *> \verbatim
126: *> ILOZ is INTEGER
127: *> \endverbatim
128: *>
129: *> \param[in] IHIZ
130: *> \verbatim
131: *> IHIZ is INTEGER
132: *> Specify the rows of Z to which transformations must be
1.21 bertrand 133: *> applied if WANTZ is .TRUE.. 1 <= ILOZ <= IHIZ <= N.
1.9 bertrand 134: *> \endverbatim
135: *>
136: *> \param[in,out] Z
137: *> \verbatim
138: *> Z is DOUBLE PRECISION array, dimension (LDZ,N)
139: *> IF WANTZ is .TRUE., then on output, the orthogonal
140: *> similarity transformation mentioned above has been
1.15 bertrand 141: *> accumulated into Z(ILOZ:IHIZ,ILOZ:IHIZ) from the right.
1.9 bertrand 142: *> If WANTZ is .FALSE., then Z is unreferenced.
143: *> \endverbatim
144: *>
145: *> \param[in] LDZ
146: *> \verbatim
1.19 bertrand 147: *> LDZ is INTEGER
1.9 bertrand 148: *> The leading dimension of Z just as declared in the
1.21 bertrand 149: *> calling subroutine. 1 <= LDZ.
1.9 bertrand 150: *> \endverbatim
151: *>
152: *> \param[out] NS
153: *> \verbatim
1.19 bertrand 154: *> NS is INTEGER
1.9 bertrand 155: *> The number of unconverged (ie approximate) eigenvalues
156: *> returned in SR and SI that may be used as shifts by the
157: *> calling subroutine.
158: *> \endverbatim
159: *>
160: *> \param[out] ND
161: *> \verbatim
1.19 bertrand 162: *> ND is INTEGER
1.9 bertrand 163: *> The number of converged eigenvalues uncovered by this
164: *> subroutine.
165: *> \endverbatim
166: *>
167: *> \param[out] SR
168: *> \verbatim
169: *> SR is DOUBLE PRECISION array, dimension (KBOT)
170: *> \endverbatim
171: *>
172: *> \param[out] SI
173: *> \verbatim
174: *> SI is DOUBLE PRECISION array, dimension (KBOT)
175: *> On output, the real and imaginary parts of approximate
176: *> eigenvalues that may be used for shifts are stored in
177: *> SR(KBOT-ND-NS+1) through SR(KBOT-ND) and
178: *> SI(KBOT-ND-NS+1) through SI(KBOT-ND), respectively.
179: *> The real and imaginary parts of converged eigenvalues
180: *> are stored in SR(KBOT-ND+1) through SR(KBOT) and
181: *> SI(KBOT-ND+1) through SI(KBOT), respectively.
182: *> \endverbatim
183: *>
184: *> \param[out] V
185: *> \verbatim
186: *> V is DOUBLE PRECISION array, dimension (LDV,NW)
187: *> An NW-by-NW work array.
188: *> \endverbatim
189: *>
190: *> \param[in] LDV
191: *> \verbatim
1.19 bertrand 192: *> LDV is INTEGER
1.9 bertrand 193: *> The leading dimension of V just as declared in the
1.21 bertrand 194: *> calling subroutine. NW <= LDV
1.9 bertrand 195: *> \endverbatim
196: *>
197: *> \param[in] NH
198: *> \verbatim
1.19 bertrand 199: *> NH is INTEGER
1.21 bertrand 200: *> The number of columns of T. NH >= NW.
1.9 bertrand 201: *> \endverbatim
202: *>
203: *> \param[out] T
204: *> \verbatim
205: *> T is DOUBLE PRECISION array, dimension (LDT,NW)
206: *> \endverbatim
207: *>
208: *> \param[in] LDT
209: *> \verbatim
1.19 bertrand 210: *> LDT is INTEGER
1.9 bertrand 211: *> The leading dimension of T just as declared in the
1.21 bertrand 212: *> calling subroutine. NW <= LDT
1.9 bertrand 213: *> \endverbatim
214: *>
215: *> \param[in] NV
216: *> \verbatim
1.19 bertrand 217: *> NV is INTEGER
1.9 bertrand 218: *> The number of rows of work array WV available for
1.21 bertrand 219: *> workspace. NV >= NW.
1.9 bertrand 220: *> \endverbatim
221: *>
222: *> \param[out] WV
223: *> \verbatim
224: *> WV is DOUBLE PRECISION array, dimension (LDWV,NW)
225: *> \endverbatim
226: *>
227: *> \param[in] LDWV
228: *> \verbatim
1.19 bertrand 229: *> LDWV is INTEGER
1.9 bertrand 230: *> The leading dimension of W just as declared in the
1.21 bertrand 231: *> calling subroutine. NW <= LDV
1.9 bertrand 232: *> \endverbatim
233: *>
234: *> \param[out] WORK
235: *> \verbatim
236: *> WORK is DOUBLE PRECISION array, dimension (LWORK)
237: *> On exit, WORK(1) is set to an estimate of the optimal value
238: *> of LWORK for the given values of N, NW, KTOP and KBOT.
239: *> \endverbatim
240: *>
241: *> \param[in] LWORK
242: *> \verbatim
1.19 bertrand 243: *> LWORK is INTEGER
1.9 bertrand 244: *> The dimension of the work array WORK. LWORK = 2*NW
245: *> suffices, but greater efficiency may result from larger
246: *> values of LWORK.
247: *>
248: *> If LWORK = -1, then a workspace query is assumed; DLAQR3
249: *> only estimates the optimal workspace size for the given
250: *> values of N, NW, KTOP and KBOT. The estimate is returned
251: *> in WORK(1). No error message related to LWORK is issued
252: *> by XERBLA. Neither H nor Z are accessed.
253: *> \endverbatim
254: *
255: * Authors:
256: * ========
257: *
1.17 bertrand 258: *> \author Univ. of Tennessee
259: *> \author Univ. of California Berkeley
260: *> \author Univ. of Colorado Denver
261: *> \author NAG Ltd.
1.9 bertrand 262: *
263: *> \ingroup doubleOTHERauxiliary
264: *
265: *> \par Contributors:
266: * ==================
267: *>
268: *> Karen Braman and Ralph Byers, Department of Mathematics,
269: *> University of Kansas, USA
270: *>
271: * =====================================================================
1.1 bertrand 272: SUBROUTINE DLAQR3( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
273: $ IHIZ, Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T,
274: $ LDT, NV, WV, LDWV, WORK, LWORK )
275: *
1.22 ! bertrand 276: * -- LAPACK auxiliary routine --
1.9 bertrand 277: * -- LAPACK is a software package provided by Univ. of Tennessee, --
278: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.1 bertrand 279: *
280: * .. Scalar Arguments ..
281: INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
282: $ LDZ, LWORK, N, ND, NH, NS, NV, NW
283: LOGICAL WANTT, WANTZ
284: * ..
285: * .. Array Arguments ..
286: DOUBLE PRECISION H( LDH, * ), SI( * ), SR( * ), T( LDT, * ),
287: $ V( LDV, * ), WORK( * ), WV( LDWV, * ),
288: $ Z( LDZ, * )
289: * ..
290: *
1.9 bertrand 291: * ================================================================
1.1 bertrand 292: * .. Parameters ..
293: DOUBLE PRECISION ZERO, ONE
294: PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0 )
295: * ..
296: * .. Local Scalars ..
297: DOUBLE PRECISION AA, BB, BETA, CC, CS, DD, EVI, EVK, FOO, S,
298: $ SAFMAX, SAFMIN, SMLNUM, SN, TAU, ULP
299: INTEGER I, IFST, ILST, INFO, INFQR, J, JW, K, KCOL,
300: $ KEND, KLN, KROW, KWTOP, LTOP, LWK1, LWK2, LWK3,
301: $ LWKOPT, NMIN
302: LOGICAL BULGE, SORTED
303: * ..
304: * .. External Functions ..
305: DOUBLE PRECISION DLAMCH
306: INTEGER ILAENV
307: EXTERNAL DLAMCH, ILAENV
308: * ..
309: * .. External Subroutines ..
310: EXTERNAL DCOPY, DGEHRD, DGEMM, DLABAD, DLACPY, DLAHQR,
311: $ DLANV2, DLAQR4, DLARF, DLARFG, DLASET, DORMHR,
312: $ DTREXC
313: * ..
314: * .. Intrinsic Functions ..
315: INTRINSIC ABS, DBLE, INT, MAX, MIN, SQRT
316: * ..
317: * .. Executable Statements ..
318: *
319: * ==== Estimate optimal workspace. ====
320: *
321: JW = MIN( NW, KBOT-KTOP+1 )
322: IF( JW.LE.2 ) THEN
323: LWKOPT = 1
324: ELSE
325: *
326: * ==== Workspace query call to DGEHRD ====
327: *
328: CALL DGEHRD( JW, 1, JW-1, T, LDT, WORK, WORK, -1, INFO )
329: LWK1 = INT( WORK( 1 ) )
330: *
331: * ==== Workspace query call to DORMHR ====
332: *
333: CALL DORMHR( 'R', 'N', JW, JW, 1, JW-1, T, LDT, WORK, V, LDV,
334: $ WORK, -1, INFO )
335: LWK2 = INT( WORK( 1 ) )
336: *
337: * ==== Workspace query call to DLAQR4 ====
338: *
339: CALL DLAQR4( .true., .true., JW, 1, JW, T, LDT, SR, SI, 1, JW,
340: $ V, LDV, WORK, -1, INFQR )
341: LWK3 = INT( WORK( 1 ) )
342: *
343: * ==== Optimal workspace ====
344: *
345: LWKOPT = MAX( JW+MAX( LWK1, LWK2 ), LWK3 )
346: END IF
347: *
348: * ==== Quick return in case of workspace query. ====
349: *
350: IF( LWORK.EQ.-1 ) THEN
351: WORK( 1 ) = DBLE( LWKOPT )
352: RETURN
353: END IF
354: *
355: * ==== Nothing to do ...
356: * ... for an empty active block ... ====
357: NS = 0
358: ND = 0
359: WORK( 1 ) = ONE
360: IF( KTOP.GT.KBOT )
361: $ RETURN
362: * ... nor for an empty deflation window. ====
363: IF( NW.LT.1 )
364: $ RETURN
365: *
366: * ==== Machine constants ====
367: *
368: SAFMIN = DLAMCH( 'SAFE MINIMUM' )
369: SAFMAX = ONE / SAFMIN
370: CALL DLABAD( SAFMIN, SAFMAX )
371: ULP = DLAMCH( 'PRECISION' )
372: SMLNUM = SAFMIN*( DBLE( N ) / ULP )
373: *
374: * ==== Setup deflation window ====
375: *
376: JW = MIN( NW, KBOT-KTOP+1 )
377: KWTOP = KBOT - JW + 1
378: IF( KWTOP.EQ.KTOP ) THEN
379: S = ZERO
380: ELSE
381: S = H( KWTOP, KWTOP-1 )
382: END IF
383: *
384: IF( KBOT.EQ.KWTOP ) THEN
385: *
386: * ==== 1-by-1 deflation window: not much to do ====
387: *
388: SR( KWTOP ) = H( KWTOP, KWTOP )
389: SI( KWTOP ) = ZERO
390: NS = 1
391: ND = 0
392: IF( ABS( S ).LE.MAX( SMLNUM, ULP*ABS( H( KWTOP, KWTOP ) ) ) )
393: $ THEN
394: NS = 0
395: ND = 1
396: IF( KWTOP.GT.KTOP )
397: $ H( KWTOP, KWTOP-1 ) = ZERO
398: END IF
399: WORK( 1 ) = ONE
400: RETURN
401: END IF
402: *
403: * ==== Convert to spike-triangular form. (In case of a
404: * . rare QR failure, this routine continues to do
405: * . aggressive early deflation using that part of
406: * . the deflation window that converged using INFQR
407: * . here and there to keep track.) ====
408: *
409: CALL DLACPY( 'U', JW, JW, H( KWTOP, KWTOP ), LDH, T, LDT )
410: CALL DCOPY( JW-1, H( KWTOP+1, KWTOP ), LDH+1, T( 2, 1 ), LDT+1 )
411: *
412: CALL DLASET( 'A', JW, JW, ZERO, ONE, V, LDV )
413: NMIN = ILAENV( 12, 'DLAQR3', 'SV', JW, 1, JW, LWORK )
414: IF( JW.GT.NMIN ) THEN
415: CALL DLAQR4( .true., .true., JW, 1, JW, T, LDT, SR( KWTOP ),
416: $ SI( KWTOP ), 1, JW, V, LDV, WORK, LWORK, INFQR )
417: ELSE
418: CALL DLAHQR( .true., .true., JW, 1, JW, T, LDT, SR( KWTOP ),
419: $ SI( KWTOP ), 1, JW, V, LDV, INFQR )
420: END IF
421: *
422: * ==== DTREXC needs a clean margin near the diagonal ====
423: *
424: DO 10 J = 1, JW - 3
425: T( J+2, J ) = ZERO
426: T( J+3, J ) = ZERO
427: 10 CONTINUE
428: IF( JW.GT.2 )
429: $ T( JW, JW-2 ) = ZERO
430: *
431: * ==== Deflation detection loop ====
432: *
433: NS = JW
434: ILST = INFQR + 1
435: 20 CONTINUE
436: IF( ILST.LE.NS ) THEN
437: IF( NS.EQ.1 ) THEN
438: BULGE = .FALSE.
439: ELSE
440: BULGE = T( NS, NS-1 ).NE.ZERO
441: END IF
442: *
443: * ==== Small spike tip test for deflation ====
444: *
1.9 bertrand 445: IF( .NOT. BULGE ) THEN
1.1 bertrand 446: *
447: * ==== Real eigenvalue ====
448: *
449: FOO = ABS( T( NS, NS ) )
450: IF( FOO.EQ.ZERO )
451: $ FOO = ABS( S )
452: IF( ABS( S*V( 1, NS ) ).LE.MAX( SMLNUM, ULP*FOO ) ) THEN
453: *
454: * ==== Deflatable ====
455: *
456: NS = NS - 1
457: ELSE
458: *
459: * ==== Undeflatable. Move it up out of the way.
460: * . (DTREXC can not fail in this case.) ====
461: *
462: IFST = NS
463: CALL DTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, WORK,
464: $ INFO )
465: ILST = ILST + 1
466: END IF
467: ELSE
468: *
469: * ==== Complex conjugate pair ====
470: *
471: FOO = ABS( T( NS, NS ) ) + SQRT( ABS( T( NS, NS-1 ) ) )*
472: $ SQRT( ABS( T( NS-1, NS ) ) )
473: IF( FOO.EQ.ZERO )
474: $ FOO = ABS( S )
475: IF( MAX( ABS( S*V( 1, NS ) ), ABS( S*V( 1, NS-1 ) ) ).LE.
476: $ MAX( SMLNUM, ULP*FOO ) ) THEN
477: *
478: * ==== Deflatable ====
479: *
480: NS = NS - 2
481: ELSE
482: *
483: * ==== Undeflatable. Move them up out of the way.
484: * . Fortunately, DTREXC does the right thing with
485: * . ILST in case of a rare exchange failure. ====
486: *
487: IFST = NS
488: CALL DTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, WORK,
489: $ INFO )
490: ILST = ILST + 2
491: END IF
492: END IF
493: *
494: * ==== End deflation detection loop ====
495: *
496: GO TO 20
497: END IF
498: *
499: * ==== Return to Hessenberg form ====
500: *
501: IF( NS.EQ.0 )
502: $ S = ZERO
503: *
504: IF( NS.LT.JW ) THEN
505: *
506: * ==== sorting diagonal blocks of T improves accuracy for
507: * . graded matrices. Bubble sort deals well with
508: * . exchange failures. ====
509: *
510: SORTED = .false.
511: I = NS + 1
512: 30 CONTINUE
513: IF( SORTED )
514: $ GO TO 50
515: SORTED = .true.
516: *
517: KEND = I - 1
518: I = INFQR + 1
519: IF( I.EQ.NS ) THEN
520: K = I + 1
521: ELSE IF( T( I+1, I ).EQ.ZERO ) THEN
522: K = I + 1
523: ELSE
524: K = I + 2
525: END IF
526: 40 CONTINUE
527: IF( K.LE.KEND ) THEN
528: IF( K.EQ.I+1 ) THEN
529: EVI = ABS( T( I, I ) )
530: ELSE
531: EVI = ABS( T( I, I ) ) + SQRT( ABS( T( I+1, I ) ) )*
532: $ SQRT( ABS( T( I, I+1 ) ) )
533: END IF
534: *
535: IF( K.EQ.KEND ) THEN
536: EVK = ABS( T( K, K ) )
537: ELSE IF( T( K+1, K ).EQ.ZERO ) THEN
538: EVK = ABS( T( K, K ) )
539: ELSE
540: EVK = ABS( T( K, K ) ) + SQRT( ABS( T( K+1, K ) ) )*
541: $ SQRT( ABS( T( K, K+1 ) ) )
542: END IF
543: *
544: IF( EVI.GE.EVK ) THEN
545: I = K
546: ELSE
547: SORTED = .false.
548: IFST = I
549: ILST = K
550: CALL DTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, WORK,
551: $ INFO )
552: IF( INFO.EQ.0 ) THEN
553: I = ILST
554: ELSE
555: I = K
556: END IF
557: END IF
558: IF( I.EQ.KEND ) THEN
559: K = I + 1
560: ELSE IF( T( I+1, I ).EQ.ZERO ) THEN
561: K = I + 1
562: ELSE
563: K = I + 2
564: END IF
565: GO TO 40
566: END IF
567: GO TO 30
568: 50 CONTINUE
569: END IF
570: *
571: * ==== Restore shift/eigenvalue array from T ====
572: *
573: I = JW
574: 60 CONTINUE
575: IF( I.GE.INFQR+1 ) THEN
576: IF( I.EQ.INFQR+1 ) THEN
577: SR( KWTOP+I-1 ) = T( I, I )
578: SI( KWTOP+I-1 ) = ZERO
579: I = I - 1
580: ELSE IF( T( I, I-1 ).EQ.ZERO ) THEN
581: SR( KWTOP+I-1 ) = T( I, I )
582: SI( KWTOP+I-1 ) = ZERO
583: I = I - 1
584: ELSE
585: AA = T( I-1, I-1 )
586: CC = T( I, I-1 )
587: BB = T( I-1, I )
588: DD = T( I, I )
589: CALL DLANV2( AA, BB, CC, DD, SR( KWTOP+I-2 ),
590: $ SI( KWTOP+I-2 ), SR( KWTOP+I-1 ),
591: $ SI( KWTOP+I-1 ), CS, SN )
592: I = I - 2
593: END IF
594: GO TO 60
595: END IF
596: *
597: IF( NS.LT.JW .OR. S.EQ.ZERO ) THEN
598: IF( NS.GT.1 .AND. S.NE.ZERO ) THEN
599: *
600: * ==== Reflect spike back into lower triangle ====
601: *
602: CALL DCOPY( NS, V, LDV, WORK, 1 )
603: BETA = WORK( 1 )
604: CALL DLARFG( NS, BETA, WORK( 2 ), 1, TAU )
605: WORK( 1 ) = ONE
606: *
607: CALL DLASET( 'L', JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), LDT )
608: *
609: CALL DLARF( 'L', NS, JW, WORK, 1, TAU, T, LDT,
610: $ WORK( JW+1 ) )
611: CALL DLARF( 'R', NS, NS, WORK, 1, TAU, T, LDT,
612: $ WORK( JW+1 ) )
613: CALL DLARF( 'R', JW, NS, WORK, 1, TAU, V, LDV,
614: $ WORK( JW+1 ) )
615: *
616: CALL DGEHRD( JW, 1, NS, T, LDT, WORK, WORK( JW+1 ),
617: $ LWORK-JW, INFO )
618: END IF
619: *
620: * ==== Copy updated reduced window into place ====
621: *
622: IF( KWTOP.GT.1 )
623: $ H( KWTOP, KWTOP-1 ) = S*V( 1, 1 )
624: CALL DLACPY( 'U', JW, JW, T, LDT, H( KWTOP, KWTOP ), LDH )
625: CALL DCOPY( JW-1, T( 2, 1 ), LDT+1, H( KWTOP+1, KWTOP ),
626: $ LDH+1 )
627: *
628: * ==== Accumulate orthogonal matrix in order update
629: * . H and Z, if requested. ====
630: *
631: IF( NS.GT.1 .AND. S.NE.ZERO )
632: $ CALL DORMHR( 'R', 'N', JW, NS, 1, NS, T, LDT, WORK, V, LDV,
633: $ WORK( JW+1 ), LWORK-JW, INFO )
634: *
635: * ==== Update vertical slab in H ====
636: *
637: IF( WANTT ) THEN
638: LTOP = 1
639: ELSE
640: LTOP = KTOP
641: END IF
642: DO 70 KROW = LTOP, KWTOP - 1, NV
643: KLN = MIN( NV, KWTOP-KROW )
644: CALL DGEMM( 'N', 'N', KLN, JW, JW, ONE, H( KROW, KWTOP ),
645: $ LDH, V, LDV, ZERO, WV, LDWV )
646: CALL DLACPY( 'A', KLN, JW, WV, LDWV, H( KROW, KWTOP ), LDH )
647: 70 CONTINUE
648: *
649: * ==== Update horizontal slab in H ====
650: *
651: IF( WANTT ) THEN
652: DO 80 KCOL = KBOT + 1, N, NH
653: KLN = MIN( NH, N-KCOL+1 )
654: CALL DGEMM( 'C', 'N', JW, KLN, JW, ONE, V, LDV,
655: $ H( KWTOP, KCOL ), LDH, ZERO, T, LDT )
656: CALL DLACPY( 'A', JW, KLN, T, LDT, H( KWTOP, KCOL ),
657: $ LDH )
658: 80 CONTINUE
659: END IF
660: *
661: * ==== Update vertical slab in Z ====
662: *
663: IF( WANTZ ) THEN
664: DO 90 KROW = ILOZ, IHIZ, NV
665: KLN = MIN( NV, IHIZ-KROW+1 )
666: CALL DGEMM( 'N', 'N', KLN, JW, JW, ONE, Z( KROW, KWTOP ),
667: $ LDZ, V, LDV, ZERO, WV, LDWV )
668: CALL DLACPY( 'A', KLN, JW, WV, LDWV, Z( KROW, KWTOP ),
669: $ LDZ )
670: 90 CONTINUE
671: END IF
672: END IF
673: *
674: * ==== Return the number of deflations ... ====
675: *
676: ND = JW - NS
677: *
678: * ==== ... and the number of shifts. (Subtracting
679: * . INFQR from the spike length takes care
680: * . of the case of a rare QR failure while
681: * . calculating eigenvalues of the deflation
682: * . window.) ====
683: *
684: NS = NS - INFQR
685: *
686: * ==== Return optimal workspace. ====
687: *
688: WORK( 1 ) = DBLE( LWKOPT )
689: *
690: * ==== End of DLAQR3 ====
691: *
692: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaqr3.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack