1 /* ./src_f77/slaexc.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 integer c__4 = 4;
12 static logical c_false = FALSE_;
13 static integer c_n1 = -1;
14 static integer c__2 = 2;
15 static integer c__3 = 3;
16
slaexc_(logical * wantq,integer * n,real * t,integer * ldt,real * q,integer * ldq,integer * j1,integer * n1,integer * n2,real * work,integer * info)17 /* Subroutine */ int slaexc_(logical *wantq, integer *n, real *t, integer *
18 ldt, real *q, integer *ldq, integer *j1, integer *n1, integer *n2,
19 real *work, integer *info)
20 {
21 /* System generated locals */
22 integer q_dim1, q_offset, t_dim1, t_offset, i__1;
23 real r__1, r__2, r__3;
24
25 /* Local variables */
26 static real d__[16] /* was [4][4] */;
27 static integer k;
28 static real u[3], x[4] /* was [2][2] */;
29 static integer j2, j3, j4;
30 static real u1[3], u2[3];
31 static integer nd;
32 static real cs, t11, t22, t33, sn, wi1, wi2, wr1, wr2, eps, tau, tau1,
33 tau2;
34 static integer ierr;
35 static real temp;
36 extern /* Subroutine */ int srot_(integer *, real *, integer *, real *,
37 integer *, real *, real *);
38 static real scale, dnorm, xnorm;
39 extern /* Subroutine */ int slanv2_(real *, real *, real *, real *, real *
40 , real *, real *, real *, real *, real *), slasy2_(logical *,
41 logical *, integer *, integer *, integer *, real *, integer *,
42 real *, integer *, real *, integer *, real *, real *, integer *,
43 real *, integer *);
44 extern doublereal slamch_(char *, ftnlen), slange_(char *, integer *,
45 integer *, real *, integer *, real *, ftnlen);
46 extern /* Subroutine */ int slarfg_(integer *, real *, real *, integer *,
47 real *), slacpy_(char *, integer *, integer *, real *, integer *,
48 real *, integer *, ftnlen), slartg_(real *, real *, real *, real *
49 , real *);
50 static real thresh;
51 extern /* Subroutine */ int slarfx_(char *, integer *, integer *, real *,
52 real *, real *, integer *, real *, ftnlen);
53 static real smlnum;
54
55
56 /* -- LAPACK auxiliary routine (version 3.0) -- */
57 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
58 /* Courant Institute, Argonne National Lab, and Rice University */
59 /* February 29, 1992 */
60
61 /* .. Scalar Arguments .. */
62 /* .. */
63 /* .. Array Arguments .. */
64 /* .. */
65
66 /* Purpose */
67 /* ======= */
68
69 /* SLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in */
70 /* an upper quasi-triangular matrix T by an orthogonal similarity */
71 /* transformation. */
72
73 /* T must be in Schur canonical form, that is, block upper triangular */
74 /* with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block */
75 /* has its diagonal elemnts equal and its off-diagonal elements of */
76 /* opposite sign. */
77
78 /* Arguments */
79 /* ========= */
80
81 /* WANTQ (input) LOGICAL */
82 /* = .TRUE. : accumulate the transformation in the matrix Q; */
83 /* = .FALSE.: do not accumulate the transformation. */
84
85 /* N (input) INTEGER */
86 /* The order of the matrix T. N >= 0. */
87
88 /* T (input/output) REAL array, dimension (LDT,N) */
89 /* On entry, the upper quasi-triangular matrix T, in Schur */
90 /* canonical form. */
91 /* On exit, the updated matrix T, again in Schur canonical form. */
92
93 /* LDT (input) INTEGER */
94 /* The leading dimension of the array T. LDT >= max(1,N). */
95
96 /* Q (input/output) REAL array, dimension (LDQ,N) */
97 /* On entry, if WANTQ is .TRUE., the orthogonal matrix Q. */
98 /* On exit, if WANTQ is .TRUE., the updated matrix Q. */
99 /* If WANTQ is .FALSE., Q is not referenced. */
100
101 /* LDQ (input) INTEGER */
102 /* The leading dimension of the array Q. */
103 /* LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N. */
104
105 /* J1 (input) INTEGER */
106 /* The index of the first row of the first block T11. */
107
108 /* N1 (input) INTEGER */
109 /* The order of the first block T11. N1 = 0, 1 or 2. */
110
111 /* N2 (input) INTEGER */
112 /* The order of the second block T22. N2 = 0, 1 or 2. */
113
114 /* WORK (workspace) REAL array, dimension (N) */
115
116 /* INFO (output) INTEGER */
117 /* = 0: successful exit */
118 /* = 1: the transformed matrix T would be too far from Schur */
119 /* form; the blocks are not swapped and T and Q are */
120 /* unchanged. */
121
122 /* ===================================================================== */
123
124 /* .. Parameters .. */
125 /* .. */
126 /* .. Local Scalars .. */
127 /* .. */
128 /* .. Local Arrays .. */
129 /* .. */
130 /* .. External Functions .. */
131 /* .. */
132 /* .. External Subroutines .. */
133 /* .. */
134 /* .. Intrinsic Functions .. */
135 /* .. */
136 /* .. Executable Statements .. */
137
138 /* Parameter adjustments */
139 t_dim1 = *ldt;
140 t_offset = 1 + t_dim1;
141 t -= t_offset;
142 q_dim1 = *ldq;
143 q_offset = 1 + q_dim1;
144 q -= q_offset;
145 --work;
146
147 /* Function Body */
148 *info = 0;
149
150 /* Quick return if possible */
151
152 if (*n == 0 || *n1 == 0 || *n2 == 0) {
153 return 0;
154 }
155 if (*j1 + *n1 > *n) {
156 return 0;
157 }
158
159 j2 = *j1 + 1;
160 j3 = *j1 + 2;
161 j4 = *j1 + 3;
162
163 if (*n1 == 1 && *n2 == 1) {
164
165 /* Swap two 1-by-1 blocks. */
166
167 t11 = t[*j1 + *j1 * t_dim1];
168 t22 = t[j2 + j2 * t_dim1];
169
170 /* Determine the transformation to perform the interchange. */
171
172 r__1 = t22 - t11;
173 slartg_(&t[*j1 + j2 * t_dim1], &r__1, &cs, &sn, &temp);
174
175 /* Apply transformation to the matrix T. */
176
177 if (j3 <= *n) {
178 i__1 = *n - *j1 - 1;
179 srot_(&i__1, &t[*j1 + j3 * t_dim1], ldt, &t[j2 + j3 * t_dim1],
180 ldt, &cs, &sn);
181 }
182 i__1 = *j1 - 1;
183 srot_(&i__1, &t[*j1 * t_dim1 + 1], &c__1, &t[j2 * t_dim1 + 1], &c__1,
184 &cs, &sn);
185
186 t[*j1 + *j1 * t_dim1] = t22;
187 t[j2 + j2 * t_dim1] = t11;
188
189 if (*wantq) {
190
191 /* Accumulate transformation in the matrix Q. */
192
193 srot_(n, &q[*j1 * q_dim1 + 1], &c__1, &q[j2 * q_dim1 + 1], &c__1,
194 &cs, &sn);
195 }
196
197 } else {
198
199 /* Swapping involves at least one 2-by-2 block. */
200
201 /* Copy the diagonal block of order N1+N2 to the local array D */
202 /* and compute its norm. */
203
204 nd = *n1 + *n2;
205 slacpy_("Full", &nd, &nd, &t[*j1 + *j1 * t_dim1], ldt, d__, &c__4, (
206 ftnlen)4);
207 dnorm = slange_("Max", &nd, &nd, d__, &c__4, &work[1], (ftnlen)3);
208
209 /* Compute machine-dependent threshold for test for accepting */
210 /* swap. */
211
212 eps = slamch_("P", (ftnlen)1);
213 smlnum = slamch_("S", (ftnlen)1) / eps;
214 /* Computing MAX */
215 r__1 = eps * 10.f * dnorm;
216 thresh = dmax(r__1,smlnum);
217
218 /* Solve T11*X - X*T22 = scale*T12 for X. */
219
220 slasy2_(&c_false, &c_false, &c_n1, n1, n2, d__, &c__4, &d__[*n1 + 1 +
221 (*n1 + 1 << 2) - 5], &c__4, &d__[(*n1 + 1 << 2) - 4], &c__4, &
222 scale, x, &c__2, &xnorm, &ierr);
223
224 /* Swap the adjacent diagonal blocks. */
225
226 k = *n1 + *n1 + *n2 - 3;
227 switch (k) {
228 case 1: goto L10;
229 case 2: goto L20;
230 case 3: goto L30;
231 }
232
233 L10:
234
235 /* N1 = 1, N2 = 2: generate elementary reflector H so that: */
236
237 /* ( scale, X11, X12 ) H = ( 0, 0, * ) */
238
239 u[0] = scale;
240 u[1] = x[0];
241 u[2] = x[2];
242 slarfg_(&c__3, &u[2], u, &c__1, &tau);
243 u[2] = 1.f;
244 t11 = t[*j1 + *j1 * t_dim1];
245
246 /* Perform swap provisionally on diagonal block in D. */
247
248 slarfx_("L", &c__3, &c__3, u, &tau, d__, &c__4, &work[1], (ftnlen)1);
249 slarfx_("R", &c__3, &c__3, u, &tau, d__, &c__4, &work[1], (ftnlen)1);
250
251 /* Test whether to reject swap. */
252
253 /* Computing MAX */
254 r__2 = dabs(d__[2]), r__3 = dabs(d__[6]), r__2 = max(r__2,r__3), r__3
255 = (r__1 = d__[10] - t11, dabs(r__1));
256 if (dmax(r__2,r__3) > thresh) {
257 goto L50;
258 }
259
260 /* Accept swap: apply transformation to the entire matrix T. */
261
262 i__1 = *n - *j1 + 1;
263 slarfx_("L", &c__3, &i__1, u, &tau, &t[*j1 + *j1 * t_dim1], ldt, &
264 work[1], (ftnlen)1);
265 slarfx_("R", &j2, &c__3, u, &tau, &t[*j1 * t_dim1 + 1], ldt, &work[1],
266 (ftnlen)1);
267
268 t[j3 + *j1 * t_dim1] = 0.f;
269 t[j3 + j2 * t_dim1] = 0.f;
270 t[j3 + j3 * t_dim1] = t11;
271
272 if (*wantq) {
273
274 /* Accumulate transformation in the matrix Q. */
275
276 slarfx_("R", n, &c__3, u, &tau, &q[*j1 * q_dim1 + 1], ldq, &work[
277 1], (ftnlen)1);
278 }
279 goto L40;
280
281 L20:
282
283 /* N1 = 2, N2 = 1: generate elementary reflector H so that: */
284
285 /* H ( -X11 ) = ( * ) */
286 /* ( -X21 ) = ( 0 ) */
287 /* ( scale ) = ( 0 ) */
288
289 u[0] = -x[0];
290 u[1] = -x[1];
291 u[2] = scale;
292 slarfg_(&c__3, u, &u[1], &c__1, &tau);
293 u[0] = 1.f;
294 t33 = t[j3 + j3 * t_dim1];
295
296 /* Perform swap provisionally on diagonal block in D. */
297
298 slarfx_("L", &c__3, &c__3, u, &tau, d__, &c__4, &work[1], (ftnlen)1);
299 slarfx_("R", &c__3, &c__3, u, &tau, d__, &c__4, &work[1], (ftnlen)1);
300
301 /* Test whether to reject swap. */
302
303 /* Computing MAX */
304 r__2 = dabs(d__[1]), r__3 = dabs(d__[2]), r__2 = max(r__2,r__3), r__3
305 = (r__1 = d__[0] - t33, dabs(r__1));
306 if (dmax(r__2,r__3) > thresh) {
307 goto L50;
308 }
309
310 /* Accept swap: apply transformation to the entire matrix T. */
311
312 slarfx_("R", &j3, &c__3, u, &tau, &t[*j1 * t_dim1 + 1], ldt, &work[1],
313 (ftnlen)1);
314 i__1 = *n - *j1;
315 slarfx_("L", &c__3, &i__1, u, &tau, &t[*j1 + j2 * t_dim1], ldt, &work[
316 1], (ftnlen)1);
317
318 t[*j1 + *j1 * t_dim1] = t33;
319 t[j2 + *j1 * t_dim1] = 0.f;
320 t[j3 + *j1 * t_dim1] = 0.f;
321
322 if (*wantq) {
323
324 /* Accumulate transformation in the matrix Q. */
325
326 slarfx_("R", n, &c__3, u, &tau, &q[*j1 * q_dim1 + 1], ldq, &work[
327 1], (ftnlen)1);
328 }
329 goto L40;
330
331 L30:
332
333 /* N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so */
334 /* that: */
335
336 /* H(2) H(1) ( -X11 -X12 ) = ( * * ) */
337 /* ( -X21 -X22 ) ( 0 * ) */
338 /* ( scale 0 ) ( 0 0 ) */
339 /* ( 0 scale ) ( 0 0 ) */
340
341 u1[0] = -x[0];
342 u1[1] = -x[1];
343 u1[2] = scale;
344 slarfg_(&c__3, u1, &u1[1], &c__1, &tau1);
345 u1[0] = 1.f;
346
347 temp = -tau1 * (x[2] + u1[1] * x[3]);
348 u2[0] = -temp * u1[1] - x[3];
349 u2[1] = -temp * u1[2];
350 u2[2] = scale;
351 slarfg_(&c__3, u2, &u2[1], &c__1, &tau2);
352 u2[0] = 1.f;
353
354 /* Perform swap provisionally on diagonal block in D. */
355
356 slarfx_("L", &c__3, &c__4, u1, &tau1, d__, &c__4, &work[1], (ftnlen)1)
357 ;
358 slarfx_("R", &c__4, &c__3, u1, &tau1, d__, &c__4, &work[1], (ftnlen)1)
359 ;
360 slarfx_("L", &c__3, &c__4, u2, &tau2, &d__[1], &c__4, &work[1], (
361 ftnlen)1);
362 slarfx_("R", &c__4, &c__3, u2, &tau2, &d__[4], &c__4, &work[1], (
363 ftnlen)1);
364
365 /* Test whether to reject swap. */
366
367 /* Computing MAX */
368 r__1 = dabs(d__[2]), r__2 = dabs(d__[6]), r__1 = max(r__1,r__2), r__2
369 = dabs(d__[3]), r__1 = max(r__1,r__2), r__2 = dabs(d__[7]);
370 if (dmax(r__1,r__2) > thresh) {
371 goto L50;
372 }
373
374 /* Accept swap: apply transformation to the entire matrix T. */
375
376 i__1 = *n - *j1 + 1;
377 slarfx_("L", &c__3, &i__1, u1, &tau1, &t[*j1 + *j1 * t_dim1], ldt, &
378 work[1], (ftnlen)1);
379 slarfx_("R", &j4, &c__3, u1, &tau1, &t[*j1 * t_dim1 + 1], ldt, &work[
380 1], (ftnlen)1);
381 i__1 = *n - *j1 + 1;
382 slarfx_("L", &c__3, &i__1, u2, &tau2, &t[j2 + *j1 * t_dim1], ldt, &
383 work[1], (ftnlen)1);
384 slarfx_("R", &j4, &c__3, u2, &tau2, &t[j2 * t_dim1 + 1], ldt, &work[1]
385 , (ftnlen)1);
386
387 t[j3 + *j1 * t_dim1] = 0.f;
388 t[j3 + j2 * t_dim1] = 0.f;
389 t[j4 + *j1 * t_dim1] = 0.f;
390 t[j4 + j2 * t_dim1] = 0.f;
391
392 if (*wantq) {
393
394 /* Accumulate transformation in the matrix Q. */
395
396 slarfx_("R", n, &c__3, u1, &tau1, &q[*j1 * q_dim1 + 1], ldq, &
397 work[1], (ftnlen)1);
398 slarfx_("R", n, &c__3, u2, &tau2, &q[j2 * q_dim1 + 1], ldq, &work[
399 1], (ftnlen)1);
400 }
401
402 L40:
403
404 if (*n2 == 2) {
405
406 /* Standardize new 2-by-2 block T11 */
407
408 slanv2_(&t[*j1 + *j1 * t_dim1], &t[*j1 + j2 * t_dim1], &t[j2 + *
409 j1 * t_dim1], &t[j2 + j2 * t_dim1], &wr1, &wi1, &wr2, &
410 wi2, &cs, &sn);
411 i__1 = *n - *j1 - 1;
412 srot_(&i__1, &t[*j1 + (*j1 + 2) * t_dim1], ldt, &t[j2 + (*j1 + 2)
413 * t_dim1], ldt, &cs, &sn);
414 i__1 = *j1 - 1;
415 srot_(&i__1, &t[*j1 * t_dim1 + 1], &c__1, &t[j2 * t_dim1 + 1], &
416 c__1, &cs, &sn);
417 if (*wantq) {
418 srot_(n, &q[*j1 * q_dim1 + 1], &c__1, &q[j2 * q_dim1 + 1], &
419 c__1, &cs, &sn);
420 }
421 }
422
423 if (*n1 == 2) {
424
425 /* Standardize new 2-by-2 block T22 */
426
427 j3 = *j1 + *n2;
428 j4 = j3 + 1;
429 slanv2_(&t[j3 + j3 * t_dim1], &t[j3 + j4 * t_dim1], &t[j4 + j3 *
430 t_dim1], &t[j4 + j4 * t_dim1], &wr1, &wi1, &wr2, &wi2, &
431 cs, &sn);
432 if (j3 + 2 <= *n) {
433 i__1 = *n - j3 - 1;
434 srot_(&i__1, &t[j3 + (j3 + 2) * t_dim1], ldt, &t[j4 + (j3 + 2)
435 * t_dim1], ldt, &cs, &sn);
436 }
437 i__1 = j3 - 1;
438 srot_(&i__1, &t[j3 * t_dim1 + 1], &c__1, &t[j4 * t_dim1 + 1], &
439 c__1, &cs, &sn);
440 if (*wantq) {
441 srot_(n, &q[j3 * q_dim1 + 1], &c__1, &q[j4 * q_dim1 + 1], &
442 c__1, &cs, &sn);
443 }
444 }
445
446 }
447 return 0;
448
449 /* Exit with INFO = 1 if swap was rejected. */
450
451 L50:
452 *info = 1;
453 return 0;
454
455 /* End of SLAEXC */
456
457 } /* slaexc_ */
458
459