|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE STGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 00002 $ LDZ, IFST, ILST, WORK, LWORK, 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 LOGICAL WANTQ, WANTZ 00011 INTEGER IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N 00012 * .. 00013 * .. Array Arguments .. 00014 REAL A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 00015 $ WORK( * ), Z( LDZ, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * STGEXC reorders the generalized real Schur decomposition of a real 00022 * matrix pair (A,B) using an orthogonal equivalence transformation 00023 * 00024 * (A, B) = Q * (A, B) * Z', 00025 * 00026 * so that the diagonal block of (A, B) with row index IFST is moved 00027 * to row ILST. 00028 * 00029 * (A, B) must be in generalized real Schur canonical form (as returned 00030 * by SGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2 00031 * diagonal blocks. B is upper triangular. 00032 * 00033 * Optionally, the matrices Q and Z of generalized Schur vectors are 00034 * updated. 00035 * 00036 * Q(in) * A(in) * Z(in)' = Q(out) * A(out) * Z(out)' 00037 * Q(in) * B(in) * Z(in)' = Q(out) * B(out) * Z(out)' 00038 * 00039 * 00040 * Arguments 00041 * ========= 00042 * 00043 * WANTQ (input) LOGICAL 00044 * .TRUE. : update the left transformation matrix Q; 00045 * .FALSE.: do not update Q. 00046 * 00047 * WANTZ (input) LOGICAL 00048 * .TRUE. : update the right transformation matrix Z; 00049 * .FALSE.: do not update Z. 00050 * 00051 * N (input) INTEGER 00052 * The order of the matrices A and B. N >= 0. 00053 * 00054 * A (input/output) REAL array, dimension (LDA,N) 00055 * On entry, the matrix A in generalized real Schur canonical 00056 * form. 00057 * On exit, the updated matrix A, again in generalized 00058 * real Schur canonical form. 00059 * 00060 * LDA (input) INTEGER 00061 * The leading dimension of the array A. LDA >= max(1,N). 00062 * 00063 * B (input/output) REAL array, dimension (LDB,N) 00064 * On entry, the matrix B in generalized real Schur canonical 00065 * form (A,B). 00066 * On exit, the updated matrix B, again in generalized 00067 * real Schur canonical form (A,B). 00068 * 00069 * LDB (input) INTEGER 00070 * The leading dimension of the array B. LDB >= max(1,N). 00071 * 00072 * Q (input/output) REAL array, dimension (LDZ,N) 00073 * On entry, if WANTQ = .TRUE., the orthogonal matrix Q. 00074 * On exit, the updated matrix Q. 00075 * If WANTQ = .FALSE., Q is not referenced. 00076 * 00077 * LDQ (input) INTEGER 00078 * The leading dimension of the array Q. LDQ >= 1. 00079 * If WANTQ = .TRUE., LDQ >= N. 00080 * 00081 * Z (input/output) REAL array, dimension (LDZ,N) 00082 * On entry, if WANTZ = .TRUE., the orthogonal matrix Z. 00083 * On exit, the updated matrix Z. 00084 * If WANTZ = .FALSE., Z is not referenced. 00085 * 00086 * LDZ (input) INTEGER 00087 * The leading dimension of the array Z. LDZ >= 1. 00088 * If WANTZ = .TRUE., LDZ >= N. 00089 * 00090 * IFST (input/output) INTEGER 00091 * ILST (input/output) INTEGER 00092 * Specify the reordering of the diagonal blocks of (A, B). 00093 * The block with row index IFST is moved to row ILST, by a 00094 * sequence of swapping between adjacent blocks. 00095 * On exit, if IFST pointed on entry to the second row of 00096 * a 2-by-2 block, it is changed to point to the first row; 00097 * ILST always points to the first row of the block in its 00098 * final position (which may differ from its input value by 00099 * +1 or -1). 1 <= IFST, ILST <= N. 00100 * 00101 * WORK (workspace/output) REAL array, dimension (MAX(1,LWORK)) 00102 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00103 * 00104 * LWORK (input) INTEGER 00105 * The dimension of the array WORK. 00106 * LWORK >= 1 when N <= 1, otherwise LWORK >= 4*N + 16. 00107 * 00108 * If LWORK = -1, then a workspace query is assumed; the routine 00109 * only calculates the optimal size of the WORK array, returns 00110 * this value as the first entry of the WORK array, and no error 00111 * message related to LWORK is issued by XERBLA. 00112 * 00113 * INFO (output) INTEGER 00114 * =0: successful exit. 00115 * <0: if INFO = -i, the i-th argument had an illegal value. 00116 * =1: The transformed matrix pair (A, B) would be too far 00117 * from generalized Schur form; the problem is ill- 00118 * conditioned. (A, B) may have been partially reordered, 00119 * and ILST points to the first row of the current 00120 * position of the block being moved. 00121 * 00122 * Further Details 00123 * =============== 00124 * 00125 * Based on contributions by 00126 * Bo Kagstrom and Peter Poromaa, Department of Computing Science, 00127 * Umea University, S-901 87 Umea, Sweden. 00128 * 00129 * [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the 00130 * Generalized Real Schur Form of a Regular Matrix Pair (A, B), in 00131 * M.S. Moonen et al (eds), Linear Algebra for Large Scale and 00132 * Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218. 00133 * 00134 * ===================================================================== 00135 * 00136 * .. Parameters .. 00137 REAL ZERO 00138 PARAMETER ( ZERO = 0.0E+0 ) 00139 * .. 00140 * .. Local Scalars .. 00141 LOGICAL LQUERY 00142 INTEGER HERE, LWMIN, NBF, NBL, NBNEXT 00143 * .. 00144 * .. External Subroutines .. 00145 EXTERNAL STGEX2, XERBLA 00146 * .. 00147 * .. Intrinsic Functions .. 00148 INTRINSIC MAX 00149 * .. 00150 * .. Executable Statements .. 00151 * 00152 * Decode and test input arguments. 00153 * 00154 INFO = 0 00155 LQUERY = ( LWORK.EQ.-1 ) 00156 IF( N.LT.0 ) THEN 00157 INFO = -3 00158 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00159 INFO = -5 00160 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00161 INFO = -7 00162 ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. ( LDQ.LT.MAX( 1, N ) ) ) THEN 00163 INFO = -9 00164 ELSE IF( LDZ.LT.1 .OR. WANTZ .AND. ( LDZ.LT.MAX( 1, N ) ) ) THEN 00165 INFO = -11 00166 ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN 00167 INFO = -12 00168 ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN 00169 INFO = -13 00170 END IF 00171 * 00172 IF( INFO.EQ.0 ) THEN 00173 IF( N.LE.1 ) THEN 00174 LWMIN = 1 00175 ELSE 00176 LWMIN = 4*N + 16 00177 END IF 00178 WORK(1) = LWMIN 00179 * 00180 IF (LWORK.LT.LWMIN .AND. .NOT.LQUERY) THEN 00181 INFO = -15 00182 END IF 00183 END IF 00184 * 00185 IF( INFO.NE.0 ) THEN 00186 CALL XERBLA( 'STGEXC', -INFO ) 00187 RETURN 00188 ELSE IF( LQUERY ) THEN 00189 RETURN 00190 END IF 00191 * 00192 * Quick return if possible 00193 * 00194 IF( N.LE.1 ) 00195 $ RETURN 00196 * 00197 * Determine the first row of the specified block and find out 00198 * if it is 1-by-1 or 2-by-2. 00199 * 00200 IF( IFST.GT.1 ) THEN 00201 IF( A( IFST, IFST-1 ).NE.ZERO ) 00202 $ IFST = IFST - 1 00203 END IF 00204 NBF = 1 00205 IF( IFST.LT.N ) THEN 00206 IF( A( IFST+1, IFST ).NE.ZERO ) 00207 $ NBF = 2 00208 END IF 00209 * 00210 * Determine the first row of the final block 00211 * and find out if it is 1-by-1 or 2-by-2. 00212 * 00213 IF( ILST.GT.1 ) THEN 00214 IF( A( ILST, ILST-1 ).NE.ZERO ) 00215 $ ILST = ILST - 1 00216 END IF 00217 NBL = 1 00218 IF( ILST.LT.N ) THEN 00219 IF( A( ILST+1, ILST ).NE.ZERO ) 00220 $ NBL = 2 00221 END IF 00222 IF( IFST.EQ.ILST ) 00223 $ RETURN 00224 * 00225 IF( IFST.LT.ILST ) THEN 00226 * 00227 * Update ILST. 00228 * 00229 IF( NBF.EQ.2 .AND. NBL.EQ.1 ) 00230 $ ILST = ILST - 1 00231 IF( NBF.EQ.1 .AND. NBL.EQ.2 ) 00232 $ ILST = ILST + 1 00233 * 00234 HERE = IFST 00235 * 00236 10 CONTINUE 00237 * 00238 * Swap with next one below. 00239 * 00240 IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN 00241 * 00242 * Current block either 1-by-1 or 2-by-2. 00243 * 00244 NBNEXT = 1 00245 IF( HERE+NBF+1.LE.N ) THEN 00246 IF( A( HERE+NBF+1, HERE+NBF ).NE.ZERO ) 00247 $ NBNEXT = 2 00248 END IF 00249 CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 00250 $ LDZ, HERE, NBF, NBNEXT, WORK, LWORK, INFO ) 00251 IF( INFO.NE.0 ) THEN 00252 ILST = HERE 00253 RETURN 00254 END IF 00255 HERE = HERE + NBNEXT 00256 * 00257 * Test if 2-by-2 block breaks into two 1-by-1 blocks. 00258 * 00259 IF( NBF.EQ.2 ) THEN 00260 IF( A( HERE+1, HERE ).EQ.ZERO ) 00261 $ NBF = 3 00262 END IF 00263 * 00264 ELSE 00265 * 00266 * Current block consists of two 1-by-1 blocks, each of which 00267 * must be swapped individually. 00268 * 00269 NBNEXT = 1 00270 IF( HERE+3.LE.N ) THEN 00271 IF( A( HERE+3, HERE+2 ).NE.ZERO ) 00272 $ NBNEXT = 2 00273 END IF 00274 CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 00275 $ LDZ, HERE+1, 1, NBNEXT, WORK, LWORK, INFO ) 00276 IF( INFO.NE.0 ) THEN 00277 ILST = HERE 00278 RETURN 00279 END IF 00280 IF( NBNEXT.EQ.1 ) THEN 00281 * 00282 * Swap two 1-by-1 blocks. 00283 * 00284 CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 00285 $ LDZ, HERE, 1, 1, WORK, LWORK, INFO ) 00286 IF( INFO.NE.0 ) THEN 00287 ILST = HERE 00288 RETURN 00289 END IF 00290 HERE = HERE + 1 00291 * 00292 ELSE 00293 * 00294 * Recompute NBNEXT in case of 2-by-2 split. 00295 * 00296 IF( A( HERE+2, HERE+1 ).EQ.ZERO ) 00297 $ NBNEXT = 1 00298 IF( NBNEXT.EQ.2 ) THEN 00299 * 00300 * 2-by-2 block did not split. 00301 * 00302 CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, 00303 $ Z, LDZ, HERE, 1, NBNEXT, WORK, LWORK, 00304 $ INFO ) 00305 IF( INFO.NE.0 ) THEN 00306 ILST = HERE 00307 RETURN 00308 END IF 00309 HERE = HERE + 2 00310 ELSE 00311 * 00312 * 2-by-2 block did split. 00313 * 00314 CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, 00315 $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO ) 00316 IF( INFO.NE.0 ) THEN 00317 ILST = HERE 00318 RETURN 00319 END IF 00320 HERE = HERE + 1 00321 CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, 00322 $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO ) 00323 IF( INFO.NE.0 ) THEN 00324 ILST = HERE 00325 RETURN 00326 END IF 00327 HERE = HERE + 1 00328 END IF 00329 * 00330 END IF 00331 END IF 00332 IF( HERE.LT.ILST ) 00333 $ GO TO 10 00334 ELSE 00335 HERE = IFST 00336 * 00337 20 CONTINUE 00338 * 00339 * Swap with next one below. 00340 * 00341 IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN 00342 * 00343 * Current block either 1-by-1 or 2-by-2. 00344 * 00345 NBNEXT = 1 00346 IF( HERE.GE.3 ) THEN 00347 IF( A( HERE-1, HERE-2 ).NE.ZERO ) 00348 $ NBNEXT = 2 00349 END IF 00350 CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 00351 $ LDZ, HERE-NBNEXT, NBNEXT, NBF, WORK, LWORK, 00352 $ INFO ) 00353 IF( INFO.NE.0 ) THEN 00354 ILST = HERE 00355 RETURN 00356 END IF 00357 HERE = HERE - NBNEXT 00358 * 00359 * Test if 2-by-2 block breaks into two 1-by-1 blocks. 00360 * 00361 IF( NBF.EQ.2 ) THEN 00362 IF( A( HERE+1, HERE ).EQ.ZERO ) 00363 $ NBF = 3 00364 END IF 00365 * 00366 ELSE 00367 * 00368 * Current block consists of two 1-by-1 blocks, each of which 00369 * must be swapped individually. 00370 * 00371 NBNEXT = 1 00372 IF( HERE.GE.3 ) THEN 00373 IF( A( HERE-1, HERE-2 ).NE.ZERO ) 00374 $ NBNEXT = 2 00375 END IF 00376 CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 00377 $ LDZ, HERE-NBNEXT, NBNEXT, 1, WORK, LWORK, 00378 $ INFO ) 00379 IF( INFO.NE.0 ) THEN 00380 ILST = HERE 00381 RETURN 00382 END IF 00383 IF( NBNEXT.EQ.1 ) THEN 00384 * 00385 * Swap two 1-by-1 blocks. 00386 * 00387 CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 00388 $ LDZ, HERE, NBNEXT, 1, WORK, LWORK, INFO ) 00389 IF( INFO.NE.0 ) THEN 00390 ILST = HERE 00391 RETURN 00392 END IF 00393 HERE = HERE - 1 00394 ELSE 00395 * 00396 * Recompute NBNEXT in case of 2-by-2 split. 00397 * 00398 IF( A( HERE, HERE-1 ).EQ.ZERO ) 00399 $ NBNEXT = 1 00400 IF( NBNEXT.EQ.2 ) THEN 00401 * 00402 * 2-by-2 block did not split. 00403 * 00404 CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, 00405 $ Z, LDZ, HERE-1, 2, 1, WORK, LWORK, INFO ) 00406 IF( INFO.NE.0 ) THEN 00407 ILST = HERE 00408 RETURN 00409 END IF 00410 HERE = HERE - 2 00411 ELSE 00412 * 00413 * 2-by-2 block did split. 00414 * 00415 CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, 00416 $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO ) 00417 IF( INFO.NE.0 ) THEN 00418 ILST = HERE 00419 RETURN 00420 END IF 00421 HERE = HERE - 1 00422 CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, 00423 $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO ) 00424 IF( INFO.NE.0 ) THEN 00425 ILST = HERE 00426 RETURN 00427 END IF 00428 HERE = HERE - 1 00429 END IF 00430 END IF 00431 END IF 00432 IF( HERE.GT.ILST ) 00433 $ GO TO 20 00434 END IF 00435 ILST = HERE 00436 WORK( 1 ) = LWMIN 00437 RETURN 00438 * 00439 * End of STGEXC 00440 * 00441 END
1.7.3 
