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