001:       SUBROUTINE STGEXC( 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: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       LOGICAL            WANTQ, WANTZ
011:       INTEGER            IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N
012: *     ..
013: *     .. Array Arguments ..
014:       REAL               A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
015:      $                   WORK( * ), Z( LDZ, * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  STGEXC reorders the generalized real Schur decomposition of a real
022: *  matrix pair (A,B) using an orthogonal equivalence transformation
023: *
024: *                 (A, B) = Q * (A, B) * Z',
025: *
026: *  so that the diagonal block of (A, B) with row index IFST is moved
027: *  to row ILST.
028: *
029: *  (A, B) must be in generalized real Schur canonical form (as returned
030: *  by SGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2
031: *  diagonal blocks. B is upper triangular.
032: *
033: *  Optionally, the matrices Q and Z of generalized Schur vectors are
034: *  updated.
035: *
036: *         Q(in) * A(in) * Z(in)' = Q(out) * A(out) * Z(out)'
037: *         Q(in) * B(in) * Z(in)' = Q(out) * B(out) * Z(out)'
038: *
039: *
040: *  Arguments
041: *  =========
042: *
043: *  WANTQ   (input) LOGICAL
044: *          .TRUE. : update the left transformation matrix Q;
045: *          .FALSE.: do not update Q.
046: *
047: *  WANTZ   (input) LOGICAL
048: *          .TRUE. : update the right transformation matrix Z;
049: *          .FALSE.: do not update Z.
050: *
051: *  N       (input) INTEGER
052: *          The order of the matrices A and B. N >= 0.
053: *
054: *  A       (input/output) REAL array, dimension (LDA,N)
055: *          On entry, the matrix A in generalized real Schur canonical
056: *          form.
057: *          On exit, the updated matrix A, again in generalized
058: *          real Schur canonical form.
059: *
060: *  LDA     (input)  INTEGER
061: *          The leading dimension of the array A. LDA >= max(1,N).
062: *
063: *  B       (input/output) REAL array, dimension (LDB,N)
064: *          On entry, the matrix B in generalized real Schur canonical
065: *          form (A,B).
066: *          On exit, the updated matrix B, again in generalized
067: *          real Schur canonical form (A,B).
068: *
069: *  LDB     (input)  INTEGER
070: *          The leading dimension of the array B. LDB >= max(1,N).
071: *
072: *  Q       (input/output) REAL array, dimension (LDZ,N)
073: *          On entry, if WANTQ = .TRUE., the orthogonal matrix Q.
074: *          On exit, the updated matrix Q.
075: *          If WANTQ = .FALSE., Q is not referenced.
076: *
077: *  LDQ     (input) INTEGER
078: *          The leading dimension of the array Q. LDQ >= 1.
079: *          If WANTQ = .TRUE., LDQ >= N.
080: *
081: *  Z       (input/output) REAL array, dimension (LDZ,N)
082: *          On entry, if WANTZ = .TRUE., the orthogonal matrix Z.
083: *          On exit, the updated matrix Z.
084: *          If WANTZ = .FALSE., Z is not referenced.
085: *
086: *  LDZ     (input) INTEGER
087: *          The leading dimension of the array Z. LDZ >= 1.
088: *          If WANTZ = .TRUE., LDZ >= N.
089: *
090: *  IFST    (input/output) INTEGER
091: *  ILST    (input/output) INTEGER
092: *          Specify the reordering of the diagonal blocks of (A, B).
093: *          The block with row index IFST is moved to row ILST, by a
094: *          sequence of swapping between adjacent blocks.
095: *          On exit, if IFST pointed on entry to the second row of
096: *          a 2-by-2 block, it is changed to point to the first row;
097: *          ILST always points to the first row of the block in its
098: *          final position (which may differ from its input value by
099: *          +1 or -1). 1 <= IFST, ILST <= N.
100: *
101: *  WORK    (workspace/output) REAL 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:       REAL               ZERO
138:       PARAMETER          ( ZERO = 0.0E+0 )
139: *     ..
140: *     .. Local Scalars ..
141:       LOGICAL            LQUERY
142:       INTEGER            HERE, LWMIN, NBF, NBL, NBNEXT
143: *     ..
144: *     .. External Subroutines ..
145:       EXTERNAL           STGEX2, 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( 'STGEXC', -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 STGEX2( 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 STGEX2( 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 STGEX2( 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 STGEX2( 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 STGEX2( 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 STGEX2( 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 STGEX2( 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 STGEX2( 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 STGEX2( 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 STGEX2( 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 STGEX2( 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 STGEX2( 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 STGEXC
440: *
441:       END
442: