LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlaqr2.f
Go to the documentation of this file.
00001 *> \brief \b DLAQR2
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DLAQR2 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqr2.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqr2.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqr2.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DLAQR2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
00022 *                          IHIZ, Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T,
00023 *                          LDT, NV, WV, LDWV, WORK, LWORK )
00024 * 
00025 *       .. Scalar Arguments ..
00026 *       INTEGER            IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
00027 *      $                   LDZ, LWORK, N, ND, NH, NS, NV, NW
00028 *       LOGICAL            WANTT, WANTZ
00029 *       ..
00030 *       .. Array Arguments ..
00031 *       DOUBLE PRECISION   H( LDH, * ), SI( * ), SR( * ), T( LDT, * ),
00032 *      $                   V( LDV, * ), WORK( * ), WV( LDWV, * ),
00033 *      $                   Z( LDZ, * )
00034 *       ..
00035 *  
00036 *
00037 *> \par Purpose:
00038 *  =============
00039 *>
00040 *> \verbatim
00041 *>
00042 *>    DLAQR2 is identical to DLAQR3 except that it avoids
00043 *>    recursion by calling DLAHQR instead of DLAQR4.
00044 *>
00045 *>    Aggressive early deflation:
00046 *>
00047 *>    This subroutine accepts as input an upper Hessenberg matrix
00048 *>    H and performs an orthogonal similarity transformation
00049 *>    designed to detect and deflate fully converged eigenvalues from
00050 *>    a trailing principal submatrix.  On output H has been over-
00051 *>    written by a new Hessenberg matrix that is a perturbation of
00052 *>    an orthogonal similarity transformation of H.  It is to be
00053 *>    hoped that the final version of H has many zero subdiagonal
00054 *>    entries.
00055 *> \endverbatim
00056 *
00057 *  Arguments:
00058 *  ==========
00059 *
00060 *> \param[in] WANTT
00061 *> \verbatim
00062 *>          WANTT is LOGICAL
00063 *>          If .TRUE., then the Hessenberg matrix H is fully updated
00064 *>          so that the quasi-triangular Schur factor may be
00065 *>          computed (in cooperation with the calling subroutine).
00066 *>          If .FALSE., then only enough of H is updated to preserve
00067 *>          the eigenvalues.
00068 *> \endverbatim
00069 *>
00070 *> \param[in] WANTZ
00071 *> \verbatim
00072 *>          WANTZ is LOGICAL
00073 *>          If .TRUE., then the orthogonal matrix Z is updated so
00074 *>          so that the orthogonal Schur factor may be computed
00075 *>          (in cooperation with the calling subroutine).
00076 *>          If .FALSE., then Z is not referenced.
00077 *> \endverbatim
00078 *>
00079 *> \param[in] N
00080 *> \verbatim
00081 *>          N is INTEGER
00082 *>          The order of the matrix H and (if WANTZ is .TRUE.) the
00083 *>          order of the orthogonal matrix Z.
00084 *> \endverbatim
00085 *>
00086 *> \param[in] KTOP
00087 *> \verbatim
00088 *>          KTOP is INTEGER
00089 *>          It is assumed that either KTOP = 1 or H(KTOP,KTOP-1)=0.
00090 *>          KBOT and KTOP together determine an isolated block
00091 *>          along the diagonal of the Hessenberg matrix.
00092 *> \endverbatim
00093 *>
00094 *> \param[in] KBOT
00095 *> \verbatim
00096 *>          KBOT is INTEGER
00097 *>          It is assumed without a check that either
00098 *>          KBOT = N or H(KBOT+1,KBOT)=0.  KBOT and KTOP together
00099 *>          determine an isolated block along the diagonal of the
00100 *>          Hessenberg matrix.
00101 *> \endverbatim
00102 *>
00103 *> \param[in] NW
00104 *> \verbatim
00105 *>          NW is INTEGER
00106 *>          Deflation window size.  1 .LE. NW .LE. (KBOT-KTOP+1).
00107 *> \endverbatim
00108 *>
00109 *> \param[in,out] H
00110 *> \verbatim
00111 *>          H is DOUBLE PRECISION array, dimension (LDH,N)
00112 *>          On input the initial N-by-N section of H stores the
00113 *>          Hessenberg matrix undergoing aggressive early deflation.
00114 *>          On output H has been transformed by an orthogonal
00115 *>          similarity transformation, perturbed, and the returned
00116 *>          to Hessenberg form that (it is to be hoped) has some
00117 *>          zero subdiagonal entries.
00118 *> \endverbatim
00119 *>
00120 *> \param[in] LDH
00121 *> \verbatim
00122 *>          LDH is integer
00123 *>          Leading dimension of H just as declared in the calling
00124 *>          subroutine.  N .LE. LDH
00125 *> \endverbatim
00126 *>
00127 *> \param[in] ILOZ
00128 *> \verbatim
00129 *>          ILOZ is INTEGER
00130 *> \endverbatim
00131 *>
00132 *> \param[in] IHIZ
00133 *> \verbatim
00134 *>          IHIZ is INTEGER
00135 *>          Specify the rows of Z to which transformations must be
00136 *>          applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N.
00137 *> \endverbatim
00138 *>
00139 *> \param[in,out] Z
00140 *> \verbatim
00141 *>          Z is DOUBLE PRECISION array, dimension (LDZ,N)
00142 *>          IF WANTZ is .TRUE., then on output, the orthogonal
00143 *>          similarity transformation mentioned above has been
00144 *>          accumulated into Z(ILOZ:IHIZ,ILO:IHI) from the right.
00145 *>          If WANTZ is .FALSE., then Z is unreferenced.
00146 *> \endverbatim
00147 *>
00148 *> \param[in] LDZ
00149 *> \verbatim
00150 *>          LDZ is integer
00151 *>          The leading dimension of Z just as declared in the
00152 *>          calling subroutine.  1 .LE. LDZ.
00153 *> \endverbatim
00154 *>
00155 *> \param[out] NS
00156 *> \verbatim
00157 *>          NS is integer
00158 *>          The number of unconverged (ie approximate) eigenvalues
00159 *>          returned in SR and SI that may be used as shifts by the
00160 *>          calling subroutine.
00161 *> \endverbatim
00162 *>
00163 *> \param[out] ND
00164 *> \verbatim
00165 *>          ND is integer
00166 *>          The number of converged eigenvalues uncovered by this
00167 *>          subroutine.
00168 *> \endverbatim
00169 *>
00170 *> \param[out] SR
00171 *> \verbatim
00172 *>          SR is DOUBLE PRECISION array, dimension (KBOT)
00173 *> \endverbatim
00174 *>
00175 *> \param[out] SI
00176 *> \verbatim
00177 *>          SI is DOUBLE PRECISION array, dimension (KBOT)
00178 *>          On output, the real and imaginary parts of approximate
00179 *>          eigenvalues that may be used for shifts are stored in
00180 *>          SR(KBOT-ND-NS+1) through SR(KBOT-ND) and
00181 *>          SI(KBOT-ND-NS+1) through SI(KBOT-ND), respectively.
00182 *>          The real and imaginary parts of converged eigenvalues
00183 *>          are stored in SR(KBOT-ND+1) through SR(KBOT) and
00184 *>          SI(KBOT-ND+1) through SI(KBOT), respectively.
00185 *> \endverbatim
00186 *>
00187 *> \param[out] V
00188 *> \verbatim
00189 *>          V is DOUBLE PRECISION array, dimension (LDV,NW)
00190 *>          An NW-by-NW work array.
00191 *> \endverbatim
00192 *>
00193 *> \param[in] LDV
00194 *> \verbatim
00195 *>          LDV is integer scalar
00196 *>          The leading dimension of V just as declared in the
00197 *>          calling subroutine.  NW .LE. LDV
00198 *> \endverbatim
00199 *>
00200 *> \param[in] NH
00201 *> \verbatim
00202 *>          NH is integer scalar
00203 *>          The number of columns of T.  NH.GE.NW.
00204 *> \endverbatim
00205 *>
00206 *> \param[out] T
00207 *> \verbatim
00208 *>          T is DOUBLE PRECISION array, dimension (LDT,NW)
00209 *> \endverbatim
00210 *>
00211 *> \param[in] LDT
00212 *> \verbatim
00213 *>          LDT is integer
00214 *>          The leading dimension of T just as declared in the
00215 *>          calling subroutine.  NW .LE. LDT
00216 *> \endverbatim
00217 *>
00218 *> \param[in] NV
00219 *> \verbatim
00220 *>          NV is integer
00221 *>          The number of rows of work array WV available for
00222 *>          workspace.  NV.GE.NW.
00223 *> \endverbatim
00224 *>
00225 *> \param[out] WV
00226 *> \verbatim
00227 *>          WV is DOUBLE PRECISION array, dimension (LDWV,NW)
00228 *> \endverbatim
00229 *>
00230 *> \param[in] LDWV
00231 *> \verbatim
00232 *>          LDWV is integer
00233 *>          The leading dimension of W just as declared in the
00234 *>          calling subroutine.  NW .LE. LDV
00235 *> \endverbatim
00236 *>
00237 *> \param[out] WORK
00238 *> \verbatim
00239 *>          WORK is DOUBLE PRECISION array, dimension (LWORK)
00240 *>          On exit, WORK(1) is set to an estimate of the optimal value
00241 *>          of LWORK for the given values of N, NW, KTOP and KBOT.
00242 *> \endverbatim
00243 *>
00244 *> \param[in] LWORK
00245 *> \verbatim
00246 *>          LWORK is integer
00247 *>          The dimension of the work array WORK.  LWORK = 2*NW
00248 *>          suffices, but greater efficiency may result from larger
00249 *>          values of LWORK.
00250 *>
00251 *>          If LWORK = -1, then a workspace query is assumed; DLAQR2
00252 *>          only estimates the optimal workspace size for the given
00253 *>          values of N, NW, KTOP and KBOT.  The estimate is returned
00254 *>          in WORK(1).  No error message related to LWORK is issued
00255 *>          by XERBLA.  Neither H nor Z are accessed.
00256 *> \endverbatim
00257 *
00258 *  Authors:
00259 *  ========
00260 *
00261 *> \author Univ. of Tennessee 
00262 *> \author Univ. of California Berkeley 
00263 *> \author Univ. of Colorado Denver 
00264 *> \author NAG Ltd. 
00265 *
00266 *> \date November 2011
00267 *
00268 *> \ingroup doubleOTHERauxiliary
00269 *
00270 *> \par Contributors:
00271 *  ==================
00272 *>
00273 *>       Karen Braman and Ralph Byers, Department of Mathematics,
00274 *>       University of Kansas, USA
00275 *>
00276 *  =====================================================================
00277       SUBROUTINE DLAQR2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
00278      $                   IHIZ, Z, LDZ, NS, ND, SR, SI, V, LDV, NH, T,
00279      $                   LDT, NV, WV, LDWV, WORK, LWORK )
00280 *
00281 *  -- LAPACK auxiliary routine (version 3.4.0) --
00282 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00283 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00284 *     November 2011
00285 *
00286 *     .. Scalar Arguments ..
00287       INTEGER            IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
00288      $                   LDZ, LWORK, N, ND, NH, NS, NV, NW
00289       LOGICAL            WANTT, WANTZ
00290 *     ..
00291 *     .. Array Arguments ..
00292       DOUBLE PRECISION   H( LDH, * ), SI( * ), SR( * ), T( LDT, * ),
00293      $                   V( LDV, * ), WORK( * ), WV( LDWV, * ),
00294      $                   Z( LDZ, * )
00295 *     ..
00296 *
00297 *  ================================================================
00298 *     .. Parameters ..
00299       DOUBLE PRECISION   ZERO, ONE
00300       PARAMETER          ( ZERO = 0.0d0, ONE = 1.0d0 )
00301 *     ..
00302 *     .. Local Scalars ..
00303       DOUBLE PRECISION   AA, BB, BETA, CC, CS, DD, EVI, EVK, FOO, S,
00304      $                   SAFMAX, SAFMIN, SMLNUM, SN, TAU, ULP
00305       INTEGER            I, IFST, ILST, INFO, INFQR, J, JW, K, KCOL,
00306      $                   KEND, KLN, KROW, KWTOP, LTOP, LWK1, LWK2,
00307      $                   LWKOPT
00308       LOGICAL            BULGE, SORTED
00309 *     ..
00310 *     .. External Functions ..
00311       DOUBLE PRECISION   DLAMCH
00312       EXTERNAL           DLAMCH
00313 *     ..
00314 *     .. External Subroutines ..
00315       EXTERNAL           DCOPY, DGEHRD, DGEMM, DLABAD, DLACPY, DLAHQR,
00316      $                   DLANV2, DLARF, DLARFG, DLASET, DORMHR, DTREXC
00317 *     ..
00318 *     .. Intrinsic Functions ..
00319       INTRINSIC          ABS, DBLE, INT, MAX, MIN, SQRT
00320 *     ..
00321 *     .. Executable Statements ..
00322 *
00323 *     ==== Estimate optimal workspace. ====
00324 *
00325       JW = MIN( NW, KBOT-KTOP+1 )
00326       IF( JW.LE.2 ) THEN
00327          LWKOPT = 1
00328       ELSE
00329 *
00330 *        ==== Workspace query call to DGEHRD ====
00331 *
00332          CALL DGEHRD( JW, 1, JW-1, T, LDT, WORK, WORK, -1, INFO )
00333          LWK1 = INT( WORK( 1 ) )
00334 *
00335 *        ==== Workspace query call to DORMHR ====
00336 *
00337          CALL DORMHR( 'R', 'N', JW, JW, 1, JW-1, T, LDT, WORK, V, LDV,
00338      $                WORK, -1, INFO )
00339          LWK2 = INT( WORK( 1 ) )
00340 *
00341 *        ==== Optimal workspace ====
00342 *
00343          LWKOPT = JW + MAX( LWK1, LWK2 )
00344       END IF
00345 *
00346 *     ==== Quick return in case of workspace query. ====
00347 *
00348       IF( LWORK.EQ.-1 ) THEN
00349          WORK( 1 ) = DBLE( LWKOPT )
00350          RETURN
00351       END IF
00352 *
00353 *     ==== Nothing to do ...
00354 *     ... for an empty active block ... ====
00355       NS = 0
00356       ND = 0
00357       WORK( 1 ) = ONE
00358       IF( KTOP.GT.KBOT )
00359      $   RETURN
00360 *     ... nor for an empty deflation window. ====
00361       IF( NW.LT.1 )
00362      $   RETURN
00363 *
00364 *     ==== Machine constants ====
00365 *
00366       SAFMIN = DLAMCH( 'SAFE MINIMUM' )
00367       SAFMAX = ONE / SAFMIN
00368       CALL DLABAD( SAFMIN, SAFMAX )
00369       ULP = DLAMCH( 'PRECISION' )
00370       SMLNUM = SAFMIN*( DBLE( N ) / ULP )
00371 *
00372 *     ==== Setup deflation window ====
00373 *
00374       JW = MIN( NW, KBOT-KTOP+1 )
00375       KWTOP = KBOT - JW + 1
00376       IF( KWTOP.EQ.KTOP ) THEN
00377          S = ZERO
00378       ELSE
00379          S = H( KWTOP, KWTOP-1 )
00380       END IF
00381 *
00382       IF( KBOT.EQ.KWTOP ) THEN
00383 *
00384 *        ==== 1-by-1 deflation window: not much to do ====
00385 *
00386          SR( KWTOP ) = H( KWTOP, KWTOP )
00387          SI( KWTOP ) = ZERO
00388          NS = 1
00389          ND = 0
00390          IF( ABS( S ).LE.MAX( SMLNUM, ULP*ABS( H( KWTOP, KWTOP ) ) ) )
00391      $        THEN
00392             NS = 0
00393             ND = 1
00394             IF( KWTOP.GT.KTOP )
00395      $         H( KWTOP, KWTOP-1 ) = ZERO
00396          END IF
00397          WORK( 1 ) = ONE
00398          RETURN
00399       END IF
00400 *
00401 *     ==== Convert to spike-triangular form.  (In case of a
00402 *     .    rare QR failure, this routine continues to do
00403 *     .    aggressive early deflation using that part of
00404 *     .    the deflation window that converged using INFQR
00405 *     .    here and there to keep track.) ====
00406 *
00407       CALL DLACPY( 'U', JW, JW, H( KWTOP, KWTOP ), LDH, T, LDT )
00408       CALL DCOPY( JW-1, H( KWTOP+1, KWTOP ), LDH+1, T( 2, 1 ), LDT+1 )
00409 *
00410       CALL DLASET( 'A', JW, JW, ZERO, ONE, V, LDV )
00411       CALL DLAHQR( .true., .true., JW, 1, JW, T, LDT, SR( KWTOP ),
00412      $             SI( KWTOP ), 1, JW, V, LDV, INFQR )
00413 *
00414 *     ==== DTREXC needs a clean margin near the diagonal ====
00415 *
00416       DO 10 J = 1, JW - 3
00417          T( J+2, J ) = ZERO
00418          T( J+3, J ) = ZERO
00419    10 CONTINUE
00420       IF( JW.GT.2 )
00421      $   T( JW, JW-2 ) = ZERO
00422 *
00423 *     ==== Deflation detection loop ====
00424 *
00425       NS = JW
00426       ILST = INFQR + 1
00427    20 CONTINUE
00428       IF( ILST.LE.NS ) THEN
00429          IF( NS.EQ.1 ) THEN
00430             BULGE = .FALSE.
00431          ELSE
00432             BULGE = T( NS, NS-1 ).NE.ZERO
00433          END IF
00434 *
00435 *        ==== Small spike tip test for deflation ====
00436 *
00437          IF( .NOT.BULGE ) THEN
00438 *
00439 *           ==== Real eigenvalue ====
00440 *
00441             FOO = ABS( T( NS, NS ) )
00442             IF( FOO.EQ.ZERO )
00443      $         FOO = ABS( S )
00444             IF( ABS( S*V( 1, NS ) ).LE.MAX( SMLNUM, ULP*FOO ) ) THEN
00445 *
00446 *              ==== Deflatable ====
00447 *
00448                NS = NS - 1
00449             ELSE
00450 *
00451 *              ==== Undeflatable.   Move it up out of the way.
00452 *              .    (DTREXC can not fail in this case.) ====
00453 *
00454                IFST = NS
00455                CALL DTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, WORK,
00456      $                      INFO )
00457                ILST = ILST + 1
00458             END IF
00459          ELSE
00460 *
00461 *           ==== Complex conjugate pair ====
00462 *
00463             FOO = ABS( T( NS, NS ) ) + SQRT( ABS( T( NS, NS-1 ) ) )*
00464      $            SQRT( ABS( T( NS-1, NS ) ) )
00465             IF( FOO.EQ.ZERO )
00466      $         FOO = ABS( S )
00467             IF( MAX( ABS( S*V( 1, NS ) ), ABS( S*V( 1, NS-1 ) ) ).LE.
00468      $          MAX( SMLNUM, ULP*FOO ) ) THEN
00469 *
00470 *              ==== Deflatable ====
00471 *
00472                NS = NS - 2
00473             ELSE
00474 *
00475 *              ==== Undeflatable. Move them up out of the way.
00476 *              .    Fortunately, DTREXC does the right thing with
00477 *              .    ILST in case of a rare exchange failure. ====
00478 *
00479                IFST = NS
00480                CALL DTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, WORK,
00481      $                      INFO )
00482                ILST = ILST + 2
00483             END IF
00484          END IF
00485 *
00486 *        ==== End deflation detection loop ====
00487 *
00488          GO TO 20
00489       END IF
00490 *
00491 *        ==== Return to Hessenberg form ====
00492 *
00493       IF( NS.EQ.0 )
00494      $   S = ZERO
00495 *
00496       IF( NS.LT.JW ) THEN
00497 *
00498 *        ==== sorting diagonal blocks of T improves accuracy for
00499 *        .    graded matrices.  Bubble sort deals well with
00500 *        .    exchange failures. ====
00501 *
00502          SORTED = .false.
00503          I = NS + 1
00504    30    CONTINUE
00505          IF( SORTED )
00506      $      GO TO 50
00507          SORTED = .true.
00508 *
00509          KEND = I - 1
00510          I = INFQR + 1
00511          IF( I.EQ.NS ) THEN
00512             K = I + 1
00513          ELSE IF( T( I+1, I ).EQ.ZERO ) THEN
00514             K = I + 1
00515          ELSE
00516             K = I + 2
00517          END IF
00518    40    CONTINUE
00519          IF( K.LE.KEND ) THEN
00520             IF( K.EQ.I+1 ) THEN
00521                EVI = ABS( T( I, I ) )
00522             ELSE
00523                EVI = ABS( T( I, I ) ) + SQRT( ABS( T( I+1, I ) ) )*
00524      $               SQRT( ABS( T( I, I+1 ) ) )
00525             END IF
00526 *
00527             IF( K.EQ.KEND ) THEN
00528                EVK = ABS( T( K, K ) )
00529             ELSE IF( T( K+1, K ).EQ.ZERO ) THEN
00530                EVK = ABS( T( K, K ) )
00531             ELSE
00532                EVK = ABS( T( K, K ) ) + SQRT( ABS( T( K+1, K ) ) )*
00533      $               SQRT( ABS( T( K, K+1 ) ) )
00534             END IF
00535 *
00536             IF( EVI.GE.EVK ) THEN
00537                I = K
00538             ELSE
00539                SORTED = .false.
00540                IFST = I
00541                ILST = K
00542                CALL DTREXC( 'V', JW, T, LDT, V, LDV, IFST, ILST, WORK,
00543      $                      INFO )
00544                IF( INFO.EQ.0 ) THEN
00545                   I = ILST
00546                ELSE
00547                   I = K
00548                END IF
00549             END IF
00550             IF( I.EQ.KEND ) THEN
00551                K = I + 1
00552             ELSE IF( T( I+1, I ).EQ.ZERO ) THEN
00553                K = I + 1
00554             ELSE
00555                K = I + 2
00556             END IF
00557             GO TO 40
00558          END IF
00559          GO TO 30
00560    50    CONTINUE
00561       END IF
00562 *
00563 *     ==== Restore shift/eigenvalue array from T ====
00564 *
00565       I = JW
00566    60 CONTINUE
00567       IF( I.GE.INFQR+1 ) THEN
00568          IF( I.EQ.INFQR+1 ) THEN
00569             SR( KWTOP+I-1 ) = T( I, I )
00570             SI( KWTOP+I-1 ) = ZERO
00571             I = I - 1
00572          ELSE IF( T( I, I-1 ).EQ.ZERO ) THEN
00573             SR( KWTOP+I-1 ) = T( I, I )
00574             SI( KWTOP+I-1 ) = ZERO
00575             I = I - 1
00576          ELSE
00577             AA = T( I-1, I-1 )
00578             CC = T( I, I-1 )
00579             BB = T( I-1, I )
00580             DD = T( I, I )
00581             CALL DLANV2( AA, BB, CC, DD, SR( KWTOP+I-2 ),
00582      $                   SI( KWTOP+I-2 ), SR( KWTOP+I-1 ),
00583      $                   SI( KWTOP+I-1 ), CS, SN )
00584             I = I - 2
00585          END IF
00586          GO TO 60
00587       END IF
00588 *
00589       IF( NS.LT.JW .OR. S.EQ.ZERO ) THEN
00590          IF( NS.GT.1 .AND. S.NE.ZERO ) THEN
00591 *
00592 *           ==== Reflect spike back into lower triangle ====
00593 *
00594             CALL DCOPY( NS, V, LDV, WORK, 1 )
00595             BETA = WORK( 1 )
00596             CALL DLARFG( NS, BETA, WORK( 2 ), 1, TAU )
00597             WORK( 1 ) = ONE
00598 *
00599             CALL DLASET( 'L', JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), LDT )
00600 *
00601             CALL DLARF( 'L', NS, JW, WORK, 1, TAU, T, LDT,
00602      $                  WORK( JW+1 ) )
00603             CALL DLARF( 'R', NS, NS, WORK, 1, TAU, T, LDT,
00604      $                  WORK( JW+1 ) )
00605             CALL DLARF( 'R', JW, NS, WORK, 1, TAU, V, LDV,
00606      $                  WORK( JW+1 ) )
00607 *
00608             CALL DGEHRD( JW, 1, NS, T, LDT, WORK, WORK( JW+1 ),
00609      $                   LWORK-JW, INFO )
00610          END IF
00611 *
00612 *        ==== Copy updated reduced window into place ====
00613 *
00614          IF( KWTOP.GT.1 )
00615      $      H( KWTOP, KWTOP-1 ) = S*V( 1, 1 )
00616          CALL DLACPY( 'U', JW, JW, T, LDT, H( KWTOP, KWTOP ), LDH )
00617          CALL DCOPY( JW-1, T( 2, 1 ), LDT+1, H( KWTOP+1, KWTOP ),
00618      $               LDH+1 )
00619 *
00620 *        ==== Accumulate orthogonal matrix in order update
00621 *        .    H and Z, if requested.  ====
00622 *
00623          IF( NS.GT.1 .AND. S.NE.ZERO )
00624      $      CALL DORMHR( 'R', 'N', JW, NS, 1, NS, T, LDT, WORK, V, LDV,
00625      $                   WORK( JW+1 ), LWORK-JW, INFO )
00626 *
00627 *        ==== Update vertical slab in H ====
00628 *
00629          IF( WANTT ) THEN
00630             LTOP = 1
00631          ELSE
00632             LTOP = KTOP
00633          END IF
00634          DO 70 KROW = LTOP, KWTOP - 1, NV
00635             KLN = MIN( NV, KWTOP-KROW )
00636             CALL DGEMM( 'N', 'N', KLN, JW, JW, ONE, H( KROW, KWTOP ),
00637      $                  LDH, V, LDV, ZERO, WV, LDWV )
00638             CALL DLACPY( 'A', KLN, JW, WV, LDWV, H( KROW, KWTOP ), LDH )
00639    70    CONTINUE
00640 *
00641 *        ==== Update horizontal slab in H ====
00642 *
00643          IF( WANTT ) THEN
00644             DO 80 KCOL = KBOT + 1, N, NH
00645                KLN = MIN( NH, N-KCOL+1 )
00646                CALL DGEMM( 'C', 'N', JW, KLN, JW, ONE, V, LDV,
00647      $                     H( KWTOP, KCOL ), LDH, ZERO, T, LDT )
00648                CALL DLACPY( 'A', JW, KLN, T, LDT, H( KWTOP, KCOL ),
00649      $                      LDH )
00650    80       CONTINUE
00651          END IF
00652 *
00653 *        ==== Update vertical slab in Z ====
00654 *
00655          IF( WANTZ ) THEN
00656             DO 90 KROW = ILOZ, IHIZ, NV
00657                KLN = MIN( NV, IHIZ-KROW+1 )
00658                CALL DGEMM( 'N', 'N', KLN, JW, JW, ONE, Z( KROW, KWTOP ),
00659      $                     LDZ, V, LDV, ZERO, WV, LDWV )
00660                CALL DLACPY( 'A', KLN, JW, WV, LDWV, Z( KROW, KWTOP ),
00661      $                      LDZ )
00662    90       CONTINUE
00663          END IF
00664       END IF
00665 *
00666 *     ==== Return the number of deflations ... ====
00667 *
00668       ND = JW - NS
00669 *
00670 *     ==== ... and the number of shifts. (Subtracting
00671 *     .    INFQR from the spike length takes care
00672 *     .    of the case of a rare QR failure while
00673 *     .    calculating eigenvalues of the deflation
00674 *     .    window.)  ====
00675 *
00676       NS = NS - INFQR
00677 *
00678 *      ==== Return optimal workspace. ====
00679 *
00680       WORK( 1 ) = DBLE( LWKOPT )
00681 *
00682 *     ==== End of DLAQR2 ====
00683 *
00684       END
 All Files Functions