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