1 /* ../netlib/dtgexc.f -- translated by f2c (version 20100827). You must link the resulting object file with libf2c: on Microsoft Windows system, link with libf2c.lib;
2 on Linux or Unix systems, link with .../path/to/libf2c.a -lm or, if you install libf2c.a in a standard place, with -lf2c -lm -- in that order, at the end of the command line, as in cc *.o -lf2c -lm Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., http://www.netlib.org/f2c/libf2c.zip */
3 #include "FLA_f2c.h" /* Table of constant values */
4 static integer c__1 = 1;
5 static integer c__2 = 2;
6 /* > \brief \b DTGEXC */
7 /* =========== DOCUMENTATION =========== */
8 /* Online html documentation available at */
9 /* http://www.netlib.org/lapack/explore-html/ */
10 /* > \htmlonly */
11 /* > Download DTGEXC + dependencies */
12 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgexc. f"> */
13 /* > [TGZ]</a> */
14 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgexc. f"> */
15 /* > [ZIP]</a> */
16 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgexc. f"> */
17 /* > [TXT]</a> */
18 /* > \endhtmlonly */
19 /* Definition: */
20 /* =========== */
21 /* SUBROUTINE DTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, */
22 /* LDZ, IFST, ILST, WORK, LWORK, INFO ) */
23 /* .. Scalar Arguments .. */
24 /* LOGICAL WANTQ, WANTZ */
25 /* INTEGER IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N */
26 /* .. */
27 /* .. Array Arguments .. */
28 /* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ), */
29 /* $ WORK( * ), Z( LDZ, * ) */
30 /* .. */
31 /* > \par Purpose: */
32 /* ============= */
33 /* > */
34 /* > \verbatim */
35 /* > */
36 /* > DTGEXC reorders the generalized real Schur decomposition of a real */
37 /* > matrix pair (A,B) using an orthogonal equivalence transformation */
38 /* > */
39 /* > (A, B) = Q * (A, B) * Z**T, */
40 /* > */
41 /* > so that the diagonal block of (A, B) with row index IFST is moved */
42 /* > to row ILST. */
43 /* > */
44 /* > (A, B) must be in generalized real Schur canonical form (as returned */
45 /* > by DGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2 */
46 /* > diagonal blocks. B is upper triangular. */
47 /* > */
48 /* > Optionally, the matrices Q and Z of generalized Schur vectors are */
49 /* > updated. */
50 /* > */
51 /* > Q(in) * A(in) * Z(in)**T = Q(out) * A(out) * Z(out)**T */
52 /* > Q(in) * B(in) * Z(in)**T = Q(out) * B(out) * Z(out)**T */
53 /* > */
54 /* > \endverbatim */
55 /* Arguments: */
56 /* ========== */
57 /* > \param[in] WANTQ */
58 /* > \verbatim */
59 /* > WANTQ is LOGICAL */
60 /* > .TRUE. : update the left transformation matrix Q;
61 */
62 /* > .FALSE.: do not update Q. */
63 /* > \endverbatim */
64 /* > */
65 /* > \param[in] WANTZ */
66 /* > \verbatim */
67 /* > WANTZ is LOGICAL */
68 /* > .TRUE. : update the right transformation matrix Z;
69 */
70 /* > .FALSE.: do not update Z. */
71 /* > \endverbatim */
72 /* > */
73 /* > \param[in] N */
74 /* > \verbatim */
75 /* > N is INTEGER */
76 /* > The order of the matrices A and B. N >= 0. */
77 /* > \endverbatim */
78 /* > */
79 /* > \param[in,out] A */
80 /* > \verbatim */
81 /* > A is DOUBLE PRECISION array, dimension (LDA,N) */
82 /* > On entry, the matrix A in generalized real Schur canonical */
83 /* > form. */
84 /* > On exit, the updated matrix A, again in generalized */
85 /* > real Schur canonical form. */
86 /* > \endverbatim */
87 /* > */
88 /* > \param[in] LDA */
89 /* > \verbatim */
90 /* > LDA is INTEGER */
91 /* > The leading dimension of the array A. LDA >= max(1,N). */
92 /* > \endverbatim */
93 /* > */
94 /* > \param[in,out] B */
95 /* > \verbatim */
96 /* > B is DOUBLE PRECISION array, dimension (LDB,N) */
97 /* > On entry, the matrix B in generalized real Schur canonical */
98 /* > form (A,B). */
99 /* > On exit, the updated matrix B, again in generalized */
100 /* > real Schur canonical form (A,B). */
101 /* > \endverbatim */
102 /* > */
103 /* > \param[in] LDB */
104 /* > \verbatim */
105 /* > LDB is INTEGER */
106 /* > The leading dimension of the array B. LDB >= max(1,N). */
107 /* > \endverbatim */
108 /* > */
109 /* > \param[in,out] Q */
110 /* > \verbatim */
111 /* > Q is DOUBLE PRECISION array, dimension (LDQ,N) */
112 /* > On entry, if WANTQ = .TRUE., the orthogonal matrix Q. */
113 /* > On exit, the updated matrix Q. */
114 /* > If WANTQ = .FALSE., Q is not referenced. */
115 /* > \endverbatim */
116 /* > */
117 /* > \param[in] LDQ */
118 /* > \verbatim */
119 /* > LDQ is INTEGER */
120 /* > The leading dimension of the array Q. LDQ >= 1. */
121 /* > If WANTQ = .TRUE., LDQ >= N. */
122 /* > \endverbatim */
123 /* > */
124 /* > \param[in,out] Z */
125 /* > \verbatim */
126 /* > Z is DOUBLE PRECISION array, dimension (LDZ,N) */
127 /* > On entry, if WANTZ = .TRUE., the orthogonal matrix Z. */
128 /* > On exit, the updated matrix Z. */
129 /* > If WANTZ = .FALSE., Z is not referenced. */
130 /* > \endverbatim */
131 /* > */
132 /* > \param[in] LDZ */
133 /* > \verbatim */
134 /* > LDZ is INTEGER */
135 /* > The leading dimension of the array Z. LDZ >= 1. */
136 /* > If WANTZ = .TRUE., LDZ >= N. */
137 /* > \endverbatim */
138 /* > */
139 /* > \param[in,out] IFST */
140 /* > \verbatim */
141 /* > IFST is INTEGER */
142 /* > \endverbatim */
143 /* > */
144 /* > \param[in,out] ILST */
145 /* > \verbatim */
146 /* > ILST is INTEGER */
147 /* > Specify the reordering of the diagonal blocks of (A, B). */
148 /* > The block with row index IFST is moved to row ILST, by a */
149 /* > sequence of swapping between adjacent blocks. */
150 /* > On exit, if IFST pointed on entry to the second row of */
151 /* > a 2-by-2 block, it is changed to point to the first row;
152 */
153 /* > ILST always points to the first row of the block in its */
154 /* > final position (which may differ from its input value by */
155 /* > +1 or -1). 1 <= IFST, ILST <= N. */
156 /* > \endverbatim */
157 /* > */
158 /* > \param[out] WORK */
159 /* > \verbatim */
160 /* > WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */
161 /* > On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
162 /* > \endverbatim */
163 /* > */
164 /* > \param[in] LWORK */
165 /* > \verbatim */
166 /* > LWORK is INTEGER */
167 /* > The dimension of the array WORK. */
168 /* > LWORK >= 1 when N <= 1, otherwise LWORK >= 4*N + 16. */
169 /* > */
170 /* > If LWORK = -1, then a workspace query is assumed;
171 the routine */
172 /* > only calculates the optimal size of the WORK array, returns */
173 /* > this value as the first entry of the WORK array, and no error */
174 /* > message related to LWORK is issued by XERBLA. */
175 /* > \endverbatim */
176 /* > */
177 /* > \param[out] INFO */
178 /* > \verbatim */
179 /* > INFO is INTEGER */
180 /* > =0: successful exit. */
181 /* > <0: if INFO = -i, the i-th argument had an illegal value. */
182 /* > =1: The transformed matrix pair (A, B) would be too far */
183 /* > from generalized Schur form;
184 the problem is ill- */
185 /* > conditioned. (A, B) may have been partially reordered, */
186 /* > and ILST points to the first row of the current */
187 /* > position of the block being moved. */
188 /* > \endverbatim */
189 /* Authors: */
190 /* ======== */
191 /* > \author Univ. of Tennessee */
192 /* > \author Univ. of California Berkeley */
193 /* > \author Univ. of Colorado Denver */
194 /* > \author NAG Ltd. */
195 /* > \date November 2011 */
196 /* > \ingroup doubleGEcomputational */
197 /* > \par Contributors: */
198 /* ================== */
199 /* > */
200 /* > Bo Kagstrom and Peter Poromaa, Department of Computing Science, */
201 /* > Umea University, S-901 87 Umea, Sweden. */
202 /* > \par References: */
203 /* ================ */
204 /* > */
205 /* > \verbatim */
206 /* > */
207 /* > [1] B. Kagstrom;
208 A Direct Method for Reordering Eigenvalues in the */
209 /* > Generalized Real Schur Form of a Regular Matrix Pair (A, B), in */
210 /* > M.S. Moonen et al (eds), Linear Algebra for Large Scale and */
211 /* > Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218. */
212 /* > \endverbatim */
213 /* > */
214 /* ===================================================================== */
215 /* Subroutine */
dtgexc_(logical * wantq,logical * wantz,integer * n,doublereal * a,integer * lda,doublereal * b,integer * ldb,doublereal * q,integer * ldq,doublereal * z__,integer * ldz,integer * ifst,integer * ilst,doublereal * work,integer * lwork,integer * info)216 int dtgexc_(logical *wantq, logical *wantz, integer *n, doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal * q, integer *ldq, doublereal *z__, integer *ldz, integer *ifst, integer *ilst, doublereal *work, integer *lwork, integer *info)
217 {
218 /* System generated locals */
219 integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, z_dim1, z_offset, i__1;
220 /* Local variables */
221 integer nbf, nbl, here, lwmin;
222 extern /* Subroutine */
223 int dtgex2_(logical *, logical *, integer *, doublereal *, integer *, doublereal *, integer *, doublereal *, integer *, doublereal *, integer *, integer *, integer *, integer *, doublereal *, integer *, integer *), xerbla_(char *, integer *);
224 integer nbnext;
225 logical lquery;
226 /* -- LAPACK computational routine (version 3.4.0) -- */
227 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
228 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
229 /* November 2011 */
230 /* .. Scalar Arguments .. */
231 /* .. */
232 /* .. Array Arguments .. */
233 /* .. */
234 /* ===================================================================== */
235 /* .. Parameters .. */
236 /* .. */
237 /* .. Local Scalars .. */
238 /* .. */
239 /* .. External Subroutines .. */
240 /* .. */
241 /* .. Intrinsic Functions .. */
242 /* .. */
243 /* .. Executable Statements .. */
244 /* Decode and test input arguments. */
245 /* Parameter adjustments */
246 a_dim1 = *lda;
247 a_offset = 1 + a_dim1;
248 a -= a_offset;
249 b_dim1 = *ldb;
250 b_offset = 1 + b_dim1;
251 b -= b_offset;
252 q_dim1 = *ldq;
253 q_offset = 1 + q_dim1;
254 q -= q_offset;
255 z_dim1 = *ldz;
256 z_offset = 1 + z_dim1;
257 z__ -= z_offset;
258 --work;
259 /* Function Body */
260 *info = 0;
261 lquery = *lwork == -1;
262 if (*n < 0)
263 {
264 *info = -3;
265 }
266 else if (*lda < max(1,*n))
267 {
268 *info = -5;
269 }
270 else if (*ldb < max(1,*n))
271 {
272 *info = -7;
273 }
274 else if (*ldq < 1 || *wantq && *ldq < max(1,*n))
275 {
276 *info = -9;
277 }
278 else if (*ldz < 1 || *wantz && *ldz < max(1,*n))
279 {
280 *info = -11;
281 }
282 else if (*ifst < 1 || *ifst > *n)
283 {
284 *info = -12;
285 }
286 else if (*ilst < 1 || *ilst > *n)
287 {
288 *info = -13;
289 }
290 if (*info == 0)
291 {
292 if (*n <= 1)
293 {
294 lwmin = 1;
295 }
296 else
297 {
298 lwmin = (*n << 2) + 16;
299 }
300 work[1] = (doublereal) lwmin;
301 if (*lwork < lwmin && ! lquery)
302 {
303 *info = -15;
304 }
305 }
306 if (*info != 0)
307 {
308 i__1 = -(*info);
309 xerbla_("DTGEXC", &i__1);
310 return 0;
311 }
312 else if (lquery)
313 {
314 return 0;
315 }
316 /* Quick return if possible */
317 if (*n <= 1)
318 {
319 return 0;
320 }
321 /* Determine the first row of the specified block and find out */
322 /* if it is 1-by-1 or 2-by-2. */
323 if (*ifst > 1)
324 {
325 if (a[*ifst + (*ifst - 1) * a_dim1] != 0.)
326 {
327 --(*ifst);
328 }
329 }
330 nbf = 1;
331 if (*ifst < *n)
332 {
333 if (a[*ifst + 1 + *ifst * a_dim1] != 0.)
334 {
335 nbf = 2;
336 }
337 }
338 /* Determine the first row of the final block */
339 /* and find out if it is 1-by-1 or 2-by-2. */
340 if (*ilst > 1)
341 {
342 if (a[*ilst + (*ilst - 1) * a_dim1] != 0.)
343 {
344 --(*ilst);
345 }
346 }
347 nbl = 1;
348 if (*ilst < *n)
349 {
350 if (a[*ilst + 1 + *ilst * a_dim1] != 0.)
351 {
352 nbl = 2;
353 }
354 }
355 if (*ifst == *ilst)
356 {
357 return 0;
358 }
359 if (*ifst < *ilst)
360 {
361 /* Update ILST. */
362 if (nbf == 2 && nbl == 1)
363 {
364 --(*ilst);
365 }
366 if (nbf == 1 && nbl == 2)
367 {
368 ++(*ilst);
369 }
370 here = *ifst;
371 L10: /* Swap with next one below. */
372 if (nbf == 1 || nbf == 2)
373 {
374 /* Current block either 1-by-1 or 2-by-2. */
375 nbnext = 1;
376 if (here + nbf + 1 <= *n)
377 {
378 if (a[here + nbf + 1 + (here + nbf) * a_dim1] != 0.)
379 {
380 nbnext = 2;
381 }
382 }
383 dtgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[ q_offset], ldq, &z__[z_offset], ldz, &here, &nbf, &nbnext, &work[1], lwork, info);
384 if (*info != 0)
385 {
386 *ilst = here;
387 return 0;
388 }
389 here += nbnext;
390 /* Test if 2-by-2 block breaks into two 1-by-1 blocks. */
391 if (nbf == 2)
392 {
393 if (a[here + 1 + here * a_dim1] == 0.)
394 {
395 nbf = 3;
396 }
397 }
398 }
399 else
400 {
401 /* Current block consists of two 1-by-1 blocks, each of which */
402 /* must be swapped individually. */
403 nbnext = 1;
404 if (here + 3 <= *n)
405 {
406 if (a[here + 3 + (here + 2) * a_dim1] != 0.)
407 {
408 nbnext = 2;
409 }
410 }
411 i__1 = here + 1;
412 dtgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[ q_offset], ldq, &z__[z_offset], ldz, &i__1, &c__1, & nbnext, &work[1], lwork, info);
413 if (*info != 0)
414 {
415 *ilst = here;
416 return 0;
417 }
418 if (nbnext == 1)
419 {
420 /* Swap two 1-by-1 blocks. */
421 dtgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[q_offset], ldq, &z__[z_offset], ldz, &here, &c__1, &c__1, &work[1], lwork, info);
422 if (*info != 0)
423 {
424 *ilst = here;
425 return 0;
426 }
427 ++here;
428 }
429 else
430 {
431 /* Recompute NBNEXT in case of 2-by-2 split. */
432 if (a[here + 2 + (here + 1) * a_dim1] == 0.)
433 {
434 nbnext = 1;
435 }
436 if (nbnext == 2)
437 {
438 /* 2-by-2 block did not split. */
439 dtgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[q_offset], ldq, &z__[z_offset], ldz, & here, &c__1, &nbnext, &work[1], lwork, info);
440 if (*info != 0)
441 {
442 *ilst = here;
443 return 0;
444 }
445 here += 2;
446 }
447 else
448 {
449 /* 2-by-2 block did split. */
450 dtgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[q_offset], ldq, &z__[z_offset], ldz, & here, &c__1, &c__1, &work[1], lwork, info);
451 if (*info != 0)
452 {
453 *ilst = here;
454 return 0;
455 }
456 ++here;
457 dtgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[q_offset], ldq, &z__[z_offset], ldz, & here, &c__1, &c__1, &work[1], lwork, info);
458 if (*info != 0)
459 {
460 *ilst = here;
461 return 0;
462 }
463 ++here;
464 }
465 }
466 }
467 if (here < *ilst)
468 {
469 goto L10;
470 }
471 }
472 else
473 {
474 here = *ifst;
475 L20: /* Swap with next one below. */
476 if (nbf == 1 || nbf == 2)
477 {
478 /* Current block either 1-by-1 or 2-by-2. */
479 nbnext = 1;
480 if (here >= 3)
481 {
482 if (a[here - 1 + (here - 2) * a_dim1] != 0.)
483 {
484 nbnext = 2;
485 }
486 }
487 i__1 = here - nbnext;
488 dtgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[ q_offset], ldq, &z__[z_offset], ldz, &i__1, &nbnext, &nbf, &work[1], lwork, info);
489 if (*info != 0)
490 {
491 *ilst = here;
492 return 0;
493 }
494 here -= nbnext;
495 /* Test if 2-by-2 block breaks into two 1-by-1 blocks. */
496 if (nbf == 2)
497 {
498 if (a[here + 1 + here * a_dim1] == 0.)
499 {
500 nbf = 3;
501 }
502 }
503 }
504 else
505 {
506 /* Current block consists of two 1-by-1 blocks, each of which */
507 /* must be swapped individually. */
508 nbnext = 1;
509 if (here >= 3)
510 {
511 if (a[here - 1 + (here - 2) * a_dim1] != 0.)
512 {
513 nbnext = 2;
514 }
515 }
516 i__1 = here - nbnext;
517 dtgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[ q_offset], ldq, &z__[z_offset], ldz, &i__1, &nbnext, & c__1, &work[1], lwork, info);
518 if (*info != 0)
519 {
520 *ilst = here;
521 return 0;
522 }
523 if (nbnext == 1)
524 {
525 /* Swap two 1-by-1 blocks. */
526 dtgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[q_offset], ldq, &z__[z_offset], ldz, &here, & nbnext, &c__1, &work[1], lwork, info);
527 if (*info != 0)
528 {
529 *ilst = here;
530 return 0;
531 }
532 --here;
533 }
534 else
535 {
536 /* Recompute NBNEXT in case of 2-by-2 split. */
537 if (a[here + (here - 1) * a_dim1] == 0.)
538 {
539 nbnext = 1;
540 }
541 if (nbnext == 2)
542 {
543 /* 2-by-2 block did not split. */
544 i__1 = here - 1;
545 dtgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[q_offset], ldq, &z__[z_offset], ldz, & i__1, &c__2, &c__1, &work[1], lwork, info);
546 if (*info != 0)
547 {
548 *ilst = here;
549 return 0;
550 }
551 here += -2;
552 }
553 else
554 {
555 /* 2-by-2 block did split. */
556 dtgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[q_offset], ldq, &z__[z_offset], ldz, & here, &c__1, &c__1, &work[1], lwork, info);
557 if (*info != 0)
558 {
559 *ilst = here;
560 return 0;
561 }
562 --here;
563 dtgex2_(wantq, wantz, n, &a[a_offset], lda, &b[b_offset], ldb, &q[q_offset], ldq, &z__[z_offset], ldz, & here, &c__1, &c__1, &work[1], lwork, info);
564 if (*info != 0)
565 {
566 *ilst = here;
567 return 0;
568 }
569 --here;
570 }
571 }
572 }
573 if (here > *ilst)
574 {
575 goto L20;
576 }
577 }
578 *ilst = here;
579 work[1] = (doublereal) lwmin;
580 return 0;
581 /* End of DTGEXC */
582 }
583 /* dtgexc_ */
584