![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SLAQR2 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SLAQR2 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaqr2.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaqr2.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaqr2.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SLAQR2( 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 *> SLAQR2 is identical to SLAQR3 except that it avoids 00043 *> recursion by calling SLAHQR instead of SLAQR4. 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 REAL 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 REAL 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 REAL array, dimension KBOT 00173 *> \endverbatim 00174 *> 00175 *> \param[out] SI 00176 *> \verbatim 00177 *> SI is REAL 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 REAL 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 REAL 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 REAL 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 REAL 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; SLAQR2 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 realOTHERauxiliary 00269 * 00270 *> \par Contributors: 00271 * ================== 00272 *> 00273 *> Karen Braman and Ralph Byers, Department of Mathematics, 00274 *> University of Kansas, USA 00275 *> 00276 * ===================================================================== 00277 SUBROUTINE SLAQR2( 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 REAL H( LDH, * ), SI( * ), SR( * ), T( LDT, * ), 00293 $ V( LDV, * ), WORK( * ), WV( LDWV, * ), 00294 $ Z( LDZ, * ) 00295 * .. 00296 * 00297 * ================================================================ 00298 * .. Parameters .. 00299 REAL ZERO, ONE 00300 PARAMETER ( ZERO = 0.0e0, ONE = 1.0e0 ) 00301 * .. 00302 * .. Local Scalars .. 00303 REAL 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 REAL SLAMCH 00312 EXTERNAL SLAMCH 00313 * .. 00314 * .. External Subroutines .. 00315 EXTERNAL SCOPY, SGEHRD, SGEMM, SLABAD, SLACPY, SLAHQR, 00316 $ SLANV2, SLARF, SLARFG, SLASET, SORMHR, STREXC 00317 * .. 00318 * .. Intrinsic Functions .. 00319 INTRINSIC ABS, INT, MAX, MIN, REAL, 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 SGEHRD ==== 00331 * 00332 CALL SGEHRD( JW, 1, JW-1, T, LDT, WORK, WORK, -1, INFO ) 00333 LWK1 = INT( WORK( 1 ) ) 00334 * 00335 * ==== Workspace query call to SORMHR ==== 00336 * 00337 CALL SORMHR( '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 ) = REAL( 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 = SLAMCH( 'SAFE MINIMUM' ) 00367 SAFMAX = ONE / SAFMIN 00368 CALL SLABAD( SAFMIN, SAFMAX ) 00369 ULP = SLAMCH( 'PRECISION' ) 00370 SMLNUM = SAFMIN*( REAL( 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 SLACPY( 'U', JW, JW, H( KWTOP, KWTOP ), LDH, T, LDT ) 00408 CALL SCOPY( JW-1, H( KWTOP+1, KWTOP ), LDH+1, T( 2, 1 ), LDT+1 ) 00409 * 00410 CALL SLASET( 'A', JW, JW, ZERO, ONE, V, LDV ) 00411 CALL SLAHQR( .true., .true., JW, 1, JW, T, LDT, SR( KWTOP ), 00412 $ SI( KWTOP ), 1, JW, V, LDV, INFQR ) 00413 * 00414 * ==== STREXC 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 * . (STREXC can not fail in this case.) ==== 00453 * 00454 IFST = NS 00455 CALL STREXC( '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, STREXC does the right thing with 00477 * . ILST in case of a rare exchange failure. ==== 00478 * 00479 IFST = NS 00480 CALL STREXC( '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 STREXC( '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 SLANV2( 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 SCOPY( NS, V, LDV, WORK, 1 ) 00595 BETA = WORK( 1 ) 00596 CALL SLARFG( NS, BETA, WORK( 2 ), 1, TAU ) 00597 WORK( 1 ) = ONE 00598 * 00599 CALL SLASET( 'L', JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), LDT ) 00600 * 00601 CALL SLARF( 'L', NS, JW, WORK, 1, TAU, T, LDT, 00602 $ WORK( JW+1 ) ) 00603 CALL SLARF( 'R', NS, NS, WORK, 1, TAU, T, LDT, 00604 $ WORK( JW+1 ) ) 00605 CALL SLARF( 'R', JW, NS, WORK, 1, TAU, V, LDV, 00606 $ WORK( JW+1 ) ) 00607 * 00608 CALL SGEHRD( 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 SLACPY( 'U', JW, JW, T, LDT, H( KWTOP, KWTOP ), LDH ) 00617 CALL SCOPY( 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 SORMHR( '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 SGEMM( 'N', 'N', KLN, JW, JW, ONE, H( KROW, KWTOP ), 00637 $ LDH, V, LDV, ZERO, WV, LDWV ) 00638 CALL SLACPY( '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 SGEMM( 'C', 'N', JW, KLN, JW, ONE, V, LDV, 00647 $ H( KWTOP, KCOL ), LDH, ZERO, T, LDT ) 00648 CALL SLACPY( '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 SGEMM( 'N', 'N', KLN, JW, JW, ONE, Z( KROW, KWTOP ), 00659 $ LDZ, V, LDV, ZERO, WV, LDWV ) 00660 CALL SLACPY( '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 ) = REAL( LWKOPT ) 00681 * 00682 * ==== End of SLAQR2 ==== 00683 * 00684 END