![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CLAQR4 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CLAQR4 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/claqr4.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/claqr4.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/claqr4.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CLAQR4( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, 00022 * IHIZ, Z, LDZ, WORK, LWORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N 00026 * LOGICAL WANTT, WANTZ 00027 * .. 00028 * .. Array Arguments .. 00029 * COMPLEX H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * ) 00030 * .. 00031 * 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> CLAQR4 implements one level of recursion for CLAQR0. 00040 *> It is a complete implementation of the small bulge multi-shift 00041 *> QR algorithm. It may be called by CLAQR0 and, for large enough 00042 *> deflation window size, it may be called by CLAQR3. This 00043 *> subroutine is identical to CLAQR0 except that it calls CLAQR2 00044 *> instead of CLAQR3. 00045 *> 00046 *> CLAQR4 computes the eigenvalues of a Hessenberg matrix H 00047 *> and, optionally, the matrices T and Z from the Schur decomposition 00048 *> H = Z T Z**H, where T is an upper triangular matrix (the 00049 *> Schur form), and Z is the unitary matrix of Schur vectors. 00050 *> 00051 *> Optionally Z may be postmultiplied into an input unitary 00052 *> matrix Q so that this routine can give the Schur factorization 00053 *> of a matrix A which has been reduced to the Hessenberg form H 00054 *> by the unitary matrix Q: A = Q*H*Q**H = (QZ)*H*(QZ)**H. 00055 *> \endverbatim 00056 * 00057 * Arguments: 00058 * ========== 00059 * 00060 *> \param[in] WANTT 00061 *> \verbatim 00062 *> WANTT is LOGICAL 00063 *> = .TRUE. : the full Schur form T is required; 00064 *> = .FALSE.: only eigenvalues are required. 00065 *> \endverbatim 00066 *> 00067 *> \param[in] WANTZ 00068 *> \verbatim 00069 *> WANTZ is LOGICAL 00070 *> = .TRUE. : the matrix of Schur vectors Z is required; 00071 *> = .FALSE.: Schur vectors are not required. 00072 *> \endverbatim 00073 *> 00074 *> \param[in] N 00075 *> \verbatim 00076 *> N is INTEGER 00077 *> The order of the matrix H. N .GE. 0. 00078 *> \endverbatim 00079 *> 00080 *> \param[in] ILO 00081 *> \verbatim 00082 *> ILO is INTEGER 00083 *> \endverbatim 00084 *> 00085 *> \param[in] IHI 00086 *> \verbatim 00087 *> IHI is INTEGER 00088 *> It is assumed that H is already upper triangular in rows 00089 *> and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1, 00090 *> H(ILO,ILO-1) is zero. ILO and IHI are normally set by a 00091 *> previous call to CGEBAL, and then passed to CGEHRD when the 00092 *> matrix output by CGEBAL is reduced to Hessenberg form. 00093 *> Otherwise, ILO and IHI should be set to 1 and N, 00094 *> respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N. 00095 *> If N = 0, then ILO = 1 and IHI = 0. 00096 *> \endverbatim 00097 *> 00098 *> \param[in,out] H 00099 *> \verbatim 00100 *> H is COMPLEX array, dimension (LDH,N) 00101 *> On entry, the upper Hessenberg matrix H. 00102 *> On exit, if INFO = 0 and WANTT is .TRUE., then H 00103 *> contains the upper triangular matrix T from the Schur 00104 *> decomposition (the Schur form). If INFO = 0 and WANT is 00105 *> .FALSE., then the contents of H are unspecified on exit. 00106 *> (The output value of H when INFO.GT.0 is given under the 00107 *> description of INFO below.) 00108 *> 00109 *> This subroutine may explicitly set H(i,j) = 0 for i.GT.j and 00110 *> j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N. 00111 *> \endverbatim 00112 *> 00113 *> \param[in] LDH 00114 *> \verbatim 00115 *> LDH is INTEGER 00116 *> The leading dimension of the array H. LDH .GE. max(1,N). 00117 *> \endverbatim 00118 *> 00119 *> \param[out] W 00120 *> \verbatim 00121 *> W is COMPLEX array, dimension (N) 00122 *> The computed eigenvalues of H(ILO:IHI,ILO:IHI) are stored 00123 *> in W(ILO:IHI). If WANTT is .TRUE., then the eigenvalues are 00124 *> stored in the same order as on the diagonal of the Schur 00125 *> form returned in H, with W(i) = H(i,i). 00126 *> \endverbatim 00127 *> 00128 *> \param[in] ILOZ 00129 *> \verbatim 00130 *> ILOZ is INTEGER 00131 *> \endverbatim 00132 *> 00133 *> \param[in] IHIZ 00134 *> \verbatim 00135 *> IHIZ is INTEGER 00136 *> Specify the rows of Z to which transformations must be 00137 *> applied if WANTZ is .TRUE.. 00138 *> 1 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N. 00139 *> \endverbatim 00140 *> 00141 *> \param[in,out] Z 00142 *> \verbatim 00143 *> Z is COMPLEX array, dimension (LDZ,IHI) 00144 *> If WANTZ is .FALSE., then Z is not referenced. 00145 *> If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is 00146 *> replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the 00147 *> orthogonal Schur factor of H(ILO:IHI,ILO:IHI). 00148 *> (The output value of Z when INFO.GT.0 is given under 00149 *> the description of INFO below.) 00150 *> \endverbatim 00151 *> 00152 *> \param[in] LDZ 00153 *> \verbatim 00154 *> LDZ is INTEGER 00155 *> The leading dimension of the array Z. if WANTZ is .TRUE. 00156 *> then LDZ.GE.MAX(1,IHIZ). Otherwize, LDZ.GE.1. 00157 *> \endverbatim 00158 *> 00159 *> \param[out] WORK 00160 *> \verbatim 00161 *> WORK is COMPLEX array, dimension LWORK 00162 *> On exit, if LWORK = -1, WORK(1) returns an estimate of 00163 *> the optimal value for LWORK. 00164 *> \endverbatim 00165 *> 00166 *> \param[in] LWORK 00167 *> \verbatim 00168 *> LWORK is INTEGER 00169 *> The dimension of the array WORK. LWORK .GE. max(1,N) 00170 *> is sufficient, but LWORK typically as large as 6*N may 00171 *> be required for optimal performance. A workspace query 00172 *> to determine the optimal workspace size is recommended. 00173 *> 00174 *> If LWORK = -1, then CLAQR4 does a workspace query. 00175 *> In this case, CLAQR4 checks the input parameters and 00176 *> estimates the optimal workspace size for the given 00177 *> values of N, ILO and IHI. The estimate is returned 00178 *> in WORK(1). No error message related to LWORK is 00179 *> issued by XERBLA. Neither H nor Z are accessed. 00180 *> \endverbatim 00181 *> 00182 *> \param[out] INFO 00183 *> \verbatim 00184 *> \verbatim 00185 *> INFO is INTEGER 00186 *> = 0: successful exit 00187 *> .GT. 0: if INFO = i, CLAQR4 failed to compute all of 00188 *> the eigenvalues. Elements 1:ilo-1 and i+1:n of WR 00189 *> and WI contain those eigenvalues which have been 00190 *> successfully computed. (Failures are rare.) 00191 *> 00192 *> If INFO .GT. 0 and WANT is .FALSE., then on exit, 00193 *> the remaining unconverged eigenvalues are the eigen- 00194 *> values of the upper Hessenberg matrix rows and 00195 *> columns ILO through INFO of the final, output 00196 *> value of H. 00197 *> 00198 *> If INFO .GT. 0 and WANTT is .TRUE., then on exit 00199 *> 00200 *> (*) (initial value of H)*U = U*(final value of H) 00201 *> 00202 *> where U is a unitary matrix. The final 00203 *> value of H is upper Hessenberg and triangular in 00204 *> rows and columns INFO+1 through IHI. 00205 *> 00206 *> If INFO .GT. 0 and WANTZ is .TRUE., then on exit 00207 *> 00208 *> (final value of Z(ILO:IHI,ILOZ:IHIZ) 00209 *> = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U 00210 *> 00211 *> where U is the unitary matrix in (*) (regard- 00212 *> less of the value of WANTT.) 00213 *> 00214 *> If INFO .GT. 0 and WANTZ is .FALSE., then Z is not 00215 *> accessed. 00216 *> \endverbatim 00217 * 00218 * Authors: 00219 * ======== 00220 * 00221 *> \author Univ. of Tennessee 00222 *> \author Univ. of California Berkeley 00223 *> \author Univ. of Colorado Denver 00224 *> \author NAG Ltd. 00225 * 00226 *> \date November 2011 00227 * 00228 *> \ingroup complexOTHERauxiliary 00229 * 00230 *> \par Contributors: 00231 * ================== 00232 *> 00233 *> Karen Braman and Ralph Byers, Department of Mathematics, 00234 *> University of Kansas, USA 00235 * 00236 *> \par References: 00237 * ================ 00238 *> 00239 *> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR 00240 *> Algorithm Part I: Maintaining Well Focused Shifts, and Level 3 00241 *> Performance, SIAM Journal of Matrix Analysis, volume 23, pages 00242 *> 929--947, 2002. 00243 *> \n 00244 *> K. Braman, R. Byers and R. Mathias, The Multi-Shift QR 00245 *> Algorithm Part II: Aggressive Early Deflation, SIAM Journal 00246 *> of Matrix Analysis, volume 23, pages 948--973, 2002. 00247 *> 00248 * ===================================================================== 00249 SUBROUTINE CLAQR4( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, 00250 $ IHIZ, Z, LDZ, WORK, LWORK, INFO ) 00251 * 00252 * -- LAPACK auxiliary routine (version 3.4.0) -- 00253 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00254 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00255 * November 2011 00256 * 00257 * .. Scalar Arguments .. 00258 INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDH, LDZ, LWORK, N 00259 LOGICAL WANTT, WANTZ 00260 * .. 00261 * .. Array Arguments .. 00262 COMPLEX H( LDH, * ), W( * ), WORK( * ), Z( LDZ, * ) 00263 * .. 00264 * 00265 * 00266 * ================================================================ 00267 * 00268 * .. Parameters .. 00269 * 00270 * ==== Matrices of order NTINY or smaller must be processed by 00271 * . CLAHQR because of insufficient subdiagonal scratch space. 00272 * . (This is a hard limit.) ==== 00273 INTEGER NTINY 00274 PARAMETER ( NTINY = 11 ) 00275 * 00276 * ==== Exceptional deflation windows: try to cure rare 00277 * . slow convergence by varying the size of the 00278 * . deflation window after KEXNW iterations. ==== 00279 INTEGER KEXNW 00280 PARAMETER ( KEXNW = 5 ) 00281 * 00282 * ==== Exceptional shifts: try to cure rare slow convergence 00283 * . with ad-hoc exceptional shifts every KEXSH iterations. 00284 * . ==== 00285 INTEGER KEXSH 00286 PARAMETER ( KEXSH = 6 ) 00287 * 00288 * ==== The constant WILK1 is used to form the exceptional 00289 * . shifts. ==== 00290 REAL WILK1 00291 PARAMETER ( WILK1 = 0.75e0 ) 00292 COMPLEX ZERO, ONE 00293 PARAMETER ( ZERO = ( 0.0e0, 0.0e0 ), 00294 $ ONE = ( 1.0e0, 0.0e0 ) ) 00295 REAL TWO 00296 PARAMETER ( TWO = 2.0e0 ) 00297 * .. 00298 * .. Local Scalars .. 00299 COMPLEX AA, BB, CC, CDUM, DD, DET, RTDISC, SWAP, TR2 00300 REAL S 00301 INTEGER I, INF, IT, ITMAX, K, KACC22, KBOT, KDU, KS, 00302 $ KT, KTOP, KU, KV, KWH, KWTOP, KWV, LD, LS, 00303 $ LWKOPT, NDEC, NDFL, NH, NHO, NIBBLE, NMIN, NS, 00304 $ NSMAX, NSR, NVE, NW, NWMAX, NWR, NWUPBD 00305 LOGICAL SORTED 00306 CHARACTER JBCMPZ*2 00307 * .. 00308 * .. External Functions .. 00309 INTEGER ILAENV 00310 EXTERNAL ILAENV 00311 * .. 00312 * .. Local Arrays .. 00313 COMPLEX ZDUM( 1, 1 ) 00314 * .. 00315 * .. External Subroutines .. 00316 EXTERNAL CLACPY, CLAHQR, CLAQR2, CLAQR5 00317 * .. 00318 * .. Intrinsic Functions .. 00319 INTRINSIC ABS, AIMAG, CMPLX, INT, MAX, MIN, MOD, REAL, 00320 $ SQRT 00321 * .. 00322 * .. Statement Functions .. 00323 REAL CABS1 00324 * .. 00325 * .. Statement Function definitions .. 00326 CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) ) 00327 * .. 00328 * .. Executable Statements .. 00329 INFO = 0 00330 * 00331 * ==== Quick return for N = 0: nothing to do. ==== 00332 * 00333 IF( N.EQ.0 ) THEN 00334 WORK( 1 ) = ONE 00335 RETURN 00336 END IF 00337 * 00338 IF( N.LE.NTINY ) THEN 00339 * 00340 * ==== Tiny matrices must use CLAHQR. ==== 00341 * 00342 LWKOPT = 1 00343 IF( LWORK.NE.-1 ) 00344 $ CALL CLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, 00345 $ IHIZ, Z, LDZ, INFO ) 00346 ELSE 00347 * 00348 * ==== Use small bulge multi-shift QR with aggressive early 00349 * . deflation on larger-than-tiny matrices. ==== 00350 * 00351 * ==== Hope for the best. ==== 00352 * 00353 INFO = 0 00354 * 00355 * ==== Set up job flags for ILAENV. ==== 00356 * 00357 IF( WANTT ) THEN 00358 JBCMPZ( 1: 1 ) = 'S' 00359 ELSE 00360 JBCMPZ( 1: 1 ) = 'E' 00361 END IF 00362 IF( WANTZ ) THEN 00363 JBCMPZ( 2: 2 ) = 'V' 00364 ELSE 00365 JBCMPZ( 2: 2 ) = 'N' 00366 END IF 00367 * 00368 * ==== NWR = recommended deflation window size. At this 00369 * . point, N .GT. NTINY = 11, so there is enough 00370 * . subdiagonal workspace for NWR.GE.2 as required. 00371 * . (In fact, there is enough subdiagonal space for 00372 * . NWR.GE.3.) ==== 00373 * 00374 NWR = ILAENV( 13, 'CLAQR4', JBCMPZ, N, ILO, IHI, LWORK ) 00375 NWR = MAX( 2, NWR ) 00376 NWR = MIN( IHI-ILO+1, ( N-1 ) / 3, NWR ) 00377 * 00378 * ==== NSR = recommended number of simultaneous shifts. 00379 * . At this point N .GT. NTINY = 11, so there is at 00380 * . enough subdiagonal workspace for NSR to be even 00381 * . and greater than or equal to two as required. ==== 00382 * 00383 NSR = ILAENV( 15, 'CLAQR4', JBCMPZ, N, ILO, IHI, LWORK ) 00384 NSR = MIN( NSR, ( N+6 ) / 9, IHI-ILO ) 00385 NSR = MAX( 2, NSR-MOD( NSR, 2 ) ) 00386 * 00387 * ==== Estimate optimal workspace ==== 00388 * 00389 * ==== Workspace query call to CLAQR2 ==== 00390 * 00391 CALL CLAQR2( WANTT, WANTZ, N, ILO, IHI, NWR+1, H, LDH, ILOZ, 00392 $ IHIZ, Z, LDZ, LS, LD, W, H, LDH, N, H, LDH, N, H, 00393 $ LDH, WORK, -1 ) 00394 * 00395 * ==== Optimal workspace = MAX(CLAQR5, CLAQR2) ==== 00396 * 00397 LWKOPT = MAX( 3*NSR / 2, INT( WORK( 1 ) ) ) 00398 * 00399 * ==== Quick return in case of workspace query. ==== 00400 * 00401 IF( LWORK.EQ.-1 ) THEN 00402 WORK( 1 ) = CMPLX( LWKOPT, 0 ) 00403 RETURN 00404 END IF 00405 * 00406 * ==== CLAHQR/CLAQR0 crossover point ==== 00407 * 00408 NMIN = ILAENV( 12, 'CLAQR4', JBCMPZ, N, ILO, IHI, LWORK ) 00409 NMIN = MAX( NTINY, NMIN ) 00410 * 00411 * ==== Nibble crossover point ==== 00412 * 00413 NIBBLE = ILAENV( 14, 'CLAQR4', JBCMPZ, N, ILO, IHI, LWORK ) 00414 NIBBLE = MAX( 0, NIBBLE ) 00415 * 00416 * ==== Accumulate reflections during ttswp? Use block 00417 * . 2-by-2 structure during matrix-matrix multiply? ==== 00418 * 00419 KACC22 = ILAENV( 16, 'CLAQR4', JBCMPZ, N, ILO, IHI, LWORK ) 00420 KACC22 = MAX( 0, KACC22 ) 00421 KACC22 = MIN( 2, KACC22 ) 00422 * 00423 * ==== NWMAX = the largest possible deflation window for 00424 * . which there is sufficient workspace. ==== 00425 * 00426 NWMAX = MIN( ( N-1 ) / 3, LWORK / 2 ) 00427 NW = NWMAX 00428 * 00429 * ==== NSMAX = the Largest number of simultaneous shifts 00430 * . for which there is sufficient workspace. ==== 00431 * 00432 NSMAX = MIN( ( N+6 ) / 9, 2*LWORK / 3 ) 00433 NSMAX = NSMAX - MOD( NSMAX, 2 ) 00434 * 00435 * ==== NDFL: an iteration count restarted at deflation. ==== 00436 * 00437 NDFL = 1 00438 * 00439 * ==== ITMAX = iteration limit ==== 00440 * 00441 ITMAX = MAX( 30, 2*KEXSH )*MAX( 10, ( IHI-ILO+1 ) ) 00442 * 00443 * ==== Last row and column in the active block ==== 00444 * 00445 KBOT = IHI 00446 * 00447 * ==== Main Loop ==== 00448 * 00449 DO 70 IT = 1, ITMAX 00450 * 00451 * ==== Done when KBOT falls below ILO ==== 00452 * 00453 IF( KBOT.LT.ILO ) 00454 $ GO TO 80 00455 * 00456 * ==== Locate active block ==== 00457 * 00458 DO 10 K = KBOT, ILO + 1, -1 00459 IF( H( K, K-1 ).EQ.ZERO ) 00460 $ GO TO 20 00461 10 CONTINUE 00462 K = ILO 00463 20 CONTINUE 00464 KTOP = K 00465 * 00466 * ==== Select deflation window size: 00467 * . Typical Case: 00468 * . If possible and advisable, nibble the entire 00469 * . active block. If not, use size MIN(NWR,NWMAX) 00470 * . or MIN(NWR+1,NWMAX) depending upon which has 00471 * . the smaller corresponding subdiagonal entry 00472 * . (a heuristic). 00473 * . 00474 * . Exceptional Case: 00475 * . If there have been no deflations in KEXNW or 00476 * . more iterations, then vary the deflation window 00477 * . size. At first, because, larger windows are, 00478 * . in general, more powerful than smaller ones, 00479 * . rapidly increase the window to the maximum possible. 00480 * . Then, gradually reduce the window size. ==== 00481 * 00482 NH = KBOT - KTOP + 1 00483 NWUPBD = MIN( NH, NWMAX ) 00484 IF( NDFL.LT.KEXNW ) THEN 00485 NW = MIN( NWUPBD, NWR ) 00486 ELSE 00487 NW = MIN( NWUPBD, 2*NW ) 00488 END IF 00489 IF( NW.LT.NWMAX ) THEN 00490 IF( NW.GE.NH-1 ) THEN 00491 NW = NH 00492 ELSE 00493 KWTOP = KBOT - NW + 1 00494 IF( CABS1( H( KWTOP, KWTOP-1 ) ).GT. 00495 $ CABS1( H( KWTOP-1, KWTOP-2 ) ) )NW = NW + 1 00496 END IF 00497 END IF 00498 IF( NDFL.LT.KEXNW ) THEN 00499 NDEC = -1 00500 ELSE IF( NDEC.GE.0 .OR. NW.GE.NWUPBD ) THEN 00501 NDEC = NDEC + 1 00502 IF( NW-NDEC.LT.2 ) 00503 $ NDEC = 0 00504 NW = NW - NDEC 00505 END IF 00506 * 00507 * ==== Aggressive early deflation: 00508 * . split workspace under the subdiagonal into 00509 * . - an nw-by-nw work array V in the lower 00510 * . left-hand-corner, 00511 * . - an NW-by-at-least-NW-but-more-is-better 00512 * . (NW-by-NHO) horizontal work array along 00513 * . the bottom edge, 00514 * . - an at-least-NW-but-more-is-better (NHV-by-NW) 00515 * . vertical work array along the left-hand-edge. 00516 * . ==== 00517 * 00518 KV = N - NW + 1 00519 KT = NW + 1 00520 NHO = ( N-NW-1 ) - KT + 1 00521 KWV = NW + 2 00522 NVE = ( N-NW ) - KWV + 1 00523 * 00524 * ==== Aggressive early deflation ==== 00525 * 00526 CALL CLAQR2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, 00527 $ IHIZ, Z, LDZ, LS, LD, W, H( KV, 1 ), LDH, NHO, 00528 $ H( KV, KT ), LDH, NVE, H( KWV, 1 ), LDH, WORK, 00529 $ LWORK ) 00530 * 00531 * ==== Adjust KBOT accounting for new deflations. ==== 00532 * 00533 KBOT = KBOT - LD 00534 * 00535 * ==== KS points to the shifts. ==== 00536 * 00537 KS = KBOT - LS + 1 00538 * 00539 * ==== Skip an expensive QR sweep if there is a (partly 00540 * . heuristic) reason to expect that many eigenvalues 00541 * . will deflate without it. Here, the QR sweep is 00542 * . skipped if many eigenvalues have just been deflated 00543 * . or if the remaining active block is small. 00544 * 00545 IF( ( LD.EQ.0 ) .OR. ( ( 100*LD.LE.NW*NIBBLE ) .AND. ( KBOT- 00546 $ KTOP+1.GT.MIN( NMIN, NWMAX ) ) ) ) THEN 00547 * 00548 * ==== NS = nominal number of simultaneous shifts. 00549 * . This may be lowered (slightly) if CLAQR2 00550 * . did not provide that many shifts. ==== 00551 * 00552 NS = MIN( NSMAX, NSR, MAX( 2, KBOT-KTOP ) ) 00553 NS = NS - MOD( NS, 2 ) 00554 * 00555 * ==== If there have been no deflations 00556 * . in a multiple of KEXSH iterations, 00557 * . then try exceptional shifts. 00558 * . Otherwise use shifts provided by 00559 * . CLAQR2 above or from the eigenvalues 00560 * . of a trailing principal submatrix. ==== 00561 * 00562 IF( MOD( NDFL, KEXSH ).EQ.0 ) THEN 00563 KS = KBOT - NS + 1 00564 DO 30 I = KBOT, KS + 1, -2 00565 W( I ) = H( I, I ) + WILK1*CABS1( H( I, I-1 ) ) 00566 W( I-1 ) = W( I ) 00567 30 CONTINUE 00568 ELSE 00569 * 00570 * ==== Got NS/2 or fewer shifts? Use CLAHQR 00571 * . on a trailing principal submatrix to 00572 * . get more. (Since NS.LE.NSMAX.LE.(N+6)/9, 00573 * . there is enough space below the subdiagonal 00574 * . to fit an NS-by-NS scratch array.) ==== 00575 * 00576 IF( KBOT-KS+1.LE.NS / 2 ) THEN 00577 KS = KBOT - NS + 1 00578 KT = N - NS + 1 00579 CALL CLACPY( 'A', NS, NS, H( KS, KS ), LDH, 00580 $ H( KT, 1 ), LDH ) 00581 CALL CLAHQR( .false., .false., NS, 1, NS, 00582 $ H( KT, 1 ), LDH, W( KS ), 1, 1, ZDUM, 00583 $ 1, INF ) 00584 KS = KS + INF 00585 * 00586 * ==== In case of a rare QR failure use 00587 * . eigenvalues of the trailing 2-by-2 00588 * . principal submatrix. Scale to avoid 00589 * . overflows, underflows and subnormals. 00590 * . (The scale factor S can not be zero, 00591 * . because H(KBOT,KBOT-1) is nonzero.) ==== 00592 * 00593 IF( KS.GE.KBOT ) THEN 00594 S = CABS1( H( KBOT-1, KBOT-1 ) ) + 00595 $ CABS1( H( KBOT, KBOT-1 ) ) + 00596 $ CABS1( H( KBOT-1, KBOT ) ) + 00597 $ CABS1( H( KBOT, KBOT ) ) 00598 AA = H( KBOT-1, KBOT-1 ) / S 00599 CC = H( KBOT, KBOT-1 ) / S 00600 BB = H( KBOT-1, KBOT ) / S 00601 DD = H( KBOT, KBOT ) / S 00602 TR2 = ( AA+DD ) / TWO 00603 DET = ( AA-TR2 )*( DD-TR2 ) - BB*CC 00604 RTDISC = SQRT( -DET ) 00605 W( KBOT-1 ) = ( TR2+RTDISC )*S 00606 W( KBOT ) = ( TR2-RTDISC )*S 00607 * 00608 KS = KBOT - 1 00609 END IF 00610 END IF 00611 * 00612 IF( KBOT-KS+1.GT.NS ) THEN 00613 * 00614 * ==== Sort the shifts (Helps a little) ==== 00615 * 00616 SORTED = .false. 00617 DO 50 K = KBOT, KS + 1, -1 00618 IF( SORTED ) 00619 $ GO TO 60 00620 SORTED = .true. 00621 DO 40 I = KS, K - 1 00622 IF( CABS1( W( I ) ).LT.CABS1( W( I+1 ) ) ) 00623 $ THEN 00624 SORTED = .false. 00625 SWAP = W( I ) 00626 W( I ) = W( I+1 ) 00627 W( I+1 ) = SWAP 00628 END IF 00629 40 CONTINUE 00630 50 CONTINUE 00631 60 CONTINUE 00632 END IF 00633 END IF 00634 * 00635 * ==== If there are only two shifts, then use 00636 * . only one. ==== 00637 * 00638 IF( KBOT-KS+1.EQ.2 ) THEN 00639 IF( CABS1( W( KBOT )-H( KBOT, KBOT ) ).LT. 00640 $ CABS1( W( KBOT-1 )-H( KBOT, KBOT ) ) ) THEN 00641 W( KBOT-1 ) = W( KBOT ) 00642 ELSE 00643 W( KBOT ) = W( KBOT-1 ) 00644 END IF 00645 END IF 00646 * 00647 * ==== Use up to NS of the the smallest magnatiude 00648 * . shifts. If there aren't NS shifts available, 00649 * . then use them all, possibly dropping one to 00650 * . make the number of shifts even. ==== 00651 * 00652 NS = MIN( NS, KBOT-KS+1 ) 00653 NS = NS - MOD( NS, 2 ) 00654 KS = KBOT - NS + 1 00655 * 00656 * ==== Small-bulge multi-shift QR sweep: 00657 * . split workspace under the subdiagonal into 00658 * . - a KDU-by-KDU work array U in the lower 00659 * . left-hand-corner, 00660 * . - a KDU-by-at-least-KDU-but-more-is-better 00661 * . (KDU-by-NHo) horizontal work array WH along 00662 * . the bottom edge, 00663 * . - and an at-least-KDU-but-more-is-better-by-KDU 00664 * . (NVE-by-KDU) vertical work WV arrow along 00665 * . the left-hand-edge. ==== 00666 * 00667 KDU = 3*NS - 3 00668 KU = N - KDU + 1 00669 KWH = KDU + 1 00670 NHO = ( N-KDU+1-4 ) - ( KDU+1 ) + 1 00671 KWV = KDU + 4 00672 NVE = N - KDU - KWV + 1 00673 * 00674 * ==== Small-bulge multi-shift QR sweep ==== 00675 * 00676 CALL CLAQR5( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NS, 00677 $ W( KS ), H, LDH, ILOZ, IHIZ, Z, LDZ, WORK, 00678 $ 3, H( KU, 1 ), LDH, NVE, H( KWV, 1 ), LDH, 00679 $ NHO, H( KU, KWH ), LDH ) 00680 END IF 00681 * 00682 * ==== Note progress (or the lack of it). ==== 00683 * 00684 IF( LD.GT.0 ) THEN 00685 NDFL = 1 00686 ELSE 00687 NDFL = NDFL + 1 00688 END IF 00689 * 00690 * ==== End of main loop ==== 00691 70 CONTINUE 00692 * 00693 * ==== Iteration limit exceeded. Set INFO to show where 00694 * . the problem occurred and exit. ==== 00695 * 00696 INFO = KBOT 00697 80 CONTINUE 00698 END IF 00699 * 00700 * ==== Return the optimal value of LWORK. ==== 00701 * 00702 WORK( 1 ) = CMPLX( LWKOPT, 0 ) 00703 * 00704 * ==== End of CLAQR4 ==== 00705 * 00706 END