00001 /* chsein.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 logical c_false = FALSE_; 00019 static logical c_true = TRUE_; 00020 00021 /* Subroutine */ int chsein_(char *side, char *eigsrc, char *initv, logical * 00022 select, integer *n, complex *h__, integer *ldh, complex *w, complex * 00023 vl, integer *ldvl, complex *vr, integer *ldvr, integer *mm, integer * 00024 m, complex *work, real *rwork, integer *ifaill, integer *ifailr, 00025 integer *info) 00026 { 00027 /* System generated locals */ 00028 integer h_dim1, h_offset, vl_dim1, vl_offset, vr_dim1, vr_offset, i__1, 00029 i__2, i__3; 00030 real r__1, r__2; 00031 complex q__1, q__2; 00032 00033 /* Builtin functions */ 00034 double r_imag(complex *); 00035 00036 /* Local variables */ 00037 integer i__, k, kl, kr, ks; 00038 complex wk; 00039 integer kln; 00040 real ulp, eps3, unfl; 00041 extern logical lsame_(char *, char *); 00042 integer iinfo; 00043 logical leftv, bothv; 00044 real hnorm; 00045 extern /* Subroutine */ int claein_(logical *, logical *, integer *, 00046 complex *, integer *, complex *, complex *, complex *, integer *, 00047 real *, real *, real *, integer *); 00048 extern doublereal slamch_(char *), clanhs_(char *, integer *, 00049 complex *, integer *, real *); 00050 extern /* Subroutine */ int xerbla_(char *, integer *); 00051 logical noinit; 00052 integer ldwork; 00053 logical rightv, fromqr; 00054 real smlnum; 00055 00056 00057 /* -- LAPACK routine (version 3.2) -- */ 00058 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00059 /* November 2006 */ 00060 00061 /* .. Scalar Arguments .. */ 00062 /* .. */ 00063 /* .. Array Arguments .. */ 00064 /* .. */ 00065 00066 /* Purpose */ 00067 /* ======= */ 00068 00069 /* CHSEIN uses inverse iteration to find specified right and/or left */ 00070 /* eigenvectors of a complex upper Hessenberg matrix H. */ 00071 00072 /* The right eigenvector x and the left eigenvector y of the matrix H */ 00073 /* corresponding to an eigenvalue w are defined by: */ 00074 00075 /* H * x = w * x, y**h * H = w * y**h */ 00076 00077 /* where y**h denotes the conjugate transpose of the vector y. */ 00078 00079 /* Arguments */ 00080 /* ========= */ 00081 00082 /* SIDE (input) CHARACTER*1 */ 00083 /* = 'R': compute right eigenvectors only; */ 00084 /* = 'L': compute left eigenvectors only; */ 00085 /* = 'B': compute both right and left eigenvectors. */ 00086 00087 /* EIGSRC (input) CHARACTER*1 */ 00088 /* Specifies the source of eigenvalues supplied in W: */ 00089 /* = 'Q': the eigenvalues were found using CHSEQR; thus, if */ 00090 /* H has zero subdiagonal elements, and so is */ 00091 /* block-triangular, then the j-th eigenvalue can be */ 00092 /* assumed to be an eigenvalue of the block containing */ 00093 /* the j-th row/column. This property allows CHSEIN to */ 00094 /* perform inverse iteration on just one diagonal block. */ 00095 /* = 'N': no assumptions are made on the correspondence */ 00096 /* between eigenvalues and diagonal blocks. In this */ 00097 /* case, CHSEIN must always perform inverse iteration */ 00098 /* using the whole matrix H. */ 00099 00100 /* INITV (input) CHARACTER*1 */ 00101 /* = 'N': no initial vectors are supplied; */ 00102 /* = 'U': user-supplied initial vectors are stored in the arrays */ 00103 /* VL and/or VR. */ 00104 00105 /* SELECT (input) LOGICAL array, dimension (N) */ 00106 /* Specifies the eigenvectors to be computed. To select the */ 00107 /* eigenvector corresponding to the eigenvalue W(j), */ 00108 /* SELECT(j) must be set to .TRUE.. */ 00109 00110 /* N (input) INTEGER */ 00111 /* The order of the matrix H. N >= 0. */ 00112 00113 /* H (input) COMPLEX array, dimension (LDH,N) */ 00114 /* The upper Hessenberg matrix H. */ 00115 00116 /* LDH (input) INTEGER */ 00117 /* The leading dimension of the array H. LDH >= max(1,N). */ 00118 00119 /* W (input/output) COMPLEX array, dimension (N) */ 00120 /* On entry, the eigenvalues of H. */ 00121 /* On exit, the real parts of W may have been altered since */ 00122 /* close eigenvalues are perturbed slightly in searching for */ 00123 /* independent eigenvectors. */ 00124 00125 /* VL (input/output) COMPLEX array, dimension (LDVL,MM) */ 00126 /* On entry, if INITV = 'U' and SIDE = 'L' or 'B', VL must */ 00127 /* contain starting vectors for the inverse iteration for the */ 00128 /* left eigenvectors; the starting vector for each eigenvector */ 00129 /* must be in the same column in which the eigenvector will be */ 00130 /* stored. */ 00131 /* On exit, if SIDE = 'L' or 'B', the left eigenvectors */ 00132 /* specified by SELECT will be stored consecutively in the */ 00133 /* columns of VL, in the same order as their eigenvalues. */ 00134 /* If SIDE = 'R', VL is not referenced. */ 00135 00136 /* LDVL (input) INTEGER */ 00137 /* The leading dimension of the array VL. */ 00138 /* LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise. */ 00139 00140 /* VR (input/output) COMPLEX array, dimension (LDVR,MM) */ 00141 /* On entry, if INITV = 'U' and SIDE = 'R' or 'B', VR must */ 00142 /* contain starting vectors for the inverse iteration for the */ 00143 /* right eigenvectors; the starting vector for each eigenvector */ 00144 /* must be in the same column in which the eigenvector will be */ 00145 /* stored. */ 00146 /* On exit, if SIDE = 'R' or 'B', the right eigenvectors */ 00147 /* specified by SELECT will be stored consecutively in the */ 00148 /* columns of VR, in the same order as their eigenvalues. */ 00149 /* If SIDE = 'L', VR is not referenced. */ 00150 00151 /* LDVR (input) INTEGER */ 00152 /* The leading dimension of the array VR. */ 00153 /* LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise. */ 00154 00155 /* MM (input) INTEGER */ 00156 /* The number of columns in the arrays VL and/or VR. MM >= M. */ 00157 00158 /* M (output) INTEGER */ 00159 /* The number of columns in the arrays VL and/or VR required to */ 00160 /* store the eigenvectors (= the number of .TRUE. elements in */ 00161 /* SELECT). */ 00162 00163 /* WORK (workspace) COMPLEX array, dimension (N*N) */ 00164 00165 /* RWORK (workspace) REAL array, dimension (N) */ 00166 00167 /* IFAILL (output) INTEGER array, dimension (MM) */ 00168 /* If SIDE = 'L' or 'B', IFAILL(i) = j > 0 if the left */ 00169 /* eigenvector in the i-th column of VL (corresponding to the */ 00170 /* eigenvalue w(j)) failed to converge; IFAILL(i) = 0 if the */ 00171 /* eigenvector converged satisfactorily. */ 00172 /* If SIDE = 'R', IFAILL is not referenced. */ 00173 00174 /* IFAILR (output) INTEGER array, dimension (MM) */ 00175 /* If SIDE = 'R' or 'B', IFAILR(i) = j > 0 if the right */ 00176 /* eigenvector in the i-th column of VR (corresponding to the */ 00177 /* eigenvalue w(j)) failed to converge; IFAILR(i) = 0 if the */ 00178 /* eigenvector converged satisfactorily. */ 00179 /* If SIDE = 'L', IFAILR is not referenced. */ 00180 00181 /* INFO (output) INTEGER */ 00182 /* = 0: successful exit */ 00183 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00184 /* > 0: if INFO = i, i is the number of eigenvectors which */ 00185 /* failed to converge; see IFAILL and IFAILR for further */ 00186 /* details. */ 00187 00188 /* Further Details */ 00189 /* =============== */ 00190 00191 /* Each eigenvector is normalized so that the element of largest */ 00192 /* magnitude has magnitude 1; here the magnitude of a complex number */ 00193 /* (x,y) is taken to be |x|+|y|. */ 00194 00195 /* ===================================================================== */ 00196 00197 /* .. Parameters .. */ 00198 /* .. */ 00199 /* .. Local Scalars .. */ 00200 /* .. */ 00201 /* .. External Functions .. */ 00202 /* .. */ 00203 /* .. External Subroutines .. */ 00204 /* .. */ 00205 /* .. Intrinsic Functions .. */ 00206 /* .. */ 00207 /* .. Statement Functions .. */ 00208 /* .. */ 00209 /* .. Statement Function definitions .. */ 00210 /* .. */ 00211 /* .. Executable Statements .. */ 00212 00213 /* Decode and test the input parameters. */ 00214 00215 /* Parameter adjustments */ 00216 --select; 00217 h_dim1 = *ldh; 00218 h_offset = 1 + h_dim1; 00219 h__ -= h_offset; 00220 --w; 00221 vl_dim1 = *ldvl; 00222 vl_offset = 1 + vl_dim1; 00223 vl -= vl_offset; 00224 vr_dim1 = *ldvr; 00225 vr_offset = 1 + vr_dim1; 00226 vr -= vr_offset; 00227 --work; 00228 --rwork; 00229 --ifaill; 00230 --ifailr; 00231 00232 /* Function Body */ 00233 bothv = lsame_(side, "B"); 00234 rightv = lsame_(side, "R") || bothv; 00235 leftv = lsame_(side, "L") || bothv; 00236 00237 fromqr = lsame_(eigsrc, "Q"); 00238 00239 noinit = lsame_(initv, "N"); 00240 00241 /* Set M to the number of columns required to store the selected */ 00242 /* eigenvectors. */ 00243 00244 *m = 0; 00245 i__1 = *n; 00246 for (k = 1; k <= i__1; ++k) { 00247 if (select[k]) { 00248 ++(*m); 00249 } 00250 /* L10: */ 00251 } 00252 00253 *info = 0; 00254 if (! rightv && ! leftv) { 00255 *info = -1; 00256 } else if (! fromqr && ! lsame_(eigsrc, "N")) { 00257 *info = -2; 00258 } else if (! noinit && ! lsame_(initv, "U")) { 00259 *info = -3; 00260 } else if (*n < 0) { 00261 *info = -5; 00262 } else if (*ldh < max(1,*n)) { 00263 *info = -7; 00264 } else if (*ldvl < 1 || leftv && *ldvl < *n) { 00265 *info = -10; 00266 } else if (*ldvr < 1 || rightv && *ldvr < *n) { 00267 *info = -12; 00268 } else if (*mm < *m) { 00269 *info = -13; 00270 } 00271 if (*info != 0) { 00272 i__1 = -(*info); 00273 xerbla_("CHSEIN", &i__1); 00274 return 0; 00275 } 00276 00277 /* Quick return if possible. */ 00278 00279 if (*n == 0) { 00280 return 0; 00281 } 00282 00283 /* Set machine-dependent constants. */ 00284 00285 unfl = slamch_("Safe minimum"); 00286 ulp = slamch_("Precision"); 00287 smlnum = unfl * (*n / ulp); 00288 00289 ldwork = *n; 00290 00291 kl = 1; 00292 kln = 0; 00293 if (fromqr) { 00294 kr = 0; 00295 } else { 00296 kr = *n; 00297 } 00298 ks = 1; 00299 00300 i__1 = *n; 00301 for (k = 1; k <= i__1; ++k) { 00302 if (select[k]) { 00303 00304 /* Compute eigenvector(s) corresponding to W(K). */ 00305 00306 if (fromqr) { 00307 00308 /* If affiliation of eigenvalues is known, check whether */ 00309 /* the matrix splits. */ 00310 00311 /* Determine KL and KR such that 1 <= KL <= K <= KR <= N */ 00312 /* and H(KL,KL-1) and H(KR+1,KR) are zero (or KL = 1 or */ 00313 /* KR = N). */ 00314 00315 /* Then inverse iteration can be performed with the */ 00316 /* submatrix H(KL:N,KL:N) for a left eigenvector, and with */ 00317 /* the submatrix H(1:KR,1:KR) for a right eigenvector. */ 00318 00319 i__2 = kl + 1; 00320 for (i__ = k; i__ >= i__2; --i__) { 00321 i__3 = i__ + (i__ - 1) * h_dim1; 00322 if (h__[i__3].r == 0.f && h__[i__3].i == 0.f) { 00323 goto L30; 00324 } 00325 /* L20: */ 00326 } 00327 L30: 00328 kl = i__; 00329 if (k > kr) { 00330 i__2 = *n - 1; 00331 for (i__ = k; i__ <= i__2; ++i__) { 00332 i__3 = i__ + 1 + i__ * h_dim1; 00333 if (h__[i__3].r == 0.f && h__[i__3].i == 0.f) { 00334 goto L50; 00335 } 00336 /* L40: */ 00337 } 00338 L50: 00339 kr = i__; 00340 } 00341 } 00342 00343 if (kl != kln) { 00344 kln = kl; 00345 00346 /* Compute infinity-norm of submatrix H(KL:KR,KL:KR) if it */ 00347 /* has not ben computed before. */ 00348 00349 i__2 = kr - kl + 1; 00350 hnorm = clanhs_("I", &i__2, &h__[kl + kl * h_dim1], ldh, & 00351 rwork[1]); 00352 if (hnorm > 0.f) { 00353 eps3 = hnorm * ulp; 00354 } else { 00355 eps3 = smlnum; 00356 } 00357 } 00358 00359 /* Perturb eigenvalue if it is close to any previous */ 00360 /* selected eigenvalues affiliated to the submatrix */ 00361 /* H(KL:KR,KL:KR). Close roots are modified by EPS3. */ 00362 00363 i__2 = k; 00364 wk.r = w[i__2].r, wk.i = w[i__2].i; 00365 L60: 00366 i__2 = kl; 00367 for (i__ = k - 1; i__ >= i__2; --i__) { 00368 i__3 = i__; 00369 q__2.r = w[i__3].r - wk.r, q__2.i = w[i__3].i - wk.i; 00370 q__1.r = q__2.r, q__1.i = q__2.i; 00371 if (select[i__] && (r__1 = q__1.r, dabs(r__1)) + (r__2 = 00372 r_imag(&q__1), dabs(r__2)) < eps3) { 00373 q__1.r = wk.r + eps3, q__1.i = wk.i; 00374 wk.r = q__1.r, wk.i = q__1.i; 00375 goto L60; 00376 } 00377 /* L70: */ 00378 } 00379 i__2 = k; 00380 w[i__2].r = wk.r, w[i__2].i = wk.i; 00381 00382 if (leftv) { 00383 00384 /* Compute left eigenvector. */ 00385 00386 i__2 = *n - kl + 1; 00387 claein_(&c_false, &noinit, &i__2, &h__[kl + kl * h_dim1], ldh, 00388 &wk, &vl[kl + ks * vl_dim1], &work[1], &ldwork, & 00389 rwork[1], &eps3, &smlnum, &iinfo); 00390 if (iinfo > 0) { 00391 ++(*info); 00392 ifaill[ks] = k; 00393 } else { 00394 ifaill[ks] = 0; 00395 } 00396 i__2 = kl - 1; 00397 for (i__ = 1; i__ <= i__2; ++i__) { 00398 i__3 = i__ + ks * vl_dim1; 00399 vl[i__3].r = 0.f, vl[i__3].i = 0.f; 00400 /* L80: */ 00401 } 00402 } 00403 if (rightv) { 00404 00405 /* Compute right eigenvector. */ 00406 00407 claein_(&c_true, &noinit, &kr, &h__[h_offset], ldh, &wk, &vr[ 00408 ks * vr_dim1 + 1], &work[1], &ldwork, &rwork[1], & 00409 eps3, &smlnum, &iinfo); 00410 if (iinfo > 0) { 00411 ++(*info); 00412 ifailr[ks] = k; 00413 } else { 00414 ifailr[ks] = 0; 00415 } 00416 i__2 = *n; 00417 for (i__ = kr + 1; i__ <= i__2; ++i__) { 00418 i__3 = i__ + ks * vr_dim1; 00419 vr[i__3].r = 0.f, vr[i__3].i = 0.f; 00420 /* L90: */ 00421 } 00422 } 00423 ++ks; 00424 } 00425 /* L100: */ 00426 } 00427 00428 return 0; 00429 00430 /* End of CHSEIN */ 00431 00432 } /* chsein_ */