1 /* ./src_f77/slagv2.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
slagv2_(real * a,integer * lda,real * b,integer * ldb,real * alphar,real * alphai,real * beta,real * csl,real * snl,real * csr,real * snr)13 /* Subroutine */ int slagv2_(real *a, integer *lda, real *b, integer *ldb,
14 real *alphar, real *alphai, real *beta, real *csl, real *snl, real *
15 csr, real *snr)
16 {
17 /* System generated locals */
18 integer a_dim1, a_offset, b_dim1, b_offset;
19 real r__1, r__2, r__3, r__4, r__5, r__6;
20
21 /* Local variables */
22 static real r__, t, h1, h2, h3, wi, qq, rr, wr1, wr2, ulp;
23 extern /* Subroutine */ int srot_(integer *, real *, integer *, real *,
24 integer *, real *, real *), slag2_(real *, integer *, real *,
25 integer *, real *, real *, real *, real *, real *, real *);
26 static real anorm, bnorm, scale1, scale2;
27 extern /* Subroutine */ int slasv2_(real *, real *, real *, real *, real *
28 , real *, real *, real *, real *);
29 extern doublereal slapy2_(real *, real *);
30 static real ascale, bscale;
31 extern doublereal slamch_(char *, ftnlen);
32 static real safmin;
33 extern /* Subroutine */ int slartg_(real *, real *, real *, real *, real *
34 );
35
36
37 /* -- LAPACK auxiliary routine (version 3.0) -- */
38 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
39 /* Courant Institute, Argonne National Lab, and Rice University */
40 /* June 30, 1999 */
41
42 /* .. Scalar Arguments .. */
43 /* .. */
44 /* .. Array Arguments .. */
45 /* .. */
46
47 /* Purpose */
48 /* ======= */
49
50 /* SLAGV2 computes the Generalized Schur factorization of a real 2-by-2 */
51 /* matrix pencil (A,B) where B is upper triangular. This routine */
52 /* computes orthogonal (rotation) matrices given by CSL, SNL and CSR, */
53 /* SNR such that */
54
55 /* 1) if the pencil (A,B) has two real eigenvalues (include 0/0 or 1/0 */
56 /* types), then */
57
58 /* [ a11 a12 ] := [ CSL SNL ] [ a11 a12 ] [ CSR -SNR ] */
59 /* [ 0 a22 ] [ -SNL CSL ] [ a21 a22 ] [ SNR CSR ] */
60
61 /* [ b11 b12 ] := [ CSL SNL ] [ b11 b12 ] [ CSR -SNR ] */
62 /* [ 0 b22 ] [ -SNL CSL ] [ 0 b22 ] [ SNR CSR ], */
63
64 /* 2) if the pencil (A,B) has a pair of complex conjugate eigenvalues, */
65 /* then */
66
67 /* [ a11 a12 ] := [ CSL SNL ] [ a11 a12 ] [ CSR -SNR ] */
68 /* [ a21 a22 ] [ -SNL CSL ] [ a21 a22 ] [ SNR CSR ] */
69
70 /* [ b11 0 ] := [ CSL SNL ] [ b11 b12 ] [ CSR -SNR ] */
71 /* [ 0 b22 ] [ -SNL CSL ] [ 0 b22 ] [ SNR CSR ] */
72
73 /* where b11 >= b22 > 0. */
74
75
76 /* Arguments */
77 /* ========= */
78
79 /* A (input/output) REAL array, dimension (LDA, 2) */
80 /* On entry, the 2 x 2 matrix A. */
81 /* On exit, A is overwritten by the ``A-part'' of the */
82 /* generalized Schur form. */
83
84 /* LDA (input) INTEGER */
85 /* THe leading dimension of the array A. LDA >= 2. */
86
87 /* B (input/output) REAL array, dimension (LDB, 2) */
88 /* On entry, the upper triangular 2 x 2 matrix B. */
89 /* On exit, B is overwritten by the ``B-part'' of the */
90 /* generalized Schur form. */
91
92 /* LDB (input) INTEGER */
93 /* THe leading dimension of the array B. LDB >= 2. */
94
95 /* ALPHAR (output) REAL array, dimension (2) */
96 /* ALPHAI (output) REAL array, dimension (2) */
97 /* BETA (output) REAL array, dimension (2) */
98 /* (ALPHAR(k)+i*ALPHAI(k))/BETA(k) are the eigenvalues of the */
99 /* pencil (A,B), k=1,2, i = sqrt(-1). Note that BETA(k) may */
100 /* be zero. */
101
102 /* CSL (output) REAL */
103 /* The cosine of the left rotation matrix. */
104
105 /* SNL (output) REAL */
106 /* The sine of the left rotation matrix. */
107
108 /* CSR (output) REAL */
109 /* The cosine of the right rotation matrix. */
110
111 /* SNR (output) REAL */
112 /* The sine of the right rotation matrix. */
113
114 /* Further Details */
115 /* =============== */
116
117 /* Based on contributions by */
118 /* Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA */
119
120 /* ===================================================================== */
121
122 /* .. Parameters .. */
123 /* .. */
124 /* .. Local Scalars .. */
125 /* .. */
126 /* .. External Subroutines .. */
127 /* .. */
128 /* .. External Functions .. */
129 /* .. */
130 /* .. Intrinsic Functions .. */
131 /* .. */
132 /* .. Executable Statements .. */
133
134 /* Parameter adjustments */
135 a_dim1 = *lda;
136 a_offset = 1 + a_dim1;
137 a -= a_offset;
138 b_dim1 = *ldb;
139 b_offset = 1 + b_dim1;
140 b -= b_offset;
141 --alphar;
142 --alphai;
143 --beta;
144
145 /* Function Body */
146 safmin = slamch_("S", (ftnlen)1);
147 ulp = slamch_("P", (ftnlen)1);
148
149 /* Scale A */
150
151 /* Computing MAX */
152 r__5 = (r__1 = a[a_dim1 + 1], dabs(r__1)) + (r__2 = a[a_dim1 + 2], dabs(
153 r__2)), r__6 = (r__3 = a[(a_dim1 << 1) + 1], dabs(r__3)) + (r__4 =
154 a[(a_dim1 << 1) + 2], dabs(r__4)), r__5 = max(r__5,r__6);
155 anorm = dmax(r__5,safmin);
156 ascale = 1.f / anorm;
157 a[a_dim1 + 1] = ascale * a[a_dim1 + 1];
158 a[(a_dim1 << 1) + 1] = ascale * a[(a_dim1 << 1) + 1];
159 a[a_dim1 + 2] = ascale * a[a_dim1 + 2];
160 a[(a_dim1 << 1) + 2] = ascale * a[(a_dim1 << 1) + 2];
161
162 /* Scale B */
163
164 /* Computing MAX */
165 r__4 = (r__3 = b[b_dim1 + 1], dabs(r__3)), r__5 = (r__1 = b[(b_dim1 << 1)
166 + 1], dabs(r__1)) + (r__2 = b[(b_dim1 << 1) + 2], dabs(r__2)),
167 r__4 = max(r__4,r__5);
168 bnorm = dmax(r__4,safmin);
169 bscale = 1.f / bnorm;
170 b[b_dim1 + 1] = bscale * b[b_dim1 + 1];
171 b[(b_dim1 << 1) + 1] = bscale * b[(b_dim1 << 1) + 1];
172 b[(b_dim1 << 1) + 2] = bscale * b[(b_dim1 << 1) + 2];
173
174 /* Check if A can be deflated */
175
176 if ((r__1 = a[a_dim1 + 2], dabs(r__1)) <= ulp) {
177 *csl = 1.f;
178 *snl = 0.f;
179 *csr = 1.f;
180 *snr = 0.f;
181 a[a_dim1 + 2] = 0.f;
182 b[b_dim1 + 2] = 0.f;
183
184 /* Check if B is singular */
185
186 } else if ((r__1 = b[b_dim1 + 1], dabs(r__1)) <= ulp) {
187 slartg_(&a[a_dim1 + 1], &a[a_dim1 + 2], csl, snl, &r__);
188 *csr = 1.f;
189 *snr = 0.f;
190 srot_(&c__2, &a[a_dim1 + 1], lda, &a[a_dim1 + 2], lda, csl, snl);
191 srot_(&c__2, &b[b_dim1 + 1], ldb, &b[b_dim1 + 2], ldb, csl, snl);
192 a[a_dim1 + 2] = 0.f;
193 b[b_dim1 + 1] = 0.f;
194 b[b_dim1 + 2] = 0.f;
195
196 } else if ((r__1 = b[(b_dim1 << 1) + 2], dabs(r__1)) <= ulp) {
197 slartg_(&a[(a_dim1 << 1) + 2], &a[a_dim1 + 2], csr, snr, &t);
198 *snr = -(*snr);
199 srot_(&c__2, &a[a_dim1 + 1], &c__1, &a[(a_dim1 << 1) + 1], &c__1, csr,
200 snr);
201 srot_(&c__2, &b[b_dim1 + 1], &c__1, &b[(b_dim1 << 1) + 1], &c__1, csr,
202 snr);
203 *csl = 1.f;
204 *snl = 0.f;
205 a[a_dim1 + 2] = 0.f;
206 b[b_dim1 + 2] = 0.f;
207 b[(b_dim1 << 1) + 2] = 0.f;
208
209 } else {
210
211 /* B is nonsingular, first compute the eigenvalues of (A,B) */
212
213 slag2_(&a[a_offset], lda, &b[b_offset], ldb, &safmin, &scale1, &
214 scale2, &wr1, &wr2, &wi);
215
216 if (wi == 0.f) {
217
218 /* two real eigenvalues, compute s*A-w*B */
219
220 h1 = scale1 * a[a_dim1 + 1] - wr1 * b[b_dim1 + 1];
221 h2 = scale1 * a[(a_dim1 << 1) + 1] - wr1 * b[(b_dim1 << 1) + 1];
222 h3 = scale1 * a[(a_dim1 << 1) + 2] - wr1 * b[(b_dim1 << 1) + 2];
223
224 rr = slapy2_(&h1, &h2);
225 r__1 = scale1 * a[a_dim1 + 2];
226 qq = slapy2_(&r__1, &h3);
227
228 if (rr > qq) {
229
230 /* find right rotation matrix to zero 1,1 element of */
231 /* (sA - wB) */
232
233 slartg_(&h2, &h1, csr, snr, &t);
234
235 } else {
236
237 /* find right rotation matrix to zero 2,1 element of */
238 /* (sA - wB) */
239
240 r__1 = scale1 * a[a_dim1 + 2];
241 slartg_(&h3, &r__1, csr, snr, &t);
242
243 }
244
245 *snr = -(*snr);
246 srot_(&c__2, &a[a_dim1 + 1], &c__1, &a[(a_dim1 << 1) + 1], &c__1,
247 csr, snr);
248 srot_(&c__2, &b[b_dim1 + 1], &c__1, &b[(b_dim1 << 1) + 1], &c__1,
249 csr, snr);
250
251 /* compute inf norms of A and B */
252
253 /* Computing MAX */
254 r__5 = (r__1 = a[a_dim1 + 1], dabs(r__1)) + (r__2 = a[(a_dim1 <<
255 1) + 1], dabs(r__2)), r__6 = (r__3 = a[a_dim1 + 2], dabs(
256 r__3)) + (r__4 = a[(a_dim1 << 1) + 2], dabs(r__4));
257 h1 = dmax(r__5,r__6);
258 /* Computing MAX */
259 r__5 = (r__1 = b[b_dim1 + 1], dabs(r__1)) + (r__2 = b[(b_dim1 <<
260 1) + 1], dabs(r__2)), r__6 = (r__3 = b[b_dim1 + 2], dabs(
261 r__3)) + (r__4 = b[(b_dim1 << 1) + 2], dabs(r__4));
262 h2 = dmax(r__5,r__6);
263
264 if (scale1 * h1 >= dabs(wr1) * h2) {
265
266 /* find left rotation matrix Q to zero out B(2,1) */
267
268 slartg_(&b[b_dim1 + 1], &b[b_dim1 + 2], csl, snl, &r__);
269
270 } else {
271
272 /* find left rotation matrix Q to zero out A(2,1) */
273
274 slartg_(&a[a_dim1 + 1], &a[a_dim1 + 2], csl, snl, &r__);
275
276 }
277
278 srot_(&c__2, &a[a_dim1 + 1], lda, &a[a_dim1 + 2], lda, csl, snl);
279 srot_(&c__2, &b[b_dim1 + 1], ldb, &b[b_dim1 + 2], ldb, csl, snl);
280
281 a[a_dim1 + 2] = 0.f;
282 b[b_dim1 + 2] = 0.f;
283
284 } else {
285
286 /* a pair of complex conjugate eigenvalues */
287 /* first compute the SVD of the matrix B */
288
289 slasv2_(&b[b_dim1 + 1], &b[(b_dim1 << 1) + 1], &b[(b_dim1 << 1) +
290 2], &r__, &t, snr, csr, snl, csl);
291
292 /* Form (A,B) := Q(A,B)Z' where Q is left rotation matrix and */
293 /* Z is right rotation matrix computed from SLASV2 */
294
295 srot_(&c__2, &a[a_dim1 + 1], lda, &a[a_dim1 + 2], lda, csl, snl);
296 srot_(&c__2, &b[b_dim1 + 1], ldb, &b[b_dim1 + 2], ldb, csl, snl);
297 srot_(&c__2, &a[a_dim1 + 1], &c__1, &a[(a_dim1 << 1) + 1], &c__1,
298 csr, snr);
299 srot_(&c__2, &b[b_dim1 + 1], &c__1, &b[(b_dim1 << 1) + 1], &c__1,
300 csr, snr);
301
302 b[b_dim1 + 2] = 0.f;
303 b[(b_dim1 << 1) + 1] = 0.f;
304
305 }
306
307 }
308
309 /* Unscaling */
310
311 a[a_dim1 + 1] = anorm * a[a_dim1 + 1];
312 a[a_dim1 + 2] = anorm * a[a_dim1 + 2];
313 a[(a_dim1 << 1) + 1] = anorm * a[(a_dim1 << 1) + 1];
314 a[(a_dim1 << 1) + 2] = anorm * a[(a_dim1 << 1) + 2];
315 b[b_dim1 + 1] = bnorm * b[b_dim1 + 1];
316 b[b_dim1 + 2] = bnorm * b[b_dim1 + 2];
317 b[(b_dim1 << 1) + 1] = bnorm * b[(b_dim1 << 1) + 1];
318 b[(b_dim1 << 1) + 2] = bnorm * b[(b_dim1 << 1) + 2];
319
320 if (wi == 0.f) {
321 alphar[1] = a[a_dim1 + 1];
322 alphar[2] = a[(a_dim1 << 1) + 2];
323 alphai[1] = 0.f;
324 alphai[2] = 0.f;
325 beta[1] = b[b_dim1 + 1];
326 beta[2] = b[(b_dim1 << 1) + 2];
327 } else {
328 alphar[1] = anorm * wr1 / scale1 / bnorm;
329 alphai[1] = anorm * wi / scale1 / bnorm;
330 alphar[2] = alphar[1];
331 alphai[2] = -alphai[1];
332 beta[1] = 1.f;
333 beta[2] = 1.f;
334 }
335
336 /* L10: */
337
338 return 0;
339
340 /* End of SLAGV2 */
341
342 } /* slagv2_ */
343
344