|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE STREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK, 00002 $ INFO ) 00003 * 00004 * -- LAPACK routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER COMPQ 00011 INTEGER IFST, ILST, INFO, LDQ, LDT, N 00012 * .. 00013 * .. Array Arguments .. 00014 REAL Q( LDQ, * ), T( LDT, * ), WORK( * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * STREXC reorders the real Schur factorization of a real matrix 00021 * A = Q*T*Q**T, so that the diagonal block of T with row index IFST is 00022 * moved to row ILST. 00023 * 00024 * The real Schur form T is reordered by an orthogonal similarity 00025 * transformation Z**T*T*Z, and optionally the matrix Q of Schur vectors 00026 * is updated by postmultiplying it with Z. 00027 * 00028 * T must be in Schur canonical form (as returned by SHSEQR), that is, 00029 * block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each 00030 * 2-by-2 diagonal block has its diagonal elements equal and its 00031 * off-diagonal elements of opposite sign. 00032 * 00033 * Arguments 00034 * ========= 00035 * 00036 * COMPQ (input) CHARACTER*1 00037 * = 'V': update the matrix Q of Schur vectors; 00038 * = 'N': do not update Q. 00039 * 00040 * N (input) INTEGER 00041 * The order of the matrix T. N >= 0. 00042 * 00043 * T (input/output) REAL array, dimension (LDT,N) 00044 * On entry, the upper quasi-triangular matrix T, in Schur 00045 * Schur canonical form. 00046 * On exit, the reordered upper quasi-triangular matrix, again 00047 * in Schur canonical form. 00048 * 00049 * LDT (input) INTEGER 00050 * The leading dimension of the array T. LDT >= max(1,N). 00051 * 00052 * Q (input/output) REAL array, dimension (LDQ,N) 00053 * On entry, if COMPQ = 'V', the matrix Q of Schur vectors. 00054 * On exit, if COMPQ = 'V', Q has been postmultiplied by the 00055 * orthogonal transformation matrix Z which reorders T. 00056 * If COMPQ = 'N', Q is not referenced. 00057 * 00058 * LDQ (input) INTEGER 00059 * The leading dimension of the array Q. LDQ >= max(1,N). 00060 * 00061 * IFST (input/output) INTEGER 00062 * ILST (input/output) INTEGER 00063 * Specify the reordering of the diagonal blocks of T. 00064 * The block with row index IFST is moved to row ILST, by a 00065 * sequence of transpositions between adjacent blocks. 00066 * On exit, if IFST pointed on entry to the second row of a 00067 * 2-by-2 block, it is changed to point to the first row; ILST 00068 * always points to the first row of the block in its final 00069 * position (which may differ from its input value by +1 or -1). 00070 * 1 <= IFST <= N; 1 <= ILST <= N. 00071 * 00072 * WORK (workspace) REAL array, dimension (N) 00073 * 00074 * INFO (output) INTEGER 00075 * = 0: successful exit 00076 * < 0: if INFO = -i, the i-th argument had an illegal value 00077 * = 1: two adjacent blocks were too close to swap (the problem 00078 * is very ill-conditioned); T may have been partially 00079 * reordered, and ILST points to the first row of the 00080 * current position of the block being moved. 00081 * 00082 * ===================================================================== 00083 * 00084 * .. Parameters .. 00085 REAL ZERO 00086 PARAMETER ( ZERO = 0.0E+0 ) 00087 * .. 00088 * .. Local Scalars .. 00089 LOGICAL WANTQ 00090 INTEGER HERE, NBF, NBL, NBNEXT 00091 * .. 00092 * .. External Functions .. 00093 LOGICAL LSAME 00094 EXTERNAL LSAME 00095 * .. 00096 * .. External Subroutines .. 00097 EXTERNAL SLAEXC, XERBLA 00098 * .. 00099 * .. Intrinsic Functions .. 00100 INTRINSIC MAX 00101 * .. 00102 * .. Executable Statements .. 00103 * 00104 * Decode and test the input arguments. 00105 * 00106 INFO = 0 00107 WANTQ = LSAME( COMPQ, 'V' ) 00108 IF( .NOT.WANTQ .AND. .NOT.LSAME( COMPQ, 'N' ) ) THEN 00109 INFO = -1 00110 ELSE IF( N.LT.0 ) THEN 00111 INFO = -2 00112 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN 00113 INFO = -4 00114 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.MAX( 1, N ) ) ) THEN 00115 INFO = -6 00116 ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN 00117 INFO = -7 00118 ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN 00119 INFO = -8 00120 END IF 00121 IF( INFO.NE.0 ) THEN 00122 CALL XERBLA( 'STREXC', -INFO ) 00123 RETURN 00124 END IF 00125 * 00126 * Quick return if possible 00127 * 00128 IF( N.LE.1 ) 00129 $ RETURN 00130 * 00131 * Determine the first row of specified block 00132 * and find out it is 1 by 1 or 2 by 2. 00133 * 00134 IF( IFST.GT.1 ) THEN 00135 IF( T( IFST, IFST-1 ).NE.ZERO ) 00136 $ IFST = IFST - 1 00137 END IF 00138 NBF = 1 00139 IF( IFST.LT.N ) THEN 00140 IF( T( IFST+1, IFST ).NE.ZERO ) 00141 $ NBF = 2 00142 END IF 00143 * 00144 * Determine the first row of the final block 00145 * and find out it is 1 by 1 or 2 by 2. 00146 * 00147 IF( ILST.GT.1 ) THEN 00148 IF( T( ILST, ILST-1 ).NE.ZERO ) 00149 $ ILST = ILST - 1 00150 END IF 00151 NBL = 1 00152 IF( ILST.LT.N ) THEN 00153 IF( T( ILST+1, ILST ).NE.ZERO ) 00154 $ NBL = 2 00155 END IF 00156 * 00157 IF( IFST.EQ.ILST ) 00158 $ RETURN 00159 * 00160 IF( IFST.LT.ILST ) THEN 00161 * 00162 * Update ILST 00163 * 00164 IF( NBF.EQ.2 .AND. NBL.EQ.1 ) 00165 $ ILST = ILST - 1 00166 IF( NBF.EQ.1 .AND. NBL.EQ.2 ) 00167 $ ILST = ILST + 1 00168 * 00169 HERE = IFST 00170 * 00171 10 CONTINUE 00172 * 00173 * Swap block with next one below 00174 * 00175 IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN 00176 * 00177 * Current block either 1 by 1 or 2 by 2 00178 * 00179 NBNEXT = 1 00180 IF( HERE+NBF+1.LE.N ) THEN 00181 IF( T( HERE+NBF+1, HERE+NBF ).NE.ZERO ) 00182 $ NBNEXT = 2 00183 END IF 00184 CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBF, NBNEXT, 00185 $ WORK, INFO ) 00186 IF( INFO.NE.0 ) THEN 00187 ILST = HERE 00188 RETURN 00189 END IF 00190 HERE = HERE + NBNEXT 00191 * 00192 * Test if 2 by 2 block breaks into two 1 by 1 blocks 00193 * 00194 IF( NBF.EQ.2 ) THEN 00195 IF( T( HERE+1, HERE ).EQ.ZERO ) 00196 $ NBF = 3 00197 END IF 00198 * 00199 ELSE 00200 * 00201 * Current block consists of two 1 by 1 blocks each of which 00202 * must be swapped individually 00203 * 00204 NBNEXT = 1 00205 IF( HERE+3.LE.N ) THEN 00206 IF( T( HERE+3, HERE+2 ).NE.ZERO ) 00207 $ NBNEXT = 2 00208 END IF 00209 CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, NBNEXT, 00210 $ WORK, INFO ) 00211 IF( INFO.NE.0 ) THEN 00212 ILST = HERE 00213 RETURN 00214 END IF 00215 IF( NBNEXT.EQ.1 ) THEN 00216 * 00217 * Swap two 1 by 1 blocks, no problems possible 00218 * 00219 CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, NBNEXT, 00220 $ WORK, INFO ) 00221 HERE = HERE + 1 00222 ELSE 00223 * 00224 * Recompute NBNEXT in case 2 by 2 split 00225 * 00226 IF( T( HERE+2, HERE+1 ).EQ.ZERO ) 00227 $ NBNEXT = 1 00228 IF( NBNEXT.EQ.2 ) THEN 00229 * 00230 * 2 by 2 Block did not split 00231 * 00232 CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 00233 $ NBNEXT, WORK, INFO ) 00234 IF( INFO.NE.0 ) THEN 00235 ILST = HERE 00236 RETURN 00237 END IF 00238 HERE = HERE + 2 00239 ELSE 00240 * 00241 * 2 by 2 Block did split 00242 * 00243 CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1, 00244 $ WORK, INFO ) 00245 CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, 1, 00246 $ WORK, INFO ) 00247 HERE = HERE + 2 00248 END IF 00249 END IF 00250 END IF 00251 IF( HERE.LT.ILST ) 00252 $ GO TO 10 00253 * 00254 ELSE 00255 * 00256 HERE = IFST 00257 20 CONTINUE 00258 * 00259 * Swap block with next one above 00260 * 00261 IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN 00262 * 00263 * Current block either 1 by 1 or 2 by 2 00264 * 00265 NBNEXT = 1 00266 IF( HERE.GE.3 ) THEN 00267 IF( T( HERE-1, HERE-2 ).NE.ZERO ) 00268 $ NBNEXT = 2 00269 END IF 00270 CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT, 00271 $ NBF, WORK, INFO ) 00272 IF( INFO.NE.0 ) THEN 00273 ILST = HERE 00274 RETURN 00275 END IF 00276 HERE = HERE - NBNEXT 00277 * 00278 * Test if 2 by 2 block breaks into two 1 by 1 blocks 00279 * 00280 IF( NBF.EQ.2 ) THEN 00281 IF( T( HERE+1, HERE ).EQ.ZERO ) 00282 $ NBF = 3 00283 END IF 00284 * 00285 ELSE 00286 * 00287 * Current block consists of two 1 by 1 blocks each of which 00288 * must be swapped individually 00289 * 00290 NBNEXT = 1 00291 IF( HERE.GE.3 ) THEN 00292 IF( T( HERE-1, HERE-2 ).NE.ZERO ) 00293 $ NBNEXT = 2 00294 END IF 00295 CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT, 00296 $ 1, WORK, INFO ) 00297 IF( INFO.NE.0 ) THEN 00298 ILST = HERE 00299 RETURN 00300 END IF 00301 IF( NBNEXT.EQ.1 ) THEN 00302 * 00303 * Swap two 1 by 1 blocks, no problems possible 00304 * 00305 CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBNEXT, 1, 00306 $ WORK, INFO ) 00307 HERE = HERE - 1 00308 ELSE 00309 * 00310 * Recompute NBNEXT in case 2 by 2 split 00311 * 00312 IF( T( HERE, HERE-1 ).EQ.ZERO ) 00313 $ NBNEXT = 1 00314 IF( NBNEXT.EQ.2 ) THEN 00315 * 00316 * 2 by 2 Block did not split 00317 * 00318 CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 2, 1, 00319 $ WORK, INFO ) 00320 IF( INFO.NE.0 ) THEN 00321 ILST = HERE 00322 RETURN 00323 END IF 00324 HERE = HERE - 2 00325 ELSE 00326 * 00327 * 2 by 2 Block did split 00328 * 00329 CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1, 00330 $ WORK, INFO ) 00331 CALL SLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 1, 1, 00332 $ WORK, INFO ) 00333 HERE = HERE - 2 00334 END IF 00335 END IF 00336 END IF 00337 IF( HERE.GT.ILST ) 00338 $ GO TO 20 00339 END IF 00340 ILST = HERE 00341 * 00342 RETURN 00343 * 00344 * End of STREXC 00345 * 00346 END
1.7.3 
