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