1 /* dtrsen.f -- translated by f2c (version 20061008).
2 You must link the resulting object file with libf2c:
3 on Microsoft Windows system, link with libf2c.lib;
4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5 or, if you install libf2c.a in a standard place, with -lf2c -lm
6 -- in that order, at the end of the command line, as in
7 cc *.o -lf2c -lm
8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
9
10 http://www.netlib.org/f2c/libf2c.zip
11 */
12
13 #include "f2c.h"
14 #include "blaswrap.h"
15
16 /* Table of constant values */
17
18 static integer c_n1 = -1;
19
dtrsen_(char * job,char * compq,logical * select,integer * n,doublereal * t,integer * ldt,doublereal * q,integer * ldq,doublereal * wr,doublereal * wi,integer * m,doublereal * s,doublereal * sep,doublereal * work,integer * lwork,integer * iwork,integer * liwork,integer * info)20 /* Subroutine */ int dtrsen_(char *job, char *compq, logical *select, integer
21 *n, doublereal *t, integer *ldt, doublereal *q, integer *ldq,
22 doublereal *wr, doublereal *wi, integer *m, doublereal *s, doublereal
23 *sep, doublereal *work, integer *lwork, integer *iwork, integer *
24 liwork, integer *info)
25 {
26 /* System generated locals */
27 integer q_dim1, q_offset, t_dim1, t_offset, i__1, i__2;
28 doublereal d__1, d__2;
29
30 /* Builtin functions */
31 double sqrt(doublereal);
32
33 /* Local variables */
34 integer k, n1, n2, kk, nn, ks;
35 doublereal est;
36 integer kase;
37 logical pair;
38 integer ierr;
39 logical swap;
40 doublereal scale;
41 extern logical lsame_(char *, char *);
42 integer isave[3], lwmin;
43 logical wantq, wants;
44 doublereal rnorm;
45 extern /* Subroutine */ int dlacn2_(integer *, doublereal *, doublereal *,
46 integer *, doublereal *, integer *, integer *);
47 extern doublereal dlange_(char *, integer *, integer *, doublereal *,
48 integer *, doublereal *);
49 extern /* Subroutine */ int dlacpy_(char *, integer *, integer *,
50 doublereal *, integer *, doublereal *, integer *),
51 xerbla_(char *, integer *);
52 logical wantbh;
53 extern /* Subroutine */ int dtrexc_(char *, integer *, doublereal *,
54 integer *, doublereal *, integer *, integer *, integer *,
55 doublereal *, integer *);
56 integer liwmin;
57 logical wantsp, lquery;
58 extern /* Subroutine */ int dtrsyl_(char *, char *, integer *, integer *,
59 integer *, doublereal *, integer *, doublereal *, integer *,
60 doublereal *, integer *, doublereal *, integer *);
61
62
63 /* -- LAPACK routine (version 3.2) -- */
64 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
65 /* November 2006 */
66
67 /* .. Scalar Arguments .. */
68 /* .. */
69 /* .. Array Arguments .. */
70 /* .. */
71
72 /* Purpose */
73 /* ======= */
74
75 /* DTRSEN reorders the real Schur factorization of a real matrix */
76 /* A = Q*T*Q**T, so that a selected cluster of eigenvalues appears in */
77 /* the leading diagonal blocks of the upper quasi-triangular matrix T, */
78 /* and the leading columns of Q form an orthonormal basis of the */
79 /* corresponding right invariant subspace. */
80
81 /* Optionally the routine computes the reciprocal condition numbers of */
82 /* the cluster of eigenvalues and/or the invariant subspace. */
83
84 /* T must be in Schur canonical form (as returned by DHSEQR), that is, */
85 /* block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each */
86 /* 2-by-2 diagonal block has its diagonal elemnts equal and its */
87 /* off-diagonal elements of opposite sign. */
88
89 /* Arguments */
90 /* ========= */
91
92 /* JOB (input) CHARACTER*1 */
93 /* Specifies whether condition numbers are required for the */
94 /* cluster of eigenvalues (S) or the invariant subspace (SEP): */
95 /* = 'N': none; */
96 /* = 'E': for eigenvalues only (S); */
97 /* = 'V': for invariant subspace only (SEP); */
98 /* = 'B': for both eigenvalues and invariant subspace (S and */
99 /* SEP). */
100
101 /* COMPQ (input) CHARACTER*1 */
102 /* = 'V': update the matrix Q of Schur vectors; */
103 /* = 'N': do not update Q. */
104
105 /* SELECT (input) LOGICAL array, dimension (N) */
106 /* SELECT specifies the eigenvalues in the selected cluster. To */
107 /* select a real eigenvalue w(j), SELECT(j) must be set to */
108 /* .TRUE.. To select a complex conjugate pair of eigenvalues */
109 /* w(j) and w(j+1), corresponding to a 2-by-2 diagonal block, */
110 /* either SELECT(j) or SELECT(j+1) or both must be set to */
111 /* .TRUE.; a complex conjugate pair of eigenvalues must be */
112 /* either both included in the cluster or both excluded. */
113
114 /* N (input) INTEGER */
115 /* The order of the matrix T. N >= 0. */
116
117 /* T (input/output) DOUBLE PRECISION array, dimension (LDT,N) */
118 /* On entry, the upper quasi-triangular matrix T, in Schur */
119 /* canonical form. */
120 /* On exit, T is overwritten by the reordered matrix T, again in */
121 /* Schur canonical form, with the selected eigenvalues in the */
122 /* leading diagonal blocks. */
123
124 /* LDT (input) INTEGER */
125 /* The leading dimension of the array T. LDT >= max(1,N). */
126
127 /* Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N) */
128 /* On entry, if COMPQ = 'V', the matrix Q of Schur vectors. */
129 /* On exit, if COMPQ = 'V', Q has been postmultiplied by the */
130 /* orthogonal transformation matrix which reorders T; the */
131 /* leading M columns of Q form an orthonormal basis for the */
132 /* specified invariant subspace. */
133 /* If COMPQ = 'N', Q is not referenced. */
134
135 /* LDQ (input) INTEGER */
136 /* The leading dimension of the array Q. */
137 /* LDQ >= 1; and if COMPQ = 'V', LDQ >= N. */
138
139 /* WR (output) DOUBLE PRECISION array, dimension (N) */
140 /* WI (output) DOUBLE PRECISION array, dimension (N) */
141 /* The real and imaginary parts, respectively, of the reordered */
142 /* eigenvalues of T. The eigenvalues are stored in the same */
143 /* order as on the diagonal of T, with WR(i) = T(i,i) and, if */
144 /* T(i:i+1,i:i+1) is a 2-by-2 diagonal block, WI(i) > 0 and */
145 /* WI(i+1) = -WI(i). Note that if a complex eigenvalue is */
146 /* sufficiently ill-conditioned, then its value may differ */
147 /* significantly from its value before reordering. */
148
149 /* M (output) INTEGER */
150 /* The dimension of the specified invariant subspace. */
151 /* 0 < = M <= N. */
152
153 /* S (output) DOUBLE PRECISION */
154 /* If JOB = 'E' or 'B', S is a lower bound on the reciprocal */
155 /* condition number for the selected cluster of eigenvalues. */
156 /* S cannot underestimate the true reciprocal condition number */
157 /* by more than a factor of sqrt(N). If M = 0 or N, S = 1. */
158 /* If JOB = 'N' or 'V', S is not referenced. */
159
160 /* SEP (output) DOUBLE PRECISION */
161 /* If JOB = 'V' or 'B', SEP is the estimated reciprocal */
162 /* condition number of the specified invariant subspace. If */
163 /* M = 0 or N, SEP = norm(T). */
164 /* If JOB = 'N' or 'E', SEP is not referenced. */
165
166 /* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */
167 /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
168
169 /* LWORK (input) INTEGER */
170 /* The dimension of the array WORK. */
171 /* If JOB = 'N', LWORK >= max(1,N); */
172 /* if JOB = 'E', LWORK >= max(1,M*(N-M)); */
173 /* if JOB = 'V' or 'B', LWORK >= max(1,2*M*(N-M)). */
174
175 /* If LWORK = -1, then a workspace query is assumed; the routine */
176 /* only calculates the optimal size of the WORK array, returns */
177 /* this value as the first entry of the WORK array, and no error */
178 /* message related to LWORK is issued by XERBLA. */
179
180 /* IWORK (workspace) INTEGER array, dimension (MAX(1,LIWORK)) */
181 /* On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. */
182
183 /* LIWORK (input) INTEGER */
184 /* The dimension of the array IWORK. */
185 /* If JOB = 'N' or 'E', LIWORK >= 1; */
186 /* if JOB = 'V' or 'B', LIWORK >= max(1,M*(N-M)). */
187
188 /* If LIWORK = -1, then a workspace query is assumed; the */
189 /* routine only calculates the optimal size of the IWORK array, */
190 /* returns this value as the first entry of the IWORK array, and */
191 /* no error message related to LIWORK is issued by XERBLA. */
192
193 /* INFO (output) INTEGER */
194 /* = 0: successful exit */
195 /* < 0: if INFO = -i, the i-th argument had an illegal value */
196 /* = 1: reordering of T failed because some eigenvalues are too */
197 /* close to separate (the problem is very ill-conditioned); */
198 /* T may have been partially reordered, and WR and WI */
199 /* contain the eigenvalues in the same order as in T; S and */
200 /* SEP (if requested) are set to zero. */
201
202 /* Further Details */
203 /* =============== */
204
205 /* DTRSEN first collects the selected eigenvalues by computing an */
206 /* orthogonal transformation Z to move them to the top left corner of T. */
207 /* In other words, the selected eigenvalues are the eigenvalues of T11 */
208 /* in: */
209
210 /* Z'*T*Z = ( T11 T12 ) n1 */
211 /* ( 0 T22 ) n2 */
212 /* n1 n2 */
213
214 /* where N = n1+n2 and Z' means the transpose of Z. The first n1 columns */
215 /* of Z span the specified invariant subspace of T. */
216
217 /* If T has been obtained from the real Schur factorization of a matrix */
218 /* A = Q*T*Q', then the reordered real Schur factorization of A is given */
219 /* by A = (Q*Z)*(Z'*T*Z)*(Q*Z)', and the first n1 columns of Q*Z span */
220 /* the corresponding invariant subspace of A. */
221
222 /* The reciprocal condition number of the average of the eigenvalues of */
223 /* T11 may be returned in S. S lies between 0 (very badly conditioned) */
224 /* and 1 (very well conditioned). It is computed as follows. First we */
225 /* compute R so that */
226
227 /* P = ( I R ) n1 */
228 /* ( 0 0 ) n2 */
229 /* n1 n2 */
230
231 /* is the projector on the invariant subspace associated with T11. */
232 /* R is the solution of the Sylvester equation: */
233
234 /* T11*R - R*T22 = T12. */
235
236 /* Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote */
237 /* the two-norm of M. Then S is computed as the lower bound */
238
239 /* (1 + F-norm(R)**2)**(-1/2) */
240
241 /* on the reciprocal of 2-norm(P), the true reciprocal condition number. */
242 /* S cannot underestimate 1 / 2-norm(P) by more than a factor of */
243 /* sqrt(N). */
244
245 /* An approximate error bound for the computed average of the */
246 /* eigenvalues of T11 is */
247
248 /* EPS * norm(T) / S */
249
250 /* where EPS is the machine precision. */
251
252 /* The reciprocal condition number of the right invariant subspace */
253 /* spanned by the first n1 columns of Z (or of Q*Z) is returned in SEP. */
254 /* SEP is defined as the separation of T11 and T22: */
255
256 /* sep( T11, T22 ) = sigma-min( C ) */
257
258 /* where sigma-min(C) is the smallest singular value of the */
259 /* n1*n2-by-n1*n2 matrix */
260
261 /* C = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) ) */
262
263 /* I(m) is an m by m identity matrix, and kprod denotes the Kronecker */
264 /* product. We estimate sigma-min(C) by the reciprocal of an estimate of */
265 /* the 1-norm of inverse(C). The true reciprocal 1-norm of inverse(C) */
266 /* cannot differ from sigma-min(C) by more than a factor of sqrt(n1*n2). */
267
268 /* When SEP is small, small changes in T can cause large changes in */
269 /* the invariant subspace. An approximate bound on the maximum angular */
270 /* error in the computed right invariant subspace is */
271
272 /* EPS * norm(T) / SEP */
273
274 /* ===================================================================== */
275
276 /* .. Parameters .. */
277 /* .. */
278 /* .. Local Scalars .. */
279 /* .. */
280 /* .. Local Arrays .. */
281 /* .. */
282 /* .. External Functions .. */
283 /* .. */
284 /* .. External Subroutines .. */
285 /* .. */
286 /* .. Intrinsic Functions .. */
287 /* .. */
288 /* .. Executable Statements .. */
289
290 /* Decode and test the input parameters */
291
292 /* Parameter adjustments */
293 --select;
294 t_dim1 = *ldt;
295 t_offset = 1 + t_dim1;
296 t -= t_offset;
297 q_dim1 = *ldq;
298 q_offset = 1 + q_dim1;
299 q -= q_offset;
300 --wr;
301 --wi;
302 --work;
303 --iwork;
304
305 /* Function Body */
306 wantbh = lsame_(job, "B");
307 wants = lsame_(job, "E") || wantbh;
308 wantsp = lsame_(job, "V") || wantbh;
309 wantq = lsame_(compq, "V");
310
311 *info = 0;
312 lquery = *lwork == -1;
313 if (! lsame_(job, "N") && ! wants && ! wantsp) {
314 *info = -1;
315 } else if (! lsame_(compq, "N") && ! wantq) {
316 *info = -2;
317 } else if (*n < 0) {
318 *info = -4;
319 } else if (*ldt < max(1,*n)) {
320 *info = -6;
321 } else if (*ldq < 1 || wantq && *ldq < *n) {
322 *info = -8;
323 } else {
324
325 /* Set M to the dimension of the specified invariant subspace, */
326 /* and test LWORK and LIWORK. */
327
328 *m = 0;
329 pair = FALSE_;
330 i__1 = *n;
331 for (k = 1; k <= i__1; ++k) {
332 if (pair) {
333 pair = FALSE_;
334 } else {
335 if (k < *n) {
336 if (t[k + 1 + k * t_dim1] == 0.) {
337 if (select[k]) {
338 ++(*m);
339 }
340 } else {
341 pair = TRUE_;
342 if (select[k] || select[k + 1]) {
343 *m += 2;
344 }
345 }
346 } else {
347 if (select[*n]) {
348 ++(*m);
349 }
350 }
351 }
352 /* L10: */
353 }
354
355 n1 = *m;
356 n2 = *n - *m;
357 nn = n1 * n2;
358
359 if (wantsp) {
360 /* Computing MAX */
361 i__1 = 1, i__2 = nn << 1;
362 lwmin = max(i__1,i__2);
363 liwmin = max(1,nn);
364 } else if (lsame_(job, "N")) {
365 lwmin = max(1,*n);
366 liwmin = 1;
367 } else if (lsame_(job, "E")) {
368 lwmin = max(1,nn);
369 liwmin = 1;
370 }
371
372 if (*lwork < lwmin && ! lquery) {
373 *info = -15;
374 } else if (*liwork < liwmin && ! lquery) {
375 *info = -17;
376 }
377 }
378
379 if (*info == 0) {
380 work[1] = (doublereal) lwmin;
381 iwork[1] = liwmin;
382 }
383
384 if (*info != 0) {
385 i__1 = -(*info);
386 xerbla_("DTRSEN", &i__1);
387 return 0;
388 } else if (lquery) {
389 return 0;
390 }
391
392 /* Quick return if possible. */
393
394 if (*m == *n || *m == 0) {
395 if (wants) {
396 *s = 1.;
397 }
398 if (wantsp) {
399 *sep = dlange_("1", n, n, &t[t_offset], ldt, &work[1]);
400 }
401 goto L40;
402 }
403
404 /* Collect the selected blocks at the top-left corner of T. */
405
406 ks = 0;
407 pair = FALSE_;
408 i__1 = *n;
409 for (k = 1; k <= i__1; ++k) {
410 if (pair) {
411 pair = FALSE_;
412 } else {
413 swap = select[k];
414 if (k < *n) {
415 if (t[k + 1 + k * t_dim1] != 0.) {
416 pair = TRUE_;
417 swap = swap || select[k + 1];
418 }
419 }
420 if (swap) {
421 ++ks;
422
423 /* Swap the K-th block to position KS. */
424
425 ierr = 0;
426 kk = k;
427 if (k != ks) {
428 dtrexc_(compq, n, &t[t_offset], ldt, &q[q_offset], ldq, &
429 kk, &ks, &work[1], &ierr);
430 }
431 if (ierr == 1 || ierr == 2) {
432
433 /* Blocks too close to swap: exit. */
434
435 *info = 1;
436 if (wants) {
437 *s = 0.;
438 }
439 if (wantsp) {
440 *sep = 0.;
441 }
442 goto L40;
443 }
444 if (pair) {
445 ++ks;
446 }
447 }
448 }
449 /* L20: */
450 }
451
452 if (wants) {
453
454 /* Solve Sylvester equation for R: */
455
456 /* T11*R - R*T22 = scale*T12 */
457
458 dlacpy_("F", &n1, &n2, &t[(n1 + 1) * t_dim1 + 1], ldt, &work[1], &n1);
459 dtrsyl_("N", "N", &c_n1, &n1, &n2, &t[t_offset], ldt, &t[n1 + 1 + (n1
460 + 1) * t_dim1], ldt, &work[1], &n1, &scale, &ierr);
461
462 /* Estimate the reciprocal of the condition number of the cluster */
463 /* of eigenvalues. */
464
465 rnorm = dlange_("F", &n1, &n2, &work[1], &n1, &work[1]);
466 if (rnorm == 0.) {
467 *s = 1.;
468 } else {
469 *s = scale / (sqrt(scale * scale / rnorm + rnorm) * sqrt(rnorm));
470 }
471 }
472
473 if (wantsp) {
474
475 /* Estimate sep(T11,T22). */
476
477 est = 0.;
478 kase = 0;
479 L30:
480 dlacn2_(&nn, &work[nn + 1], &work[1], &iwork[1], &est, &kase, isave);
481 if (kase != 0) {
482 if (kase == 1) {
483
484 /* Solve T11*R - R*T22 = scale*X. */
485
486 dtrsyl_("N", "N", &c_n1, &n1, &n2, &t[t_offset], ldt, &t[n1 +
487 1 + (n1 + 1) * t_dim1], ldt, &work[1], &n1, &scale, &
488 ierr);
489 } else {
490
491 /* Solve T11'*R - R*T22' = scale*X. */
492
493 dtrsyl_("T", "T", &c_n1, &n1, &n2, &t[t_offset], ldt, &t[n1 +
494 1 + (n1 + 1) * t_dim1], ldt, &work[1], &n1, &scale, &
495 ierr);
496 }
497 goto L30;
498 }
499
500 *sep = scale / est;
501 }
502
503 L40:
504
505 /* Store the output eigenvalues in WR and WI. */
506
507 i__1 = *n;
508 for (k = 1; k <= i__1; ++k) {
509 wr[k] = t[k + k * t_dim1];
510 wi[k] = 0.;
511 /* L50: */
512 }
513 i__1 = *n - 1;
514 for (k = 1; k <= i__1; ++k) {
515 if (t[k + 1 + k * t_dim1] != 0.) {
516 wi[k] = sqrt((d__1 = t[k + (k + 1) * t_dim1], abs(d__1))) * sqrt((
517 d__2 = t[k + 1 + k * t_dim1], abs(d__2)));
518 wi[k + 1] = -wi[k];
519 }
520 /* L60: */
521 }
522
523 work[1] = (doublereal) lwmin;
524 iwork[1] = liwmin;
525
526 return 0;
527
528 /* End of DTRSEN */
529
530 } /* dtrsen_ */
531