1 SUBROUTINE DTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
2 $ LDZ, IFST, ILST, WORK, LWORK, INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 LOGICAL WANTQ, WANTZ
11 INTEGER IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
15 $ WORK( * ), Z( LDZ, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DTGEXC reorders the generalized real Schur decomposition of a real
22 * matrix pair (A,B) using an orthogonal equivalence transformation
23 *
24 * (A, B) = Q * (A, B) * Z**T,
25 *
26 * so that the diagonal block of (A, B) with row index IFST is moved
27 * to row ILST.
28 *
29 * (A, B) must be in generalized real Schur canonical form (as returned
30 * by DGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2
31 * diagonal blocks. B is upper triangular.
32 *
33 * Optionally, the matrices Q and Z of generalized Schur vectors are
34 * updated.
35 *
36 * Q(in) * A(in) * Z(in)**T = Q(out) * A(out) * Z(out)**T
37 * Q(in) * B(in) * Z(in)**T = Q(out) * B(out) * Z(out)**T
38 *
39 *
40 * Arguments
41 * =========
42 *
43 * WANTQ (input) LOGICAL
44 * .TRUE. : update the left transformation matrix Q;
45 * .FALSE.: do not update Q.
46 *
47 * WANTZ (input) LOGICAL
48 * .TRUE. : update the right transformation matrix Z;
49 * .FALSE.: do not update Z.
50 *
51 * N (input) INTEGER
52 * The order of the matrices A and B. N >= 0.
53 *
54 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
55 * On entry, the matrix A in generalized real Schur canonical
56 * form.
57 * On exit, the updated matrix A, again in generalized
58 * real Schur canonical form.
59 *
60 * LDA (input) INTEGER
61 * The leading dimension of the array A. LDA >= max(1,N).
62 *
63 * B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
64 * On entry, the matrix B in generalized real Schur canonical
65 * form (A,B).
66 * On exit, the updated matrix B, again in generalized
67 * real Schur canonical form (A,B).
68 *
69 * LDB (input) INTEGER
70 * The leading dimension of the array B. LDB >= max(1,N).
71 *
72 * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
73 * On entry, if WANTQ = .TRUE., the orthogonal matrix Q.
74 * On exit, the updated matrix Q.
75 * If WANTQ = .FALSE., Q is not referenced.
76 *
77 * LDQ (input) INTEGER
78 * The leading dimension of the array Q. LDQ >= 1.
79 * If WANTQ = .TRUE., LDQ >= N.
80 *
81 * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
82 * On entry, if WANTZ = .TRUE., the orthogonal matrix Z.
83 * On exit, the updated matrix Z.
84 * If WANTZ = .FALSE., Z is not referenced.
85 *
86 * LDZ (input) INTEGER
87 * The leading dimension of the array Z. LDZ >= 1.
88 * If WANTZ = .TRUE., LDZ >= N.
89 *
90 * IFST (input/output) INTEGER
91 * ILST (input/output) INTEGER
92 * Specify the reordering of the diagonal blocks of (A, B).
93 * The block with row index IFST is moved to row ILST, by a
94 * sequence of swapping between adjacent blocks.
95 * On exit, if IFST pointed on entry to the second row of
96 * a 2-by-2 block, it is changed to point to the first row;
97 * ILST always points to the first row of the block in its
98 * final position (which may differ from its input value by
99 * +1 or -1). 1 <= IFST, ILST <= N.
100 *
101 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
102 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
103 *
104 * LWORK (input) INTEGER
105 * The dimension of the array WORK.
106 * LWORK >= 1 when N <= 1, otherwise LWORK >= 4*N + 16.
107 *
108 * If LWORK = -1, then a workspace query is assumed; the routine
109 * only calculates the optimal size of the WORK array, returns
110 * this value as the first entry of the WORK array, and no error
111 * message related to LWORK is issued by XERBLA.
112 *
113 * INFO (output) INTEGER
114 * =0: successful exit.
115 * <0: if INFO = -i, the i-th argument had an illegal value.
116 * =1: The transformed matrix pair (A, B) would be too far
117 * from generalized Schur form; the problem is ill-
118 * conditioned. (A, B) may have been partially reordered,
119 * and ILST points to the first row of the current
120 * position of the block being moved.
121 *
122 * Further Details
123 * ===============
124 *
125 * Based on contributions by
126 * Bo Kagstrom and Peter Poromaa, Department of Computing Science,
127 * Umea University, S-901 87 Umea, Sweden.
128 *
129 * [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
130 * Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
131 * M.S. Moonen et al (eds), Linear Algebra for Large Scale and
132 * Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
133 *
134 * =====================================================================
135 *
136 * .. Parameters ..
137 DOUBLE PRECISION ZERO
138 PARAMETER ( ZERO = 0.0D+0 )
139 * ..
140 * .. Local Scalars ..
141 LOGICAL LQUERY
142 INTEGER HERE, LWMIN, NBF, NBL, NBNEXT
143 * ..
144 * .. External Subroutines ..
145 EXTERNAL DTGEX2, XERBLA
146 * ..
147 * .. Intrinsic Functions ..
148 INTRINSIC MAX
149 * ..
150 * .. Executable Statements ..
151 *
152 * Decode and test input arguments.
153 *
154 INFO = 0
155 LQUERY = ( LWORK.EQ.-1 )
156 IF( N.LT.0 ) THEN
157 INFO = -3
158 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
159 INFO = -5
160 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
161 INFO = -7
162 ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. ( LDQ.LT.MAX( 1, N ) ) ) THEN
163 INFO = -9
164 ELSE IF( LDZ.LT.1 .OR. WANTZ .AND. ( LDZ.LT.MAX( 1, N ) ) ) THEN
165 INFO = -11
166 ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
167 INFO = -12
168 ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
169 INFO = -13
170 END IF
171 *
172 IF( INFO.EQ.0 ) THEN
173 IF( N.LE.1 ) THEN
174 LWMIN = 1
175 ELSE
176 LWMIN = 4*N + 16
177 END IF
178 WORK(1) = LWMIN
179 *
180 IF (LWORK.LT.LWMIN .AND. .NOT.LQUERY) THEN
181 INFO = -15
182 END IF
183 END IF
184 *
185 IF( INFO.NE.0 ) THEN
186 CALL XERBLA( 'DTGEXC', -INFO )
187 RETURN
188 ELSE IF( LQUERY ) THEN
189 RETURN
190 END IF
191 *
192 * Quick return if possible
193 *
194 IF( N.LE.1 )
195 $ RETURN
196 *
197 * Determine the first row of the specified block and find out
198 * if it is 1-by-1 or 2-by-2.
199 *
200 IF( IFST.GT.1 ) THEN
201 IF( A( IFST, IFST-1 ).NE.ZERO )
202 $ IFST = IFST - 1
203 END IF
204 NBF = 1
205 IF( IFST.LT.N ) THEN
206 IF( A( IFST+1, IFST ).NE.ZERO )
207 $ NBF = 2
208 END IF
209 *
210 * Determine the first row of the final block
211 * and find out if it is 1-by-1 or 2-by-2.
212 *
213 IF( ILST.GT.1 ) THEN
214 IF( A( ILST, ILST-1 ).NE.ZERO )
215 $ ILST = ILST - 1
216 END IF
217 NBL = 1
218 IF( ILST.LT.N ) THEN
219 IF( A( ILST+1, ILST ).NE.ZERO )
220 $ NBL = 2
221 END IF
222 IF( IFST.EQ.ILST )
223 $ RETURN
224 *
225 IF( IFST.LT.ILST ) THEN
226 *
227 * Update ILST.
228 *
229 IF( NBF.EQ.2 .AND. NBL.EQ.1 )
230 $ ILST = ILST - 1
231 IF( NBF.EQ.1 .AND. NBL.EQ.2 )
232 $ ILST = ILST + 1
233 *
234 HERE = IFST
235 *
236 10 CONTINUE
237 *
238 * Swap with next one below.
239 *
240 IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
241 *
242 * Current block either 1-by-1 or 2-by-2.
243 *
244 NBNEXT = 1
245 IF( HERE+NBF+1.LE.N ) THEN
246 IF( A( HERE+NBF+1, HERE+NBF ).NE.ZERO )
247 $ NBNEXT = 2
248 END IF
249 CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
250 $ LDZ, HERE, NBF, NBNEXT, WORK, LWORK, INFO )
251 IF( INFO.NE.0 ) THEN
252 ILST = HERE
253 RETURN
254 END IF
255 HERE = HERE + NBNEXT
256 *
257 * Test if 2-by-2 block breaks into two 1-by-1 blocks.
258 *
259 IF( NBF.EQ.2 ) THEN
260 IF( A( HERE+1, HERE ).EQ.ZERO )
261 $ NBF = 3
262 END IF
263 *
264 ELSE
265 *
266 * Current block consists of two 1-by-1 blocks, each of which
267 * must be swapped individually.
268 *
269 NBNEXT = 1
270 IF( HERE+3.LE.N ) THEN
271 IF( A( HERE+3, HERE+2 ).NE.ZERO )
272 $ NBNEXT = 2
273 END IF
274 CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
275 $ LDZ, HERE+1, 1, NBNEXT, WORK, LWORK, INFO )
276 IF( INFO.NE.0 ) THEN
277 ILST = HERE
278 RETURN
279 END IF
280 IF( NBNEXT.EQ.1 ) THEN
281 *
282 * Swap two 1-by-1 blocks.
283 *
284 CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
285 $ LDZ, HERE, 1, 1, WORK, LWORK, INFO )
286 IF( INFO.NE.0 ) THEN
287 ILST = HERE
288 RETURN
289 END IF
290 HERE = HERE + 1
291 *
292 ELSE
293 *
294 * Recompute NBNEXT in case of 2-by-2 split.
295 *
296 IF( A( HERE+2, HERE+1 ).EQ.ZERO )
297 $ NBNEXT = 1
298 IF( NBNEXT.EQ.2 ) THEN
299 *
300 * 2-by-2 block did not split.
301 *
302 CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
303 $ Z, LDZ, HERE, 1, NBNEXT, WORK, LWORK,
304 $ INFO )
305 IF( INFO.NE.0 ) THEN
306 ILST = HERE
307 RETURN
308 END IF
309 HERE = HERE + 2
310 ELSE
311 *
312 * 2-by-2 block did split.
313 *
314 CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
315 $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
316 IF( INFO.NE.0 ) THEN
317 ILST = HERE
318 RETURN
319 END IF
320 HERE = HERE + 1
321 CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
322 $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
323 IF( INFO.NE.0 ) THEN
324 ILST = HERE
325 RETURN
326 END IF
327 HERE = HERE + 1
328 END IF
329 *
330 END IF
331 END IF
332 IF( HERE.LT.ILST )
333 $ GO TO 10
334 ELSE
335 HERE = IFST
336 *
337 20 CONTINUE
338 *
339 * Swap with next one below.
340 *
341 IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
342 *
343 * Current block either 1-by-1 or 2-by-2.
344 *
345 NBNEXT = 1
346 IF( HERE.GE.3 ) THEN
347 IF( A( HERE-1, HERE-2 ).NE.ZERO )
348 $ NBNEXT = 2
349 END IF
350 CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
351 $ LDZ, HERE-NBNEXT, NBNEXT, NBF, WORK, LWORK,
352 $ INFO )
353 IF( INFO.NE.0 ) THEN
354 ILST = HERE
355 RETURN
356 END IF
357 HERE = HERE - NBNEXT
358 *
359 * Test if 2-by-2 block breaks into two 1-by-1 blocks.
360 *
361 IF( NBF.EQ.2 ) THEN
362 IF( A( HERE+1, HERE ).EQ.ZERO )
363 $ NBF = 3
364 END IF
365 *
366 ELSE
367 *
368 * Current block consists of two 1-by-1 blocks, each of which
369 * must be swapped individually.
370 *
371 NBNEXT = 1
372 IF( HERE.GE.3 ) THEN
373 IF( A( HERE-1, HERE-2 ).NE.ZERO )
374 $ NBNEXT = 2
375 END IF
376 CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
377 $ LDZ, HERE-NBNEXT, NBNEXT, 1, WORK, LWORK,
378 $ INFO )
379 IF( INFO.NE.0 ) THEN
380 ILST = HERE
381 RETURN
382 END IF
383 IF( NBNEXT.EQ.1 ) THEN
384 *
385 * Swap two 1-by-1 blocks.
386 *
387 CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
388 $ LDZ, HERE, NBNEXT, 1, WORK, LWORK, INFO )
389 IF( INFO.NE.0 ) THEN
390 ILST = HERE
391 RETURN
392 END IF
393 HERE = HERE - 1
394 ELSE
395 *
396 * Recompute NBNEXT in case of 2-by-2 split.
397 *
398 IF( A( HERE, HERE-1 ).EQ.ZERO )
399 $ NBNEXT = 1
400 IF( NBNEXT.EQ.2 ) THEN
401 *
402 * 2-by-2 block did not split.
403 *
404 CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
405 $ Z, LDZ, HERE-1, 2, 1, WORK, LWORK, INFO )
406 IF( INFO.NE.0 ) THEN
407 ILST = HERE
408 RETURN
409 END IF
410 HERE = HERE - 2
411 ELSE
412 *
413 * 2-by-2 block did split.
414 *
415 CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
416 $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
417 IF( INFO.NE.0 ) THEN
418 ILST = HERE
419 RETURN
420 END IF
421 HERE = HERE - 1
422 CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
423 $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
424 IF( INFO.NE.0 ) THEN
425 ILST = HERE
426 RETURN
427 END IF
428 HERE = HERE - 1
429 END IF
430 END IF
431 END IF
432 IF( HERE.GT.ILST )
433 $ GO TO 20
434 END IF
435 ILST = HERE
436 WORK( 1 ) = LWMIN
437 RETURN
438 *
439 * End of DTGEXC
440 *
441 END
2 $ LDZ, IFST, ILST, WORK, LWORK, INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 LOGICAL WANTQ, WANTZ
11 INTEGER IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
15 $ WORK( * ), Z( LDZ, * )
16 * ..
17 *
18 * Purpose
19 * =======
20 *
21 * DTGEXC reorders the generalized real Schur decomposition of a real
22 * matrix pair (A,B) using an orthogonal equivalence transformation
23 *
24 * (A, B) = Q * (A, B) * Z**T,
25 *
26 * so that the diagonal block of (A, B) with row index IFST is moved
27 * to row ILST.
28 *
29 * (A, B) must be in generalized real Schur canonical form (as returned
30 * by DGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2
31 * diagonal blocks. B is upper triangular.
32 *
33 * Optionally, the matrices Q and Z of generalized Schur vectors are
34 * updated.
35 *
36 * Q(in) * A(in) * Z(in)**T = Q(out) * A(out) * Z(out)**T
37 * Q(in) * B(in) * Z(in)**T = Q(out) * B(out) * Z(out)**T
38 *
39 *
40 * Arguments
41 * =========
42 *
43 * WANTQ (input) LOGICAL
44 * .TRUE. : update the left transformation matrix Q;
45 * .FALSE.: do not update Q.
46 *
47 * WANTZ (input) LOGICAL
48 * .TRUE. : update the right transformation matrix Z;
49 * .FALSE.: do not update Z.
50 *
51 * N (input) INTEGER
52 * The order of the matrices A and B. N >= 0.
53 *
54 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
55 * On entry, the matrix A in generalized real Schur canonical
56 * form.
57 * On exit, the updated matrix A, again in generalized
58 * real Schur canonical form.
59 *
60 * LDA (input) INTEGER
61 * The leading dimension of the array A. LDA >= max(1,N).
62 *
63 * B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
64 * On entry, the matrix B in generalized real Schur canonical
65 * form (A,B).
66 * On exit, the updated matrix B, again in generalized
67 * real Schur canonical form (A,B).
68 *
69 * LDB (input) INTEGER
70 * The leading dimension of the array B. LDB >= max(1,N).
71 *
72 * Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
73 * On entry, if WANTQ = .TRUE., the orthogonal matrix Q.
74 * On exit, the updated matrix Q.
75 * If WANTQ = .FALSE., Q is not referenced.
76 *
77 * LDQ (input) INTEGER
78 * The leading dimension of the array Q. LDQ >= 1.
79 * If WANTQ = .TRUE., LDQ >= N.
80 *
81 * Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
82 * On entry, if WANTZ = .TRUE., the orthogonal matrix Z.
83 * On exit, the updated matrix Z.
84 * If WANTZ = .FALSE., Z is not referenced.
85 *
86 * LDZ (input) INTEGER
87 * The leading dimension of the array Z. LDZ >= 1.
88 * If WANTZ = .TRUE., LDZ >= N.
89 *
90 * IFST (input/output) INTEGER
91 * ILST (input/output) INTEGER
92 * Specify the reordering of the diagonal blocks of (A, B).
93 * The block with row index IFST is moved to row ILST, by a
94 * sequence of swapping between adjacent blocks.
95 * On exit, if IFST pointed on entry to the second row of
96 * a 2-by-2 block, it is changed to point to the first row;
97 * ILST always points to the first row of the block in its
98 * final position (which may differ from its input value by
99 * +1 or -1). 1 <= IFST, ILST <= N.
100 *
101 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
102 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
103 *
104 * LWORK (input) INTEGER
105 * The dimension of the array WORK.
106 * LWORK >= 1 when N <= 1, otherwise LWORK >= 4*N + 16.
107 *
108 * If LWORK = -1, then a workspace query is assumed; the routine
109 * only calculates the optimal size of the WORK array, returns
110 * this value as the first entry of the WORK array, and no error
111 * message related to LWORK is issued by XERBLA.
112 *
113 * INFO (output) INTEGER
114 * =0: successful exit.
115 * <0: if INFO = -i, the i-th argument had an illegal value.
116 * =1: The transformed matrix pair (A, B) would be too far
117 * from generalized Schur form; the problem is ill-
118 * conditioned. (A, B) may have been partially reordered,
119 * and ILST points to the first row of the current
120 * position of the block being moved.
121 *
122 * Further Details
123 * ===============
124 *
125 * Based on contributions by
126 * Bo Kagstrom and Peter Poromaa, Department of Computing Science,
127 * Umea University, S-901 87 Umea, Sweden.
128 *
129 * [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
130 * Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
131 * M.S. Moonen et al (eds), Linear Algebra for Large Scale and
132 * Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
133 *
134 * =====================================================================
135 *
136 * .. Parameters ..
137 DOUBLE PRECISION ZERO
138 PARAMETER ( ZERO = 0.0D+0 )
139 * ..
140 * .. Local Scalars ..
141 LOGICAL LQUERY
142 INTEGER HERE, LWMIN, NBF, NBL, NBNEXT
143 * ..
144 * .. External Subroutines ..
145 EXTERNAL DTGEX2, XERBLA
146 * ..
147 * .. Intrinsic Functions ..
148 INTRINSIC MAX
149 * ..
150 * .. Executable Statements ..
151 *
152 * Decode and test input arguments.
153 *
154 INFO = 0
155 LQUERY = ( LWORK.EQ.-1 )
156 IF( N.LT.0 ) THEN
157 INFO = -3
158 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
159 INFO = -5
160 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
161 INFO = -7
162 ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. ( LDQ.LT.MAX( 1, N ) ) ) THEN
163 INFO = -9
164 ELSE IF( LDZ.LT.1 .OR. WANTZ .AND. ( LDZ.LT.MAX( 1, N ) ) ) THEN
165 INFO = -11
166 ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
167 INFO = -12
168 ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
169 INFO = -13
170 END IF
171 *
172 IF( INFO.EQ.0 ) THEN
173 IF( N.LE.1 ) THEN
174 LWMIN = 1
175 ELSE
176 LWMIN = 4*N + 16
177 END IF
178 WORK(1) = LWMIN
179 *
180 IF (LWORK.LT.LWMIN .AND. .NOT.LQUERY) THEN
181 INFO = -15
182 END IF
183 END IF
184 *
185 IF( INFO.NE.0 ) THEN
186 CALL XERBLA( 'DTGEXC', -INFO )
187 RETURN
188 ELSE IF( LQUERY ) THEN
189 RETURN
190 END IF
191 *
192 * Quick return if possible
193 *
194 IF( N.LE.1 )
195 $ RETURN
196 *
197 * Determine the first row of the specified block and find out
198 * if it is 1-by-1 or 2-by-2.
199 *
200 IF( IFST.GT.1 ) THEN
201 IF( A( IFST, IFST-1 ).NE.ZERO )
202 $ IFST = IFST - 1
203 END IF
204 NBF = 1
205 IF( IFST.LT.N ) THEN
206 IF( A( IFST+1, IFST ).NE.ZERO )
207 $ NBF = 2
208 END IF
209 *
210 * Determine the first row of the final block
211 * and find out if it is 1-by-1 or 2-by-2.
212 *
213 IF( ILST.GT.1 ) THEN
214 IF( A( ILST, ILST-1 ).NE.ZERO )
215 $ ILST = ILST - 1
216 END IF
217 NBL = 1
218 IF( ILST.LT.N ) THEN
219 IF( A( ILST+1, ILST ).NE.ZERO )
220 $ NBL = 2
221 END IF
222 IF( IFST.EQ.ILST )
223 $ RETURN
224 *
225 IF( IFST.LT.ILST ) THEN
226 *
227 * Update ILST.
228 *
229 IF( NBF.EQ.2 .AND. NBL.EQ.1 )
230 $ ILST = ILST - 1
231 IF( NBF.EQ.1 .AND. NBL.EQ.2 )
232 $ ILST = ILST + 1
233 *
234 HERE = IFST
235 *
236 10 CONTINUE
237 *
238 * Swap with next one below.
239 *
240 IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
241 *
242 * Current block either 1-by-1 or 2-by-2.
243 *
244 NBNEXT = 1
245 IF( HERE+NBF+1.LE.N ) THEN
246 IF( A( HERE+NBF+1, HERE+NBF ).NE.ZERO )
247 $ NBNEXT = 2
248 END IF
249 CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
250 $ LDZ, HERE, NBF, NBNEXT, WORK, LWORK, INFO )
251 IF( INFO.NE.0 ) THEN
252 ILST = HERE
253 RETURN
254 END IF
255 HERE = HERE + NBNEXT
256 *
257 * Test if 2-by-2 block breaks into two 1-by-1 blocks.
258 *
259 IF( NBF.EQ.2 ) THEN
260 IF( A( HERE+1, HERE ).EQ.ZERO )
261 $ NBF = 3
262 END IF
263 *
264 ELSE
265 *
266 * Current block consists of two 1-by-1 blocks, each of which
267 * must be swapped individually.
268 *
269 NBNEXT = 1
270 IF( HERE+3.LE.N ) THEN
271 IF( A( HERE+3, HERE+2 ).NE.ZERO )
272 $ NBNEXT = 2
273 END IF
274 CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
275 $ LDZ, HERE+1, 1, NBNEXT, WORK, LWORK, INFO )
276 IF( INFO.NE.0 ) THEN
277 ILST = HERE
278 RETURN
279 END IF
280 IF( NBNEXT.EQ.1 ) THEN
281 *
282 * Swap two 1-by-1 blocks.
283 *
284 CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
285 $ LDZ, HERE, 1, 1, WORK, LWORK, INFO )
286 IF( INFO.NE.0 ) THEN
287 ILST = HERE
288 RETURN
289 END IF
290 HERE = HERE + 1
291 *
292 ELSE
293 *
294 * Recompute NBNEXT in case of 2-by-2 split.
295 *
296 IF( A( HERE+2, HERE+1 ).EQ.ZERO )
297 $ NBNEXT = 1
298 IF( NBNEXT.EQ.2 ) THEN
299 *
300 * 2-by-2 block did not split.
301 *
302 CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
303 $ Z, LDZ, HERE, 1, NBNEXT, WORK, LWORK,
304 $ INFO )
305 IF( INFO.NE.0 ) THEN
306 ILST = HERE
307 RETURN
308 END IF
309 HERE = HERE + 2
310 ELSE
311 *
312 * 2-by-2 block did split.
313 *
314 CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
315 $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
316 IF( INFO.NE.0 ) THEN
317 ILST = HERE
318 RETURN
319 END IF
320 HERE = HERE + 1
321 CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
322 $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
323 IF( INFO.NE.0 ) THEN
324 ILST = HERE
325 RETURN
326 END IF
327 HERE = HERE + 1
328 END IF
329 *
330 END IF
331 END IF
332 IF( HERE.LT.ILST )
333 $ GO TO 10
334 ELSE
335 HERE = IFST
336 *
337 20 CONTINUE
338 *
339 * Swap with next one below.
340 *
341 IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
342 *
343 * Current block either 1-by-1 or 2-by-2.
344 *
345 NBNEXT = 1
346 IF( HERE.GE.3 ) THEN
347 IF( A( HERE-1, HERE-2 ).NE.ZERO )
348 $ NBNEXT = 2
349 END IF
350 CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
351 $ LDZ, HERE-NBNEXT, NBNEXT, NBF, WORK, LWORK,
352 $ INFO )
353 IF( INFO.NE.0 ) THEN
354 ILST = HERE
355 RETURN
356 END IF
357 HERE = HERE - NBNEXT
358 *
359 * Test if 2-by-2 block breaks into two 1-by-1 blocks.
360 *
361 IF( NBF.EQ.2 ) THEN
362 IF( A( HERE+1, HERE ).EQ.ZERO )
363 $ NBF = 3
364 END IF
365 *
366 ELSE
367 *
368 * Current block consists of two 1-by-1 blocks, each of which
369 * must be swapped individually.
370 *
371 NBNEXT = 1
372 IF( HERE.GE.3 ) THEN
373 IF( A( HERE-1, HERE-2 ).NE.ZERO )
374 $ NBNEXT = 2
375 END IF
376 CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
377 $ LDZ, HERE-NBNEXT, NBNEXT, 1, WORK, LWORK,
378 $ INFO )
379 IF( INFO.NE.0 ) THEN
380 ILST = HERE
381 RETURN
382 END IF
383 IF( NBNEXT.EQ.1 ) THEN
384 *
385 * Swap two 1-by-1 blocks.
386 *
387 CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
388 $ LDZ, HERE, NBNEXT, 1, WORK, LWORK, INFO )
389 IF( INFO.NE.0 ) THEN
390 ILST = HERE
391 RETURN
392 END IF
393 HERE = HERE - 1
394 ELSE
395 *
396 * Recompute NBNEXT in case of 2-by-2 split.
397 *
398 IF( A( HERE, HERE-1 ).EQ.ZERO )
399 $ NBNEXT = 1
400 IF( NBNEXT.EQ.2 ) THEN
401 *
402 * 2-by-2 block did not split.
403 *
404 CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
405 $ Z, LDZ, HERE-1, 2, 1, WORK, LWORK, INFO )
406 IF( INFO.NE.0 ) THEN
407 ILST = HERE
408 RETURN
409 END IF
410 HERE = HERE - 2
411 ELSE
412 *
413 * 2-by-2 block did split.
414 *
415 CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
416 $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
417 IF( INFO.NE.0 ) THEN
418 ILST = HERE
419 RETURN
420 END IF
421 HERE = HERE - 1
422 CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
423 $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
424 IF( INFO.NE.0 ) THEN
425 ILST = HERE
426 RETURN
427 END IF
428 HERE = HERE - 1
429 END IF
430 END IF
431 END IF
432 IF( HERE.GT.ILST )
433 $ GO TO 20
434 END IF
435 ILST = HERE
436 WORK( 1 ) = LWMIN
437 RETURN
438 *
439 * End of DTGEXC
440 *
441 END