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