Annotation of rpl/lapack/lapack/dlaqr2.f, revision 1.9
1.9 ! bertrand 1: *> \brief \b DLAQR2
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 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">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 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 )
! 24: *
! 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: * ..
! 35: *
! 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
! 106: *> Deflation window size. 1 .LE. NW .LE. (KBOT-KTOP+1).
! 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
! 122: *> LDH is integer
! 123: *> Leading dimension of H just as declared in the calling
! 124: *> subroutine. N .LE. LDH
! 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
! 136: *> applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.
! 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
! 144: *> accumulated into Z(ILOZ:IHIZ,ILO:IHI) from the right.
! 145: *> If WANTZ is .FALSE., then Z is unreferenced.
! 146: *> \endverbatim
! 147: *>
! 148: *> \param[in] LDZ
! 149: *> \verbatim
! 150: *> LDZ is integer
! 151: *> The leading dimension of Z just as declared in the
! 152: *> calling subroutine. 1 .LE. LDZ.
! 153: *> \endverbatim
! 154: *>
! 155: *> \param[out] NS
! 156: *> \verbatim
! 157: *> NS is integer
! 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
! 165: *> ND is integer
! 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
! 195: *> LDV is integer scalar
! 196: *> The leading dimension of V just as declared in the
! 197: *> calling subroutine. NW .LE. LDV
! 198: *> \endverbatim
! 199: *>
! 200: *> \param[in] NH
! 201: *> \verbatim
! 202: *> NH is integer scalar
! 203: *> The number of columns of T. NH.GE.NW.
! 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
! 213: *> LDT is integer
! 214: *> The leading dimension of T just as declared in the
! 215: *> calling subroutine. NW .LE. LDT
! 216: *> \endverbatim
! 217: *>
! 218: *> \param[in] NV
! 219: *> \verbatim
! 220: *> NV is integer
! 221: *> The number of rows of work array WV available for
! 222: *> workspace. NV.GE.NW.
! 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
! 232: *> LDWV is integer
! 233: *> The leading dimension of W just as declared in the
! 234: *> calling subroutine. NW .LE. LDV
! 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
! 246: *> LWORK is integer
! 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: *
! 261: *> \author Univ. of Tennessee
! 262: *> \author Univ. of California Berkeley
! 263: *> \author Univ. of Colorado Denver
! 264: *> \author NAG Ltd.
! 265: *
! 266: *> \date November 2011
! 267: *
! 268: *> \ingroup doubleOTHERauxiliary
! 269: *
! 270: *> \par Contributors:
! 271: * ==================
! 272: *>
! 273: *> Karen Braman and Ralph Byers, Department of Mathematics,
! 274: *> University of Kansas, USA
! 275: *>
! 276: * =====================================================================
1.1 bertrand 277: SUBROUTINE DLAQR2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
278: $ IHIZ, Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T,
279: $ LDT, NV, WV, LDWV, WORK, LWORK )
280: *
1.9 ! bertrand 281: * -- LAPACK auxiliary routine (version 3.4.0) --
! 282: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 283: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 284: * November 2011
1.1 bertrand 285: *
286: * .. Scalar Arguments ..
287: INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
288: $ LDZ, LWORK, N, ND, NH, NS, NV, NW
289: LOGICAL WANTT, WANTZ
290: * ..
291: * .. Array Arguments ..
292: DOUBLE PRECISION H( LDH, * ), SI( * ), SR( * ), T( LDT, * ),
293: $ V( LDV, * ), WORK( * ), WV( LDWV, * ),
294: $ Z( LDZ, * )
295: * ..
296: *
1.9 ! bertrand 297: * ================================================================
1.1 bertrand 298: * .. Parameters ..
299: DOUBLE PRECISION ZERO, ONE
300: PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0 )
301: * ..
302: * .. Local Scalars ..
303: DOUBLE PRECISION AA, BB, BETA, CC, CS, DD, EVI, EVK, FOO, S,
304: $ SAFMAX, SAFMIN, SMLNUM, SN, TAU, ULP
305: INTEGER I, IFST, ILST, INFO, INFQR, J, JW, K, KCOL,
306: $ KEND, KLN, KROW, KWTOP, LTOP, LWK1, LWK2,
307: $ LWKOPT
308: LOGICAL BULGE, SORTED
309: * ..
310: * .. External Functions ..
311: DOUBLE PRECISION DLAMCH
312: EXTERNAL DLAMCH
313: * ..
314: * .. External Subroutines ..
315: EXTERNAL DCOPY, DGEHRD, DGEMM, DLABAD, DLACPY, DLAHQR,
316: $ DLANV2, DLARF, DLARFG, DLASET, DORMHR, DTREXC
317: * ..
318: * .. Intrinsic Functions ..
319: INTRINSIC ABS, DBLE, INT, MAX, MIN, SQRT
320: * ..
321: * .. Executable Statements ..
322: *
323: * ==== Estimate optimal workspace. ====
324: *
325: JW = MIN( NW, KBOT-KTOP+1 )
326: IF( JW.LE.2 ) THEN
327: LWKOPT = 1
328: ELSE
329: *
330: * ==== Workspace query call to DGEHRD ====
331: *
332: CALL DGEHRD( JW, 1, JW-1, T, LDT, WORK, WORK, -1, INFO )
333: LWK1 = INT( WORK( 1 ) )
334: *
335: * ==== Workspace query call to DORMHR ====
336: *
337: CALL DORMHR( 'R', 'N', JW, JW, 1, JW-1, T, LDT, WORK, V, LDV,
338: $ WORK, -1, INFO )
339: LWK2 = INT( WORK( 1 ) )
340: *
341: * ==== Optimal workspace ====
342: *
343: LWKOPT = JW + MAX( LWK1, LWK2 )
344: END IF
345: *
346: * ==== Quick return in case of workspace query. ====
347: *
348: IF( LWORK.EQ.-1 ) THEN
349: WORK( 1 ) = DBLE( LWKOPT )
350: RETURN
351: END IF
352: *
353: * ==== Nothing to do ...
354: * ... for an empty active block ... ====
355: NS = 0
356: ND = 0
357: WORK( 1 ) = ONE
358: IF( KTOP.GT.KBOT )
359: $ RETURN
360: * ... nor for an empty deflation window. ====
361: IF( NW.LT.1 )
362: $ RETURN
363: *
364: * ==== Machine constants ====
365: *
366: SAFMIN = DLAMCH( 'SAFE MINIMUM' )
367: SAFMAX = ONE / SAFMIN
368: CALL DLABAD( SAFMIN, SAFMAX )
369: ULP = DLAMCH( 'PRECISION' )
370: SMLNUM = SAFMIN*( DBLE( N ) / ULP )
371: *
372: * ==== Setup deflation window ====
373: *
374: JW = MIN( NW, KBOT-KTOP+1 )
375: KWTOP = KBOT - JW + 1
376: IF( KWTOP.EQ.KTOP ) THEN
377: S = ZERO
378: ELSE
379: S = H( KWTOP, KWTOP-1 )
380: END IF
381: *
382: IF( KBOT.EQ.KWTOP ) THEN
383: *
384: * ==== 1-by-1 deflation window: not much to do ====
385: *
386: SR( KWTOP ) = H( KWTOP, KWTOP )
387: SI( KWTOP ) = ZERO
388: NS = 1
389: ND = 0
390: IF( ABS( S ).LE.MAX( SMLNUM, ULP*ABS( H( KWTOP, KWTOP ) ) ) )
391: $ THEN
392: NS = 0
393: ND = 1
394: IF( KWTOP.GT.KTOP )
395: $ H( KWTOP, KWTOP-1 ) = ZERO
396: END IF
397: WORK( 1 ) = ONE
398: RETURN
399: END IF
400: *
401: * ==== Convert to spike-triangular form. (In case of a
402: * . rare QR failure, this routine continues to do
403: * . aggressive early deflation using that part of
404: * . the deflation window that converged using INFQR
405: * . here and there to keep track.) ====
406: *
407: CALL DLACPY( 'U', JW, JW, H( KWTOP, KWTOP ), LDH, T, LDT )
408: CALL DCOPY( JW-1, H( KWTOP+1, KWTOP ), LDH+1, T( 2, 1 ), LDT+1 )
409: *
410: CALL DLASET( 'A', JW, JW, ZERO, ONE, V, LDV )
411: CALL DLAHQR( .true., .true., JW, 1, JW, T, LDT, SR( KWTOP ),
412: $ SI( KWTOP ), 1, JW, V, LDV, INFQR )
413: *
414: * ==== DTREXC needs a clean margin near the diagonal ====
415: *
416: DO 10 J = 1, JW - 3
417: T( J+2, J ) = ZERO
418: T( J+3, J ) = ZERO
419: 10 CONTINUE
420: IF( JW.GT.2 )
421: $ T( JW, JW-2 ) = ZERO
422: *
423: * ==== Deflation detection loop ====
424: *
425: NS = JW
426: ILST = INFQR + 1
427: 20 CONTINUE
428: IF( ILST.LE.NS ) THEN
429: IF( NS.EQ.1 ) THEN
430: BULGE = .FALSE.
431: ELSE
432: BULGE = T( NS, NS-1 ).NE.ZERO
433: END IF
434: *
435: * ==== Small spike tip test for deflation ====
436: *
437: IF( .NOT.BULGE ) THEN
438: *
439: * ==== Real eigenvalue ====
440: *
441: FOO = ABS( T( NS, NS ) )
442: IF( FOO.EQ.ZERO )
443: $ FOO = ABS( S )
444: IF( ABS( S*V( 1, NS ) ).LE.MAX( SMLNUM, ULP*FOO ) ) THEN
445: *
446: * ==== Deflatable ====
447: *
448: NS = NS - 1
449: ELSE
450: *
451: * ==== Undeflatable. Move it up out of the way.
452: * . (DTREXC can not fail in this case.) ====
453: *
454: IFST = NS
455: CALL DTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, WORK,
456: $ INFO )
457: ILST = ILST + 1
458: END IF
459: ELSE
460: *
461: * ==== Complex conjugate pair ====
462: *
463: FOO = ABS( T( NS, NS ) ) + SQRT( ABS( T( NS, NS-1 ) ) )*
464: $ SQRT( ABS( T( NS-1, NS ) ) )
465: IF( FOO.EQ.ZERO )
466: $ FOO = ABS( S )
467: IF( MAX( ABS( S*V( 1, NS ) ), ABS( S*V( 1, NS-1 ) ) ).LE.
468: $ MAX( SMLNUM, ULP*FOO ) ) THEN
469: *
470: * ==== Deflatable ====
471: *
472: NS = NS - 2
473: ELSE
474: *
475: * ==== Undeflatable. Move them up out of the way.
476: * . Fortunately, DTREXC does the right thing with
477: * . ILST in case of a rare exchange failure. ====
478: *
479: IFST = NS
480: CALL DTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, WORK,
481: $ INFO )
482: ILST = ILST + 2
483: END IF
484: END IF
485: *
486: * ==== End deflation detection loop ====
487: *
488: GO TO 20
489: END IF
490: *
491: * ==== Return to Hessenberg form ====
492: *
493: IF( NS.EQ.0 )
494: $ S = ZERO
495: *
496: IF( NS.LT.JW ) THEN
497: *
498: * ==== sorting diagonal blocks of T improves accuracy for
499: * . graded matrices. Bubble sort deals well with
500: * . exchange failures. ====
501: *
502: SORTED = .false.
503: I = NS + 1
504: 30 CONTINUE
505: IF( SORTED )
506: $ GO TO 50
507: SORTED = .true.
508: *
509: KEND = I - 1
510: I = INFQR + 1
511: IF( I.EQ.NS ) THEN
512: K = I + 1
513: ELSE IF( T( I+1, I ).EQ.ZERO ) THEN
514: K = I + 1
515: ELSE
516: K = I + 2
517: END IF
518: 40 CONTINUE
519: IF( K.LE.KEND ) THEN
520: IF( K.EQ.I+1 ) THEN
521: EVI = ABS( T( I, I ) )
522: ELSE
523: EVI = ABS( T( I, I ) ) + SQRT( ABS( T( I+1, I ) ) )*
524: $ SQRT( ABS( T( I, I+1 ) ) )
525: END IF
526: *
527: IF( K.EQ.KEND ) THEN
528: EVK = ABS( T( K, K ) )
529: ELSE IF( T( K+1, K ).EQ.ZERO ) THEN
530: EVK = ABS( T( K, K ) )
531: ELSE
532: EVK = ABS( T( K, K ) ) + SQRT( ABS( T( K+1, K ) ) )*
533: $ SQRT( ABS( T( K, K+1 ) ) )
534: END IF
535: *
536: IF( EVI.GE.EVK ) THEN
537: I = K
538: ELSE
539: SORTED = .false.
540: IFST = I
541: ILST = K
542: CALL DTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, WORK,
543: $ INFO )
544: IF( INFO.EQ.0 ) THEN
545: I = ILST
546: ELSE
547: I = K
548: END IF
549: END IF
550: IF( I.EQ.KEND ) THEN
551: K = I + 1
552: ELSE IF( T( I+1, I ).EQ.ZERO ) THEN
553: K = I + 1
554: ELSE
555: K = I + 2
556: END IF
557: GO TO 40
558: END IF
559: GO TO 30
560: 50 CONTINUE
561: END IF
562: *
563: * ==== Restore shift/eigenvalue array from T ====
564: *
565: I = JW
566: 60 CONTINUE
567: IF( I.GE.INFQR+1 ) THEN
568: IF( I.EQ.INFQR+1 ) THEN
569: SR( KWTOP+I-1 ) = T( I, I )
570: SI( KWTOP+I-1 ) = ZERO
571: I = I - 1
572: ELSE IF( T( I, I-1 ).EQ.ZERO ) THEN
573: SR( KWTOP+I-1 ) = T( I, I )
574: SI( KWTOP+I-1 ) = ZERO
575: I = I - 1
576: ELSE
577: AA = T( I-1, I-1 )
578: CC = T( I, I-1 )
579: BB = T( I-1, I )
580: DD = T( I, I )
581: CALL DLANV2( AA, BB, CC, DD, SR( KWTOP+I-2 ),
582: $ SI( KWTOP+I-2 ), SR( KWTOP+I-1 ),
583: $ SI( KWTOP+I-1 ), CS, SN )
584: I = I - 2
585: END IF
586: GO TO 60
587: END IF
588: *
589: IF( NS.LT.JW .OR. S.EQ.ZERO ) THEN
590: IF( NS.GT.1 .AND. S.NE.ZERO ) THEN
591: *
592: * ==== Reflect spike back into lower triangle ====
593: *
594: CALL DCOPY( NS, V, LDV, WORK, 1 )
595: BETA = WORK( 1 )
596: CALL DLARFG( NS, BETA, WORK( 2 ), 1, TAU )
597: WORK( 1 ) = ONE
598: *
599: CALL DLASET( 'L', JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), LDT )
600: *
601: CALL DLARF( 'L', NS, JW, WORK, 1, TAU, T, LDT,
602: $ WORK( JW+1 ) )
603: CALL DLARF( 'R', NS, NS, WORK, 1, TAU, T, LDT,
604: $ WORK( JW+1 ) )
605: CALL DLARF( 'R', JW, NS, WORK, 1, TAU, V, LDV,
606: $ WORK( JW+1 ) )
607: *
608: CALL DGEHRD( JW, 1, NS, T, LDT, WORK, WORK( JW+1 ),
609: $ LWORK-JW, INFO )
610: END IF
611: *
612: * ==== Copy updated reduced window into place ====
613: *
614: IF( KWTOP.GT.1 )
615: $ H( KWTOP, KWTOP-1 ) = S*V( 1, 1 )
616: CALL DLACPY( 'U', JW, JW, T, LDT, H( KWTOP, KWTOP ), LDH )
617: CALL DCOPY( JW-1, T( 2, 1 ), LDT+1, H( KWTOP+1, KWTOP ),
618: $ LDH+1 )
619: *
620: * ==== Accumulate orthogonal matrix in order update
621: * . H and Z, if requested. ====
622: *
623: IF( NS.GT.1 .AND. S.NE.ZERO )
624: $ CALL DORMHR( 'R', 'N', JW, NS, 1, NS, T, LDT, WORK, V, LDV,
625: $ WORK( JW+1 ), LWORK-JW, INFO )
626: *
627: * ==== Update vertical slab in H ====
628: *
629: IF( WANTT ) THEN
630: LTOP = 1
631: ELSE
632: LTOP = KTOP
633: END IF
634: DO 70 KROW = LTOP, KWTOP - 1, NV
635: KLN = MIN( NV, KWTOP-KROW )
636: CALL DGEMM( 'N', 'N', KLN, JW, JW, ONE, H( KROW, KWTOP ),
637: $ LDH, V, LDV, ZERO, WV, LDWV )
638: CALL DLACPY( 'A', KLN, JW, WV, LDWV, H( KROW, KWTOP ), LDH )
639: 70 CONTINUE
640: *
641: * ==== Update horizontal slab in H ====
642: *
643: IF( WANTT ) THEN
644: DO 80 KCOL = KBOT + 1, N, NH
645: KLN = MIN( NH, N-KCOL+1 )
646: CALL DGEMM( 'C', 'N', JW, KLN, JW, ONE, V, LDV,
647: $ H( KWTOP, KCOL ), LDH, ZERO, T, LDT )
648: CALL DLACPY( 'A', JW, KLN, T, LDT, H( KWTOP, KCOL ),
649: $ LDH )
650: 80 CONTINUE
651: END IF
652: *
653: * ==== Update vertical slab in Z ====
654: *
655: IF( WANTZ ) THEN
656: DO 90 KROW = ILOZ, IHIZ, NV
657: KLN = MIN( NV, IHIZ-KROW+1 )
658: CALL DGEMM( 'N', 'N', KLN, JW, JW, ONE, Z( KROW, KWTOP ),
659: $ LDZ, V, LDV, ZERO, WV, LDWV )
660: CALL DLACPY( 'A', KLN, JW, WV, LDWV, Z( KROW, KWTOP ),
661: $ LDZ )
662: 90 CONTINUE
663: END IF
664: END IF
665: *
666: * ==== Return the number of deflations ... ====
667: *
668: ND = JW - NS
669: *
670: * ==== ... and the number of shifts. (Subtracting
671: * . INFQR from the spike length takes care
672: * . of the case of a rare QR failure while
673: * . calculating eigenvalues of the deflation
674: * . window.) ====
675: *
676: NS = NS - INFQR
677: *
678: * ==== Return optimal workspace. ====
679: *
680: WORK( 1 ) = DBLE( LWKOPT )
681: *
682: * ==== End of DLAQR2 ====
683: *
684: 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