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