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