00001 /* slaqr4.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 slaqr4_(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 slaqr2_(logical *, logical *, integer *, 00053 integer *, integer *, integer *, real *, integer *, integer *, 00054 integer *, real *, integer *, integer *, integer *, real *, real * 00055 , real *, integer *, integer *, real *, integer *, integer *, 00056 real *, integer *, real *, integer *), slanv2_(real *, real *, 00057 real *, real *, real *, real *, real *, real *, real *, real *), 00058 slaqr5_(logical *, logical *, integer *, integer *, integer *, 00059 integer *, integer *, real *, real *, real *, integer *, integer * 00060 , integer *, real *, integer *, real *, integer *, real *, 00061 integer *, integer *, real *, integer *, integer *, real *, 00062 integer *); 00063 integer nibble; 00064 extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 00065 integer *, integer *); 00066 char jbcmpz[2]; 00067 extern /* Subroutine */ int slahqr_(logical *, logical *, integer *, 00068 integer *, integer *, real *, integer *, real *, real *, integer * 00069 , integer *, real *, integer *, integer *), slacpy_(char *, 00070 integer *, integer *, real *, integer *, real *, integer *); 00071 integer nwupbd; 00072 logical sorted; 00073 integer lwkopt; 00074 00075 00076 /* -- LAPACK auxiliary routine (version 3.2) -- */ 00077 /* Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd.. */ 00078 /* November 2006 */ 00079 00080 /* .. Scalar Arguments .. */ 00081 /* .. */ 00082 /* .. Array Arguments .. */ 00083 /* .. */ 00084 00085 /* This subroutine implements one level of recursion for SLAQR0. */ 00086 /* It is a complete implementation of the small bulge multi-shift */ 00087 /* QR algorithm. It may be called by SLAQR0 and, for large enough */ 00088 /* deflation window size, it may be called by SLAQR3. This */ 00089 /* subroutine is identical to SLAQR0 except that it calls SLAQR2 */ 00090 /* instead of SLAQR3. */ 00091 00092 /* Purpose */ 00093 /* ======= */ 00094 00095 /* SLAQR4 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 SGEBAL, and then passed to SGEHRD when the */ 00125 /* matrix output by SGEBAL 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) REAL 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) REAL array, dimension (IHI) */ 00149 /* WI (output) REAL 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) REAL 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) REAL 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 SLAQR4 does a workspace query. */ 00191 /* In this case, SLAQR4 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, SLAQR4 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 /* . SLAHQR 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.f; 00293 return 0; 00294 } 00295 00296 if (*n <= 11) { 00297 00298 /* ==== Tiny matrices must use SLAHQR. ==== */ 00299 00300 lwkopt = 1; 00301 if (*lwork != -1) { 00302 slahqr_(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, "SLAQR4", 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, "SLAQR4", 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 SLAQR2 ==== */ 00356 00357 i__1 = nwr + 1; 00358 slaqr2_(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(SLAQR5, SLAQR2) ==== */ 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] = (real) lwkopt; 00373 return 0; 00374 } 00375 00376 /* ==== SLAHQR/SLAQR0 crossover point ==== */ 00377 00378 nmin = ilaenv_(&c__12, "SLAQR4", jbcmpz, n, ilo, ihi, lwork); 00379 nmin = max(11,nmin); 00380 00381 /* ==== Nibble crossover point ==== */ 00382 00383 nibble = ilaenv_(&c__14, "SLAQR4", 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, "SLAQR4", 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.f) { 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 ((r__1 = h__[kwtop + (kwtop - 1) * h_dim1], dabs(r__1)) 00478 > (r__2 = h__[kwtop - 1 + (kwtop - 2) * h_dim1], 00479 dabs(r__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 slaqr2_(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 SLAQR2 */ 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 /* . SLAQR2 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 = (r__1 = h__[i__ + (i__ - 1) * h_dim1], dabs(r__1) 00560 ) + (r__2 = h__[i__ - 1 + (i__ - 2) * h_dim1], 00561 dabs(r__2)); 00562 aa = ss * .75f + h__[i__ + i__ * h_dim1]; 00563 bb = ss; 00564 cc = ss * -.4375f; 00565 dd = aa; 00566 slanv2_(&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.f; 00573 wr[ks] = wr[ks + 1]; 00574 wi[ks] = wi[ks + 1]; 00575 } 00576 } else { 00577 00578 /* ==== Got NS/2 or fewer shifts? Use SLAHQR */ 00579 /* . 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 slacpy_("A", &ns, &ns, &h__[ks + ks * h_dim1], ldh, & 00588 h__[kt + h_dim1], ldh); 00589 slahqr_(&c_false, &c_false, &ns, &c__1, &ns, &h__[kt 00590 + h_dim1], ldh, &wr[ks], &wi[ks], &c__1, & 00591 c__1, zdum, &c__1, &inf); 00592 ks += inf; 00593 00594 /* ==== In case of a rare QR failure use */ 00595 /* . eigenvalues of the trailing 2-by-2 */ 00596 /* . principal submatrix. ==== */ 00597 00598 if (ks >= kbot) { 00599 aa = h__[kbot - 1 + (kbot - 1) * h_dim1]; 00600 cc = h__[kbot + (kbot - 1) * h_dim1]; 00601 bb = h__[kbot - 1 + kbot * h_dim1]; 00602 dd = h__[kbot + kbot * h_dim1]; 00603 slanv2_(&aa, &bb, &cc, &dd, &wr[kbot - 1], &wi[ 00604 kbot - 1], &wr[kbot], &wi[kbot], &cs, &sn) 00605 ; 00606 ks = kbot - 1; 00607 } 00608 } 00609 00610 if (kbot - ks + 1 > ns) { 00611 00612 /* ==== Sort the shifts (Helps a little) */ 00613 /* . Bubble sort keeps complex conjugate */ 00614 /* . pairs together. ==== */ 00615 00616 sorted = FALSE_; 00617 i__2 = ks + 1; 00618 for (k = kbot; k >= i__2; --k) { 00619 if (sorted) { 00620 goto L60; 00621 } 00622 sorted = TRUE_; 00623 i__3 = k - 1; 00624 for (i__ = ks; i__ <= i__3; ++i__) { 00625 if ((r__1 = wr[i__], dabs(r__1)) + (r__2 = wi[ 00626 i__], dabs(r__2)) < (r__3 = wr[i__ + 00627 1], dabs(r__3)) + (r__4 = wi[i__ + 1], 00628 dabs(r__4))) { 00629 sorted = FALSE_; 00630 00631 swap = wr[i__]; 00632 wr[i__] = wr[i__ + 1]; 00633 wr[i__ + 1] = swap; 00634 00635 swap = wi[i__]; 00636 wi[i__] = wi[i__ + 1]; 00637 wi[i__ + 1] = swap; 00638 } 00639 /* L40: */ 00640 } 00641 /* L50: */ 00642 } 00643 L60: 00644 ; 00645 } 00646 00647 /* ==== Shuffle shifts into pairs of real shifts */ 00648 /* . and pairs of complex conjugate shifts */ 00649 /* . assuming complex conjugate shifts are */ 00650 /* . already adjacent to one another. (Yes, */ 00651 /* . they are.) ==== */ 00652 00653 i__2 = ks + 2; 00654 for (i__ = kbot; i__ >= i__2; i__ += -2) { 00655 if (wi[i__] != -wi[i__ - 1]) { 00656 00657 swap = wr[i__]; 00658 wr[i__] = wr[i__ - 1]; 00659 wr[i__ - 1] = wr[i__ - 2]; 00660 wr[i__ - 2] = swap; 00661 00662 swap = wi[i__]; 00663 wi[i__] = wi[i__ - 1]; 00664 wi[i__ - 1] = wi[i__ - 2]; 00665 wi[i__ - 2] = swap; 00666 } 00667 /* L70: */ 00668 } 00669 } 00670 00671 /* ==== If there are only two shifts and both are */ 00672 /* . real, then use only one. ==== */ 00673 00674 if (kbot - ks + 1 == 2) { 00675 if (wi[kbot] == 0.f) { 00676 if ((r__1 = wr[kbot] - h__[kbot + kbot * h_dim1], 00677 dabs(r__1)) < (r__2 = wr[kbot - 1] - h__[kbot 00678 + kbot * h_dim1], dabs(r__2))) { 00679 wr[kbot - 1] = wr[kbot]; 00680 } else { 00681 wr[kbot] = wr[kbot - 1]; 00682 } 00683 } 00684 } 00685 00686 /* ==== Use up to NS of the the smallest magnatiude */ 00687 /* . shifts. If there aren't NS shifts available, */ 00688 /* . then use them all, possibly dropping one to */ 00689 /* . make the number of shifts even. ==== */ 00690 00691 /* Computing MIN */ 00692 i__2 = ns, i__3 = kbot - ks + 1; 00693 ns = min(i__2,i__3); 00694 ns -= ns % 2; 00695 ks = kbot - ns + 1; 00696 00697 /* ==== Small-bulge multi-shift QR sweep: */ 00698 /* . split workspace under the subdiagonal into */ 00699 /* . - a KDU-by-KDU work array U in the lower */ 00700 /* . left-hand-corner, */ 00701 /* . - a KDU-by-at-least-KDU-but-more-is-better */ 00702 /* . (KDU-by-NHo) horizontal work array WH along */ 00703 /* . the bottom edge, */ 00704 /* . - and an at-least-KDU-but-more-is-better-by-KDU */ 00705 /* . (NVE-by-KDU) vertical work WV arrow along */ 00706 /* . the left-hand-edge. ==== */ 00707 00708 kdu = ns * 3 - 3; 00709 ku = *n - kdu + 1; 00710 kwh = kdu + 1; 00711 nho = *n - kdu - 3 - (kdu + 1) + 1; 00712 kwv = kdu + 4; 00713 nve = *n - kdu - kwv + 1; 00714 00715 /* ==== Small-bulge multi-shift QR sweep ==== */ 00716 00717 slaqr5_(wantt, wantz, &kacc22, n, &ktop, &kbot, &ns, &wr[ks], 00718 &wi[ks], &h__[h_offset], ldh, iloz, ihiz, &z__[ 00719 z_offset], ldz, &work[1], &c__3, &h__[ku + h_dim1], 00720 ldh, &nve, &h__[kwv + h_dim1], ldh, &nho, &h__[ku + 00721 kwh * h_dim1], ldh); 00722 } 00723 00724 /* ==== Note progress (or the lack of it). ==== */ 00725 00726 if (ld > 0) { 00727 ndfl = 1; 00728 } else { 00729 ++ndfl; 00730 } 00731 00732 /* ==== End of main loop ==== */ 00733 /* L80: */ 00734 } 00735 00736 /* ==== Iteration limit exceeded. Set INFO to show where */ 00737 /* . the problem occurred and exit. ==== */ 00738 00739 *info = kbot; 00740 L90: 00741 ; 00742 } 00743 00744 /* ==== Return the optimal value of LWORK. ==== */ 00745 00746 work[1] = (real) lwkopt; 00747 00748 /* ==== End of SLAQR4 ==== */ 00749 00750 return 0; 00751 } /* slaqr4_ */