![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLAQR3 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DLAQR3 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqr3.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqr3.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqr3.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DLAQR3( 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 *> Aggressive early deflation: 00043 *> 00044 *> DLAQR3 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (KBOT) 00170 *> \endverbatim 00171 *> 00172 *> \param[out] SI 00173 *> \verbatim 00174 *> SI is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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; DLAQR3 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 doubleOTHERauxiliary 00266 * 00267 *> \par Contributors: 00268 * ================== 00269 *> 00270 *> Karen Braman and Ralph Byers, Department of Mathematics, 00271 *> University of Kansas, USA 00272 *> 00273 * ===================================================================== 00274 SUBROUTINE DLAQR3( 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 DOUBLE PRECISION H( LDH, * ), SI( * ), SR( * ), T( LDT, * ), 00290 $ V( LDV, * ), WORK( * ), WV( LDWV, * ), 00291 $ Z( LDZ, * ) 00292 * .. 00293 * 00294 * ================================================================ 00295 * .. Parameters .. 00296 DOUBLE PRECISION ZERO, ONE 00297 PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0 ) 00298 * .. 00299 * .. Local Scalars .. 00300 DOUBLE PRECISION 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 DOUBLE PRECISION DLAMCH 00309 INTEGER ILAENV 00310 EXTERNAL DLAMCH, ILAENV 00311 * .. 00312 * .. External Subroutines .. 00313 EXTERNAL DCOPY, DGEHRD, DGEMM, DLABAD, DLACPY, DLAHQR, 00314 $ DLANV2, DLAQR4, DLARF, DLARFG, DLASET, DORMHR, 00315 $ DTREXC 00316 * .. 00317 * .. Intrinsic Functions .. 00318 INTRINSIC ABS, DBLE, INT, MAX, MIN, 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 DGEHRD ==== 00330 * 00331 CALL DGEHRD( JW, 1, JW-1, T, LDT, WORK, WORK, -1, INFO ) 00332 LWK1 = INT( WORK( 1 ) ) 00333 * 00334 * ==== Workspace query call to DORMHR ==== 00335 * 00336 CALL DORMHR( '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 DLAQR4 ==== 00341 * 00342 CALL DLAQR4( .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 ) = DBLE( 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 = DLAMCH( 'SAFE MINIMUM' ) 00372 SAFMAX = ONE / SAFMIN 00373 CALL DLABAD( SAFMIN, SAFMAX ) 00374 ULP = DLAMCH( 'PRECISION' ) 00375 SMLNUM = SAFMIN*( DBLE( 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 DLACPY( 'U', JW, JW, H( KWTOP, KWTOP ), LDH, T, LDT ) 00413 CALL DCOPY( JW-1, H( KWTOP+1, KWTOP ), LDH+1, T( 2, 1 ), LDT+1 ) 00414 * 00415 CALL DLASET( 'A', JW, JW, ZERO, ONE, V, LDV ) 00416 NMIN = ILAENV( 12, 'DLAQR3', 'SV', JW, 1, JW, LWORK ) 00417 IF( JW.GT.NMIN ) THEN 00418 CALL DLAQR4( .true., .true., JW, 1, JW, T, LDT, SR( KWTOP ), 00419 $ SI( KWTOP ), 1, JW, V, LDV, WORK, LWORK, INFQR ) 00420 ELSE 00421 CALL DLAHQR( .true., .true., JW, 1, JW, T, LDT, SR( KWTOP ), 00422 $ SI( KWTOP ), 1, JW, V, LDV, INFQR ) 00423 END IF 00424 * 00425 * ==== DTREXC 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 * . (DTREXC can not fail in this case.) ==== 00464 * 00465 IFST = NS 00466 CALL DTREXC( '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, DTREXC does the right thing with 00488 * . ILST in case of a rare exchange failure. ==== 00489 * 00490 IFST = NS 00491 CALL DTREXC( '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 DTREXC( '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 DLANV2( 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 DCOPY( NS, V, LDV, WORK, 1 ) 00606 BETA = WORK( 1 ) 00607 CALL DLARFG( NS, BETA, WORK( 2 ), 1, TAU ) 00608 WORK( 1 ) = ONE 00609 * 00610 CALL DLASET( 'L', JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), LDT ) 00611 * 00612 CALL DLARF( 'L', NS, JW, WORK, 1, TAU, T, LDT, 00613 $ WORK( JW+1 ) ) 00614 CALL DLARF( 'R', NS, NS, WORK, 1, TAU, T, LDT, 00615 $ WORK( JW+1 ) ) 00616 CALL DLARF( 'R', JW, NS, WORK, 1, TAU, V, LDV, 00617 $ WORK( JW+1 ) ) 00618 * 00619 CALL DGEHRD( 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 DLACPY( 'U', JW, JW, T, LDT, H( KWTOP, KWTOP ), LDH ) 00628 CALL DCOPY( 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 DORMHR( '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 DGEMM( 'N', 'N', KLN, JW, JW, ONE, H( KROW, KWTOP ), 00648 $ LDH, V, LDV, ZERO, WV, LDWV ) 00649 CALL DLACPY( '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 DGEMM( 'C', 'N', JW, KLN, JW, ONE, V, LDV, 00658 $ H( KWTOP, KCOL ), LDH, ZERO, T, LDT ) 00659 CALL DLACPY( '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 DGEMM( 'N', 'N', KLN, JW, JW, ONE, Z( KROW, KWTOP ), 00670 $ LDZ, V, LDV, ZERO, WV, LDWV ) 00671 CALL DLACPY( '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 ) = DBLE( LWKOPT ) 00692 * 00693 * ==== End of DLAQR3 ==== 00694 * 00695 END