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