001:       SUBROUTINE DTGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
002:      $                   LDZ, IFST, ILST, WORK, LWORK, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       LOGICAL            WANTQ, WANTZ
010:       INTEGER            IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
011: *     ..
012: *     .. Array Arguments ..
013:       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
014:      $                   WORK( * ), Z( LDZ, * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  DTGEXC reorders the generalized real Schur decomposition of a real
021: *  matrix pair (A,B) using an orthogonal equivalence transformation
022: *
023: *                 (A, B) = Q * (A, B) * Z',
024: *
025: *  so that the diagonal block of (A, B) with row index IFST is moved
026: *  to row ILST.
027: *
028: *  (A, B) must be in generalized real Schur canonical form (as returned
029: *  by DGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2
030: *  diagonal blocks. B is upper triangular.
031: *
032: *  Optionally, the matrices Q and Z of generalized Schur vectors are
033: *  updated.
034: *
035: *         Q(in) * A(in) * Z(in)' = Q(out) * A(out) * Z(out)'
036: *         Q(in) * B(in) * Z(in)' = Q(out) * B(out) * Z(out)'
037: *
038: *
039: *  Arguments
040: *  =========
041: *
042: *  WANTQ   (input) LOGICAL
043: *          .TRUE. : update the left transformation matrix Q;
044: *          .FALSE.: do not update Q.
045: *
046: *  WANTZ   (input) LOGICAL
047: *          .TRUE. : update the right transformation matrix Z;
048: *          .FALSE.: do not update Z.
049: *
050: *  N       (input) INTEGER
051: *          The order of the matrices A and B. N >= 0.
052: *
053: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
054: *          On entry, the matrix A in generalized real Schur canonical
055: *          form.
056: *          On exit, the updated matrix A, again in generalized
057: *          real Schur canonical form.
058: *
059: *  LDA     (input)  INTEGER
060: *          The leading dimension of the array A. LDA >= max(1,N).
061: *
062: *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,N)
063: *          On entry, the matrix B in generalized real Schur canonical
064: *          form (A,B).
065: *          On exit, the updated matrix B, again in generalized
066: *          real Schur canonical form (A,B).
067: *
068: *  LDB     (input)  INTEGER
069: *          The leading dimension of the array B. LDB >= max(1,N).
070: *
071: *  Q       (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
072: *          On entry, if WANTQ = .TRUE., the orthogonal matrix Q.
073: *          On exit, the updated matrix Q.
074: *          If WANTQ = .FALSE., Q is not referenced.
075: *
076: *  LDQ     (input) INTEGER
077: *          The leading dimension of the array Q. LDQ >= 1.
078: *          If WANTQ = .TRUE., LDQ >= N.
079: *
080: *  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
081: *          On entry, if WANTZ = .TRUE., the orthogonal matrix Z.
082: *          On exit, the updated matrix Z.
083: *          If WANTZ = .FALSE., Z is not referenced.
084: *
085: *  LDZ     (input) INTEGER
086: *          The leading dimension of the array Z. LDZ >= 1.
087: *          If WANTZ = .TRUE., LDZ >= N.
088: *
089: *  IFST    (input/output) INTEGER
090: *  ILST    (input/output) INTEGER
091: *          Specify the reordering of the diagonal blocks of (A, B).
092: *          The block with row index IFST is moved to row ILST, by a
093: *          sequence of swapping between adjacent blocks.
094: *          On exit, if IFST pointed on entry to the second row of
095: *          a 2-by-2 block, it is changed to point to the first row;
096: *          ILST always points to the first row of the block in its
097: *          final position (which may differ from its input value by
098: *          +1 or -1). 1 <= IFST, ILST <= N.
099: *
100: *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
101: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
102: *
103: *  LWORK   (input) INTEGER
104: *          The dimension of the array WORK.
105: *          LWORK >= 1 when N <= 1, otherwise LWORK >= 4*N + 16.
106: *
107: *          If LWORK = -1, then a workspace query is assumed; the routine
108: *          only calculates the optimal size of the WORK array, returns
109: *          this value as the first entry of the WORK array, and no error
110: *          message related to LWORK is issued by XERBLA.
111: *
112: *  INFO    (output) INTEGER
113: *           =0:  successful exit.
114: *           <0:  if INFO = -i, the i-th argument had an illegal value.
115: *           =1:  The transformed matrix pair (A, B) would be too far
116: *                from generalized Schur form; the problem is ill-
117: *                conditioned. (A, B) may have been partially reordered,
118: *                and ILST points to the first row of the current
119: *                position of the block being moved.
120: *
121: *  Further Details
122: *  ===============
123: *
124: *  Based on contributions by
125: *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
126: *     Umea University, S-901 87 Umea, Sweden.
127: *
128: *  [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the
129: *      Generalized Real Schur Form of a Regular Matrix Pair (A, B), in
130: *      M.S. Moonen et al (eds), Linear Algebra for Large Scale and
131: *      Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218.
132: *
133: *  =====================================================================
134: *
135: *     .. Parameters ..
136:       DOUBLE PRECISION   ZERO
137:       PARAMETER          ( ZERO = 0.0D+0 )
138: *     ..
139: *     .. Local Scalars ..
140:       LOGICAL            LQUERY
141:       INTEGER            HERE, LWMIN, NBF, NBL, NBNEXT
142: *     ..
143: *     .. External Subroutines ..
144:       EXTERNAL           DTGEX2, XERBLA
145: *     ..
146: *     .. Intrinsic Functions ..
147:       INTRINSIC          MAX
148: *     ..
149: *     .. Executable Statements ..
150: *
151: *     Decode and test input arguments.
152: *
153:       INFO = 0
154:       LQUERY = ( LWORK.EQ.-1 )
155:       IF( N.LT.0 ) THEN
156:          INFO = -3
157:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
158:          INFO = -5
159:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
160:          INFO = -7
161:       ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. ( LDQ.LT.MAX( 1, N ) ) ) THEN
162:          INFO = -9
163:       ELSE IF( LDZ.LT.1 .OR. WANTZ .AND. ( LDZ.LT.MAX( 1, N ) ) ) THEN
164:          INFO = -11
165:       ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
166:          INFO = -12
167:       ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
168:          INFO = -13
169:       END IF
170: *
171:       IF( INFO.EQ.0 ) THEN
172:          IF( N.LE.1 ) THEN
173:             LWMIN = 1
174:          ELSE
175:             LWMIN = 4*N + 16
176:          END IF
177:          WORK(1) = LWMIN
178: *
179:          IF (LWORK.LT.LWMIN .AND. .NOT.LQUERY) THEN
180:             INFO = -15
181:          END IF
182:       END IF
183: *
184:       IF( INFO.NE.0 ) THEN
185:          CALL XERBLA( 'DTGEXC', -INFO )
186:          RETURN
187:       ELSE IF( LQUERY ) THEN
188:          RETURN
189:       END IF
190: *
191: *     Quick return if possible
192: *
193:       IF( N.LE.1 )
194:      $   RETURN
195: *
196: *     Determine the first row of the specified block and find out
197: *     if it is 1-by-1 or 2-by-2.
198: *
199:       IF( IFST.GT.1 ) THEN
200:          IF( A( IFST, IFST-1 ).NE.ZERO )
201:      $      IFST = IFST - 1
202:       END IF
203:       NBF = 1
204:       IF( IFST.LT.N ) THEN
205:          IF( A( IFST+1, IFST ).NE.ZERO )
206:      $      NBF = 2
207:       END IF
208: *
209: *     Determine the first row of the final block
210: *     and find out if it is 1-by-1 or 2-by-2.
211: *
212:       IF( ILST.GT.1 ) THEN
213:          IF( A( ILST, ILST-1 ).NE.ZERO )
214:      $      ILST = ILST - 1
215:       END IF
216:       NBL = 1
217:       IF( ILST.LT.N ) THEN
218:          IF( A( ILST+1, ILST ).NE.ZERO )
219:      $      NBL = 2
220:       END IF
221:       IF( IFST.EQ.ILST )
222:      $   RETURN
223: *
224:       IF( IFST.LT.ILST ) THEN
225: *
226: *        Update ILST.
227: *
228:          IF( NBF.EQ.2 .AND. NBL.EQ.1 )
229:      $      ILST = ILST - 1
230:          IF( NBF.EQ.1 .AND. NBL.EQ.2 )
231:      $      ILST = ILST + 1
232: *
233:          HERE = IFST
234: *
235:    10    CONTINUE
236: *
237: *        Swap with next one below.
238: *
239:          IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
240: *
241: *           Current block either 1-by-1 or 2-by-2.
242: *
243:             NBNEXT = 1
244:             IF( HERE+NBF+1.LE.N ) THEN
245:                IF( A( HERE+NBF+1, HERE+NBF ).NE.ZERO )
246:      $            NBNEXT = 2
247:             END IF
248:             CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
249:      $                   LDZ, HERE, NBF, NBNEXT, WORK, LWORK, INFO )
250:             IF( INFO.NE.0 ) THEN
251:                ILST = HERE
252:                RETURN
253:             END IF
254:             HERE = HERE + NBNEXT
255: *
256: *           Test if 2-by-2 block breaks into two 1-by-1 blocks.
257: *
258:             IF( NBF.EQ.2 ) THEN
259:                IF( A( HERE+1, HERE ).EQ.ZERO )
260:      $            NBF = 3
261:             END IF
262: *
263:          ELSE
264: *
265: *           Current block consists of two 1-by-1 blocks, each of which
266: *           must be swapped individually.
267: *
268:             NBNEXT = 1
269:             IF( HERE+3.LE.N ) THEN
270:                IF( A( HERE+3, HERE+2 ).NE.ZERO )
271:      $            NBNEXT = 2
272:             END IF
273:             CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
274:      $                   LDZ, HERE+1, 1, NBNEXT, WORK, LWORK, INFO )
275:             IF( INFO.NE.0 ) THEN
276:                ILST = HERE
277:                RETURN
278:             END IF
279:             IF( NBNEXT.EQ.1 ) THEN
280: *
281: *              Swap two 1-by-1 blocks.
282: *
283:                CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
284:      $                      LDZ, HERE, 1, 1, WORK, LWORK, INFO )
285:                IF( INFO.NE.0 ) THEN
286:                   ILST = HERE
287:                   RETURN
288:                END IF
289:                HERE = HERE + 1
290: *
291:             ELSE
292: *
293: *              Recompute NBNEXT in case of 2-by-2 split.
294: *
295:                IF( A( HERE+2, HERE+1 ).EQ.ZERO )
296:      $            NBNEXT = 1
297:                IF( NBNEXT.EQ.2 ) THEN
298: *
299: *                 2-by-2 block did not split.
300: *
301:                   CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
302:      $                         Z, LDZ, HERE, 1, NBNEXT, WORK, LWORK,
303:      $                         INFO )
304:                   IF( INFO.NE.0 ) THEN
305:                      ILST = HERE
306:                      RETURN
307:                   END IF
308:                   HERE = HERE + 2
309:                ELSE
310: *
311: *                 2-by-2 block did split.
312: *
313:                   CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
314:      $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
315:                   IF( INFO.NE.0 ) THEN
316:                      ILST = HERE
317:                      RETURN
318:                   END IF
319:                   HERE = HERE + 1
320:                   CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
321:      $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
322:                   IF( INFO.NE.0 ) THEN
323:                      ILST = HERE
324:                      RETURN
325:                   END IF
326:                   HERE = HERE + 1
327:                END IF
328: *
329:             END IF
330:          END IF
331:          IF( HERE.LT.ILST )
332:      $      GO TO 10
333:       ELSE
334:          HERE = IFST
335: *
336:    20    CONTINUE
337: *
338: *        Swap with next one below.
339: *
340:          IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
341: *
342: *           Current block either 1-by-1 or 2-by-2.
343: *
344:             NBNEXT = 1
345:             IF( HERE.GE.3 ) THEN
346:                IF( A( HERE-1, HERE-2 ).NE.ZERO )
347:      $            NBNEXT = 2
348:             END IF
349:             CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
350:      $                   LDZ, HERE-NBNEXT, NBNEXT, NBF, WORK, LWORK,
351:      $                   INFO )
352:             IF( INFO.NE.0 ) THEN
353:                ILST = HERE
354:                RETURN
355:             END IF
356:             HERE = HERE - NBNEXT
357: *
358: *           Test if 2-by-2 block breaks into two 1-by-1 blocks.
359: *
360:             IF( NBF.EQ.2 ) THEN
361:                IF( A( HERE+1, HERE ).EQ.ZERO )
362:      $            NBF = 3
363:             END IF
364: *
365:          ELSE
366: *
367: *           Current block consists of two 1-by-1 blocks, each of which
368: *           must be swapped individually.
369: *
370:             NBNEXT = 1
371:             IF( HERE.GE.3 ) THEN
372:                IF( A( HERE-1, HERE-2 ).NE.ZERO )
373:      $            NBNEXT = 2
374:             END IF
375:             CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
376:      $                   LDZ, HERE-NBNEXT, NBNEXT, 1, WORK, LWORK,
377:      $                   INFO )
378:             IF( INFO.NE.0 ) THEN
379:                ILST = HERE
380:                RETURN
381:             END IF
382:             IF( NBNEXT.EQ.1 ) THEN
383: *
384: *              Swap two 1-by-1 blocks.
385: *
386:                CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z,
387:      $                      LDZ, HERE, NBNEXT, 1, WORK, LWORK, INFO )
388:                IF( INFO.NE.0 ) THEN
389:                   ILST = HERE
390:                   RETURN
391:                END IF
392:                HERE = HERE - 1
393:             ELSE
394: *
395: *             Recompute NBNEXT in case of 2-by-2 split.
396: *
397:                IF( A( HERE, HERE-1 ).EQ.ZERO )
398:      $            NBNEXT = 1
399:                IF( NBNEXT.EQ.2 ) THEN
400: *
401: *                 2-by-2 block did not split.
402: *
403:                   CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
404:      $                         Z, LDZ, HERE-1, 2, 1, WORK, LWORK, INFO )
405:                   IF( INFO.NE.0 ) THEN
406:                      ILST = HERE
407:                      RETURN
408:                   END IF
409:                   HERE = HERE - 2
410:                ELSE
411: *
412: *                 2-by-2 block did split.
413: *
414:                   CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
415:      $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
416:                   IF( INFO.NE.0 ) THEN
417:                      ILST = HERE
418:                      RETURN
419:                   END IF
420:                   HERE = HERE - 1
421:                   CALL DTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ,
422:      $                         Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO )
423:                   IF( INFO.NE.0 ) THEN
424:                      ILST = HERE
425:                      RETURN
426:                   END IF
427:                   HERE = HERE - 1
428:                END IF
429:             END IF
430:          END IF
431:          IF( HERE.GT.ILST )
432:      $      GO TO 20
433:       END IF
434:       ILST = HERE
435:       WORK( 1 ) = LWMIN
436:       RETURN
437: *
438: *     End of DTGEXC
439: *
440:       END
441: