00001 /* ctrsen.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 ctrsen_(char *job, char *compq, logical *select, integer 00021 *n, complex *t, integer *ldt, complex *q, integer *ldq, complex *w, 00022 integer *m, real *s, real *sep, complex *work, integer *lwork, 00023 integer *info) 00024 { 00025 /* System generated locals */ 00026 integer q_dim1, q_offset, t_dim1, t_offset, i__1, i__2, i__3; 00027 00028 /* Builtin functions */ 00029 double sqrt(doublereal); 00030 00031 /* Local variables */ 00032 integer k, n1, n2, nn, ks; 00033 real est; 00034 integer kase, ierr; 00035 real scale; 00036 extern logical lsame_(char *, char *); 00037 integer isave[3], lwmin; 00038 logical wantq, wants; 00039 real rnorm; 00040 extern /* Subroutine */ int clacn2_(integer *, complex *, complex *, real 00041 *, integer *, integer *); 00042 real rwork[1]; 00043 extern doublereal clange_(char *, integer *, integer *, complex *, 00044 integer *, real *); 00045 extern /* Subroutine */ int clacpy_(char *, integer *, integer *, complex 00046 *, integer *, complex *, integer *), xerbla_(char *, 00047 integer *); 00048 logical wantbh; 00049 extern /* Subroutine */ int ctrexc_(char *, integer *, complex *, integer 00050 *, complex *, integer *, integer *, integer *, integer *); 00051 logical wantsp; 00052 extern /* Subroutine */ int ctrsyl_(char *, char *, integer *, integer *, 00053 integer *, complex *, integer *, complex *, integer *, complex *, 00054 integer *, real *, integer *); 00055 logical lquery; 00056 00057 00058 /* -- LAPACK routine (version 3.2) -- */ 00059 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00060 /* November 2006 */ 00061 00062 /* Modified to call CLACN2 in place of CLACON, 10 Feb 03, SJH. */ 00063 00064 /* .. Scalar Arguments .. */ 00065 /* .. */ 00066 /* .. Array Arguments .. */ 00067 /* .. */ 00068 00069 /* Purpose */ 00070 /* ======= */ 00071 00072 /* CTRSEN reorders the Schur factorization of a complex matrix */ 00073 /* A = Q*T*Q**H, so that a selected cluster of eigenvalues appears in */ 00074 /* the leading positions on the diagonal of the upper triangular matrix */ 00075 /* T, and the leading columns of Q form an orthonormal basis of the */ 00076 /* corresponding right invariant subspace. */ 00077 00078 /* Optionally the routine computes the reciprocal condition numbers of */ 00079 /* the cluster of eigenvalues and/or the invariant subspace. */ 00080 00081 /* Arguments */ 00082 /* ========= */ 00083 00084 /* JOB (input) CHARACTER*1 */ 00085 /* Specifies whether condition numbers are required for the */ 00086 /* cluster of eigenvalues (S) or the invariant subspace (SEP): */ 00087 /* = 'N': none; */ 00088 /* = 'E': for eigenvalues only (S); */ 00089 /* = 'V': for invariant subspace only (SEP); */ 00090 /* = 'B': for both eigenvalues and invariant subspace (S and */ 00091 /* SEP). */ 00092 00093 /* COMPQ (input) CHARACTER*1 */ 00094 /* = 'V': update the matrix Q of Schur vectors; */ 00095 /* = 'N': do not update Q. */ 00096 00097 /* SELECT (input) LOGICAL array, dimension (N) */ 00098 /* SELECT specifies the eigenvalues in the selected cluster. To */ 00099 /* select the j-th eigenvalue, SELECT(j) must be set to .TRUE.. */ 00100 00101 /* N (input) INTEGER */ 00102 /* The order of the matrix T. N >= 0. */ 00103 00104 /* T (input/output) COMPLEX array, dimension (LDT,N) */ 00105 /* On entry, the upper triangular matrix T. */ 00106 /* On exit, T is overwritten by the reordered matrix T, with the */ 00107 /* selected eigenvalues as the leading diagonal elements. */ 00108 00109 /* LDT (input) INTEGER */ 00110 /* The leading dimension of the array T. LDT >= max(1,N). */ 00111 00112 /* Q (input/output) COMPLEX array, dimension (LDQ,N) */ 00113 /* On entry, if COMPQ = 'V', the matrix Q of Schur vectors. */ 00114 /* On exit, if COMPQ = 'V', Q has been postmultiplied by the */ 00115 /* unitary transformation matrix which reorders T; the leading M */ 00116 /* columns of Q form an orthonormal basis for the specified */ 00117 /* invariant subspace. */ 00118 /* If COMPQ = 'N', Q is not referenced. */ 00119 00120 /* LDQ (input) INTEGER */ 00121 /* The leading dimension of the array Q. */ 00122 /* LDQ >= 1; and if COMPQ = 'V', LDQ >= N. */ 00123 00124 /* W (output) COMPLEX array, dimension (N) */ 00125 /* The reordered eigenvalues of T, in the same order as they */ 00126 /* appear on the diagonal of T. */ 00127 00128 /* M (output) INTEGER */ 00129 /* The dimension of the specified invariant subspace. */ 00130 /* 0 <= M <= N. */ 00131 00132 /* S (output) REAL */ 00133 /* If JOB = 'E' or 'B', S is a lower bound on the reciprocal */ 00134 /* condition number for the selected cluster of eigenvalues. */ 00135 /* S cannot underestimate the true reciprocal condition number */ 00136 /* by more than a factor of sqrt(N). If M = 0 or N, S = 1. */ 00137 /* If JOB = 'N' or 'V', S is not referenced. */ 00138 00139 /* SEP (output) REAL */ 00140 /* If JOB = 'V' or 'B', SEP is the estimated reciprocal */ 00141 /* condition number of the specified invariant subspace. If */ 00142 /* M = 0 or N, SEP = norm(T). */ 00143 /* If JOB = 'N' or 'E', SEP is not referenced. */ 00144 00145 /* WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK)) */ 00146 /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */ 00147 00148 /* LWORK (input) INTEGER */ 00149 /* The dimension of the array WORK. */ 00150 /* If JOB = 'N', LWORK >= 1; */ 00151 /* if JOB = 'E', LWORK = max(1,M*(N-M)); */ 00152 /* if JOB = 'V' or 'B', LWORK >= max(1,2*M*(N-M)). */ 00153 00154 /* If LWORK = -1, then a workspace query is assumed; the routine */ 00155 /* only calculates the optimal size of the WORK array, returns */ 00156 /* this value as the first entry of the WORK array, and no error */ 00157 /* message related to LWORK is issued by XERBLA. */ 00158 00159 /* INFO (output) INTEGER */ 00160 /* = 0: successful exit */ 00161 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00162 00163 /* Further Details */ 00164 /* =============== */ 00165 00166 /* CTRSEN first collects the selected eigenvalues by computing a unitary */ 00167 /* transformation Z to move them to the top left corner of T. In other */ 00168 /* words, the selected eigenvalues are the eigenvalues of T11 in: */ 00169 00170 /* Z'*T*Z = ( T11 T12 ) n1 */ 00171 /* ( 0 T22 ) n2 */ 00172 /* n1 n2 */ 00173 00174 /* where N = n1+n2 and Z' means the conjugate transpose of Z. The first */ 00175 /* n1 columns of Z span the specified invariant subspace of T. */ 00176 00177 /* If T has been obtained from the Schur factorization of a matrix */ 00178 /* A = Q*T*Q', then the reordered Schur factorization of A is given by */ 00179 /* A = (Q*Z)*(Z'*T*Z)*(Q*Z)', and the first n1 columns of Q*Z span the */ 00180 /* corresponding invariant subspace of A. */ 00181 00182 /* The reciprocal condition number of the average of the eigenvalues of */ 00183 /* T11 may be returned in S. S lies between 0 (very badly conditioned) */ 00184 /* and 1 (very well conditioned). It is computed as follows. First we */ 00185 /* compute R so that */ 00186 00187 /* P = ( I R ) n1 */ 00188 /* ( 0 0 ) n2 */ 00189 /* n1 n2 */ 00190 00191 /* is the projector on the invariant subspace associated with T11. */ 00192 /* R is the solution of the Sylvester equation: */ 00193 00194 /* T11*R - R*T22 = T12. */ 00195 00196 /* Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote */ 00197 /* the two-norm of M. Then S is computed as the lower bound */ 00198 00199 /* (1 + F-norm(R)**2)**(-1/2) */ 00200 00201 /* on the reciprocal of 2-norm(P), the true reciprocal condition number. */ 00202 /* S cannot underestimate 1 / 2-norm(P) by more than a factor of */ 00203 /* sqrt(N). */ 00204 00205 /* An approximate error bound for the computed average of the */ 00206 /* eigenvalues of T11 is */ 00207 00208 /* EPS * norm(T) / S */ 00209 00210 /* where EPS is the machine precision. */ 00211 00212 /* The reciprocal condition number of the right invariant subspace */ 00213 /* spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP. */ 00214 /* SEP is defined as the separation of T11 and T22: */ 00215 00216 /* sep( T11, T22 ) = sigma-min( C ) */ 00217 00218 /* where sigma-min(C) is the smallest singular value of the */ 00219 /* n1*n2-by-n1*n2 matrix */ 00220 00221 /* C = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) ) */ 00222 00223 /* I(m) is an m by m identity matrix, and kprod denotes the Kronecker */ 00224 /* product. We estimate sigma-min(C) by the reciprocal of an estimate of */ 00225 /* the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C) */ 00226 /* cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2). */ 00227 00228 /* When SEP is small, small changes in T can cause large changes in */ 00229 /* the invariant subspace. An approximate bound on the maximum angular */ 00230 /* error in the computed right invariant subspace is */ 00231 00232 /* EPS * norm(T) / SEP */ 00233 00234 /* ===================================================================== */ 00235 00236 /* .. Parameters .. */ 00237 /* .. */ 00238 /* .. Local Scalars .. */ 00239 /* .. */ 00240 /* .. Local Arrays .. */ 00241 /* .. */ 00242 /* .. External Functions .. */ 00243 /* .. */ 00244 /* .. External Subroutines .. */ 00245 /* .. */ 00246 /* .. Intrinsic Functions .. */ 00247 /* .. */ 00248 /* .. Executable Statements .. */ 00249 00250 /* Decode and test the input parameters. */ 00251 00252 /* Parameter adjustments */ 00253 --select; 00254 t_dim1 = *ldt; 00255 t_offset = 1 + t_dim1; 00256 t -= t_offset; 00257 q_dim1 = *ldq; 00258 q_offset = 1 + q_dim1; 00259 q -= q_offset; 00260 --w; 00261 --work; 00262 00263 /* Function Body */ 00264 wantbh = lsame_(job, "B"); 00265 wants = lsame_(job, "E") || wantbh; 00266 wantsp = lsame_(job, "V") || wantbh; 00267 wantq = lsame_(compq, "V"); 00268 00269 /* Set M to the number of selected eigenvalues. */ 00270 00271 *m = 0; 00272 i__1 = *n; 00273 for (k = 1; k <= i__1; ++k) { 00274 if (select[k]) { 00275 ++(*m); 00276 } 00277 /* L10: */ 00278 } 00279 00280 n1 = *m; 00281 n2 = *n - *m; 00282 nn = n1 * n2; 00283 00284 *info = 0; 00285 lquery = *lwork == -1; 00286 00287 if (wantsp) { 00288 /* Computing MAX */ 00289 i__1 = 1, i__2 = nn << 1; 00290 lwmin = max(i__1,i__2); 00291 } else if (lsame_(job, "N")) { 00292 lwmin = 1; 00293 } else if (lsame_(job, "E")) { 00294 lwmin = max(1,nn); 00295 } 00296 00297 if (! lsame_(job, "N") && ! wants && ! wantsp) { 00298 *info = -1; 00299 } else if (! lsame_(compq, "N") && ! wantq) { 00300 *info = -2; 00301 } else if (*n < 0) { 00302 *info = -4; 00303 } else if (*ldt < max(1,*n)) { 00304 *info = -6; 00305 } else if (*ldq < 1 || wantq && *ldq < *n) { 00306 *info = -8; 00307 } else if (*lwork < lwmin && ! lquery) { 00308 *info = -14; 00309 } 00310 00311 if (*info == 0) { 00312 work[1].r = (real) lwmin, work[1].i = 0.f; 00313 } 00314 00315 if (*info != 0) { 00316 i__1 = -(*info); 00317 xerbla_("CTRSEN", &i__1); 00318 return 0; 00319 } else if (lquery) { 00320 return 0; 00321 } 00322 00323 /* Quick return if possible */ 00324 00325 if (*m == *n || *m == 0) { 00326 if (wants) { 00327 *s = 1.f; 00328 } 00329 if (wantsp) { 00330 *sep = clange_("1", n, n, &t[t_offset], ldt, rwork); 00331 } 00332 goto L40; 00333 } 00334 00335 /* Collect the selected eigenvalues at the top left corner of T. */ 00336 00337 ks = 0; 00338 i__1 = *n; 00339 for (k = 1; k <= i__1; ++k) { 00340 if (select[k]) { 00341 ++ks; 00342 00343 /* Swap the K-th eigenvalue to position KS. */ 00344 00345 if (k != ks) { 00346 ctrexc_(compq, n, &t[t_offset], ldt, &q[q_offset], ldq, &k, & 00347 ks, &ierr); 00348 } 00349 } 00350 /* L20: */ 00351 } 00352 00353 if (wants) { 00354 00355 /* Solve the Sylvester equation for R: */ 00356 00357 /* T11*R - R*T22 = scale*T12 */ 00358 00359 clacpy_("F", &n1, &n2, &t[(n1 + 1) * t_dim1 + 1], ldt, &work[1], &n1); 00360 ctrsyl_("N", "N", &c_n1, &n1, &n2, &t[t_offset], ldt, &t[n1 + 1 + (n1 00361 + 1) * t_dim1], ldt, &work[1], &n1, &scale, &ierr); 00362 00363 /* Estimate the reciprocal of the condition number of the cluster */ 00364 /* of eigenvalues. */ 00365 00366 rnorm = clange_("F", &n1, &n2, &work[1], &n1, rwork); 00367 if (rnorm == 0.f) { 00368 *s = 1.f; 00369 } else { 00370 *s = scale / (sqrt(scale * scale / rnorm + rnorm) * sqrt(rnorm)); 00371 } 00372 } 00373 00374 if (wantsp) { 00375 00376 /* Estimate sep(T11,T22). */ 00377 00378 est = 0.f; 00379 kase = 0; 00380 L30: 00381 clacn2_(&nn, &work[nn + 1], &work[1], &est, &kase, isave); 00382 if (kase != 0) { 00383 if (kase == 1) { 00384 00385 /* Solve T11*R - R*T22 = scale*X. */ 00386 00387 ctrsyl_("N", "N", &c_n1, &n1, &n2, &t[t_offset], ldt, &t[n1 + 00388 1 + (n1 + 1) * t_dim1], ldt, &work[1], &n1, &scale, & 00389 ierr); 00390 } else { 00391 00392 /* Solve T11'*R - R*T22' = scale*X. */ 00393 00394 ctrsyl_("C", "C", &c_n1, &n1, &n2, &t[t_offset], ldt, &t[n1 + 00395 1 + (n1 + 1) * t_dim1], ldt, &work[1], &n1, &scale, & 00396 ierr); 00397 } 00398 goto L30; 00399 } 00400 00401 *sep = scale / est; 00402 } 00403 00404 L40: 00405 00406 /* Copy reordered eigenvalues to W. */ 00407 00408 i__1 = *n; 00409 for (k = 1; k <= i__1; ++k) { 00410 i__2 = k; 00411 i__3 = k + k * t_dim1; 00412 w[i__2].r = t[i__3].r, w[i__2].i = t[i__3].i; 00413 /* L50: */ 00414 } 00415 00416 work[1].r = (real) lwmin, work[1].i = 0.f; 00417 00418 return 0; 00419 00420 /* End of CTRSEN */ 00421 00422 } /* ctrsen_ */