3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download STREXC + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/strexc.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/strexc.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/strexc.f">
21 * SUBROUTINE STREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
24 * .. Scalar Arguments ..
26 * INTEGER IFST, ILST, INFO, LDQ, LDT, N
28 * .. Array Arguments ..
29 * REAL Q( LDQ, * ), T( LDT, * ), WORK( * )
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
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.
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.
57 *> COMPQ is CHARACTER*1
58 *> = 'V': update the matrix Q of Schur vectors;
59 *> = 'N': do not update Q.
65 *> The order of the matrix T. N >= 0.
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.
80 *> The leading dimension of the array T. LDT >= max(1,N).
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.
95 *> The leading dimension of the array Q. LDQ >= max(1,N).
98 *> \param[in,out] IFST
103 *> \param[in,out] ILST
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.
119 *> WORK is REAL array, dimension (N)
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.
136 *> \author Univ. of Tennessee
137 *> \author Univ. of California Berkeley
138 *> \author Univ. of Colorado Denver
141 *> \date November 2011
143 *> \ingroup realOTHERcomputational
145 * =====================================================================
146 SUBROUTINE STREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK,
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..--
154 * .. Scalar Arguments ..
156 INTEGER IFST, ILST, INFO, LDQ, LDT, N
158 * .. Array Arguments ..
159 REAL Q( LDQ, * ), T( LDT, * ), WORK( * )
162 * =====================================================================
166 PARAMETER ( ZERO = 0.0E+0 )
168 * .. Local Scalars ..
170 INTEGER HERE, NBF, NBL, NBNEXT
172 * .. External Functions ..
176 * .. External Subroutines ..
177 EXTERNAL SLAEXC, XERBLA
179 * .. Intrinsic Functions ..
182 * .. Executable Statements ..
184 * Decode and test the input arguments.
187 WANTQ = LSAME( COMPQ, 'V' )
188 IF( .NOT.WANTQ .AND. .NOT.LSAME( COMPQ, 'N' ) ) THEN
190 ELSE IF( N.LT.0 ) THEN
192 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
194 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.MAX( 1, N ) ) ) THEN
196 ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN
198 ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN
202 CALL XERBLA( 'STREXC', -INFO )
206 * Quick return if possible
211 * Determine the first row of specified block
212 * and find out it is 1 by 1 or 2 by 2.
215 IF( T( IFST, IFST-1 ).NE.ZERO )
220 IF( T( IFST+1, IFST ).NE.ZERO )
224 * Determine the first row of the final block
225 * and find out it is 1 by 1 or 2 by 2.
228 IF( T( ILST, ILST-1 ).NE.ZERO )
233 IF( T( ILST+1, ILST ).NE.ZERO )
240 IF( IFST.LT.ILST ) THEN
244 IF( NBF.EQ.2 .AND. NBL.EQ.1 )
246 IF( NBF.EQ.1 .AND. NBL.EQ.2 )
253 * Swap block with next one below
255 IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
257 * Current block either 1 by 1 or 2 by 2
260 IF( HERE+NBF+1.LE.N ) THEN
261 IF( T( HERE+NBF+1, HERE+NBF ).NE.ZERO )
264 CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBF, NBNEXT,
272 * Test if 2 by 2 block breaks into two 1 by 1 blocks
275 IF( T( HERE+1, HERE ).EQ.ZERO )
281 * Current block consists of two 1 by 1 blocks each of which
282 * must be swapped individually
285 IF( HERE+3.LE.N ) THEN
286 IF( T( HERE+3, HERE+2 ).NE.ZERO )
289 CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, NBNEXT,
295 IF( NBNEXT.EQ.1 ) THEN
297 * Swap two 1 by 1 blocks, no problems possible
299 CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, NBNEXT,
304 * Recompute NBNEXT in case 2 by 2 split
306 IF( T( HERE+2, HERE+1 ).EQ.ZERO )
308 IF( NBNEXT.EQ.2 ) THEN
310 * 2 by 2 Block did not split
312 CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1,
313 $ NBNEXT, WORK, INFO )
321 * 2 by 2 Block did split
323 CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
325 CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, 1,
339 * Swap block with next one above
341 IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN
343 * Current block either 1 by 1 or 2 by 2
347 IF( T( HERE-1, HERE-2 ).NE.ZERO )
350 CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
358 * Test if 2 by 2 block breaks into two 1 by 1 blocks
361 IF( T( HERE+1, HERE ).EQ.ZERO )
367 * Current block consists of two 1 by 1 blocks each of which
368 * must be swapped individually
372 IF( T( HERE-1, HERE-2 ).NE.ZERO )
375 CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT,
381 IF( NBNEXT.EQ.1 ) THEN
383 * Swap two 1 by 1 blocks, no problems possible
385 CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBNEXT, 1,
390 * Recompute NBNEXT in case 2 by 2 split
392 IF( T( HERE, HERE-1 ).EQ.ZERO )
394 IF( NBNEXT.EQ.2 ) THEN
396 * 2 by 2 Block did not split
398 CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 2, 1,
407 * 2 by 2 Block did split
409 CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1,
411 CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 1, 1,