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