 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DLAQR0( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, 00002 $ ILOZ, 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 DOUBLE PRECISION H( LDH, * ), WI( * ), WORK( * ), WR( * ), 00014 $ Z( LDZ, * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * DLAQR0 computes the eigenvalues of a Hessenberg matrix H 00021 * and, optionally, the matrices T and Z from the Schur decomposition 00022 * H = Z T Z**T, where T is an upper quasi-triangular matrix (the 00023 * Schur form), and Z is the orthogonal matrix of Schur vectors. 00024 * 00025 * Optionally Z may be postmultiplied into an input orthogonal 00026 * matrix Q so that this routine can give the Schur factorization 00027 * of a matrix A which has been reduced to the Hessenberg form H 00028 * by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T. 00029 * 00030 * Arguments 00031 * ========= 00032 * 00033 * WANTT (input) LOGICAL 00034 * = .TRUE. : the full Schur form T is required; 00035 * = .FALSE.: only eigenvalues are required. 00036 * 00037 * WANTZ (input) LOGICAL 00038 * = .TRUE. : the matrix of Schur vectors Z is required; 00039 * = .FALSE.: Schur vectors are not required. 00040 * 00041 * N (input) INTEGER 00042 * The order of the matrix H. N .GE. 0. 00043 * 00044 * ILO (input) INTEGER 00045 * IHI (input) INTEGER 00046 * It is assumed that H is already upper triangular in rows 00047 * and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1, 00048 * H(ILO,ILO-1) is zero. ILO and IHI are normally set by a 00049 * previous call to DGEBAL, and then passed to DGEHRD when the 00050 * matrix output by DGEBAL is reduced to Hessenberg form. 00051 * Otherwise, ILO and IHI should be set to 1 and N, 00052 * respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N. 00053 * If N = 0, then ILO = 1 and IHI = 0. 00054 * 00055 * H (input/output) DOUBLE PRECISION array, dimension (LDH,N) 00056 * On entry, the upper Hessenberg matrix H. 00057 * On exit, if INFO = 0 and WANTT is .TRUE., then H contains 00058 * the upper quasi-triangular matrix T from the Schur 00059 * decomposition (the Schur form); 2-by-2 diagonal blocks 00060 * (corresponding to complex conjugate pairs of eigenvalues) 00061 * are returned in standard form, with H(i,i) = H(i+1,i+1) 00062 * and H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and WANTT is 00063 * .FALSE., then the contents of H are unspecified on exit. 00064 * (The output value of H when INFO.GT.0 is given under the 00065 * description of INFO below.) 00066 * 00067 * This subroutine may explicitly set H(i,j) = 0 for i.GT.j and 00068 * j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N. 00069 * 00070 * LDH (input) INTEGER 00071 * The leading dimension of the array H. LDH .GE. max(1,N). 00072 * 00073 * WR (output) DOUBLE PRECISION array, dimension (IHI) 00074 * WI (output) DOUBLE PRECISION array, dimension (IHI) 00075 * The real and imaginary parts, respectively, of the computed 00076 * eigenvalues of H(ILO:IHI,ILO:IHI) are stored in WR(ILO:IHI) 00077 * and WI(ILO:IHI). If two eigenvalues are computed as a 00078 * complex conjugate pair, they are stored in consecutive 00079 * elements of WR and WI, say the i-th and (i+1)th, with 00080 * WI(i) .GT. 0 and WI(i+1) .LT. 0. If WANTT is .TRUE., then 00081 * the eigenvalues are stored in the same order as on the 00082 * diagonal of the Schur form returned in H, with 00083 * WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal 00084 * block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and 00085 * WI(i+1) = -WI(i). 00086 * 00087 * ILOZ (input) INTEGER 00088 * IHIZ (input) INTEGER 00089 * Specify the rows of Z to which transformations must be 00090 * applied if WANTZ is .TRUE.. 00091 * 1 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N. 00092 * 00093 * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,IHI) 00094 * If WANTZ is .FALSE., then Z is not referenced. 00095 * If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is 00096 * replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the 00097 * orthogonal Schur factor of H(ILO:IHI,ILO:IHI). 00098 * (The output value of Z when INFO.GT.0 is given under 00099 * the description of INFO below.) 00100 * 00101 * LDZ (input) INTEGER 00102 * The leading dimension of the array Z. if WANTZ is .TRUE. 00103 * then LDZ.GE.MAX(1,IHIZ). Otherwize, LDZ.GE.1. 00104 * 00105 * WORK (workspace/output) DOUBLE PRECISION array, dimension LWORK 00106 * On exit, if LWORK = -1, WORK(1) returns an estimate of 00107 * the optimal value for LWORK. 00108 * 00109 * LWORK (input) INTEGER 00110 * The dimension of the array WORK. LWORK .GE. max(1,N) 00111 * is sufficient, but LWORK typically as large as 6*N may 00112 * be required for optimal performance. A workspace query 00113 * to determine the optimal workspace size is recommended. 00114 * 00115 * If LWORK = -1, then DLAQR0 does a workspace query. 00116 * In this case, DLAQR0 checks the input parameters and 00117 * estimates the optimal workspace size for the given 00118 * values of N, ILO and IHI. The estimate is returned 00119 * in WORK(1). No error message related to LWORK is 00120 * issued by XERBLA. Neither H nor Z are accessed. 00121 * 00122 * 00123 * INFO (output) INTEGER 00124 * = 0: successful exit 00125 * .GT. 0: if INFO = i, DLAQR0 failed to compute all of 00126 * the eigenvalues. Elements 1:ilo-1 and i+1:n of WR 00127 * and WI contain those eigenvalues which have been 00128 * successfully computed. (Failures are rare.) 00129 * 00130 * If INFO .GT. 0 and WANT is .FALSE., then on exit, 00131 * the remaining unconverged eigenvalues are the eigen- 00132 * values of the upper Hessenberg matrix rows and 00133 * columns ILO through INFO of the final, output 00134 * value of H. 00135 * 00136 * If INFO .GT. 0 and WANTT is .TRUE., then on exit 00137 * 00138 * (*) (initial value of H)*U = U*(final value of H) 00139 * 00140 * where U is an orthogonal matrix. The final 00141 * value of H is upper Hessenberg and quasi-triangular 00142 * in rows and columns INFO+1 through IHI. 00143 * 00144 * If INFO .GT. 0 and WANTZ is .TRUE., then on exit 00145 * 00146 * (final value of Z(ILO:IHI,ILOZ:IHIZ) 00147 * = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U 00148 * 00149 * where U is the orthogonal matrix in (*) (regard- 00150 * less of the value of WANTT.) 00151 * 00152 * If INFO .GT. 0 and WANTZ is .FALSE., then Z is not 00153 * accessed. 00154 * 00155 * ================================================================ 00156 * Based on contributions by 00157 * Karen Braman and Ralph Byers, Department of Mathematics, 00158 * University of Kansas, USA 00159 * 00160 * ================================================================ 00161 * References: 00162 * K. Braman, R. Byers and R. Mathias, The Multi-Shift QR 00163 * Algorithm Part I: Maintaining Well Focused Shifts, and Level 3 00164 * Performance, SIAM Journal of Matrix Analysis, volume 23, pages 00165 * 929--947, 2002. 00166 * 00167 * K. Braman, R. Byers and R. Mathias, The Multi-Shift QR 00168 * Algorithm Part II: Aggressive Early Deflation, SIAM Journal 00169 * of Matrix Analysis, volume 23, pages 948--973, 2002. 00170 * 00171 * ================================================================ 00172 * .. Parameters .. 00173 * 00174 * ==== Matrices of order NTINY or smaller must be processed by 00175 * . DLAHQR because of insufficient subdiagonal scratch space. 00176 * . (This is a hard limit.) ==== 00177 INTEGER NTINY 00178 PARAMETER ( NTINY = 11 ) 00179 * 00180 * ==== Exceptional deflation windows: try to cure rare 00181 * . slow convergence by varying the size of the 00182 * . deflation window after KEXNW iterations. ==== 00183 INTEGER KEXNW 00184 PARAMETER ( KEXNW = 5 ) 00185 * 00186 * ==== Exceptional shifts: try to cure rare slow convergence 00187 * . with ad-hoc exceptional shifts every KEXSH iterations. 00188 * . ==== 00189 INTEGER KEXSH 00190 PARAMETER ( KEXSH = 6 ) 00191 * 00192 * ==== The constants WILK1 and WILK2 are used to form the 00193 * . exceptional shifts. ==== 00194 DOUBLE PRECISION WILK1, WILK2 00195 PARAMETER ( WILK1 = 0.75d0, WILK2 = -0.4375d0 ) 00196 DOUBLE PRECISION ZERO, ONE 00197 PARAMETER ( ZERO = 0.0d0, ONE = 1.0d0 ) 00198 * .. 00199 * .. Local Scalars .. 00200 DOUBLE PRECISION AA, BB, CC, CS, DD, SN, SS, SWAP 00201 INTEGER I, INF, IT, ITMAX, K, KACC22, KBOT, KDU, KS, 00202 $ KT, KTOP, KU, KV, KWH, KWTOP, KWV, LD, LS, 00203 $ LWKOPT, NDEC, NDFL, NH, NHO, NIBBLE, NMIN, NS, 00204 $ NSMAX, NSR, NVE, NW, NWMAX, NWR, NWUPBD 00205 LOGICAL SORTED 00206 CHARACTER JBCMPZ*2 00207 * .. 00208 * .. External Functions .. 00209 INTEGER ILAENV 00210 EXTERNAL ILAENV 00211 * .. 00212 * .. Local Arrays .. 00213 DOUBLE PRECISION ZDUM( 1, 1 ) 00214 * .. 00215 * .. External Subroutines .. 00216 EXTERNAL DLACPY, DLAHQR, DLANV2, DLAQR3, DLAQR4, DLAQR5 00217 * .. 00218 * .. Intrinsic Functions .. 00219 INTRINSIC ABS, DBLE, INT, MAX, MIN, MOD 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 DLAHQR. ==== 00234 * 00235 LWKOPT = 1 00236 IF( LWORK.NE.-1 ) 00237 $ CALL DLAHQR( WANTT, WANTZ, N, ILO, IHI, H, LDH, WR, WI, 00238 $ ILOZ, 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, 'DLAQR0', 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, 'DLAQR0', 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 DLAQR3 ==== 00283 * 00284 CALL DLAQR3( WANTT, WANTZ, N, ILO, IHI, NWR+1, H, LDH, ILOZ, 00285 $ IHIZ, Z, LDZ, LS, LD, WR, WI, H, LDH, N, H, LDH, 00286 $ N, H, LDH, WORK, -1 ) 00287 * 00288 * ==== Optimal workspace = MAX(DLAQR5, DLAQR3) ==== 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 ) = DBLE( LWKOPT ) 00296 RETURN 00297 END IF 00298 * 00299 * ==== DLAHQR/DLAQR0 crossover point ==== 00300 * 00301 NMIN = ILAENV( 12, 'DLAQR0', JBCMPZ, N, ILO, IHI, LWORK ) 00302 NMIN = MAX( NTINY, NMIN ) 00303 * 00304 * ==== Nibble crossover point ==== 00305 * 00306 NIBBLE = ILAENV( 14, 'DLAQR0', 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, 'DLAQR0', 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 80 IT = 1, ITMAX 00343 * 00344 * ==== Done when KBOT falls below ILO ==== 00345 * 00346 IF( KBOT.LT.ILO ) 00347 $ GO TO 90 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( ABS( H( KWTOP, KWTOP-1 ) ).GT. 00388 $ ABS( 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 DLAQR3( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, 00420 $ IHIZ, Z, LDZ, LS, LD, WR, WI, H( KV, 1 ), LDH, 00421 $ NHO, H( KV, KT ), LDH, NVE, H( KWV, 1 ), LDH, 00422 $ WORK, 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 DLAQR3 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 * . DLAQR3 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, MAX( KS+1, KTOP+2 ), -2 00458 SS = ABS( H( I, I-1 ) ) + ABS( H( I-1, I-2 ) ) 00459 AA = WILK1*SS + H( I, I ) 00460 BB = SS 00461 CC = WILK2*SS 00462 DD = AA 00463 CALL DLANV2( AA, BB, CC, DD, WR( I-1 ), WI( I-1 ), 00464 $ WR( I ), WI( I ), CS, SN ) 00465 30 CONTINUE 00466 IF( KS.EQ.KTOP ) THEN 00467 WR( KS+1 ) = H( KS+1, KS+1 ) 00468 WI( KS+1 ) = ZERO 00469 WR( KS ) = WR( KS+1 ) 00470 WI( KS ) = WI( KS+1 ) 00471 END IF 00472 ELSE 00473 * 00474 * ==== Got NS/2 or fewer shifts? Use DLAQR4 or 00475 * . DLAHQR on a trailing principal submatrix to 00476 * . get more. (Since NS.LE.NSMAX.LE.(N+6)/9, 00477 * . there is enough space below the subdiagonal 00478 * . to fit an NS-by-NS scratch array.) ==== 00479 * 00480 IF( KBOT-KS+1.LE.NS / 2 ) THEN 00481 KS = KBOT - NS + 1 00482 KT = N - NS + 1 00483 CALL DLACPY( 'A', NS, NS, H( KS, KS ), LDH, 00484 $ H( KT, 1 ), LDH ) 00485 IF( NS.GT.NMIN ) THEN 00486 CALL DLAQR4( .false., .false., NS, 1, NS, 00487 $ H( KT, 1 ), LDH, WR( KS ), 00488 $ WI( KS ), 1, 1, ZDUM, 1, WORK, 00489 $ LWORK, INF ) 00490 ELSE 00491 CALL DLAHQR( .false., .false., NS, 1, NS, 00492 $ H( KT, 1 ), LDH, WR( KS ), 00493 $ WI( KS ), 1, 1, ZDUM, 1, INF ) 00494 END IF 00495 KS = KS + INF 00496 * 00497 * ==== In case of a rare QR failure use 00498 * . eigenvalues of the trailing 2-by-2 00499 * . principal submatrix. ==== 00500 * 00501 IF( KS.GE.KBOT ) THEN 00502 AA = H( KBOT-1, KBOT-1 ) 00503 CC = H( KBOT, KBOT-1 ) 00504 BB = H( KBOT-1, KBOT ) 00505 DD = H( KBOT, KBOT ) 00506 CALL DLANV2( AA, BB, CC, DD, WR( KBOT-1 ), 00507 $ WI( KBOT-1 ), WR( KBOT ), 00508 $ WI( KBOT ), CS, SN ) 00509 KS = KBOT - 1 00510 END IF 00511 END IF 00512 * 00513 IF( KBOT-KS+1.GT.NS ) THEN 00514 * 00515 * ==== Sort the shifts (Helps a little) 00516 * . Bubble sort keeps complex conjugate 00517 * . pairs together. ==== 00518 * 00519 SORTED = .false. 00520 DO 50 K = KBOT, KS + 1, -1 00521 IF( SORTED ) 00522 $ GO TO 60 00523 SORTED = .true. 00524 DO 40 I = KS, K - 1 00525 IF( ABS( WR( I ) )+ABS( WI( I ) ).LT. 00526 $ ABS( WR( I+1 ) )+ABS( WI( I+1 ) ) ) THEN 00527 SORTED = .false. 00528 * 00529 SWAP = WR( I ) 00530 WR( I ) = WR( I+1 ) 00531 WR( I+1 ) = SWAP 00532 * 00533 SWAP = WI( I ) 00534 WI( I ) = WI( I+1 ) 00535 WI( I+1 ) = SWAP 00536 END IF 00537 40 CONTINUE 00538 50 CONTINUE 00539 60 CONTINUE 00540 END IF 00541 * 00542 * ==== Shuffle shifts into pairs of real shifts 00543 * . and pairs of complex conjugate shifts 00544 * . assuming complex conjugate shifts are 00545 * . already adjacent to one another. (Yes, 00546 * . they are.) ==== 00547 * 00548 DO 70 I = KBOT, KS + 2, -2 00549 IF( WI( I ).NE.-WI( I-1 ) ) THEN 00550 * 00551 SWAP = WR( I ) 00552 WR( I ) = WR( I-1 ) 00553 WR( I-1 ) = WR( I-2 ) 00554 WR( I-2 ) = SWAP 00555 * 00556 SWAP = WI( I ) 00557 WI( I ) = WI( I-1 ) 00558 WI( I-1 ) = WI( I-2 ) 00559 WI( I-2 ) = SWAP 00560 END IF 00561 70 CONTINUE 00562 END IF 00563 * 00564 * ==== If there are only two shifts and both are 00565 * . real, then use only one. ==== 00566 * 00567 IF( KBOT-KS+1.EQ.2 ) THEN 00568 IF( WI( KBOT ).EQ.ZERO ) THEN 00569 IF( ABS( WR( KBOT )-H( KBOT, KBOT ) ).LT. 00570 $ ABS( WR( KBOT-1 )-H( KBOT, KBOT ) ) ) THEN 00571 WR( KBOT-1 ) = WR( KBOT ) 00572 ELSE 00573 WR( KBOT ) = WR( KBOT-1 ) 00574 END IF 00575 END IF 00576 END IF 00577 * 00578 * ==== Use up to NS of the the smallest magnatiude 00579 * . shifts. If there aren't NS shifts available, 00580 * . then use them all, possibly dropping one to 00581 * . make the number of shifts even. ==== 00582 * 00583 NS = MIN( NS, KBOT-KS+1 ) 00584 NS = NS - MOD( NS, 2 ) 00585 KS = KBOT - NS + 1 00586 * 00587 * ==== Small-bulge multi-shift QR sweep: 00588 * . split workspace under the subdiagonal into 00589 * . - a KDU-by-KDU work array U in the lower 00590 * . left-hand-corner, 00591 * . - a KDU-by-at-least-KDU-but-more-is-better 00592 * . (KDU-by-NHo) horizontal work array WH along 00593 * . the bottom edge, 00594 * . - and an at-least-KDU-but-more-is-better-by-KDU 00595 * . (NVE-by-KDU) vertical work WV arrow along 00596 * . the left-hand-edge. ==== 00597 * 00598 KDU = 3*NS - 3 00599 KU = N - KDU + 1 00600 KWH = KDU + 1 00601 NHO = ( N-KDU+1-4 ) - ( KDU+1 ) + 1 00602 KWV = KDU + 4 00603 NVE = N - KDU - KWV + 1 00604 * 00605 * ==== Small-bulge multi-shift QR sweep ==== 00606 * 00607 CALL DLAQR5( WANTT, WANTZ, KACC22, N, KTOP, KBOT, NS, 00608 $ WR( KS ), WI( KS ), H, LDH, ILOZ, IHIZ, Z, 00609 $ LDZ, WORK, 3, H( KU, 1 ), LDH, NVE, 00610 $ H( KWV, 1 ), LDH, NHO, H( KU, KWH ), LDH ) 00611 END IF 00612 * 00613 * ==== Note progress (or the lack of it). ==== 00614 * 00615 IF( LD.GT.0 ) THEN 00616 NDFL = 1 00617 ELSE 00618 NDFL = NDFL + 1 00619 END IF 00620 * 00621 * ==== End of main loop ==== 00622 80 CONTINUE 00623 * 00624 * ==== Iteration limit exceeded. Set INFO to show where 00625 * . the problem occurred and exit. ==== 00626 * 00627 INFO = KBOT 00628 90 CONTINUE 00629 END IF 00630 * 00631 * ==== Return the optimal value of LWORK. ==== 00632 * 00633 WORK( 1 ) = DBLE( LWKOPT ) 00634 * 00635 * ==== End of DLAQR0 ==== 00636 * 00637 END
1.7.3 
