00001 /* dhseqr.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 doublereal c_b11 = 0.; 00019 static doublereal c_b12 = 1.; 00020 static integer c__12 = 12; 00021 static integer c__2 = 2; 00022 static integer c__49 = 49; 00023 00024 /* Subroutine */ int dhseqr_(char *job, char *compz, integer *n, integer *ilo, 00025 integer *ihi, doublereal *h__, integer *ldh, doublereal *wr, 00026 doublereal *wi, doublereal *z__, integer *ldz, doublereal *work, 00027 integer *lwork, integer *info) 00028 { 00029 /* System generated locals */ 00030 address a__1[2]; 00031 integer h_dim1, h_offset, z_dim1, z_offset, i__1, i__2[2], i__3; 00032 doublereal d__1; 00033 char ch__1[2]; 00034 00035 /* Builtin functions */ 00036 /* Subroutine */ int s_cat(char *, char **, integer *, integer *, ftnlen); 00037 00038 /* Local variables */ 00039 integer i__; 00040 doublereal hl[2401] /* was [49][49] */; 00041 integer kbot, nmin; 00042 extern logical lsame_(char *, char *); 00043 logical initz; 00044 doublereal workl[49]; 00045 logical wantt, wantz; 00046 extern /* Subroutine */ int dlaqr0_(logical *, logical *, integer *, 00047 integer *, integer *, doublereal *, integer *, doublereal *, 00048 doublereal *, integer *, integer *, doublereal *, integer *, 00049 doublereal *, integer *, integer *), dlahqr_(logical *, logical *, 00050 integer *, integer *, integer *, doublereal *, integer *, 00051 doublereal *, doublereal *, integer *, integer *, doublereal *, 00052 integer *, integer *), dlacpy_(char *, integer *, integer *, 00053 doublereal *, integer *, doublereal *, integer *), 00054 dlaset_(char *, integer *, integer *, doublereal *, doublereal *, 00055 doublereal *, integer *); 00056 extern integer ilaenv_(integer *, char *, char *, integer *, integer *, 00057 integer *, integer *); 00058 extern /* Subroutine */ int xerbla_(char *, integer *); 00059 logical lquery; 00060 00061 00062 /* -- LAPACK driver routine (version 3.2) -- */ 00063 /* Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd.. */ 00064 /* November 2006 */ 00065 00066 /* .. Scalar Arguments .. */ 00067 /* .. */ 00068 /* .. Array Arguments .. */ 00069 /* .. */ 00070 /* Purpose */ 00071 /* ======= */ 00072 00073 /* DHSEQR computes the eigenvalues of a Hessenberg matrix H */ 00074 /* and, optionally, the matrices T and Z from the Schur decomposition */ 00075 /* H = Z T Z**T, where T is an upper quasi-triangular matrix (the */ 00076 /* Schur form), and Z is the orthogonal matrix of Schur vectors. */ 00077 00078 /* Optionally Z may be postmultiplied into an input orthogonal */ 00079 /* matrix Q so that this routine can give the Schur factorization */ 00080 /* of a matrix A which has been reduced to the Hessenberg form H */ 00081 /* by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T. */ 00082 00083 /* Arguments */ 00084 /* ========= */ 00085 00086 /* JOB (input) CHARACTER*1 */ 00087 /* = 'E': compute eigenvalues only; */ 00088 /* = 'S': compute eigenvalues and the Schur form T. */ 00089 00090 /* COMPZ (input) CHARACTER*1 */ 00091 /* = 'N': no Schur vectors are computed; */ 00092 /* = 'I': Z is initialized to the unit matrix and the matrix Z */ 00093 /* of Schur vectors of H is returned; */ 00094 /* = 'V': Z must contain an orthogonal matrix Q on entry, and */ 00095 /* the product Q*Z is returned. */ 00096 00097 /* N (input) INTEGER */ 00098 /* The order of the matrix H. N .GE. 0. */ 00099 00100 /* ILO (input) INTEGER */ 00101 /* IHI (input) INTEGER */ 00102 /* It is assumed that H is already upper triangular in rows */ 00103 /* and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally */ 00104 /* set by a previous call to DGEBAL, and then passed to DGEHRD */ 00105 /* when the matrix output by DGEBAL is reduced to Hessenberg */ 00106 /* form. Otherwise ILO and IHI should be set to 1 and N */ 00107 /* respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N. */ 00108 /* If N = 0, then ILO = 1 and IHI = 0. */ 00109 00110 /* H (input/output) DOUBLE PRECISION array, dimension (LDH,N) */ 00111 /* On entry, the upper Hessenberg matrix H. */ 00112 /* On exit, if INFO = 0 and JOB = 'S', then H contains the */ 00113 /* upper quasi-triangular matrix T from the Schur decomposition */ 00114 /* (the Schur form); 2-by-2 diagonal blocks (corresponding to */ 00115 /* complex conjugate pairs of eigenvalues) are returned in */ 00116 /* standard form, with H(i,i) = H(i+1,i+1) and */ 00117 /* H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and JOB = 'E', the */ 00118 /* contents of H are unspecified on exit. (The output value of */ 00119 /* H when INFO.GT.0 is given under the description of INFO */ 00120 /* below.) */ 00121 00122 /* Unlike earlier versions of DHSEQR, this subroutine may */ 00123 /* explicitly H(i,j) = 0 for i.GT.j and j = 1, 2, ... ILO-1 */ 00124 /* or j = IHI+1, IHI+2, ... N. */ 00125 00126 /* LDH (input) INTEGER */ 00127 /* The leading dimension of the array H. LDH .GE. max(1,N). */ 00128 00129 /* WR (output) DOUBLE PRECISION array, dimension (N) */ 00130 /* WI (output) DOUBLE PRECISION array, dimension (N) */ 00131 /* The real and imaginary parts, respectively, of the computed */ 00132 /* eigenvalues. If two eigenvalues are computed as a complex */ 00133 /* conjugate pair, they are stored in consecutive elements of */ 00134 /* WR and WI, say the i-th and (i+1)th, with WI(i) .GT. 0 and */ 00135 /* WI(i+1) .LT. 0. If JOB = 'S', the eigenvalues are stored in */ 00136 /* the same order as on the diagonal of the Schur form returned */ 00137 /* in H, with WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 */ 00138 /* diagonal block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and */ 00139 /* WI(i+1) = -WI(i). */ 00140 00141 /* Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N) */ 00142 /* If COMPZ = 'N', Z is not referenced. */ 00143 /* If COMPZ = 'I', on entry Z need not be set and on exit, */ 00144 /* if INFO = 0, Z contains the orthogonal matrix Z of the Schur */ 00145 /* vectors of H. If COMPZ = 'V', on entry Z must contain an */ 00146 /* N-by-N matrix Q, which is assumed to be equal to the unit */ 00147 /* matrix except for the submatrix Z(ILO:IHI,ILO:IHI). On exit, */ 00148 /* if INFO = 0, Z contains Q*Z. */ 00149 /* Normally Q is the orthogonal matrix generated by DORGHR */ 00150 /* after the call to DGEHRD which formed the Hessenberg matrix */ 00151 /* H. (The output value of Z when INFO.GT.0 is given under */ 00152 /* the description of INFO below.) */ 00153 00154 /* LDZ (input) INTEGER */ 00155 /* The leading dimension of the array Z. if COMPZ = 'I' or */ 00156 /* COMPZ = 'V', then LDZ.GE.MAX(1,N). Otherwize, LDZ.GE.1. */ 00157 00158 /* WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) */ 00159 /* On exit, if INFO = 0, WORK(1) returns an estimate of */ 00160 /* the optimal value for LWORK. */ 00161 00162 /* LWORK (input) INTEGER */ 00163 /* The dimension of the array WORK. LWORK .GE. max(1,N) */ 00164 /* is sufficient and delivers very good and sometimes */ 00165 /* optimal performance. However, LWORK as large as 11*N */ 00166 /* may be required for optimal performance. A workspace */ 00167 /* query is recommended to determine the optimal workspace */ 00168 /* size. */ 00169 00170 /* If LWORK = -1, then DHSEQR does a workspace query. */ 00171 /* In this case, DHSEQR checks the input parameters and */ 00172 /* estimates the optimal workspace size for the given */ 00173 /* values of N, ILO and IHI. The estimate is returned */ 00174 /* in WORK(1). No error message related to LWORK is */ 00175 /* issued by XERBLA. Neither H nor Z are accessed. */ 00176 00177 00178 /* INFO (output) INTEGER */ 00179 /* = 0: successful exit */ 00180 /* .LT. 0: if INFO = -i, the i-th argument had an illegal */ 00181 /* value */ 00182 /* .GT. 0: if INFO = i, DHSEQR failed to compute all of */ 00183 /* the eigenvalues. Elements 1:ilo-1 and i+1:n of WR */ 00184 /* and WI contain those eigenvalues which have been */ 00185 /* successfully computed. (Failures are rare.) */ 00186 00187 /* If INFO .GT. 0 and JOB = 'E', then on exit, the */ 00188 /* remaining unconverged eigenvalues are the eigen- */ 00189 /* values of the upper Hessenberg matrix rows and */ 00190 /* columns ILO through INFO of the final, output */ 00191 /* value of H. */ 00192 00193 /* If INFO .GT. 0 and JOB = 'S', then on exit */ 00194 00195 /* (*) (initial value of H)*U = U*(final value of H) */ 00196 00197 /* where U is an orthogonal matrix. The final */ 00198 /* value of H is upper Hessenberg and quasi-triangular */ 00199 /* in rows and columns INFO+1 through IHI. */ 00200 00201 /* If INFO .GT. 0 and COMPZ = 'V', then on exit */ 00202 00203 /* (final value of Z) = (initial value of Z)*U */ 00204 00205 /* where U is the orthogonal matrix in (*) (regard- */ 00206 /* less of the value of JOB.) */ 00207 00208 /* If INFO .GT. 0 and COMPZ = 'I', then on exit */ 00209 /* (final value of Z) = U */ 00210 /* where U is the orthogonal matrix in (*) (regard- */ 00211 /* less of the value of JOB.) */ 00212 00213 /* If INFO .GT. 0 and COMPZ = 'N', then Z is not */ 00214 /* accessed. */ 00215 00216 /* ================================================================ */ 00217 /* Default values supplied by */ 00218 /* ILAENV(ISPEC,'DHSEQR',JOB(:1)//COMPZ(:1),N,ILO,IHI,LWORK). */ 00219 /* It is suggested that these defaults be adjusted in order */ 00220 /* to attain best performance in each particular */ 00221 /* computational environment. */ 00222 00223 /* ISPEC=12: The DLAHQR vs DLAQR0 crossover point. */ 00224 /* Default: 75. (Must be at least 11.) */ 00225 00226 /* ISPEC=13: Recommended deflation window size. */ 00227 /* This depends on ILO, IHI and NS. NS is the */ 00228 /* number of simultaneous shifts returned */ 00229 /* by ILAENV(ISPEC=15). (See ISPEC=15 below.) */ 00230 /* The default for (IHI-ILO+1).LE.500 is NS. */ 00231 /* The default for (IHI-ILO+1).GT.500 is 3*NS/2. */ 00232 00233 /* ISPEC=14: Nibble crossover point. (See IPARMQ for */ 00234 /* details.) Default: 14% of deflation window */ 00235 /* size. */ 00236 00237 /* ISPEC=15: Number of simultaneous shifts in a multishift */ 00238 /* QR iteration. */ 00239 00240 /* If IHI-ILO+1 is ... */ 00241 00242 /* greater than ...but less ... the */ 00243 /* or equal to ... than default is */ 00244 00245 /* 1 30 NS = 2(+) */ 00246 /* 30 60 NS = 4(+) */ 00247 /* 60 150 NS = 10(+) */ 00248 /* 150 590 NS = ** */ 00249 /* 590 3000 NS = 64 */ 00250 /* 3000 6000 NS = 128 */ 00251 /* 6000 infinity NS = 256 */ 00252 00253 /* (+) By default some or all matrices of this order */ 00254 /* are passed to the implicit double shift routine */ 00255 /* DLAHQR and this parameter is ignored. See */ 00256 /* ISPEC=12 above and comments in IPARMQ for */ 00257 /* details. */ 00258 00259 /* (**) The asterisks (**) indicate an ad-hoc */ 00260 /* function of N increasing from 10 to 64. */ 00261 00262 /* ISPEC=16: Select structured matrix multiply. */ 00263 /* If the number of simultaneous shifts (specified */ 00264 /* by ISPEC=15) is less than 14, then the default */ 00265 /* for ISPEC=16 is 0. Otherwise the default for */ 00266 /* ISPEC=16 is 2. */ 00267 00268 /* ================================================================ */ 00269 /* Based on contributions by */ 00270 /* Karen Braman and Ralph Byers, Department of Mathematics, */ 00271 /* University of Kansas, USA */ 00272 00273 /* ================================================================ */ 00274 /* References: */ 00275 /* K. Braman, R. Byers and R. Mathias, The Multi-Shift QR */ 00276 /* Algorithm Part I: Maintaining Well Focused Shifts, and Level 3 */ 00277 /* Performance, SIAM Journal of Matrix Analysis, volume 23, pages */ 00278 /* 929--947, 2002. */ 00279 00280 /* K. Braman, R. Byers and R. Mathias, The Multi-Shift QR */ 00281 /* Algorithm Part II: Aggressive Early Deflation, SIAM Journal */ 00282 /* of Matrix Analysis, volume 23, pages 948--973, 2002. */ 00283 00284 /* ================================================================ */ 00285 /* .. Parameters .. */ 00286 00287 /* ==== Matrices of order NTINY or smaller must be processed by */ 00288 /* . DLAHQR because of insufficient subdiagonal scratch space. */ 00289 /* . (This is a hard limit.) ==== */ 00290 00291 /* ==== NL allocates some local workspace to help small matrices */ 00292 /* . through a rare DLAHQR failure. NL .GT. NTINY = 11 is */ 00293 /* . required and NL .LE. NMIN = ILAENV(ISPEC=12,...) is recom- */ 00294 /* . mended. (The default value of NMIN is 75.) Using NL = 49 */ 00295 /* . allows up to six simultaneous shifts and a 16-by-16 */ 00296 /* . deflation window. ==== */ 00297 /* .. */ 00298 /* .. Local Arrays .. */ 00299 /* .. */ 00300 /* .. Local Scalars .. */ 00301 /* .. */ 00302 /* .. External Functions .. */ 00303 /* .. */ 00304 /* .. External Subroutines .. */ 00305 /* .. */ 00306 /* .. Intrinsic Functions .. */ 00307 /* .. */ 00308 /* .. Executable Statements .. */ 00309 00310 /* ==== Decode and check the input parameters. ==== */ 00311 00312 /* Parameter adjustments */ 00313 h_dim1 = *ldh; 00314 h_offset = 1 + h_dim1; 00315 h__ -= h_offset; 00316 --wr; 00317 --wi; 00318 z_dim1 = *ldz; 00319 z_offset = 1 + z_dim1; 00320 z__ -= z_offset; 00321 --work; 00322 00323 /* Function Body */ 00324 wantt = lsame_(job, "S"); 00325 initz = lsame_(compz, "I"); 00326 wantz = initz || lsame_(compz, "V"); 00327 work[1] = (doublereal) max(1,*n); 00328 lquery = *lwork == -1; 00329 00330 *info = 0; 00331 if (! lsame_(job, "E") && ! wantt) { 00332 *info = -1; 00333 } else if (! lsame_(compz, "N") && ! wantz) { 00334 *info = -2; 00335 } else if (*n < 0) { 00336 *info = -3; 00337 } else if (*ilo < 1 || *ilo > max(1,*n)) { 00338 *info = -4; 00339 } else if (*ihi < min(*ilo,*n) || *ihi > *n) { 00340 *info = -5; 00341 } else if (*ldh < max(1,*n)) { 00342 *info = -7; 00343 } else if (*ldz < 1 || wantz && *ldz < max(1,*n)) { 00344 *info = -11; 00345 } else if (*lwork < max(1,*n) && ! lquery) { 00346 *info = -13; 00347 } 00348 00349 if (*info != 0) { 00350 00351 /* ==== Quick return in case of invalid argument. ==== */ 00352 00353 i__1 = -(*info); 00354 xerbla_("DHSEQR", &i__1); 00355 return 0; 00356 00357 } else if (*n == 0) { 00358 00359 /* ==== Quick return in case N = 0; nothing to do. ==== */ 00360 00361 return 0; 00362 00363 } else if (lquery) { 00364 00365 /* ==== Quick return in case of a workspace query ==== */ 00366 00367 dlaqr0_(&wantt, &wantz, n, ilo, ihi, &h__[h_offset], ldh, &wr[1], &wi[ 00368 1], ilo, ihi, &z__[z_offset], ldz, &work[1], lwork, info); 00369 /* ==== Ensure reported workspace size is backward-compatible with */ 00370 /* . previous LAPACK versions. ==== */ 00371 /* Computing MAX */ 00372 d__1 = (doublereal) max(1,*n); 00373 work[1] = max(d__1,work[1]); 00374 return 0; 00375 00376 } else { 00377 00378 /* ==== copy eigenvalues isolated by DGEBAL ==== */ 00379 00380 i__1 = *ilo - 1; 00381 for (i__ = 1; i__ <= i__1; ++i__) { 00382 wr[i__] = h__[i__ + i__ * h_dim1]; 00383 wi[i__] = 0.; 00384 /* L10: */ 00385 } 00386 i__1 = *n; 00387 for (i__ = *ihi + 1; i__ <= i__1; ++i__) { 00388 wr[i__] = h__[i__ + i__ * h_dim1]; 00389 wi[i__] = 0.; 00390 /* L20: */ 00391 } 00392 00393 /* ==== Initialize Z, if requested ==== */ 00394 00395 if (initz) { 00396 dlaset_("A", n, n, &c_b11, &c_b12, &z__[z_offset], ldz) 00397 ; 00398 } 00399 00400 /* ==== Quick return if possible ==== */ 00401 00402 if (*ilo == *ihi) { 00403 wr[*ilo] = h__[*ilo + *ilo * h_dim1]; 00404 wi[*ilo] = 0.; 00405 return 0; 00406 } 00407 00408 /* ==== DLAHQR/DLAQR0 crossover point ==== */ 00409 00410 /* Writing concatenation */ 00411 i__2[0] = 1, a__1[0] = job; 00412 i__2[1] = 1, a__1[1] = compz; 00413 s_cat(ch__1, a__1, i__2, &c__2, (ftnlen)2); 00414 nmin = ilaenv_(&c__12, "DHSEQR", ch__1, n, ilo, ihi, lwork); 00415 nmin = max(11,nmin); 00416 00417 /* ==== DLAQR0 for big matrices; DLAHQR for small ones ==== */ 00418 00419 if (*n > nmin) { 00420 dlaqr0_(&wantt, &wantz, n, ilo, ihi, &h__[h_offset], ldh, &wr[1], 00421 &wi[1], ilo, ihi, &z__[z_offset], ldz, &work[1], lwork, 00422 info); 00423 } else { 00424 00425 /* ==== Small matrix ==== */ 00426 00427 dlahqr_(&wantt, &wantz, n, ilo, ihi, &h__[h_offset], ldh, &wr[1], 00428 &wi[1], ilo, ihi, &z__[z_offset], ldz, info); 00429 00430 if (*info > 0) { 00431 00432 /* ==== A rare DLAHQR failure! DLAQR0 sometimes succeeds */ 00433 /* . when DLAHQR fails. ==== */ 00434 00435 kbot = *info; 00436 00437 if (*n >= 49) { 00438 00439 /* ==== Larger matrices have enough subdiagonal scratch */ 00440 /* . space to call DLAQR0 directly. ==== */ 00441 00442 dlaqr0_(&wantt, &wantz, n, ilo, &kbot, &h__[h_offset], 00443 ldh, &wr[1], &wi[1], ilo, ihi, &z__[z_offset], 00444 ldz, &work[1], lwork, info); 00445 00446 } else { 00447 00448 /* ==== Tiny matrices don't have enough subdiagonal */ 00449 /* . scratch space to benefit from DLAQR0. Hence, */ 00450 /* . tiny matrices must be copied into a larger */ 00451 /* . array before calling DLAQR0. ==== */ 00452 00453 dlacpy_("A", n, n, &h__[h_offset], ldh, hl, &c__49); 00454 hl[*n + 1 + *n * 49 - 50] = 0.; 00455 i__1 = 49 - *n; 00456 dlaset_("A", &c__49, &i__1, &c_b11, &c_b11, &hl[(*n + 1) * 00457 49 - 49], &c__49); 00458 dlaqr0_(&wantt, &wantz, &c__49, ilo, &kbot, hl, &c__49, & 00459 wr[1], &wi[1], ilo, ihi, &z__[z_offset], ldz, 00460 workl, &c__49, info); 00461 if (wantt || *info != 0) { 00462 dlacpy_("A", n, n, hl, &c__49, &h__[h_offset], ldh); 00463 } 00464 } 00465 } 00466 } 00467 00468 /* ==== Clear out the trash, if necessary. ==== */ 00469 00470 if ((wantt || *info != 0) && *n > 2) { 00471 i__1 = *n - 2; 00472 i__3 = *n - 2; 00473 dlaset_("L", &i__1, &i__3, &c_b11, &c_b11, &h__[h_dim1 + 3], ldh); 00474 } 00475 00476 /* ==== Ensure reported workspace size is backward-compatible with */ 00477 /* . previous LAPACK versions. ==== */ 00478 00479 /* Computing MAX */ 00480 d__1 = (doublereal) max(1,*n); 00481 work[1] = max(d__1,work[1]); 00482 } 00483 00484 /* ==== End of DHSEQR ==== */ 00485 00486 return 0; 00487 } /* dhseqr_ */