1*> \brief \b DTREXC
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DTREXC + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrexc.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrexc.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrexc.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DTREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
22*                          INFO )
23*
24*       .. Scalar Arguments ..
25*       CHARACTER          COMPQ
26*       INTEGER            IFST, ILST, INFO, LDQ, LDT, N
27*       ..
28*       .. Array Arguments ..
29*       DOUBLE PRECISION   Q( LDQ, * ), T( LDT, * ), WORK( * )
30*       ..
31*
32*
33*> \par Purpose:
34*  =============
35*>
36*> \verbatim
37*>
38*> DTREXC reorders the real Schur factorization of a real matrix
39*> A = Q*T*Q**T, so that the diagonal block of T with row index IFST is
40*> moved to row ILST.
41*>
42*> The real Schur form T is reordered by an orthogonal similarity
43*> transformation Z**T*T*Z, and optionally the matrix Q of Schur vectors
44*> is updated by postmultiplying it with Z.
45*>
46*> T must be in Schur canonical form (as returned by DHSEQR), that is,
47*> block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each
48*> 2-by-2 diagonal block has its diagonal elements equal and its
49*> off-diagonal elements of opposite sign.
50*> \endverbatim
51*
52*  Arguments:
53*  ==========
54*
55*> \param[in] COMPQ
56*> \verbatim
57*>          COMPQ is CHARACTER*1
58*>          = 'V':  update the matrix Q of Schur vectors;
59*>          = 'N':  do not update Q.
60*> \endverbatim
61*>
62*> \param[in] N
63*> \verbatim
64*>          N is INTEGER
65*>          The order of the matrix T. N >= 0.
66*>          If N == 0 arguments ILST and IFST may be any value.
67*> \endverbatim
68*>
69*> \param[in,out] T
70*> \verbatim
71*>          T is DOUBLE PRECISION array, dimension (LDT,N)
72*>          On entry, the upper quasi-triangular matrix T, in Schur
73*>          Schur canonical form.
74*>          On exit, the reordered upper quasi-triangular matrix, again
75*>          in Schur canonical form.
76*> \endverbatim
77*>
78*> \param[in] LDT
79*> \verbatim
80*>          LDT is INTEGER
81*>          The leading dimension of the array T. LDT >= max(1,N).
82*> \endverbatim
83*>
84*> \param[in,out] Q
85*> \verbatim
86*>          Q is DOUBLE PRECISION array, dimension (LDQ,N)
87*>          On entry, if COMPQ = 'V', the matrix Q of Schur vectors.
88*>          On exit, if COMPQ = 'V', Q has been postmultiplied by the
89*>          orthogonal transformation matrix Z which reorders T.
90*>          If COMPQ = 'N', Q is not referenced.
91*> \endverbatim
92*>
93*> \param[in] LDQ
94*> \verbatim
95*>          LDQ is INTEGER
96*>          The leading dimension of the array Q.  LDQ >= 1, and if
97*>          COMPQ = 'V', LDQ >= max(1,N).
98*> \endverbatim
99*>
100*> \param[in,out] IFST
101*> \verbatim
102*>          IFST is INTEGER
103*> \endverbatim
104*>
105*> \param[in,out] ILST
106*> \verbatim
107*>          ILST is INTEGER
108*>
109*>          Specify the reordering of the diagonal blocks of T.
110*>          The block with row index IFST is moved to row ILST, by a
111*>          sequence of transpositions between adjacent blocks.
112*>          On exit, if IFST pointed on entry to the second row of a
113*>          2-by-2 block, it is changed to point to the first row; ILST
114*>          always points to the first row of the block in its final
115*>          position (which may differ from its input value by +1 or -1).
116*>          1 <= IFST <= N; 1 <= ILST <= N.
117*> \endverbatim
118*>
119*> \param[out] WORK
120*> \verbatim
121*>          WORK is DOUBLE PRECISION array, dimension (N)
122*> \endverbatim
123*>
124*> \param[out] INFO
125*> \verbatim
126*>          INFO is INTEGER
127*>          = 0:  successful exit
128*>          < 0:  if INFO = -i, the i-th argument had an illegal value
129*>          = 1:  two adjacent blocks were too close to swap (the problem
130*>                is very ill-conditioned); T may have been partially
131*>                reordered, and ILST points to the first row of the
132*>                current position of the block being moved.
133*> \endverbatim
134*
135*  Authors:
136*  ========
137*
138*> \author Univ. of Tennessee
139*> \author Univ. of California Berkeley
140*> \author Univ. of Colorado Denver
141*> \author NAG Ltd.
142*
143*> \ingroup doubleOTHERcomputational
144*
145*  =====================================================================
146      SUBROUTINE DTREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
147     $                   INFO )
148*
149*  -- LAPACK computational routine --
150*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
151*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
152*
153*     .. Scalar Arguments ..
154      CHARACTER          COMPQ
155      INTEGER            IFST, ILST, INFO, LDQ, LDT, N
156*     ..
157*     .. Array Arguments ..
158      DOUBLE PRECISION   Q( LDQ, * ), T( LDT, * ), WORK( * )
159*     ..
160*
161*  =====================================================================
162*
163*     .. Parameters ..
164      DOUBLE PRECISION   ZERO
165      PARAMETER          ( ZERO = 0.0D+0 )
166*     ..
167*     .. Local Scalars ..
168      LOGICAL            WANTQ
169      INTEGER            HERE, NBF, NBL, NBNEXT
170*     ..
171*     .. External Functions ..
172      LOGICAL            LSAME
173      EXTERNAL           LSAME
174*     ..
175*     .. External Subroutines ..
176      EXTERNAL           DLAEXC, XERBLA
177*     ..
178*     .. Intrinsic Functions ..
179      INTRINSIC          MAX
180*     ..
181*     .. Executable Statements ..
182*
183*     Decode and test the input arguments.
184*
185      INFO = 0
186      WANTQ = LSAME( COMPQ, 'V' )
187      IF( .NOT.WANTQ .AND. .NOT.LSAME( COMPQ, 'N' ) ) THEN
188         INFO = -1
189      ELSE IF( N.LT.0 ) THEN
190         INFO = -2
191      ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
192         INFO = -4
193      ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.MAX( 1, N ) ) ) THEN
194         INFO = -6
195      ELSE IF(( IFST.LT.1 .OR. IFST.GT.N ).AND.( N.GT.0 )) THEN
196         INFO = -7
197      ELSE IF(( ILST.LT.1 .OR. ILST.GT.N ).AND.( N.GT.0 )) THEN
198         INFO = -8
199      END IF
200      IF( INFO.NE.0 ) THEN
201         CALL XERBLA( 'DTREXC', -INFO )
202         RETURN
203      END IF
204*
205*     Quick return if possible
206*
207      IF( N.LE.1 )
208     $   RETURN
209*
210*     Determine the first row of specified block
211*     and find out it is 1 by 1 or 2 by 2.
212*
213      IF( IFST.GT.1 ) THEN
214         IF( T( IFST, IFST-1 ).NE.ZERO )
215     $      IFST = IFST - 1
216      END IF
217      NBF = 1
218      IF( IFST.LT.N ) THEN
219         IF( T( IFST+1, IFST ).NE.ZERO )
220     $      NBF = 2
221      END IF
222*
223*     Determine the first row of the final block
224*     and find out it is 1 by 1 or 2 by 2.
225*
226      IF( ILST.GT.1 ) THEN
227         IF( T( ILST, ILST-1 ).NE.ZERO )
228     $      ILST = ILST - 1
229      END IF
230      NBL = 1
231      IF( ILST.LT.N ) THEN
232         IF( T( ILST+1, ILST ).NE.ZERO )
233     $      NBL = 2
234      END IF
235*
236      IF( IFST.EQ.ILST )
237     $   RETURN
238*
239      IF( IFST.LT.ILST ) THEN
240*
241*        Update ILST
242*
243         IF( NBF.EQ.2 .AND. NBL.EQ.1 )
244     $      ILST = ILST - 1
245         IF( NBF.EQ.1 .AND. NBL.EQ.2 )
246     $      ILST = ILST + 1
247*
248         HERE = IFST
249*
250   10    CONTINUE
251*
252*        Swap block with next one below
253*
254         IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
255*
256*           Current block either 1 by 1 or 2 by 2
257*
258            NBNEXT = 1
259            IF( HERE+NBF+1.LE.N ) THEN
260               IF( T( HERE+NBF+1, HERE+NBF ).NE.ZERO )
261     $            NBNEXT = 2
262            END IF
263            CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBF, NBNEXT,
264     $                   WORK, INFO )
265            IF( INFO.NE.0 ) THEN
266               ILST = HERE
267               RETURN
268            END IF
269            HERE = HERE + NBNEXT
270*
271*           Test if 2 by 2 block breaks into two 1 by 1 blocks
272*
273            IF( NBF.EQ.2 ) THEN
274               IF( T( HERE+1, HERE ).EQ.ZERO )
275     $            NBF = 3
276            END IF
277*
278         ELSE
279*
280*           Current block consists of two 1 by 1 blocks each of which
281*           must be swapped individually
282*
283            NBNEXT = 1
284            IF( HERE+3.LE.N ) THEN
285               IF( T( HERE+3, HERE+2 ).NE.ZERO )
286     $            NBNEXT = 2
287            END IF
288            CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, NBNEXT,
289     $                   WORK, INFO )
290            IF( INFO.NE.0 ) THEN
291               ILST = HERE
292               RETURN
293            END IF
294            IF( NBNEXT.EQ.1 ) THEN
295*
296*              Swap two 1 by 1 blocks, no problems possible
297*
298               CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, NBNEXT,
299     $                      WORK, INFO )
300               HERE = HERE + 1
301            ELSE
302*
303*              Recompute NBNEXT in case 2 by 2 split
304*
305               IF( T( HERE+2, HERE+1 ).EQ.ZERO )
306     $            NBNEXT = 1
307               IF( NBNEXT.EQ.2 ) THEN
308*
309*                 2 by 2 Block did not split
310*
311                  CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1,
312     $                         NBNEXT, WORK, INFO )
313                  IF( INFO.NE.0 ) THEN
314                     ILST = HERE
315                     RETURN
316                  END IF
317                  HERE = HERE + 2
318               ELSE
319*
320*                 2 by 2 Block did split
321*
322                  CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
323     $                         WORK, INFO )
324                  CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, 1,
325     $                         WORK, INFO )
326                  HERE = HERE + 2
327               END IF
328            END IF
329         END IF
330         IF( HERE.LT.ILST )
331     $      GO TO 10
332*
333      ELSE
334*
335         HERE = IFST
336   20    CONTINUE
337*
338*        Swap block with next one above
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( T( HERE-1, HERE-2 ).NE.ZERO )
347     $            NBNEXT = 2
348            END IF
349            CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
350     $                   NBF, WORK, INFO )
351            IF( INFO.NE.0 ) THEN
352               ILST = HERE
353               RETURN
354            END IF
355            HERE = HERE - NBNEXT
356*
357*           Test if 2 by 2 block breaks into two 1 by 1 blocks
358*
359            IF( NBF.EQ.2 ) THEN
360               IF( T( HERE+1, HERE ).EQ.ZERO )
361     $            NBF = 3
362            END IF
363*
364         ELSE
365*
366*           Current block consists of two 1 by 1 blocks each of which
367*           must be swapped individually
368*
369            NBNEXT = 1
370            IF( HERE.GE.3 ) THEN
371               IF( T( HERE-1, HERE-2 ).NE.ZERO )
372     $            NBNEXT = 2
373            END IF
374            CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
375     $                   1, WORK, INFO )
376            IF( INFO.NE.0 ) THEN
377               ILST = HERE
378               RETURN
379            END IF
380            IF( NBNEXT.EQ.1 ) THEN
381*
382*              Swap two 1 by 1 blocks, no problems possible
383*
384               CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBNEXT, 1,
385     $                      WORK, INFO )
386               HERE = HERE - 1
387            ELSE
388*
389*              Recompute NBNEXT in case 2 by 2 split
390*
391               IF( T( HERE, HERE-1 ).EQ.ZERO )
392     $            NBNEXT = 1
393               IF( NBNEXT.EQ.2 ) THEN
394*
395*                 2 by 2 Block did not split
396*
397                  CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 2, 1,
398     $                         WORK, INFO )
399                  IF( INFO.NE.0 ) THEN
400                     ILST = HERE
401                     RETURN
402                  END IF
403                  HERE = HERE - 2
404               ELSE
405*
406*                 2 by 2 Block did split
407*
408                  CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
409     $                         WORK, INFO )
410                  CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 1, 1,
411     $                         WORK, INFO )
412                  HERE = HERE - 2
413               END IF
414            END IF
415         END IF
416         IF( HERE.GT.ILST )
417     $      GO TO 20
418      END IF
419      ILST = HERE
420*
421      RETURN
422*
423*     End of DTREXC
424*
425      END
426