1 /* ./src_f77/ctgsy2.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__2 = 2;
11 static integer c__1 = 1;
12
ctgsy2_(char * trans,integer * ijob,integer * m,integer * n,complex * a,integer * lda,complex * b,integer * ldb,complex * c__,integer * ldc,complex * d__,integer * ldd,complex * e,integer * lde,complex * f,integer * ldf,real * scale,real * rdsum,real * rdscal,integer * info,ftnlen trans_len)13 /* Subroutine */ int ctgsy2_(char *trans, integer *ijob, integer *m, integer *
14 n, complex *a, integer *lda, complex *b, integer *ldb, complex *c__,
15 integer *ldc, complex *d__, integer *ldd, complex *e, integer *lde,
16 complex *f, integer *ldf, real *scale, real *rdsum, real *rdscal,
17 integer *info, ftnlen trans_len)
18 {
19 /* System generated locals */
20 integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, d_dim1,
21 d_offset, e_dim1, e_offset, f_dim1, f_offset, i__1, i__2, i__3,
22 i__4;
23 complex q__1, q__2, q__3, q__4, q__5, q__6;
24
25 /* Builtin functions */
26 void r_cnjg(complex *, complex *);
27
28 /* Local variables */
29 static integer i__, j, k;
30 static complex z__[4] /* was [2][2] */, rhs[2];
31 static integer ierr, ipiv[2], jpiv[2];
32 static complex alpha;
33 extern /* Subroutine */ int cscal_(integer *, complex *, complex *,
34 integer *);
35 extern logical lsame_(char *, char *, ftnlen, ftnlen);
36 extern /* Subroutine */ int caxpy_(integer *, complex *, complex *,
37 integer *, complex *, integer *), cgesc2_(integer *, complex *,
38 integer *, complex *, integer *, integer *, real *), cgetc2_(
39 integer *, complex *, integer *, integer *, integer *, integer *),
40 clatdf_(integer *, integer *, complex *, integer *, complex *,
41 real *, real *, integer *, integer *);
42 static real scaloc;
43 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
44 static logical notran;
45
46
47 /* -- LAPACK auxiliary routine (version 3.0) -- */
48 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
49 /* Courant Institute, Argonne National Lab, and Rice University */
50 /* June 30, 1999 */
51
52 /* .. Scalar Arguments .. */
53 /* .. */
54 /* .. Array Arguments .. */
55 /* .. */
56
57 /* Purpose */
58 /* ======= */
59
60 /* CTGSY2 solves the generalized Sylvester equation */
61
62 /* A * R - L * B = scale * C (1) */
63 /* D * R - L * E = scale * F */
64
65 /* using Level 1 and 2 BLAS, where R and L are unknown M-by-N matrices, */
66 /* (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M, */
67 /* N-by-N and M-by-N, respectively. A, B, D and E are upper triangular */
68 /* (i.e., (A,D) and (B,E) in generalized Schur form). */
69
70 /* The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an output */
71 /* scaling factor chosen to avoid overflow. */
72
73 /* In matrix notation solving equation (1) corresponds to solve */
74 /* Zx = scale * b, where Z is defined as */
75
76 /* Z = [ kron(In, A) -kron(B', Im) ] (2) */
77 /* [ kron(In, D) -kron(E', Im) ], */
78
79 /* Ik is the identity matrix of size k and X' is the transpose of X. */
80 /* kron(X, Y) is the Kronecker product between the matrices X and Y. */
81
82 /* If TRANS = 'C', y in the conjugate transposed system Z'y = scale*b */
83 /* is solved for, which is equivalent to solve for R and L in */
84
85 /* A' * R + D' * L = scale * C (3) */
86 /* R * B' + L * E' = scale * -F */
87
88 /* This case is used to compute an estimate of Dif[(A, D), (B, E)] = */
89 /* = sigma_min(Z) using reverse communicaton with CLACON. */
90
91 /* CTGSY2 also (IJOB >= 1) contributes to the computation in CTGSYL */
92 /* of an upper bound on the separation between to matrix pairs. Then */
93 /* the input (A, D), (B, E) are sub-pencils of two matrix pairs in */
94 /* CTGSYL. */
95
96 /* Arguments */
97 /* ========= */
98
99 /* TRANS (input) CHARACTER */
100 /* = 'N', solve the generalized Sylvester equation (1). */
101 /* = 'T': solve the 'transposed' system (3). */
102
103 /* IJOB (input) INTEGER */
104 /* Specifies what kind of functionality to be performed. */
105 /* =0: solve (1) only. */
106 /* =1: A contribution from this subsystem to a Frobenius */
107 /* norm-based estimate of the separation between two matrix */
108 /* pairs is computed. (look ahead strategy is used). */
109 /* =2: A contribution from this subsystem to a Frobenius */
110 /* norm-based estimate of the separation between two matrix */
111 /* pairs is computed. (SGECON on sub-systems is used.) */
112 /* Not referenced if TRANS = 'T'. */
113
114 /* M (input) INTEGER */
115 /* On entry, M specifies the order of A and D, and the row */
116 /* dimension of C, F, R and L. */
117
118 /* N (input) INTEGER */
119 /* On entry, N specifies the order of B and E, and the column */
120 /* dimension of C, F, R and L. */
121
122 /* A (input) COMPLEX array, dimension (LDA, M) */
123 /* On entry, A contains an upper triangular matrix. */
124
125 /* LDA (input) INTEGER */
126 /* The leading dimension of the matrix A. LDA >= max(1, M). */
127
128 /* B (input) COMPLEX array, dimension (LDB, N) */
129 /* On entry, B contains an upper triangular matrix. */
130
131 /* LDB (input) INTEGER */
132 /* The leading dimension of the matrix B. LDB >= max(1, N). */
133
134 /* C (input/ output) COMPLEX array, dimension (LDC, N) */
135 /* On entry, C contains the right-hand-side of the first matrix */
136 /* equation in (1). */
137 /* On exit, if IJOB = 0, C has been overwritten by the solution */
138 /* R. */
139
140 /* LDC (input) INTEGER */
141 /* The leading dimension of the matrix C. LDC >= max(1, M). */
142
143 /* D (input) COMPLEX array, dimension (LDD, M) */
144 /* On entry, D contains an upper triangular matrix. */
145
146 /* LDD (input) INTEGER */
147 /* The leading dimension of the matrix D. LDD >= max(1, M). */
148
149 /* E (input) COMPLEX array, dimension (LDE, N) */
150 /* On entry, E contains an upper triangular matrix. */
151
152 /* LDE (input) INTEGER */
153 /* The leading dimension of the matrix E. LDE >= max(1, N). */
154
155 /* F (input/ output) COMPLEX array, dimension (LDF, N) */
156 /* On entry, F contains the right-hand-side of the second matrix */
157 /* equation in (1). */
158 /* On exit, if IJOB = 0, F has been overwritten by the solution */
159 /* L. */
160
161 /* LDF (input) INTEGER */
162 /* The leading dimension of the matrix F. LDF >= max(1, M). */
163
164 /* SCALE (output) REAL */
165 /* On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions */
166 /* R and L (C and F on entry) will hold the solutions to a */
167 /* slightly perturbed system but the input matrices A, B, D and */
168 /* E have not been changed. If SCALE = 0, R and L will hold the */
169 /* solutions to the homogeneous system with C = F = 0. */
170 /* Normally, SCALE = 1. */
171
172 /* RDSUM (input/output) REAL */
173 /* On entry, the sum of squares of computed contributions to */
174 /* the Dif-estimate under computation by CTGSYL, where the */
175 /* scaling factor RDSCAL (see below) has been factored out. */
176 /* On exit, the corresponding sum of squares updated with the */
177 /* contributions from the current sub-system. */
178 /* If TRANS = 'T' RDSUM is not touched. */
179 /* NOTE: RDSUM only makes sense when CTGSY2 is called by */
180 /* CTGSYL. */
181
182 /* RDSCAL (input/output) REAL */
183 /* On entry, scaling factor used to prevent overflow in RDSUM. */
184 /* On exit, RDSCAL is updated w.r.t. the current contributions */
185 /* in RDSUM. */
186 /* If TRANS = 'T', RDSCAL is not touched. */
187 /* NOTE: RDSCAL only makes sense when CTGSY2 is called by */
188 /* CTGSYL. */
189
190 /* INFO (output) INTEGER */
191 /* On exit, if INFO is set to */
192 /* =0: Successful exit */
193 /* <0: If INFO = -i, input argument number i is illegal. */
194 /* >0: The matrix pairs (A, D) and (B, E) have common or very */
195 /* close eigenvalues. */
196
197 /* Further Details */
198 /* =============== */
199
200 /* Based on contributions by */
201 /* Bo Kagstrom and Peter Poromaa, Department of Computing Science, */
202 /* Umea University, S-901 87 Umea, Sweden. */
203
204 /* ===================================================================== */
205
206 /* .. Parameters .. */
207 /* .. */
208 /* .. Local Scalars .. */
209 /* .. */
210 /* .. Local Arrays .. */
211 /* .. */
212 /* .. External Functions .. */
213 /* .. */
214 /* .. External Subroutines .. */
215 /* .. */
216 /* .. Intrinsic Functions .. */
217 /* .. */
218 /* .. Executable Statements .. */
219
220 /* Decode and test input parameters */
221
222 /* Parameter adjustments */
223 a_dim1 = *lda;
224 a_offset = 1 + a_dim1;
225 a -= a_offset;
226 b_dim1 = *ldb;
227 b_offset = 1 + b_dim1;
228 b -= b_offset;
229 c_dim1 = *ldc;
230 c_offset = 1 + c_dim1;
231 c__ -= c_offset;
232 d_dim1 = *ldd;
233 d_offset = 1 + d_dim1;
234 d__ -= d_offset;
235 e_dim1 = *lde;
236 e_offset = 1 + e_dim1;
237 e -= e_offset;
238 f_dim1 = *ldf;
239 f_offset = 1 + f_dim1;
240 f -= f_offset;
241
242 /* Function Body */
243 *info = 0;
244 ierr = 0;
245 notran = lsame_(trans, "N", (ftnlen)1, (ftnlen)1);
246 if (! notran && ! lsame_(trans, "C", (ftnlen)1, (ftnlen)1)) {
247 *info = -1;
248 } else if (*ijob < 0 || *ijob > 2) {
249 *info = -2;
250 } else if (*m <= 0) {
251 *info = -3;
252 } else if (*n <= 0) {
253 *info = -4;
254 } else if (*lda < max(1,*m)) {
255 *info = -5;
256 } else if (*ldb < max(1,*n)) {
257 *info = -8;
258 } else if (*ldc < max(1,*m)) {
259 *info = -10;
260 } else if (*ldd < max(1,*m)) {
261 *info = -12;
262 } else if (*lde < max(1,*n)) {
263 *info = -14;
264 } else if (*ldf < max(1,*m)) {
265 *info = -16;
266 }
267 if (*info != 0) {
268 i__1 = -(*info);
269 xerbla_("CTGSY2", &i__1, (ftnlen)6);
270 return 0;
271 }
272
273 if (notran) {
274
275 /* Solve (I, J) - system */
276 /* A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J) */
277 /* D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J) */
278 /* for I = M, M - 1, ..., 1; J = 1, 2, ..., N */
279
280 *scale = 1.f;
281 scaloc = 1.f;
282 i__1 = *n;
283 for (j = 1; j <= i__1; ++j) {
284 for (i__ = *m; i__ >= 1; --i__) {
285
286 /* Build 2 by 2 system */
287
288 i__2 = i__ + i__ * a_dim1;
289 z__[0].r = a[i__2].r, z__[0].i = a[i__2].i;
290 i__2 = i__ + i__ * d_dim1;
291 z__[1].r = d__[i__2].r, z__[1].i = d__[i__2].i;
292 i__2 = j + j * b_dim1;
293 q__1.r = -b[i__2].r, q__1.i = -b[i__2].i;
294 z__[2].r = q__1.r, z__[2].i = q__1.i;
295 i__2 = j + j * e_dim1;
296 q__1.r = -e[i__2].r, q__1.i = -e[i__2].i;
297 z__[3].r = q__1.r, z__[3].i = q__1.i;
298
299 /* Set up right hand side(s) */
300
301 i__2 = i__ + j * c_dim1;
302 rhs[0].r = c__[i__2].r, rhs[0].i = c__[i__2].i;
303 i__2 = i__ + j * f_dim1;
304 rhs[1].r = f[i__2].r, rhs[1].i = f[i__2].i;
305
306 /* Solve Z * x = RHS */
307
308 cgetc2_(&c__2, z__, &c__2, ipiv, jpiv, &ierr);
309 if (ierr > 0) {
310 *info = ierr;
311 }
312 if (*ijob == 0) {
313 cgesc2_(&c__2, z__, &c__2, rhs, ipiv, jpiv, &scaloc);
314 if (scaloc != 1.f) {
315 i__2 = *n;
316 for (k = 1; k <= i__2; ++k) {
317 q__1.r = scaloc, q__1.i = 0.f;
318 cscal_(m, &q__1, &c__[k * c_dim1 + 1], &c__1);
319 q__1.r = scaloc, q__1.i = 0.f;
320 cscal_(m, &q__1, &f[k * f_dim1 + 1], &c__1);
321 /* L10: */
322 }
323 *scale *= scaloc;
324 }
325 } else {
326 clatdf_(ijob, &c__2, z__, &c__2, rhs, rdsum, rdscal, ipiv,
327 jpiv);
328 }
329
330 /* Unpack solution vector(s) */
331
332 i__2 = i__ + j * c_dim1;
333 c__[i__2].r = rhs[0].r, c__[i__2].i = rhs[0].i;
334 i__2 = i__ + j * f_dim1;
335 f[i__2].r = rhs[1].r, f[i__2].i = rhs[1].i;
336
337 /* Substitute R(I, J) and L(I, J) into remaining equation. */
338
339 if (i__ > 1) {
340 q__1.r = -rhs[0].r, q__1.i = -rhs[0].i;
341 alpha.r = q__1.r, alpha.i = q__1.i;
342 i__2 = i__ - 1;
343 caxpy_(&i__2, &alpha, &a[i__ * a_dim1 + 1], &c__1, &c__[j
344 * c_dim1 + 1], &c__1);
345 i__2 = i__ - 1;
346 caxpy_(&i__2, &alpha, &d__[i__ * d_dim1 + 1], &c__1, &f[j
347 * f_dim1 + 1], &c__1);
348 }
349 if (j < *n) {
350 i__2 = *n - j;
351 caxpy_(&i__2, &rhs[1], &b[j + (j + 1) * b_dim1], ldb, &
352 c__[i__ + (j + 1) * c_dim1], ldc);
353 i__2 = *n - j;
354 caxpy_(&i__2, &rhs[1], &e[j + (j + 1) * e_dim1], lde, &f[
355 i__ + (j + 1) * f_dim1], ldf);
356 }
357
358 /* L20: */
359 }
360 /* L30: */
361 }
362 } else {
363
364 /* Solve transposed (I, J) - system: */
365 /* A(I, I)' * R(I, J) + D(I, I)' * L(J, J) = C(I, J) */
366 /* R(I, I) * B(J, J) + L(I, J) * E(J, J) = -F(I, J) */
367 /* for I = 1, 2, ..., M, J = N, N - 1, ..., 1 */
368
369 *scale = 1.f;
370 scaloc = 1.f;
371 i__1 = *m;
372 for (i__ = 1; i__ <= i__1; ++i__) {
373 for (j = *n; j >= 1; --j) {
374
375 /* Build 2 by 2 system Z' */
376
377 r_cnjg(&q__1, &a[i__ + i__ * a_dim1]);
378 z__[0].r = q__1.r, z__[0].i = q__1.i;
379 r_cnjg(&q__2, &b[j + j * b_dim1]);
380 q__1.r = -q__2.r, q__1.i = -q__2.i;
381 z__[1].r = q__1.r, z__[1].i = q__1.i;
382 r_cnjg(&q__1, &d__[i__ + i__ * d_dim1]);
383 z__[2].r = q__1.r, z__[2].i = q__1.i;
384 r_cnjg(&q__2, &e[j + j * e_dim1]);
385 q__1.r = -q__2.r, q__1.i = -q__2.i;
386 z__[3].r = q__1.r, z__[3].i = q__1.i;
387
388
389 /* Set up right hand side(s) */
390
391 i__2 = i__ + j * c_dim1;
392 rhs[0].r = c__[i__2].r, rhs[0].i = c__[i__2].i;
393 i__2 = i__ + j * f_dim1;
394 rhs[1].r = f[i__2].r, rhs[1].i = f[i__2].i;
395
396 /* Solve Z' * x = RHS */
397
398 cgetc2_(&c__2, z__, &c__2, ipiv, jpiv, &ierr);
399 if (ierr > 0) {
400 *info = ierr;
401 }
402 cgesc2_(&c__2, z__, &c__2, rhs, ipiv, jpiv, &scaloc);
403 if (scaloc != 1.f) {
404 i__2 = *n;
405 for (k = 1; k <= i__2; ++k) {
406 q__1.r = scaloc, q__1.i = 0.f;
407 cscal_(m, &q__1, &c__[k * c_dim1 + 1], &c__1);
408 q__1.r = scaloc, q__1.i = 0.f;
409 cscal_(m, &q__1, &f[k * f_dim1 + 1], &c__1);
410 /* L40: */
411 }
412 *scale *= scaloc;
413 }
414
415 /* Unpack solution vector(s) */
416
417 i__2 = i__ + j * c_dim1;
418 c__[i__2].r = rhs[0].r, c__[i__2].i = rhs[0].i;
419 i__2 = i__ + j * f_dim1;
420 f[i__2].r = rhs[1].r, f[i__2].i = rhs[1].i;
421
422 /* Substitute R(I, J) and L(I, J) into remaining equation. */
423
424 i__2 = j - 1;
425 for (k = 1; k <= i__2; ++k) {
426 i__3 = i__ + k * f_dim1;
427 i__4 = i__ + k * f_dim1;
428 r_cnjg(&q__4, &b[k + j * b_dim1]);
429 q__3.r = rhs[0].r * q__4.r - rhs[0].i * q__4.i, q__3.i =
430 rhs[0].r * q__4.i + rhs[0].i * q__4.r;
431 q__2.r = f[i__4].r + q__3.r, q__2.i = f[i__4].i + q__3.i;
432 r_cnjg(&q__6, &e[k + j * e_dim1]);
433 q__5.r = rhs[1].r * q__6.r - rhs[1].i * q__6.i, q__5.i =
434 rhs[1].r * q__6.i + rhs[1].i * q__6.r;
435 q__1.r = q__2.r + q__5.r, q__1.i = q__2.i + q__5.i;
436 f[i__3].r = q__1.r, f[i__3].i = q__1.i;
437 /* L50: */
438 }
439 i__2 = *m;
440 for (k = i__ + 1; k <= i__2; ++k) {
441 i__3 = k + j * c_dim1;
442 i__4 = k + j * c_dim1;
443 r_cnjg(&q__4, &a[i__ + k * a_dim1]);
444 q__3.r = q__4.r * rhs[0].r - q__4.i * rhs[0].i, q__3.i =
445 q__4.r * rhs[0].i + q__4.i * rhs[0].r;
446 q__2.r = c__[i__4].r - q__3.r, q__2.i = c__[i__4].i -
447 q__3.i;
448 r_cnjg(&q__6, &d__[i__ + k * d_dim1]);
449 q__5.r = q__6.r * rhs[1].r - q__6.i * rhs[1].i, q__5.i =
450 q__6.r * rhs[1].i + q__6.i * rhs[1].r;
451 q__1.r = q__2.r - q__5.r, q__1.i = q__2.i - q__5.i;
452 c__[i__3].r = q__1.r, c__[i__3].i = q__1.i;
453 /* L60: */
454 }
455
456 /* L70: */
457 }
458 /* L80: */
459 }
460 }
461 return 0;
462
463 /* End of CTGSY2 */
464
465 } /* ctgsy2_ */
466
467