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