1       SUBROUTINE DTREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
  2      $                   INFO )
  3 *
  4 *  -- LAPACK routine (version 3.2) --
  5 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  6 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  7 *     November 2006
  8 *
  9 *     .. Scalar Arguments ..
 10       CHARACTER          COMPQ
 11       INTEGER            IFST, ILST, INFO, LDQ, LDT, N
 12 *     ..
 13 *     .. Array Arguments ..
 14       DOUBLE PRECISION   Q( LDQ, * ), T( LDT, * ), WORK( * )
 15 *     ..
 16 *
 17 *  Purpose
 18 *  =======
 19 *
 20 *  DTREXC reorders the real Schur factorization of a real matrix
 21 *  A = Q*T*Q**T, so that the diagonal block of T with row index IFST is
 22 *  moved to row ILST.
 23 *
 24 *  The real Schur form T is reordered by an orthogonal similarity
 25 *  transformation Z**T*T*Z, and optionally the matrix Q of Schur vectors
 26 *  is updated by postmultiplying it with Z.
 27 *
 28 *  T must be in Schur canonical form (as returned by DHSEQR), that is,
 29 *  block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
 30 *  2-by-2 diagonal block has its diagonal elements equal and its
 31 *  off-diagonal elements of opposite sign.
 32 *
 33 *  Arguments
 34 *  =========
 35 *
 36 *  COMPQ   (input) CHARACTER*1
 37 *          = 'V':  update the matrix Q of Schur vectors;
 38 *          = 'N':  do not update Q.
 39 *
 40 *  N       (input) INTEGER
 41 *          The order of the matrix T. N >= 0.
 42 *
 43 *  T       (input/output) DOUBLE PRECISION array, dimension (LDT,N)
 44 *          On entry, the upper quasi-triangular matrix T, in Schur
 45 *          Schur canonical form.
 46 *          On exit, the reordered upper quasi-triangular matrix, again
 47 *          in Schur canonical form.
 48 *
 49 *  LDT     (input) INTEGER
 50 *          The leading dimension of the array T. LDT >= max(1,N).
 51 *
 52 *  Q       (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
 53 *          On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
 54 *          On exit, if COMPQ = 'V', Q has been postmultiplied by the
 55 *          orthogonal transformation matrix Z which reorders T.
 56 *          If COMPQ = 'N', Q is not referenced.
 57 *
 58 *  LDQ     (input) INTEGER
 59 *          The leading dimension of the array Q.  LDQ >= max(1,N).
 60 *
 61 *  IFST    (input/output) INTEGER
 62 *  ILST    (input/output) INTEGER
 63 *          Specify the reordering of the diagonal blocks of T.
 64 *          The block with row index IFST is moved to row ILST, by a
 65 *          sequence of transpositions between adjacent blocks.
 66 *          On exit, if IFST pointed on entry to the second row of a
 67 *          2-by-2 block, it is changed to point to the first row; ILST
 68 *          always points to the first row of the block in its final
 69 *          position (which may differ from its input value by +1 or -1).
 70 *          1 <= IFST <= N; 1 <= ILST <= N.
 71 *
 72 *  WORK    (workspace) DOUBLE PRECISION array, dimension (N)
 73 *
 74 *  INFO    (output) INTEGER
 75 *          = 0:  successful exit
 76 *          < 0:  if INFO = -i, the i-th argument had an illegal value
 77 *          = 1:  two adjacent blocks were too close to swap (the problem
 78 *                is very ill-conditioned); T may have been partially
 79 *                reordered, and ILST points to the first row of the
 80 *                current position of the block being moved.
 81 *
 82 *  =====================================================================
 83 *
 84 *     .. Parameters ..
 85       DOUBLE PRECISION   ZERO
 86       PARAMETER          ( ZERO = 0.0D+0 )
 87 *     ..
 88 *     .. Local Scalars ..
 89       LOGICAL            WANTQ
 90       INTEGER            HERE, NBF, NBL, NBNEXT
 91 *     ..
 92 *     .. External Functions ..
 93       LOGICAL            LSAME
 94       EXTERNAL           LSAME
 95 *     ..
 96 *     .. External Subroutines ..
 97       EXTERNAL           DLAEXC, XERBLA
 98 *     ..
 99 *     .. 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.MAX1, N ) ) THEN
113          INFO = -4
114       ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.MAX1, 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+11, 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, 11,
244      $                         WORK, INFO )
245                   CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+111,
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-121,
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, 11,
330      $                         WORK, INFO )
331                   CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-111,
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