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