00001 /* dlaqr0.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Table of constant values */ 00017 00018 static integer c__13 = 13; 00019 static integer c__15 = 15; 00020 static integer c_n1 = -1; 00021 static integer c__12 = 12; 00022 static integer c__14 = 14; 00023 static integer c__16 = 16; 00024 static logical c_false = FALSE_; 00025 static integer c__1 = 1; 00026 static integer c__3 = 3; 00027 00028 /* Subroutine */ int dlaqr0_(logical *wantt, logical *wantz, integer *n, 00029 integer *ilo, integer *ihi, doublereal *h__, integer *ldh, doublereal 00030 *wr, doublereal *wi, integer *iloz, integer *ihiz, doublereal *z__, 00031 integer *ldz, doublereal *work, integer *lwork, integer *info) 00032 { 00033 /* System generated locals */ 00034 integer h_dim1, h_offset, z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5; 00035 doublereal d__1, d__2, d__3, d__4; 00036 00037 /* Local variables */ 00038 integer i__, k; 00039 doublereal aa, bb, cc, dd; 00040 integer ld; 00041 doublereal cs; 00042 integer nh, it, ks, kt; 00043 doublereal sn; 00044 integer ku, kv, ls, ns; 00045 doublereal ss; 00046 integer nw, inf, kdu, nho, nve, kwh, nsr, nwr, kwv, ndec, ndfl, kbot, 00047 nmin; 00048 doublereal swap; 00049 integer ktop; 00050 doublereal zdum[1] /* was [1][1] */; 00051 integer kacc22, itmax, nsmax, nwmax, kwtop; 00052 extern /* Subroutine */ int dlanv2_(doublereal *, doublereal *, 00053 doublereal *, doublereal *, doublereal *, doublereal *, 00054 doublereal *, doublereal *, doublereal *, doublereal *), dlaqr3_( 00055 logical *, logical *, integer *, integer *, integer *, integer *, 00056 doublereal *, integer *, integer *, integer *, doublereal *, 00057 integer *, integer *, integer *, doublereal *, doublereal *, 00058 doublereal *, integer *, integer *, doublereal *, integer *, 00059 integer *, doublereal *, integer *, doublereal *, integer *), 00060 dlaqr4_(logical *, logical *, integer *, integer *, integer *, 00061 doublereal *, integer *, doublereal *, doublereal *, integer *, 00062 integer *, doublereal *, integer *, doublereal *, integer *, 00063 integer *), dlaqr5_(logical *, logical *, integer *, integer *, 00064 integer *, integer *, integer *, doublereal *, doublereal *, 00065 doublereal *, integer *, integer *, integer *, doublereal *, 00066 integer *, doublereal *, integer *, doublereal *, integer *, 00067 integer *, doublereal *, integer *, integer *, doublereal *, 00068 integer *); 00069 integer nibble; 00070 extern /* Subroutine */ int dlahqr_(logical *, logical *, integer *, 00071 integer *, integer *, doublereal *, integer *, doublereal *, 00072 doublereal *, integer *, integer *, doublereal *, integer *, 00073 integer *), dlacpy_(char *, integer *, integer *, doublereal *, 00074 integer *, doublereal *, integer *); 00075 extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 00076 integer *, integer *); 00077 char jbcmpz[2]; 00078 integer nwupbd; 00079 logical sorted; 00080 integer lwkopt; 00081 00082 00083 /* -- LAPACK auxiliary routine (version 3.2) -- */ 00084 /* Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd.. */ 00085 /* November 2006 */ 00086 00087 /* .. Scalar Arguments .. */ 00088 /* .. */ 00089 /* .. Array Arguments .. */ 00090 /* .. */ 00091 00092 /* Purpose */ 00093 /* ======= */ 00094 00095 /* DLAQR0 computes the eigenvalues of a Hessenberg matrix H */ 00096 /* and, optionally, the matrices T and Z from the Schur decomposition */ 00097 /* H = Z T Z**T, where T is an upper quasi-triangular matrix (the */ 00098 /* Schur form), and Z is the orthogonal matrix of Schur vectors. */ 00099 00100 /* Optionally Z may be postmultiplied into an input orthogonal */ 00101 /* matrix Q so that this routine can give the Schur factorization */ 00102 /* of a matrix A which has been reduced to the Hessenberg form H */ 00103 /* by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T. */ 00104 00105 /* Arguments */ 00106 /* ========= */ 00107 00108 /* WANTT (input) LOGICAL */ 00109 /* = .TRUE. : the full Schur form T is required; */ 00110 /* = .FALSE.: only eigenvalues are required. */ 00111 00112 /* WANTZ (input) LOGICAL */ 00113 /* = .TRUE. : the matrix of Schur vectors Z is required; */ 00114 /* = .FALSE.: Schur vectors are not required. */ 00115 00116 /* N (input) INTEGER */ 00117 /* The order of the matrix H. N .GE. 0. */ 00118 00119 /* ILO (input) INTEGER */ 00120 /* IHI (input) INTEGER */ 00121 /* It is assumed that H is already upper triangular in rows */ 00122 /* and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1, */ 00123 /* H(ILO,ILO-1) is zero. ILO and IHI are normally set by a */ 00124 /* previous call to DGEBAL, and then passed to DGEHRD when the */ 00125 /* matrix output by DGEBAL is reduced to Hessenberg form. */ 00126 /* Otherwise, ILO and IHI should be set to 1 and N, */ 00127 /* respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N. */ 00128 /* If N = 0, then ILO = 1 and IHI = 0. */ 00129 00130 /* H (input/output) DOUBLE PRECISION array, dimension (LDH,N) */ 00131 /* On entry, the upper Hessenberg matrix H. */ 00132 /* On exit, if INFO = 0 and WANTT is .TRUE., then H contains */ 00133 /* the upper quasi-triangular matrix T from the Schur */ 00134 /* decomposition (the Schur form); 2-by-2 diagonal blocks */ 00135 /* (corresponding to complex conjugate pairs of eigenvalues) */ 00136 /* are returned in standard form, with H(i,i) = H(i+1,i+1) */ 00137 /* and H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and WANTT is */ 00138 /* .FALSE., then the contents of H are unspecified on exit. */ 00139 /* (The output value of H when INFO.GT.0 is given under the */ 00140 /* description of INFO below.) */ 00141 00142 /* This subroutine may explicitly set H(i,j) = 0 for i.GT.j and */ 00143 /* j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N. */ 00144 00145 /* LDH (input) INTEGER */ 00146 /* The leading dimension of the array H. LDH .GE. max(1,N). */ 00147 00148 /* WR (output) DOUBLE PRECISION array, dimension (IHI) */ 00149 /* WI (output) DOUBLE PRECISION array, dimension (IHI) */ 00150 /* The real and imaginary parts, respectively, of the computed */ 00151 /* eigenvalues of H(ILO:IHI,ILO:IHI) are stored in WR(ILO:IHI) */ 00152 /* and WI(ILO:IHI). If two eigenvalues are computed as a */ 00153 /* complex conjugate pair, they are stored in consecutive */ 00154 /* elements of WR and WI, say the i-th and (i+1)th, with */ 00155 /* WI(i) .GT. 0 and WI(i+1) .LT. 0. If WANTT is .TRUE., then */ 00156 /* the eigenvalues are stored in the same order as on the */ 00157 /* diagonal of the Schur form returned in H, with */ 00158 /* WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal */ 00159 /* block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and */ 00160 /* WI(i+1) = -WI(i). */ 00161 00162 /* ILOZ (input) INTEGER */ 00163 /* IHIZ (input) INTEGER */ 00164 /* Specify the rows of Z to which transformations must be */ 00165 /* applied if WANTZ is .TRUE.. */ 00166 /* 1 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N. */ 00167 00168 /* Z (input/output) DOUBLE PRECISION array, dimension (LDZ,IHI) */ 00169 /* If WANTZ is .FALSE., then Z is not referenced. */ 00170 /* If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is */ 00171 /* replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the */ 00172 /* orthogonal Schur factor of H(ILO:IHI,ILO:IHI). */ 00173 /* (The output value of Z when INFO.GT.0 is given under */ 00174 /* the description of INFO below.) */ 00175 00176 /* LDZ (input) INTEGER */ 00177 /* The leading dimension of the array Z. if WANTZ is .TRUE. */ 00178 /* then LDZ.GE.MAX(1,IHIZ). Otherwize, LDZ.GE.1. */ 00179 00180 /* WORK (workspace/output) DOUBLE PRECISION array, dimension LWORK */ 00181 /* On exit, if LWORK = -1, WORK(1) returns an estimate of */ 00182 /* the optimal value for LWORK. */ 00183 00184 /* LWORK (input) INTEGER */ 00185 /* The dimension of the array WORK. LWORK .GE. max(1,N) */ 00186 /* is sufficient, but LWORK typically as large as 6*N may */ 00187 /* be required for optimal performance. A workspace query */ 00188 /* to determine the optimal workspace size is recommended. */ 00189 00190 /* If LWORK = -1, then DLAQR0 does a workspace query. */ 00191 /* In this case, DLAQR0 checks the input parameters and */ 00192 /* estimates the optimal workspace size for the given */ 00193 /* values of N, ILO and IHI. The estimate is returned */ 00194 /* in WORK(1). No error message related to LWORK is */ 00195 /* issued by XERBLA. Neither H nor Z are accessed. */ 00196 00197 00198 /* INFO (output) INTEGER */ 00199 /* = 0: successful exit */ 00200 /* .GT. 0: if INFO = i, DLAQR0 failed to compute all of */ 00201 /* the eigenvalues. Elements 1:ilo-1 and i+1:n of WR */ 00202 /* and WI contain those eigenvalues which have been */ 00203 /* successfully computed. (Failures are rare.) */ 00204 00205 /* If INFO .GT. 0 and WANT is .FALSE., then on exit, */ 00206 /* the remaining unconverged eigenvalues are the eigen- */ 00207 /* values of the upper Hessenberg matrix rows and */ 00208 /* columns ILO through INFO of the final, output */ 00209 /* value of H. */ 00210 00211 /* If INFO .GT. 0 and WANTT is .TRUE., then on exit */ 00212 00213 /* (*) (initial value of H)*U = U*(final value of H) */ 00214 00215 /* where U is an orthogonal matrix. The final */ 00216 /* value of H is upper Hessenberg and quasi-triangular */ 00217 /* in rows and columns INFO+1 through IHI. */ 00218 00219 /* If INFO .GT. 0 and WANTZ is .TRUE., then on exit */ 00220 00221 /* (final value of Z(ILO:IHI,ILOZ:IHIZ) */ 00222 /* = (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U */ 00223 00224 /* where U is the orthogonal matrix in (*) (regard- */ 00225 /* less of the value of WANTT.) */ 00226 00227 /* If INFO .GT. 0 and WANTZ is .FALSE., then Z is not */ 00228 /* accessed. */ 00229 00230 /* ================================================================ */ 00231 /* Based on contributions by */ 00232 /* Karen Braman and Ralph Byers, Department of Mathematics, */ 00233 /* University of Kansas, USA */ 00234 00235 /* ================================================================ */ 00236 /* References: */ 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 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 /* .. Parameters .. */ 00248 00249 /* ==== Matrices of order NTINY or smaller must be processed by */ 00250 /* . DLAHQR because of insufficient subdiagonal scratch space. */ 00251 /* . (This is a hard limit.) ==== */ 00252 00253 /* ==== Exceptional deflation windows: try to cure rare */ 00254 /* . slow convergence by varying the size of the */ 00255 /* . deflation window after KEXNW iterations. ==== */ 00256 00257 /* ==== Exceptional shifts: try to cure rare slow convergence */ 00258 /* . with ad-hoc exceptional shifts every KEXSH iterations. */ 00259 /* . ==== */ 00260 00261 /* ==== The constants WILK1 and WILK2 are used to form the */ 00262 /* . exceptional shifts. ==== */ 00263 /* .. */ 00264 /* .. Local Scalars .. */ 00265 /* .. */ 00266 /* .. External Functions .. */ 00267 /* .. */ 00268 /* .. Local Arrays .. */ 00269 /* .. */ 00270 /* .. External Subroutines .. */ 00271 /* .. */ 00272 /* .. Intrinsic Functions .. */ 00273 /* .. */ 00274 /* .. Executable Statements .. */ 00275 /* Parameter adjustments */ 00276 h_dim1 = *ldh; 00277 h_offset = 1 + h_dim1; 00278 h__ -= h_offset; 00279 --wr; 00280 --wi; 00281 z_dim1 = *ldz; 00282 z_offset = 1 + z_dim1; 00283 z__ -= z_offset; 00284 --work; 00285 00286 /* Function Body */ 00287 *info = 0; 00288 00289 /* ==== Quick return for N = 0: nothing to do. ==== */ 00290 00291 if (*n == 0) { 00292 work[1] = 1.; 00293 return 0; 00294 } 00295 00296 if (*n <= 11) { 00297 00298 /* ==== Tiny matrices must use DLAHQR. ==== */ 00299 00300 lwkopt = 1; 00301 if (*lwork != -1) { 00302 dlahqr_(wantt, wantz, n, ilo, ihi, &h__[h_offset], ldh, &wr[1], & 00303 wi[1], iloz, ihiz, &z__[z_offset], ldz, info); 00304 } 00305 } else { 00306 00307 /* ==== Use small bulge multi-shift QR with aggressive early */ 00308 /* . deflation on larger-than-tiny matrices. ==== */ 00309 00310 /* ==== Hope for the best. ==== */ 00311 00312 *info = 0; 00313 00314 /* ==== Set up job flags for ILAENV. ==== */ 00315 00316 if (*wantt) { 00317 *(unsigned char *)jbcmpz = 'S'; 00318 } else { 00319 *(unsigned char *)jbcmpz = 'E'; 00320 } 00321 if (*wantz) { 00322 *(unsigned char *)&jbcmpz[1] = 'V'; 00323 } else { 00324 *(unsigned char *)&jbcmpz[1] = 'N'; 00325 } 00326 00327 /* ==== NWR = recommended deflation window size. At this */ 00328 /* . point, N .GT. NTINY = 11, so there is enough */ 00329 /* . subdiagonal workspace for NWR.GE.2 as required. */ 00330 /* . (In fact, there is enough subdiagonal space for */ 00331 /* . NWR.GE.3.) ==== */ 00332 00333 nwr = ilaenv_(&c__13, "DLAQR0", jbcmpz, n, ilo, ihi, lwork); 00334 nwr = max(2,nwr); 00335 /* Computing MIN */ 00336 i__1 = *ihi - *ilo + 1, i__2 = (*n - 1) / 3, i__1 = min(i__1,i__2); 00337 nwr = min(i__1,nwr); 00338 00339 /* ==== NSR = recommended number of simultaneous shifts. */ 00340 /* . At this point N .GT. NTINY = 11, so there is at */ 00341 /* . enough subdiagonal workspace for NSR to be even */ 00342 /* . and greater than or equal to two as required. ==== */ 00343 00344 nsr = ilaenv_(&c__15, "DLAQR0", jbcmpz, n, ilo, ihi, lwork); 00345 /* Computing MIN */ 00346 i__1 = nsr, i__2 = (*n + 6) / 9, i__1 = min(i__1,i__2), i__2 = *ihi - 00347 *ilo; 00348 nsr = min(i__1,i__2); 00349 /* Computing MAX */ 00350 i__1 = 2, i__2 = nsr - nsr % 2; 00351 nsr = max(i__1,i__2); 00352 00353 /* ==== Estimate optimal workspace ==== */ 00354 00355 /* ==== Workspace query call to DLAQR3 ==== */ 00356 00357 i__1 = nwr + 1; 00358 dlaqr3_(wantt, wantz, n, ilo, ihi, &i__1, &h__[h_offset], ldh, iloz, 00359 ihiz, &z__[z_offset], ldz, &ls, &ld, &wr[1], &wi[1], &h__[ 00360 h_offset], ldh, n, &h__[h_offset], ldh, n, &h__[h_offset], 00361 ldh, &work[1], &c_n1); 00362 00363 /* ==== Optimal workspace = MAX(DLAQR5, DLAQR3) ==== */ 00364 00365 /* Computing MAX */ 00366 i__1 = nsr * 3 / 2, i__2 = (integer) work[1]; 00367 lwkopt = max(i__1,i__2); 00368 00369 /* ==== Quick return in case of workspace query. ==== */ 00370 00371 if (*lwork == -1) { 00372 work[1] = (doublereal) lwkopt; 00373 return 0; 00374 } 00375 00376 /* ==== DLAHQR/DLAQR0 crossover point ==== */ 00377 00378 nmin = ilaenv_(&c__12, "DLAQR0", jbcmpz, n, ilo, ihi, lwork); 00379 nmin = max(11,nmin); 00380 00381 /* ==== Nibble crossover point ==== */ 00382 00383 nibble = ilaenv_(&c__14, "DLAQR0", jbcmpz, n, ilo, ihi, lwork); 00384 nibble = max(0,nibble); 00385 00386 /* ==== Accumulate reflections during ttswp? Use block */ 00387 /* . 2-by-2 structure during matrix-matrix multiply? ==== */ 00388 00389 kacc22 = ilaenv_(&c__16, "DLAQR0", jbcmpz, n, ilo, ihi, lwork); 00390 kacc22 = max(0,kacc22); 00391 kacc22 = min(2,kacc22); 00392 00393 /* ==== NWMAX = the largest possible deflation window for */ 00394 /* . which there is sufficient workspace. ==== */ 00395 00396 /* Computing MIN */ 00397 i__1 = (*n - 1) / 3, i__2 = *lwork / 2; 00398 nwmax = min(i__1,i__2); 00399 nw = nwmax; 00400 00401 /* ==== NSMAX = the Largest number of simultaneous shifts */ 00402 /* . for which there is sufficient workspace. ==== */ 00403 00404 /* Computing MIN */ 00405 i__1 = (*n + 6) / 9, i__2 = (*lwork << 1) / 3; 00406 nsmax = min(i__1,i__2); 00407 nsmax -= nsmax % 2; 00408 00409 /* ==== NDFL: an iteration count restarted at deflation. ==== */ 00410 00411 ndfl = 1; 00412 00413 /* ==== ITMAX = iteration limit ==== */ 00414 00415 /* Computing MAX */ 00416 i__1 = 10, i__2 = *ihi - *ilo + 1; 00417 itmax = max(i__1,i__2) * 30; 00418 00419 /* ==== Last row and column in the active block ==== */ 00420 00421 kbot = *ihi; 00422 00423 /* ==== Main Loop ==== */ 00424 00425 i__1 = itmax; 00426 for (it = 1; it <= i__1; ++it) { 00427 00428 /* ==== Done when KBOT falls below ILO ==== */ 00429 00430 if (kbot < *ilo) { 00431 goto L90; 00432 } 00433 00434 /* ==== Locate active block ==== */ 00435 00436 i__2 = *ilo + 1; 00437 for (k = kbot; k >= i__2; --k) { 00438 if (h__[k + (k - 1) * h_dim1] == 0.) { 00439 goto L20; 00440 } 00441 /* L10: */ 00442 } 00443 k = *ilo; 00444 L20: 00445 ktop = k; 00446 00447 /* ==== Select deflation window size: */ 00448 /* . Typical Case: */ 00449 /* . If possible and advisable, nibble the entire */ 00450 /* . active block. If not, use size MIN(NWR,NWMAX) */ 00451 /* . or MIN(NWR+1,NWMAX) depending upon which has */ 00452 /* . the smaller corresponding subdiagonal entry */ 00453 /* . (a heuristic). */ 00454 /* . */ 00455 /* . Exceptional Case: */ 00456 /* . If there have been no deflations in KEXNW or */ 00457 /* . more iterations, then vary the deflation window */ 00458 /* . size. At first, because, larger windows are, */ 00459 /* . in general, more powerful than smaller ones, */ 00460 /* . rapidly increase the window to the maximum possible. */ 00461 /* . Then, gradually reduce the window size. ==== */ 00462 00463 nh = kbot - ktop + 1; 00464 nwupbd = min(nh,nwmax); 00465 if (ndfl < 5) { 00466 nw = min(nwupbd,nwr); 00467 } else { 00468 /* Computing MIN */ 00469 i__2 = nwupbd, i__3 = nw << 1; 00470 nw = min(i__2,i__3); 00471 } 00472 if (nw < nwmax) { 00473 if (nw >= nh - 1) { 00474 nw = nh; 00475 } else { 00476 kwtop = kbot - nw + 1; 00477 if ((d__1 = h__[kwtop + (kwtop - 1) * h_dim1], abs(d__1)) 00478 > (d__2 = h__[kwtop - 1 + (kwtop - 2) * h_dim1], 00479 abs(d__2))) { 00480 ++nw; 00481 } 00482 } 00483 } 00484 if (ndfl < 5) { 00485 ndec = -1; 00486 } else if (ndec >= 0 || nw >= nwupbd) { 00487 ++ndec; 00488 if (nw - ndec < 2) { 00489 ndec = 0; 00490 } 00491 nw -= ndec; 00492 } 00493 00494 /* ==== Aggressive early deflation: */ 00495 /* . split workspace under the subdiagonal into */ 00496 /* . - an nw-by-nw work array V in the lower */ 00497 /* . left-hand-corner, */ 00498 /* . - an NW-by-at-least-NW-but-more-is-better */ 00499 /* . (NW-by-NHO) horizontal work array along */ 00500 /* . the bottom edge, */ 00501 /* . - an at-least-NW-but-more-is-better (NHV-by-NW) */ 00502 /* . vertical work array along the left-hand-edge. */ 00503 /* . ==== */ 00504 00505 kv = *n - nw + 1; 00506 kt = nw + 1; 00507 nho = *n - nw - 1 - kt + 1; 00508 kwv = nw + 2; 00509 nve = *n - nw - kwv + 1; 00510 00511 /* ==== Aggressive early deflation ==== */ 00512 00513 dlaqr3_(wantt, wantz, n, &ktop, &kbot, &nw, &h__[h_offset], ldh, 00514 iloz, ihiz, &z__[z_offset], ldz, &ls, &ld, &wr[1], &wi[1], 00515 &h__[kv + h_dim1], ldh, &nho, &h__[kv + kt * h_dim1], 00516 ldh, &nve, &h__[kwv + h_dim1], ldh, &work[1], lwork); 00517 00518 /* ==== Adjust KBOT accounting for new deflations. ==== */ 00519 00520 kbot -= ld; 00521 00522 /* ==== KS points to the shifts. ==== */ 00523 00524 ks = kbot - ls + 1; 00525 00526 /* ==== Skip an expensive QR sweep if there is a (partly */ 00527 /* . heuristic) reason to expect that many eigenvalues */ 00528 /* . will deflate without it. Here, the QR sweep is */ 00529 /* . skipped if many eigenvalues have just been deflated */ 00530 /* . or if the remaining active block is small. */ 00531 00532 if (ld == 0 || ld * 100 <= nw * nibble && kbot - ktop + 1 > min( 00533 nmin,nwmax)) { 00534 00535 /* ==== NS = nominal number of simultaneous shifts. */ 00536 /* . This may be lowered (slightly) if DLAQR3 */ 00537 /* . did not provide that many shifts. ==== */ 00538 00539 /* Computing MIN */ 00540 /* Computing MAX */ 00541 i__4 = 2, i__5 = kbot - ktop; 00542 i__2 = min(nsmax,nsr), i__3 = max(i__4,i__5); 00543 ns = min(i__2,i__3); 00544 ns -= ns % 2; 00545 00546 /* ==== If there have been no deflations */ 00547 /* . in a multiple of KEXSH iterations, */ 00548 /* . then try exceptional shifts. */ 00549 /* . Otherwise use shifts provided by */ 00550 /* . DLAQR3 above or from the eigenvalues */ 00551 /* . of a trailing principal submatrix. ==== */ 00552 00553 if (ndfl % 6 == 0) { 00554 ks = kbot - ns + 1; 00555 /* Computing MAX */ 00556 i__3 = ks + 1, i__4 = ktop + 2; 00557 i__2 = max(i__3,i__4); 00558 for (i__ = kbot; i__ >= i__2; i__ += -2) { 00559 ss = (d__1 = h__[i__ + (i__ - 1) * h_dim1], abs(d__1)) 00560 + (d__2 = h__[i__ - 1 + (i__ - 2) * h_dim1], 00561 abs(d__2)); 00562 aa = ss * .75 + h__[i__ + i__ * h_dim1]; 00563 bb = ss; 00564 cc = ss * -.4375; 00565 dd = aa; 00566 dlanv2_(&aa, &bb, &cc, &dd, &wr[i__ - 1], &wi[i__ - 1] 00567 , &wr[i__], &wi[i__], &cs, &sn); 00568 /* L30: */ 00569 } 00570 if (ks == ktop) { 00571 wr[ks + 1] = h__[ks + 1 + (ks + 1) * h_dim1]; 00572 wi[ks + 1] = 0.; 00573 wr[ks] = wr[ks + 1]; 00574 wi[ks] = wi[ks + 1]; 00575 } 00576 } else { 00577 00578 /* ==== Got NS/2 or fewer shifts? Use DLAQR4 or */ 00579 /* . DLAHQR on a trailing principal submatrix to */ 00580 /* . get more. (Since NS.LE.NSMAX.LE.(N+6)/9, */ 00581 /* . there is enough space below the subdiagonal */ 00582 /* . to fit an NS-by-NS scratch array.) ==== */ 00583 00584 if (kbot - ks + 1 <= ns / 2) { 00585 ks = kbot - ns + 1; 00586 kt = *n - ns + 1; 00587 dlacpy_("A", &ns, &ns, &h__[ks + ks * h_dim1], ldh, & 00588 h__[kt + h_dim1], ldh); 00589 if (ns > nmin) { 00590 dlaqr4_(&c_false, &c_false, &ns, &c__1, &ns, &h__[ 00591 kt + h_dim1], ldh, &wr[ks], &wi[ks], & 00592 c__1, &c__1, zdum, &c__1, &work[1], lwork, 00593 &inf); 00594 } else { 00595 dlahqr_(&c_false, &c_false, &ns, &c__1, &ns, &h__[ 00596 kt + h_dim1], ldh, &wr[ks], &wi[ks], & 00597 c__1, &c__1, zdum, &c__1, &inf); 00598 } 00599 ks += inf; 00600 00601 /* ==== In case of a rare QR failure use */ 00602 /* . eigenvalues of the trailing 2-by-2 */ 00603 /* . principal submatrix. ==== */ 00604 00605 if (ks >= kbot) { 00606 aa = h__[kbot - 1 + (kbot - 1) * h_dim1]; 00607 cc = h__[kbot + (kbot - 1) * h_dim1]; 00608 bb = h__[kbot - 1 + kbot * h_dim1]; 00609 dd = h__[kbot + kbot * h_dim1]; 00610 dlanv2_(&aa, &bb, &cc, &dd, &wr[kbot - 1], &wi[ 00611 kbot - 1], &wr[kbot], &wi[kbot], &cs, &sn) 00612 ; 00613 ks = kbot - 1; 00614 } 00615 } 00616 00617 if (kbot - ks + 1 > ns) { 00618 00619 /* ==== Sort the shifts (Helps a little) */ 00620 /* . Bubble sort keeps complex conjugate */ 00621 /* . pairs together. ==== */ 00622 00623 sorted = FALSE_; 00624 i__2 = ks + 1; 00625 for (k = kbot; k >= i__2; --k) { 00626 if (sorted) { 00627 goto L60; 00628 } 00629 sorted = TRUE_; 00630 i__3 = k - 1; 00631 for (i__ = ks; i__ <= i__3; ++i__) { 00632 if ((d__1 = wr[i__], abs(d__1)) + (d__2 = wi[ 00633 i__], abs(d__2)) < (d__3 = wr[i__ + 1] 00634 , abs(d__3)) + (d__4 = wi[i__ + 1], 00635 abs(d__4))) { 00636 sorted = FALSE_; 00637 00638 swap = wr[i__]; 00639 wr[i__] = wr[i__ + 1]; 00640 wr[i__ + 1] = swap; 00641 00642 swap = wi[i__]; 00643 wi[i__] = wi[i__ + 1]; 00644 wi[i__ + 1] = swap; 00645 } 00646 /* L40: */ 00647 } 00648 /* L50: */ 00649 } 00650 L60: 00651 ; 00652 } 00653 00654 /* ==== Shuffle shifts into pairs of real shifts */ 00655 /* . and pairs of complex conjugate shifts */ 00656 /* . assuming complex conjugate shifts are */ 00657 /* . already adjacent to one another. (Yes, */ 00658 /* . they are.) ==== */ 00659 00660 i__2 = ks + 2; 00661 for (i__ = kbot; i__ >= i__2; i__ += -2) { 00662 if (wi[i__] != -wi[i__ - 1]) { 00663 00664 swap = wr[i__]; 00665 wr[i__] = wr[i__ - 1]; 00666 wr[i__ - 1] = wr[i__ - 2]; 00667 wr[i__ - 2] = swap; 00668 00669 swap = wi[i__]; 00670 wi[i__] = wi[i__ - 1]; 00671 wi[i__ - 1] = wi[i__ - 2]; 00672 wi[i__ - 2] = swap; 00673 } 00674 /* L70: */ 00675 } 00676 } 00677 00678 /* ==== If there are only two shifts and both are */ 00679 /* . real, then use only one. ==== */ 00680 00681 if (kbot - ks + 1 == 2) { 00682 if (wi[kbot] == 0.) { 00683 if ((d__1 = wr[kbot] - h__[kbot + kbot * h_dim1], abs( 00684 d__1)) < (d__2 = wr[kbot - 1] - h__[kbot + 00685 kbot * h_dim1], abs(d__2))) { 00686 wr[kbot - 1] = wr[kbot]; 00687 } else { 00688 wr[kbot] = wr[kbot - 1]; 00689 } 00690 } 00691 } 00692 00693 /* ==== Use up to NS of the the smallest magnatiude */ 00694 /* . shifts. If there aren't NS shifts available, */ 00695 /* . then use them all, possibly dropping one to */ 00696 /* . make the number of shifts even. ==== */ 00697 00698 /* Computing MIN */ 00699 i__2 = ns, i__3 = kbot - ks + 1; 00700 ns = min(i__2,i__3); 00701 ns -= ns % 2; 00702 ks = kbot - ns + 1; 00703 00704 /* ==== Small-bulge multi-shift QR sweep: */ 00705 /* . split workspace under the subdiagonal into */ 00706 /* . - a KDU-by-KDU work array U in the lower */ 00707 /* . left-hand-corner, */ 00708 /* . - a KDU-by-at-least-KDU-but-more-is-better */ 00709 /* . (KDU-by-NHo) horizontal work array WH along */ 00710 /* . the bottom edge, */ 00711 /* . - and an at-least-KDU-but-more-is-better-by-KDU */ 00712 /* . (NVE-by-KDU) vertical work WV arrow along */ 00713 /* . the left-hand-edge. ==== */ 00714 00715 kdu = ns * 3 - 3; 00716 ku = *n - kdu + 1; 00717 kwh = kdu + 1; 00718 nho = *n - kdu - 3 - (kdu + 1) + 1; 00719 kwv = kdu + 4; 00720 nve = *n - kdu - kwv + 1; 00721 00722 /* ==== Small-bulge multi-shift QR sweep ==== */ 00723 00724 dlaqr5_(wantt, wantz, &kacc22, n, &ktop, &kbot, &ns, &wr[ks], 00725 &wi[ks], &h__[h_offset], ldh, iloz, ihiz, &z__[ 00726 z_offset], ldz, &work[1], &c__3, &h__[ku + h_dim1], 00727 ldh, &nve, &h__[kwv + h_dim1], ldh, &nho, &h__[ku + 00728 kwh * h_dim1], ldh); 00729 } 00730 00731 /* ==== Note progress (or the lack of it). ==== */ 00732 00733 if (ld > 0) { 00734 ndfl = 1; 00735 } else { 00736 ++ndfl; 00737 } 00738 00739 /* ==== End of main loop ==== */ 00740 /* L80: */ 00741 } 00742 00743 /* ==== Iteration limit exceeded. Set INFO to show where */ 00744 /* . the problem occurred and exit. ==== */ 00745 00746 *info = kbot; 00747 L90: 00748 ; 00749 } 00750 00751 /* ==== Return the optimal value of LWORK. ==== */ 00752 00753 work[1] = (doublereal) lwkopt; 00754 00755 /* ==== End of DLAQR0 ==== */ 00756 00757 return 0; 00758 } /* dlaqr0_ */