001:       SUBROUTINE DTREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
002:      $                   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:       CHARACTER          COMPQ
011:       INTEGER            IFST, ILST, INFO, LDQ, LDT, N
012: *     ..
013: *     .. Array Arguments ..
014:       DOUBLE PRECISION   Q( LDQ, * ), T( LDT, * ), WORK( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  DTREXC reorders the real Schur factorization of a real matrix
021: *  A = Q*T*Q**T, so that the diagonal block of T with row index IFST is
022: *  moved to row ILST.
023: *
024: *  The real Schur form T is reordered by an orthogonal similarity
025: *  transformation Z**T*T*Z, and optionally the matrix Q of Schur vectors
026: *  is updated by postmultiplying it with Z.
027: *
028: *  T must be in Schur canonical form (as returned by DHSEQR), that is,
029: *  block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
030: *  2-by-2 diagonal block has its diagonal elements equal and its
031: *  off-diagonal elements of opposite sign.
032: *
033: *  Arguments
034: *  =========
035: *
036: *  COMPQ   (input) CHARACTER*1
037: *          = 'V':  update the matrix Q of Schur vectors;
038: *          = 'N':  do not update Q.
039: *
040: *  N       (input) INTEGER
041: *          The order of the matrix T. N >= 0.
042: *
043: *  T       (input/output) DOUBLE PRECISION array, dimension (LDT,N)
044: *          On entry, the upper quasi-triangular matrix T, in Schur
045: *          Schur canonical form.
046: *          On exit, the reordered upper quasi-triangular matrix, again
047: *          in Schur canonical form.
048: *
049: *  LDT     (input) INTEGER
050: *          The leading dimension of the array T. LDT >= max(1,N).
051: *
052: *  Q       (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
053: *          On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
054: *          On exit, if COMPQ = 'V', Q has been postmultiplied by the
055: *          orthogonal transformation matrix Z which reorders T.
056: *          If COMPQ = 'N', Q is not referenced.
057: *
058: *  LDQ     (input) INTEGER
059: *          The leading dimension of the array Q.  LDQ >= max(1,N).
060: *
061: *  IFST    (input/output) INTEGER
062: *  ILST    (input/output) INTEGER
063: *          Specify the reordering of the diagonal blocks of T.
064: *          The block with row index IFST is moved to row ILST, by a
065: *          sequence of transpositions between adjacent blocks.
066: *          On exit, if IFST pointed on entry to the second row of a
067: *          2-by-2 block, it is changed to point to the first row; ILST
068: *          always points to the first row of the block in its final
069: *          position (which may differ from its input value by +1 or -1).
070: *          1 <= IFST <= N; 1 <= ILST <= N.
071: *
072: *  WORK    (workspace) DOUBLE PRECISION array, dimension (N)
073: *
074: *  INFO    (output) INTEGER
075: *          = 0:  successful exit
076: *          < 0:  if INFO = -i, the i-th argument had an illegal value
077: *          = 1:  two adjacent blocks were too close to swap (the problem
078: *                is very ill-conditioned); T may have been partially
079: *                reordered, and ILST points to the first row of the
080: *                current position of the block being moved.
081: *
082: *  =====================================================================
083: *
084: *     .. Parameters ..
085:       DOUBLE PRECISION   ZERO
086:       PARAMETER          ( ZERO = 0.0D+0 )
087: *     ..
088: *     .. Local Scalars ..
089:       LOGICAL            WANTQ
090:       INTEGER            HERE, NBF, NBL, NBNEXT
091: *     ..
092: *     .. External Functions ..
093:       LOGICAL            LSAME
094:       EXTERNAL           LSAME
095: *     ..
096: *     .. External Subroutines ..
097:       EXTERNAL           DLAEXC, XERBLA
098: *     ..
099: *     .. Intrinsic Functions ..
100:       INTRINSIC          MAX
101: *     ..
102: *     .. Executable Statements ..
103: *
104: *     Decode and test the input arguments.
105: *
106:       INFO = 0
107:       WANTQ = LSAME( COMPQ, 'V' )
108:       IF( .NOT.WANTQ .AND. .NOT.LSAME( COMPQ, 'N' ) ) THEN
109:          INFO = -1
110:       ELSE IF( N.LT.0 ) THEN
111:          INFO = -2
112:       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
113:          INFO = -4
114:       ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.MAX( 1, N ) ) ) THEN
115:          INFO = -6
116:       ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
117:          INFO = -7
118:       ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
119:          INFO = -8
120:       END IF
121:       IF( INFO.NE.0 ) THEN
122:          CALL XERBLA( 'DTREXC', -INFO )
123:          RETURN
124:       END IF
125: *
126: *     Quick return if possible
127: *
128:       IF( N.LE.1 )
129:      $   RETURN
130: *
131: *     Determine the first row of specified block
132: *     and find out it is 1 by 1 or 2 by 2.
133: *
134:       IF( IFST.GT.1 ) THEN
135:          IF( T( IFST, IFST-1 ).NE.ZERO )
136:      $      IFST = IFST - 1
137:       END IF
138:       NBF = 1
139:       IF( IFST.LT.N ) THEN
140:          IF( T( IFST+1, IFST ).NE.ZERO )
141:      $      NBF = 2
142:       END IF
143: *
144: *     Determine the first row of the final block
145: *     and find out it is 1 by 1 or 2 by 2.
146: *
147:       IF( ILST.GT.1 ) THEN
148:          IF( T( ILST, ILST-1 ).NE.ZERO )
149:      $      ILST = ILST - 1
150:       END IF
151:       NBL = 1
152:       IF( ILST.LT.N ) THEN
153:          IF( T( ILST+1, ILST ).NE.ZERO )
154:      $      NBL = 2
155:       END IF
156: *
157:       IF( IFST.EQ.ILST )
158:      $   RETURN
159: *
160:       IF( IFST.LT.ILST ) THEN
161: *
162: *        Update ILST
163: *
164:          IF( NBF.EQ.2 .AND. NBL.EQ.1 )
165:      $      ILST = ILST - 1
166:          IF( NBF.EQ.1 .AND. NBL.EQ.2 )
167:      $      ILST = ILST + 1
168: *
169:          HERE = IFST
170: *
171:    10    CONTINUE
172: *
173: *        Swap block with next one below
174: *
175:          IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
176: *
177: *           Current block either 1 by 1 or 2 by 2
178: *
179:             NBNEXT = 1
180:             IF( HERE+NBF+1.LE.N ) THEN
181:                IF( T( HERE+NBF+1, HERE+NBF ).NE.ZERO )
182:      $            NBNEXT = 2
183:             END IF
184:             CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBF, NBNEXT,
185:      $                   WORK, INFO )
186:             IF( INFO.NE.0 ) THEN
187:                ILST = HERE
188:                RETURN
189:             END IF
190:             HERE = HERE + NBNEXT
191: *
192: *           Test if 2 by 2 block breaks into two 1 by 1 blocks
193: *
194:             IF( NBF.EQ.2 ) THEN
195:                IF( T( HERE+1, HERE ).EQ.ZERO )
196:      $            NBF = 3
197:             END IF
198: *
199:          ELSE
200: *
201: *           Current block consists of two 1 by 1 blocks each of which
202: *           must be swapped individually
203: *
204:             NBNEXT = 1
205:             IF( HERE+3.LE.N ) THEN
206:                IF( T( HERE+3, HERE+2 ).NE.ZERO )
207:      $            NBNEXT = 2
208:             END IF
209:             CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, NBNEXT,
210:      $                   WORK, INFO )
211:             IF( INFO.NE.0 ) THEN
212:                ILST = HERE
213:                RETURN
214:             END IF
215:             IF( NBNEXT.EQ.1 ) THEN
216: *
217: *              Swap two 1 by 1 blocks, no problems possible
218: *
219:                CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, NBNEXT,
220:      $                      WORK, INFO )
221:                HERE = HERE + 1
222:             ELSE
223: *
224: *              Recompute NBNEXT in case 2 by 2 split
225: *
226:                IF( T( HERE+2, HERE+1 ).EQ.ZERO )
227:      $            NBNEXT = 1
228:                IF( NBNEXT.EQ.2 ) THEN
229: *
230: *                 2 by 2 Block did not split
231: *
232:                   CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1,
233:      $                         NBNEXT, WORK, INFO )
234:                   IF( INFO.NE.0 ) THEN
235:                      ILST = HERE
236:                      RETURN
237:                   END IF
238:                   HERE = HERE + 2
239:                ELSE
240: *
241: *                 2 by 2 Block did split
242: *
243:                   CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
244:      $                         WORK, INFO )
245:                   CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, 1,
246:      $                         WORK, INFO )
247:                   HERE = HERE + 2
248:                END IF
249:             END IF
250:          END IF
251:          IF( HERE.LT.ILST )
252:      $      GO TO 10
253: *
254:       ELSE
255: *
256:          HERE = IFST
257:    20    CONTINUE
258: *
259: *        Swap block with next one above
260: *
261:          IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
262: *
263: *           Current block either 1 by 1 or 2 by 2
264: *
265:             NBNEXT = 1
266:             IF( HERE.GE.3 ) THEN
267:                IF( T( HERE-1, HERE-2 ).NE.ZERO )
268:      $            NBNEXT = 2
269:             END IF
270:             CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
271:      $                   NBF, WORK, INFO )
272:             IF( INFO.NE.0 ) THEN
273:                ILST = HERE
274:                RETURN
275:             END IF
276:             HERE = HERE - NBNEXT
277: *
278: *           Test if 2 by 2 block breaks into two 1 by 1 blocks
279: *
280:             IF( NBF.EQ.2 ) THEN
281:                IF( T( HERE+1, HERE ).EQ.ZERO )
282:      $            NBF = 3
283:             END IF
284: *
285:          ELSE
286: *
287: *           Current block consists of two 1 by 1 blocks each of which
288: *           must be swapped individually
289: *
290:             NBNEXT = 1
291:             IF( HERE.GE.3 ) THEN
292:                IF( T( HERE-1, HERE-2 ).NE.ZERO )
293:      $            NBNEXT = 2
294:             END IF
295:             CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
296:      $                   1, WORK, INFO )
297:             IF( INFO.NE.0 ) THEN
298:                ILST = HERE
299:                RETURN
300:             END IF
301:             IF( NBNEXT.EQ.1 ) THEN
302: *
303: *              Swap two 1 by 1 blocks, no problems possible
304: *
305:                CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBNEXT, 1,
306:      $                      WORK, INFO )
307:                HERE = HERE - 1
308:             ELSE
309: *
310: *              Recompute NBNEXT in case 2 by 2 split
311: *
312:                IF( T( HERE, HERE-1 ).EQ.ZERO )
313:      $            NBNEXT = 1
314:                IF( NBNEXT.EQ.2 ) THEN
315: *
316: *                 2 by 2 Block did not split
317: *
318:                   CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 2, 1,
319:      $                         WORK, INFO )
320:                   IF( INFO.NE.0 ) THEN
321:                      ILST = HERE
322:                      RETURN
323:                   END IF
324:                   HERE = HERE - 2
325:                ELSE
326: *
327: *                 2 by 2 Block did split
328: *
329:                   CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
330:      $                         WORK, INFO )
331:                   CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 1, 1,
332:      $                         WORK, INFO )
333:                   HERE = HERE - 2
334:                END IF
335:             END IF
336:          END IF
337:          IF( HERE.GT.ILST )
338:      $      GO TO 20
339:       END IF
340:       ILST = HERE
341: *
342:       RETURN
343: *
344: *     End of DTREXC
345: *
346:       END
347: