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