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