1 /* ./src_f77/ctgsna.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
8 /* Table of constant values */
9
10 static integer c__1 = 1;
11 static complex c_b19 = {1.f,0.f};
12 static complex c_b20 = {0.f,0.f};
13 static logical c_false = FALSE_;
14 static integer c__3 = 3;
15
ctgsna_(char * job,char * howmny,logical * select,integer * n,complex * a,integer * lda,complex * b,integer * ldb,complex * vl,integer * ldvl,complex * vr,integer * ldvr,real * s,real * dif,integer * mm,integer * m,complex * work,integer * lwork,integer * iwork,integer * info,ftnlen job_len,ftnlen howmny_len)16 /* Subroutine */ int ctgsna_(char *job, char *howmny, logical *select,
17 integer *n, complex *a, integer *lda, complex *b, integer *ldb,
18 complex *vl, integer *ldvl, complex *vr, integer *ldvr, real *s, real
19 *dif, integer *mm, integer *m, complex *work, integer *lwork, integer
20 *iwork, integer *info, ftnlen job_len, ftnlen howmny_len)
21 {
22 /* System generated locals */
23 integer a_dim1, a_offset, b_dim1, b_offset, vl_dim1, vl_offset, vr_dim1,
24 vr_offset, i__1, i__2;
25 real r__1, r__2;
26 complex q__1;
27
28 /* Builtin functions */
29 double c_abs(complex *);
30
31 /* Local variables */
32 static integer i__, k, n1, n2, ks;
33 static real eps, cond;
34 static integer ierr, ifst;
35 static real lnrm;
36 static complex yhax, yhbx;
37 static integer ilst;
38 static real rnrm, scale;
39 extern /* Complex */ VOID cdotc_(complex *, integer *, complex *, integer
40 *, complex *, integer *);
41 extern logical lsame_(char *, char *, ftnlen, ftnlen);
42 extern /* Subroutine */ int cgemv_(char *, integer *, integer *, complex *
43 , complex *, integer *, complex *, integer *, complex *, complex *
44 , integer *, ftnlen);
45 static integer lwmin;
46 static logical wants;
47 static integer llwrk;
48 static complex dummy[1];
49 extern doublereal scnrm2_(integer *, complex *, integer *), slapy2_(real *
50 , real *);
51 static complex dummy1[1];
52 extern /* Subroutine */ int slabad_(real *, real *);
53 extern doublereal slamch_(char *, ftnlen);
54 extern /* Subroutine */ int clacpy_(char *, integer *, integer *, complex
55 *, integer *, complex *, integer *, ftnlen), ctgexc_(logical *,
56 logical *, integer *, complex *, integer *, complex *, integer *,
57 complex *, integer *, complex *, integer *, integer *, integer *,
58 integer *), xerbla_(char *, integer *, ftnlen);
59 static real bignum;
60 static logical wantbh, wantdf, somcon;
61 extern /* Subroutine */ int ctgsyl_(char *, integer *, integer *, integer
62 *, complex *, integer *, complex *, integer *, complex *, integer
63 *, complex *, integer *, complex *, integer *, complex *, integer
64 *, real *, real *, complex *, integer *, integer *, integer *,
65 ftnlen);
66 static real smlnum;
67 static logical lquery;
68
69
70 /* -- LAPACK routine (version 3.0) -- */
71 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
72 /* Courant Institute, Argonne National Lab, and Rice University */
73 /* June 30, 1999 */
74
75 /* .. Scalar Arguments .. */
76 /* .. */
77 /* .. Array Arguments .. */
78 /* .. */
79
80 /* Purpose */
81 /* ======= */
82
83 /* CTGSNA estimates reciprocal condition numbers for specified */
84 /* eigenvalues and/or eigenvectors of a matrix pair (A, B). */
85
86 /* (A, B) must be in generalized Schur canonical form, that is, A and */
87 /* B are both upper triangular. */
88
89 /* Arguments */
90 /* ========= */
91
92 /* JOB (input) CHARACTER*1 */
93 /* Specifies whether condition numbers are required for */
94 /* eigenvalues (S) or eigenvectors (DIF): */
95 /* = 'E': for eigenvalues only (S); */
96 /* = 'V': for eigenvectors only (DIF); */
97 /* = 'B': for both eigenvalues and eigenvectors (S and DIF). */
98
99 /* HOWMNY (input) CHARACTER*1 */
100 /* = 'A': compute condition numbers for all eigenpairs; */
101 /* = 'S': compute condition numbers for selected eigenpairs */
102 /* specified by the array SELECT. */
103
104 /* SELECT (input) LOGICAL array, dimension (N) */
105 /* If HOWMNY = 'S', SELECT specifies the eigenpairs for which */
106 /* condition numbers are required. To select condition numbers */
107 /* for the corresponding j-th eigenvalue and/or eigenvector, */
108 /* SELECT(j) must be set to .TRUE.. */
109 /* If HOWMNY = 'A', SELECT is not referenced. */
110
111 /* N (input) INTEGER */
112 /* The order of the square matrix pair (A, B). N >= 0. */
113
114 /* A (input) COMPLEX array, dimension (LDA,N) */
115 /* The upper triangular matrix A in the pair (A,B). */
116
117 /* LDA (input) INTEGER */
118 /* The leading dimension of the array A. LDA >= max(1,N). */
119
120 /* B (input) COMPLEX array, dimension (LDB,N) */
121 /* The upper triangular matrix B in the pair (A, B). */
122
123 /* LDB (input) INTEGER */
124 /* The leading dimension of the array B. LDB >= max(1,N). */
125
126 /* VL (input) COMPLEX array, dimension (LDVL,M) */
127 /* IF JOB = 'E' or 'B', VL must contain left eigenvectors of */
128 /* (A, B), corresponding to the eigenpairs specified by HOWMNY */
129 /* and SELECT. The eigenvectors must be stored in consecutive */
130 /* columns of VL, as returned by CTGEVC. */
131 /* If JOB = 'V', VL is not referenced. */
132
133 /* LDVL (input) INTEGER */
134 /* The leading dimension of the array VL. LDVL >= 1; and */
135 /* If JOB = 'E' or 'B', LDVL >= N. */
136
137 /* VR (input) COMPLEX array, dimension (LDVR,M) */
138 /* IF JOB = 'E' or 'B', VR must contain right eigenvectors of */
139 /* (A, B), corresponding to the eigenpairs specified by HOWMNY */
140 /* and SELECT. The eigenvectors must be stored in consecutive */
141 /* columns of VR, as returned by CTGEVC. */
142 /* If JOB = 'V', VR is not referenced. */
143
144 /* LDVR (input) INTEGER */
145 /* The leading dimension of the array VR. LDVR >= 1; */
146 /* If JOB = 'E' or 'B', LDVR >= N. */
147
148 /* S (output) REAL array, dimension (MM) */
149 /* If JOB = 'E' or 'B', the reciprocal condition numbers of the */
150 /* selected eigenvalues, stored in consecutive elements of the */
151 /* array. */
152 /* If JOB = 'V', S is not referenced. */
153
154 /* DIF (output) REAL array, dimension (MM) */
155 /* If JOB = 'V' or 'B', the estimated reciprocal condition */
156 /* numbers of the selected eigenvectors, stored in consecutive */
157 /* elements of the array. */
158 /* If the eigenvalues cannot be reordered to compute DIF(j), */
159 /* DIF(j) is set to 0; this can only occur when the true value */
160 /* would be very small anyway. */
161 /* For each eigenvalue/vector specified by SELECT, DIF stores */
162 /* a Frobenius norm-based estimate of Difl. */
163 /* If JOB = 'E', DIF is not referenced. */
164
165 /* MM (input) INTEGER */
166 /* The number of elements in the arrays S and DIF. MM >= M. */
167
168 /* M (output) INTEGER */
169 /* The number of elements of the arrays S and DIF used to store */
170 /* the specified condition numbers; for each selected eigenvalue */
171 /* one element is used. If HOWMNY = 'A', M is set to N. */
172
173 /* WORK (workspace/output) COMPLEX array, dimension (LWORK) */
174 /* If JOB = 'E', WORK is not referenced. Otherwise, */
175 /* on exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
176
177 /* LWORK (input) INTEGER */
178 /* The dimension of the array WORK. LWORK >= 1. */
179 /* If JOB = 'V' or 'B', LWORK >= 2*N*N. */
180
181 /* IWORK (workspace) INTEGER array, dimension (N+2) */
182 /* If JOB = 'E', IWORK is not referenced. */
183
184 /* INFO (output) INTEGER */
185 /* = 0: Successful exit */
186 /* < 0: If INFO = -i, the i-th argument had an illegal value */
187
188 /* Further Details */
189 /* =============== */
190
191 /* The reciprocal of the condition number of the i-th generalized */
192 /* eigenvalue w = (a, b) is defined as */
193
194 /* S(I) = (|v'Au|**2 + |v'Bu|**2)**(1/2) / (norm(u)*norm(v)) */
195
196 /* where u and v are the right and left eigenvectors of (A, B) */
197 /* corresponding to w; |z| denotes the absolute value of the complex */
198 /* number, and norm(u) denotes the 2-norm of the vector u. The pair */
199 /* (a, b) corresponds to an eigenvalue w = a/b (= v'Au/v'Bu) of the */
200 /* matrix pair (A, B). If both a and b equal zero, then (A,B) is */
201 /* singular and S(I) = -1 is returned. */
202
203 /* An approximate error bound on the chordal distance between the i-th */
204 /* computed generalized eigenvalue w and the corresponding exact */
205 /* eigenvalue lambda is */
206
207 /* chord(w, lambda) <= EPS * norm(A, B) / S(I), */
208
209 /* where EPS is the machine precision. */
210
211 /* The reciprocal of the condition number of the right eigenvector u */
212 /* and left eigenvector v corresponding to the generalized eigenvalue w */
213 /* is defined as follows. Suppose */
214
215 /* (A, B) = ( a * ) ( b * ) 1 */
216 /* ( 0 A22 ),( 0 B22 ) n-1 */
217 /* 1 n-1 1 n-1 */
218
219 /* Then the reciprocal condition number DIF(I) is */
220
221 /* Difl[(a, b), (A22, B22)] = sigma-min( Zl ) */
222
223 /* where sigma-min(Zl) denotes the smallest singular value of */
224
225 /* Zl = [ kron(a, In-1) -kron(1, A22) ] */
226 /* [ kron(b, In-1) -kron(1, B22) ]. */
227
228 /* Here In-1 is the identity matrix of size n-1 and X' is the conjugate */
229 /* transpose of X. kron(X, Y) is the Kronecker product between the */
230 /* matrices X and Y. */
231
232 /* We approximate the smallest singular value of Zl with an upper */
233 /* bound. This is done by CLATDF. */
234
235 /* An approximate error bound for a computed eigenvector VL(i) or */
236 /* VR(i) is given by */
237
238 /* EPS * norm(A, B) / DIF(i). */
239
240 /* See ref. [2-3] for more details and further references. */
241
242 /* Based on contributions by */
243 /* Bo Kagstrom and Peter Poromaa, Department of Computing Science, */
244 /* Umea University, S-901 87 Umea, Sweden. */
245
246 /* References */
247 /* ========== */
248
249 /* [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the */
250 /* Generalized Real Schur Form of a Regular Matrix Pair (A, B), in */
251 /* M.S. Moonen et al (eds), Linear Algebra for Large Scale and */
252 /* Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218. */
253
254 /* [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified */
255 /* Eigenvalues of a Regular Matrix Pair (A, B) and Condition */
256 /* Estimation: Theory, Algorithms and Software, Report */
257 /* UMINF - 94.04, Department of Computing Science, Umea University, */
258 /* S-901 87 Umea, Sweden, 1994. Also as LAPACK Working Note 87. */
259 /* To appear in Numerical Algorithms, 1996. */
260
261 /* [3] B. Kagstrom and P. Poromaa, LAPACK-Style Algorithms and Software */
262 /* for Solving the Generalized Sylvester Equation and Estimating the */
263 /* Separation between Regular Matrix Pairs, Report UMINF - 93.23, */
264 /* Department of Computing Science, Umea University, S-901 87 Umea, */
265 /* Sweden, December 1993, Revised April 1994, Also as LAPACK Working */
266 /* Note 75. */
267 /* To appear in ACM Trans. on Math. Software, Vol 22, No 1, 1996. */
268
269 /* ===================================================================== */
270
271 /* .. Parameters .. */
272 /* .. */
273 /* .. Local Scalars .. */
274 /* .. */
275 /* .. Local Arrays .. */
276 /* .. */
277 /* .. External Functions .. */
278 /* .. */
279 /* .. External Subroutines .. */
280 /* .. */
281 /* .. Intrinsic Functions .. */
282 /* .. */
283 /* .. Executable Statements .. */
284
285 /* Decode and test the input parameters */
286
287 /* Parameter adjustments */
288 --select;
289 a_dim1 = *lda;
290 a_offset = 1 + a_dim1;
291 a -= a_offset;
292 b_dim1 = *ldb;
293 b_offset = 1 + b_dim1;
294 b -= b_offset;
295 vl_dim1 = *ldvl;
296 vl_offset = 1 + vl_dim1;
297 vl -= vl_offset;
298 vr_dim1 = *ldvr;
299 vr_offset = 1 + vr_dim1;
300 vr -= vr_offset;
301 --s;
302 --dif;
303 --work;
304 --iwork;
305
306 /* Function Body */
307 wantbh = lsame_(job, "B", (ftnlen)1, (ftnlen)1);
308 wants = lsame_(job, "E", (ftnlen)1, (ftnlen)1) || wantbh;
309 wantdf = lsame_(job, "V", (ftnlen)1, (ftnlen)1) || wantbh;
310
311 somcon = lsame_(howmny, "S", (ftnlen)1, (ftnlen)1);
312
313 *info = 0;
314 lquery = *lwork == -1;
315
316 if (lsame_(job, "V", (ftnlen)1, (ftnlen)1) || lsame_(job, "B", (ftnlen)1,
317 (ftnlen)1)) {
318 /* Computing MAX */
319 i__1 = 1, i__2 = (*n << 1) * *n;
320 lwmin = max(i__1,i__2);
321 } else {
322 lwmin = 1;
323 }
324
325 if (! wants && ! wantdf) {
326 *info = -1;
327 } else if (! lsame_(howmny, "A", (ftnlen)1, (ftnlen)1) && ! somcon) {
328 *info = -2;
329 } else if (*n < 0) {
330 *info = -4;
331 } else if (*lda < max(1,*n)) {
332 *info = -6;
333 } else if (*ldb < max(1,*n)) {
334 *info = -8;
335 } else if (wants && *ldvl < *n) {
336 *info = -10;
337 } else if (wants && *ldvr < *n) {
338 *info = -12;
339 } else {
340
341 /* Set M to the number of eigenpairs for which condition numbers */
342 /* are required, and test MM. */
343
344 if (somcon) {
345 *m = 0;
346 i__1 = *n;
347 for (k = 1; k <= i__1; ++k) {
348 if (select[k]) {
349 ++(*m);
350 }
351 /* L10: */
352 }
353 } else {
354 *m = *n;
355 }
356
357 if (*mm < *m) {
358 *info = -15;
359 } else if (*lwork < lwmin && ! lquery) {
360 *info = -18;
361 }
362 }
363
364 if (*info == 0) {
365 work[1].r = (real) lwmin, work[1].i = 0.f;
366 }
367
368 if (*info != 0) {
369 i__1 = -(*info);
370 xerbla_("CTGSNA", &i__1, (ftnlen)6);
371 return 0;
372 } else if (lquery) {
373 return 0;
374 }
375
376 /* Quick return if possible */
377
378 if (*n == 0) {
379 return 0;
380 }
381
382 /* Get machine constants */
383
384 eps = slamch_("P", (ftnlen)1);
385 smlnum = slamch_("S", (ftnlen)1) / eps;
386 bignum = 1.f / smlnum;
387 slabad_(&smlnum, &bignum);
388 llwrk = *lwork - (*n << 1) * *n;
389 ks = 0;
390 i__1 = *n;
391 for (k = 1; k <= i__1; ++k) {
392
393 /* Determine whether condition numbers are required for the k-th */
394 /* eigenpair. */
395
396 if (somcon) {
397 if (! select[k]) {
398 goto L20;
399 }
400 }
401
402 ++ks;
403
404 if (wants) {
405
406 /* Compute the reciprocal condition number of the k-th */
407 /* eigenvalue. */
408
409 rnrm = scnrm2_(n, &vr[ks * vr_dim1 + 1], &c__1);
410 lnrm = scnrm2_(n, &vl[ks * vl_dim1 + 1], &c__1);
411 cgemv_("N", n, n, &c_b19, &a[a_offset], lda, &vr[ks * vr_dim1 + 1]
412 , &c__1, &c_b20, &work[1], &c__1, (ftnlen)1);
413 cdotc_(&q__1, n, &work[1], &c__1, &vl[ks * vl_dim1 + 1], &c__1);
414 yhax.r = q__1.r, yhax.i = q__1.i;
415 cgemv_("N", n, n, &c_b19, &b[b_offset], ldb, &vr[ks * vr_dim1 + 1]
416 , &c__1, &c_b20, &work[1], &c__1, (ftnlen)1);
417 cdotc_(&q__1, n, &work[1], &c__1, &vl[ks * vl_dim1 + 1], &c__1);
418 yhbx.r = q__1.r, yhbx.i = q__1.i;
419 r__1 = c_abs(&yhax);
420 r__2 = c_abs(&yhbx);
421 cond = slapy2_(&r__1, &r__2);
422 if (cond == 0.f) {
423 s[ks] = -1.f;
424 } else {
425 s[ks] = cond / (rnrm * lnrm);
426 }
427 }
428
429 if (wantdf) {
430 if (*n == 1) {
431 r__1 = c_abs(&a[a_dim1 + 1]);
432 r__2 = c_abs(&b[b_dim1 + 1]);
433 dif[ks] = slapy2_(&r__1, &r__2);
434 goto L20;
435 }
436
437 /* Estimate the reciprocal condition number of the k-th */
438 /* eigenvectors. */
439
440 /* Copy the matrix (A, B) to the array WORK and move the */
441 /* (k,k)th pair to the (1,1) position. */
442
443 clacpy_("Full", n, n, &a[a_offset], lda, &work[1], n, (ftnlen)4);
444 clacpy_("Full", n, n, &b[b_offset], ldb, &work[*n * *n + 1], n, (
445 ftnlen)4);
446 ifst = k;
447 ilst = 1;
448
449 ctgexc_(&c_false, &c_false, n, &work[1], n, &work[*n * *n + 1], n,
450 dummy, &c__1, dummy1, &c__1, &ifst, &ilst, &ierr);
451
452 if (ierr > 0) {
453
454 /* Ill-conditioned problem - swap rejected. */
455
456 dif[ks] = 0.f;
457 } else {
458
459 /* Reordering successful, solve generalized Sylvester */
460 /* equation for R and L, */
461 /* A22 * R - L * A11 = A12 */
462 /* B22 * R - L * B11 = B12, */
463 /* and compute estimate of Difl[(A11,B11), (A22, B22)]. */
464
465 n1 = 1;
466 n2 = *n - n1;
467 i__ = *n * *n + 1;
468 ctgsyl_("N", &c__3, &n2, &n1, &work[*n * n1 + n1 + 1], n, &
469 work[1], n, &work[n1 + 1], n, &work[*n * n1 + n1 +
470 i__], n, &work[i__], n, &work[n1 + i__], n, &scale, &
471 dif[ks], &work[(*n * *n << 1) + 1], &llwrk, &iwork[1],
472 &ierr, (ftnlen)1);
473 }
474 }
475
476 L20:
477 ;
478 }
479 work[1].r = (real) lwmin, work[1].i = 0.f;
480 return 0;
481
482 /* End of CTGSNA */
483
484 } /* ctgsna_ */
485
486