00001 /* strsen.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_n1 = -1; 00019 00020 /* Subroutine */ int strsen_(char *job, char *compq, logical *select, integer 00021 *n, real *t, integer *ldt, real *q, integer *ldq, real *wr, real *wi, 00022 integer *m, real *s, real *sep, real *work, integer *lwork, integer * 00023 iwork, integer *liwork, integer *info) 00024 { 00025 /* System generated locals */ 00026 integer q_dim1, q_offset, t_dim1, t_offset, i__1, i__2; 00027 real r__1, r__2; 00028 00029 /* Builtin functions */ 00030 double sqrt(doublereal); 00031 00032 /* Local variables */ 00033 integer k, n1, n2, kk, nn, ks; 00034 real est; 00035 integer kase; 00036 logical pair; 00037 integer ierr; 00038 logical swap; 00039 real scale; 00040 extern logical lsame_(char *, char *); 00041 integer isave[3], lwmin; 00042 logical wantq, wants; 00043 real rnorm; 00044 extern /* Subroutine */ int slacn2_(integer *, real *, real *, integer *, 00045 real *, integer *, integer *); 00046 extern doublereal slange_(char *, integer *, integer *, real *, integer *, 00047 real *); 00048 extern /* Subroutine */ int xerbla_(char *, integer *); 00049 logical wantbh; 00050 extern /* Subroutine */ int slacpy_(char *, integer *, integer *, real *, 00051 integer *, real *, integer *); 00052 integer liwmin; 00053 extern /* Subroutine */ int strexc_(char *, integer *, real *, integer *, 00054 real *, integer *, integer *, integer *, real *, integer *); 00055 logical wantsp, lquery; 00056 extern /* Subroutine */ int strsyl_(char *, char *, integer *, integer *, 00057 integer *, real *, integer *, real *, integer *, real *, integer * 00058 , real *, integer *); 00059 00060 00061 /* -- LAPACK routine (version 3.2) -- */ 00062 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00063 /* November 2006 */ 00064 00065 /* Modified to call SLACN2 in place of SLACON, 7 Feb 03, SJH. */ 00066 00067 /* .. Scalar Arguments .. */ 00068 /* .. */ 00069 /* .. Array Arguments .. */ 00070 /* .. */ 00071 00072 /* Purpose */ 00073 /* ======= */ 00074 00075 /* STRSEN reorders the real Schur factorization of a real matrix */ 00076 /* A = Q*T*Q**T, so that a selected cluster of eigenvalues appears in */ 00077 /* the leading diagonal blocks of the upper quasi-triangular matrix T, */ 00078 /* and the leading columns of Q form an orthonormal basis of the */ 00079 /* corresponding right invariant subspace. */ 00080 00081 /* Optionally the routine computes the reciprocal condition numbers of */ 00082 /* the cluster of eigenvalues and/or the invariant subspace. */ 00083 00084 /* T must be in Schur canonical form (as returned by SHSEQR), that is, */ 00085 /* block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each */ 00086 /* 2-by-2 diagonal block has its diagonal elemnts equal and its */ 00087 /* off-diagonal elements of opposite sign. */ 00088 00089 /* Arguments */ 00090 /* ========= */ 00091 00092 /* JOB (input) CHARACTER*1 */ 00093 /* Specifies whether condition numbers are required for the */ 00094 /* cluster of eigenvalues (S) or the invariant subspace (SEP): */ 00095 /* = 'N': none; */ 00096 /* = 'E': for eigenvalues only (S); */ 00097 /* = 'V': for invariant subspace only (SEP); */ 00098 /* = 'B': for both eigenvalues and invariant subspace (S and */ 00099 /* SEP). */ 00100 00101 /* COMPQ (input) CHARACTER*1 */ 00102 /* = 'V': update the matrix Q of Schur vectors; */ 00103 /* = 'N': do not update Q. */ 00104 00105 /* SELECT (input) LOGICAL array, dimension (N) */ 00106 /* SELECT specifies the eigenvalues in the selected cluster. To */ 00107 /* select a real eigenvalue w(j), SELECT(j) must be set to */ 00108 /* .TRUE.. To select a complex conjugate pair of eigenvalues */ 00109 /* w(j) and w(j+1), corresponding to a 2-by-2 diagonal block, */ 00110 /* either SELECT(j) or SELECT(j+1) or both must be set to */ 00111 /* .TRUE.; a complex conjugate pair of eigenvalues must be */ 00112 /* either both included in the cluster or both excluded. */ 00113 00114 /* N (input) INTEGER */ 00115 /* The order of the matrix T. N >= 0. */ 00116 00117 /* T (input/output) REAL array, dimension (LDT,N) */ 00118 /* On entry, the upper quasi-triangular matrix T, in Schur */ 00119 /* canonical form. */ 00120 /* On exit, T is overwritten by the reordered matrix T, again in */ 00121 /* Schur canonical form, with the selected eigenvalues in the */ 00122 /* leading diagonal blocks. */ 00123 00124 /* LDT (input) INTEGER */ 00125 /* The leading dimension of the array T. LDT >= max(1,N). */ 00126 00127 /* Q (input/output) REAL array, dimension (LDQ,N) */ 00128 /* On entry, if COMPQ = 'V', the matrix Q of Schur vectors. */ 00129 /* On exit, if COMPQ = 'V', Q has been postmultiplied by the */ 00130 /* orthogonal transformation matrix which reorders T; the */ 00131 /* leading M columns of Q form an orthonormal basis for the */ 00132 /* specified invariant subspace. */ 00133 /* If COMPQ = 'N', Q is not referenced. */ 00134 00135 /* LDQ (input) INTEGER */ 00136 /* The leading dimension of the array Q. */ 00137 /* LDQ >= 1; and if COMPQ = 'V', LDQ >= N. */ 00138 00139 /* WR (output) REAL array, dimension (N) */ 00140 /* WI (output) REAL array, dimension (N) */ 00141 /* The real and imaginary parts, respectively, of the reordered */ 00142 /* eigenvalues of T. The eigenvalues are stored in the same */ 00143 /* order as on the diagonal of T, with WR(i) = T(i,i) and, if */ 00144 /* T(i:i+1,i:i+1) is a 2-by-2 diagonal block, WI(i) > 0 and */ 00145 /* WI(i+1) = -WI(i). Note that if a complex eigenvalue is */ 00146 /* sufficiently ill-conditioned, then its value may differ */ 00147 /* significantly from its value before reordering. */ 00148 00149 /* M (output) INTEGER */ 00150 /* The dimension of the specified invariant subspace. */ 00151 /* 0 < = M <= N. */ 00152 00153 /* S (output) REAL */ 00154 /* If JOB = 'E' or 'B', S is a lower bound on the reciprocal */ 00155 /* condition number for the selected cluster of eigenvalues. */ 00156 /* S cannot underestimate the true reciprocal condition number */ 00157 /* by more than a factor of sqrt(N). If M = 0 or N, S = 1. */ 00158 /* If JOB = 'N' or 'V', S is not referenced. */ 00159 00160 /* SEP (output) REAL */ 00161 /* If JOB = 'V' or 'B', SEP is the estimated reciprocal */ 00162 /* condition number of the specified invariant subspace. If */ 00163 /* M = 0 or N, SEP = norm(T). */ 00164 /* If JOB = 'N' or 'E', SEP is not referenced. */ 00165 00166 /* WORK (workspace/output) REAL array, dimension (MAX(1,LWORK)) */ 00167 /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */ 00168 00169 /* LWORK (input) INTEGER */ 00170 /* The dimension of the array WORK. */ 00171 /* If JOB = 'N', LWORK >= max(1,N); */ 00172 /* if JOB = 'E', LWORK >= max(1,M*(N-M)); */ 00173 /* if JOB = 'V' or 'B', LWORK >= max(1,2*M*(N-M)). */ 00174 00175 /* If LWORK = -1, then a workspace query is assumed; the routine */ 00176 /* only calculates the optimal size of the WORK array, returns */ 00177 /* this value as the first entry of the WORK array, and no error */ 00178 /* message related to LWORK is issued by XERBLA. */ 00179 00180 /* IWORK (workspace) INTEGER array, dimension (MAX(1,LIWORK)) */ 00181 /* On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. */ 00182 00183 /* LIWORK (input) INTEGER */ 00184 /* The dimension of the array IWORK. */ 00185 /* If JOB = 'N' or 'E', LIWORK >= 1; */ 00186 /* if JOB = 'V' or 'B', LIWORK >= max(1,M*(N-M)). */ 00187 00188 /* If LIWORK = -1, then a workspace query is assumed; the */ 00189 /* routine only calculates the optimal size of the IWORK array, */ 00190 /* returns this value as the first entry of the IWORK array, and */ 00191 /* no error message related to LIWORK is issued by XERBLA. */ 00192 00193 /* INFO (output) INTEGER */ 00194 /* = 0: successful exit */ 00195 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00196 /* = 1: reordering of T failed because some eigenvalues are too */ 00197 /* close to separate (the problem is very ill-conditioned); */ 00198 /* T may have been partially reordered, and WR and WI */ 00199 /* contain the eigenvalues in the same order as in T; S and */ 00200 /* SEP (if requested) are set to zero. */ 00201 00202 /* Further Details */ 00203 /* =============== */ 00204 00205 /* STRSEN first collects the selected eigenvalues by computing an */ 00206 /* orthogonal transformation Z to move them to the top left corner of T. */ 00207 /* In other words, the selected eigenvalues are the eigenvalues of T11 */ 00208 /* in: */ 00209 00210 /* Z'*T*Z = ( T11 T12 ) n1 */ 00211 /* ( 0 T22 ) n2 */ 00212 /* n1 n2 */ 00213 00214 /* where N = n1+n2 and Z' means the transpose of Z. The first n1 columns */ 00215 /* of Z span the specified invariant subspace of T. */ 00216 00217 /* If T has been obtained from the real Schur factorization of a matrix */ 00218 /* A = Q*T*Q', then the reordered real Schur factorization of A is given */ 00219 /* by A = (Q*Z)*(Z'*T*Z)*(Q*Z)', and the first n1 columns of Q*Z span */ 00220 /* the corresponding invariant subspace of A. */ 00221 00222 /* The reciprocal condition number of the average of the eigenvalues of */ 00223 /* T11 may be returned in S. S lies between 0 (very badly conditioned) */ 00224 /* and 1 (very well conditioned). It is computed as follows. First we */ 00225 /* compute R so that */ 00226 00227 /* P = ( I R ) n1 */ 00228 /* ( 0 0 ) n2 */ 00229 /* n1 n2 */ 00230 00231 /* is the projector on the invariant subspace associated with T11. */ 00232 /* R is the solution of the Sylvester equation: */ 00233 00234 /* T11*R - R*T22 = T12. */ 00235 00236 /* Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote */ 00237 /* the two-norm of M. Then S is computed as the lower bound */ 00238 00239 /* (1 + F-norm(R)**2)**(-1/2) */ 00240 00241 /* on the reciprocal of 2-norm(P), the true reciprocal condition number. */ 00242 /* S cannot underestimate 1 / 2-norm(P) by more than a factor of */ 00243 /* sqrt(N). */ 00244 00245 /* An approximate error bound for the computed average of the */ 00246 /* eigenvalues of T11 is */ 00247 00248 /* EPS * norm(T) / S */ 00249 00250 /* where EPS is the machine precision. */ 00251 00252 /* The reciprocal condition number of the right invariant subspace */ 00253 /* spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP. */ 00254 /* SEP is defined as the separation of T11 and T22: */ 00255 00256 /* sep( T11, T22 ) = sigma-min( C ) */ 00257 00258 /* where sigma-min(C) is the smallest singular value of the */ 00259 /* n1*n2-by-n1*n2 matrix */ 00260 00261 /* C = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) ) */ 00262 00263 /* I(m) is an m by m identity matrix, and kprod denotes the Kronecker */ 00264 /* product. We estimate sigma-min(C) by the reciprocal of an estimate of */ 00265 /* the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C) */ 00266 /* cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2). */ 00267 00268 /* When SEP is small, small changes in T can cause large changes in */ 00269 /* the invariant subspace. An approximate bound on the maximum angular */ 00270 /* error in the computed right invariant subspace is */ 00271 00272 /* EPS * norm(T) / SEP */ 00273 00274 /* ===================================================================== */ 00275 00276 /* .. Parameters .. */ 00277 /* .. */ 00278 /* .. Local Scalars .. */ 00279 /* .. */ 00280 /* .. Local Arrays .. */ 00281 /* .. */ 00282 /* .. External Functions .. */ 00283 /* .. */ 00284 /* .. External Subroutines .. */ 00285 /* .. */ 00286 /* .. Intrinsic Functions .. */ 00287 /* .. */ 00288 /* .. Executable Statements .. */ 00289 00290 /* Decode and test the input parameters */ 00291 00292 /* Parameter adjustments */ 00293 --select; 00294 t_dim1 = *ldt; 00295 t_offset = 1 + t_dim1; 00296 t -= t_offset; 00297 q_dim1 = *ldq; 00298 q_offset = 1 + q_dim1; 00299 q -= q_offset; 00300 --wr; 00301 --wi; 00302 --work; 00303 --iwork; 00304 00305 /* Function Body */ 00306 wantbh = lsame_(job, "B"); 00307 wants = lsame_(job, "E") || wantbh; 00308 wantsp = lsame_(job, "V") || wantbh; 00309 wantq = lsame_(compq, "V"); 00310 00311 *info = 0; 00312 lquery = *lwork == -1; 00313 if (! lsame_(job, "N") && ! wants && ! wantsp) { 00314 *info = -1; 00315 } else if (! lsame_(compq, "N") && ! wantq) { 00316 *info = -2; 00317 } else if (*n < 0) { 00318 *info = -4; 00319 } else if (*ldt < max(1,*n)) { 00320 *info = -6; 00321 } else if (*ldq < 1 || wantq && *ldq < *n) { 00322 *info = -8; 00323 } else { 00324 00325 /* Set M to the dimension of the specified invariant subspace, */ 00326 /* and test LWORK and LIWORK. */ 00327 00328 *m = 0; 00329 pair = FALSE_; 00330 i__1 = *n; 00331 for (k = 1; k <= i__1; ++k) { 00332 if (pair) { 00333 pair = FALSE_; 00334 } else { 00335 if (k < *n) { 00336 if (t[k + 1 + k * t_dim1] == 0.f) { 00337 if (select[k]) { 00338 ++(*m); 00339 } 00340 } else { 00341 pair = TRUE_; 00342 if (select[k] || select[k + 1]) { 00343 *m += 2; 00344 } 00345 } 00346 } else { 00347 if (select[*n]) { 00348 ++(*m); 00349 } 00350 } 00351 } 00352 /* L10: */ 00353 } 00354 00355 n1 = *m; 00356 n2 = *n - *m; 00357 nn = n1 * n2; 00358 00359 if (wantsp) { 00360 /* Computing MAX */ 00361 i__1 = 1, i__2 = nn << 1; 00362 lwmin = max(i__1,i__2); 00363 liwmin = max(1,nn); 00364 } else if (lsame_(job, "N")) { 00365 lwmin = max(1,*n); 00366 liwmin = 1; 00367 } else if (lsame_(job, "E")) { 00368 lwmin = max(1,nn); 00369 liwmin = 1; 00370 } 00371 00372 if (*lwork < lwmin && ! lquery) { 00373 *info = -15; 00374 } else if (*liwork < liwmin && ! lquery) { 00375 *info = -17; 00376 } 00377 } 00378 00379 if (*info == 0) { 00380 work[1] = (real) lwmin; 00381 iwork[1] = liwmin; 00382 } 00383 00384 if (*info != 0) { 00385 i__1 = -(*info); 00386 xerbla_("STRSEN", &i__1); 00387 return 0; 00388 } else if (lquery) { 00389 return 0; 00390 } 00391 00392 /* Quick return if possible. */ 00393 00394 if (*m == *n || *m == 0) { 00395 if (wants) { 00396 *s = 1.f; 00397 } 00398 if (wantsp) { 00399 *sep = slange_("1", n, n, &t[t_offset], ldt, &work[1]); 00400 } 00401 goto L40; 00402 } 00403 00404 /* Collect the selected blocks at the top-left corner of T. */ 00405 00406 ks = 0; 00407 pair = FALSE_; 00408 i__1 = *n; 00409 for (k = 1; k <= i__1; ++k) { 00410 if (pair) { 00411 pair = FALSE_; 00412 } else { 00413 swap = select[k]; 00414 if (k < *n) { 00415 if (t[k + 1 + k * t_dim1] != 0.f) { 00416 pair = TRUE_; 00417 swap = swap || select[k + 1]; 00418 } 00419 } 00420 if (swap) { 00421 ++ks; 00422 00423 /* Swap the K-th block to position KS. */ 00424 00425 ierr = 0; 00426 kk = k; 00427 if (k != ks) { 00428 strexc_(compq, n, &t[t_offset], ldt, &q[q_offset], ldq, & 00429 kk, &ks, &work[1], &ierr); 00430 } 00431 if (ierr == 1 || ierr == 2) { 00432 00433 /* Blocks too close to swap: exit. */ 00434 00435 *info = 1; 00436 if (wants) { 00437 *s = 0.f; 00438 } 00439 if (wantsp) { 00440 *sep = 0.f; 00441 } 00442 goto L40; 00443 } 00444 if (pair) { 00445 ++ks; 00446 } 00447 } 00448 } 00449 /* L20: */ 00450 } 00451 00452 if (wants) { 00453 00454 /* Solve Sylvester equation for R: */ 00455 00456 /* T11*R - R*T22 = scale*T12 */ 00457 00458 slacpy_("F", &n1, &n2, &t[(n1 + 1) * t_dim1 + 1], ldt, &work[1], &n1); 00459 strsyl_("N", "N", &c_n1, &n1, &n2, &t[t_offset], ldt, &t[n1 + 1 + (n1 00460 + 1) * t_dim1], ldt, &work[1], &n1, &scale, &ierr); 00461 00462 /* Estimate the reciprocal of the condition number of the cluster */ 00463 /* of eigenvalues. */ 00464 00465 rnorm = slange_("F", &n1, &n2, &work[1], &n1, &work[1]); 00466 if (rnorm == 0.f) { 00467 *s = 1.f; 00468 } else { 00469 *s = scale / (sqrt(scale * scale / rnorm + rnorm) * sqrt(rnorm)); 00470 } 00471 } 00472 00473 if (wantsp) { 00474 00475 /* Estimate sep(T11,T22). */ 00476 00477 est = 0.f; 00478 kase = 0; 00479 L30: 00480 slacn2_(&nn, &work[nn + 1], &work[1], &iwork[1], &est, &kase, isave); 00481 if (kase != 0) { 00482 if (kase == 1) { 00483 00484 /* Solve T11*R - R*T22 = scale*X. */ 00485 00486 strsyl_("N", "N", &c_n1, &n1, &n2, &t[t_offset], ldt, &t[n1 + 00487 1 + (n1 + 1) * t_dim1], ldt, &work[1], &n1, &scale, & 00488 ierr); 00489 } else { 00490 00491 /* Solve T11'*R - R*T22' = scale*X. */ 00492 00493 strsyl_("T", "T", &c_n1, &n1, &n2, &t[t_offset], ldt, &t[n1 + 00494 1 + (n1 + 1) * t_dim1], ldt, &work[1], &n1, &scale, & 00495 ierr); 00496 } 00497 goto L30; 00498 } 00499 00500 *sep = scale / est; 00501 } 00502 00503 L40: 00504 00505 /* Store the output eigenvalues in WR and WI. */ 00506 00507 i__1 = *n; 00508 for (k = 1; k <= i__1; ++k) { 00509 wr[k] = t[k + k * t_dim1]; 00510 wi[k] = 0.f; 00511 /* L50: */ 00512 } 00513 i__1 = *n - 1; 00514 for (k = 1; k <= i__1; ++k) { 00515 if (t[k + 1 + k * t_dim1] != 0.f) { 00516 wi[k] = sqrt((r__1 = t[k + (k + 1) * t_dim1], dabs(r__1))) * sqrt( 00517 (r__2 = t[k + 1 + k * t_dim1], dabs(r__2))); 00518 wi[k + 1] = -wi[k]; 00519 } 00520 /* L60: */ 00521 } 00522 00523 work[1] = (real) lwmin; 00524 iwork[1] = liwmin; 00525 00526 return 0; 00527 00528 /* End of STRSEN */ 00529 00530 } /* strsen_ */