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