![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b STGEXC 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download STGEXC + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/stgexc.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/stgexc.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/stgexc.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE STGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 00022 * LDZ, IFST, ILST, WORK, LWORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * LOGICAL WANTQ, WANTZ 00026 * INTEGER IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N 00027 * .. 00028 * .. Array Arguments .. 00029 * REAL A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 00030 * $ WORK( * ), Z( LDZ, * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> STGEXC reorders the generalized real Schur decomposition of a real 00040 *> matrix pair (A,B) using an orthogonal equivalence transformation 00041 *> 00042 *> (A, B) = Q * (A, B) * Z**T, 00043 *> 00044 *> so that the diagonal block of (A, B) with row index IFST is moved 00045 *> to row ILST. 00046 *> 00047 *> (A, B) must be in generalized real Schur canonical form (as returned 00048 *> by SGGES), i.e. A is block upper triangular with 1-by-1 and 2-by-2 00049 *> diagonal blocks. B is upper triangular. 00050 *> 00051 *> Optionally, the matrices Q and Z of generalized Schur vectors are 00052 *> updated. 00053 *> 00054 *> Q(in) * A(in) * Z(in)**T = Q(out) * A(out) * Z(out)**T 00055 *> Q(in) * B(in) * Z(in)**T = Q(out) * B(out) * Z(out)**T 00056 *> 00057 *> \endverbatim 00058 * 00059 * Arguments: 00060 * ========== 00061 * 00062 *> \param[in] WANTQ 00063 *> \verbatim 00064 *> WANTQ is LOGICAL 00065 *> .TRUE. : update the left transformation matrix Q; 00066 *> .FALSE.: do not update Q. 00067 *> \endverbatim 00068 *> 00069 *> \param[in] WANTZ 00070 *> \verbatim 00071 *> WANTZ is LOGICAL 00072 *> .TRUE. : update the right transformation matrix Z; 00073 *> .FALSE.: do not update Z. 00074 *> \endverbatim 00075 *> 00076 *> \param[in] N 00077 *> \verbatim 00078 *> N is INTEGER 00079 *> The order of the matrices A and B. N >= 0. 00080 *> \endverbatim 00081 *> 00082 *> \param[in,out] A 00083 *> \verbatim 00084 *> A is REAL array, dimension (LDA,N) 00085 *> On entry, the matrix A in generalized real Schur canonical 00086 *> form. 00087 *> On exit, the updated matrix A, again in generalized 00088 *> real Schur canonical form. 00089 *> \endverbatim 00090 *> 00091 *> \param[in] LDA 00092 *> \verbatim 00093 *> LDA is INTEGER 00094 *> The leading dimension of the array A. LDA >= max(1,N). 00095 *> \endverbatim 00096 *> 00097 *> \param[in,out] B 00098 *> \verbatim 00099 *> B is REAL array, dimension (LDB,N) 00100 *> On entry, the matrix B in generalized real Schur canonical 00101 *> form (A,B). 00102 *> On exit, the updated matrix B, again in generalized 00103 *> real Schur canonical form (A,B). 00104 *> \endverbatim 00105 *> 00106 *> \param[in] LDB 00107 *> \verbatim 00108 *> LDB is INTEGER 00109 *> The leading dimension of the array B. LDB >= max(1,N). 00110 *> \endverbatim 00111 *> 00112 *> \param[in,out] Q 00113 *> \verbatim 00114 *> Q is REAL array, dimension (LDZ,N) 00115 *> On entry, if WANTQ = .TRUE., the orthogonal matrix Q. 00116 *> On exit, the updated matrix Q. 00117 *> If WANTQ = .FALSE., Q is not referenced. 00118 *> \endverbatim 00119 *> 00120 *> \param[in] LDQ 00121 *> \verbatim 00122 *> LDQ is INTEGER 00123 *> The leading dimension of the array Q. LDQ >= 1. 00124 *> If WANTQ = .TRUE., LDQ >= N. 00125 *> \endverbatim 00126 *> 00127 *> \param[in,out] Z 00128 *> \verbatim 00129 *> Z is REAL array, dimension (LDZ,N) 00130 *> On entry, if WANTZ = .TRUE., the orthogonal matrix Z. 00131 *> On exit, the updated matrix Z. 00132 *> If WANTZ = .FALSE., Z is not referenced. 00133 *> \endverbatim 00134 *> 00135 *> \param[in] LDZ 00136 *> \verbatim 00137 *> LDZ is INTEGER 00138 *> The leading dimension of the array Z. LDZ >= 1. 00139 *> If WANTZ = .TRUE., LDZ >= N. 00140 *> \endverbatim 00141 *> 00142 *> \param[in,out] IFST 00143 *> \verbatim 00144 *> IFST is INTEGER 00145 *> \endverbatim 00146 *> 00147 *> \param[in,out] ILST 00148 *> \verbatim 00149 *> ILST is INTEGER 00150 *> Specify the reordering of the diagonal blocks of (A, B). 00151 *> The block with row index IFST is moved to row ILST, by a 00152 *> sequence of swapping between adjacent blocks. 00153 *> On exit, if IFST pointed on entry to the second row of 00154 *> a 2-by-2 block, it is changed to point to the first row; 00155 *> ILST always points to the first row of the block in its 00156 *> final position (which may differ from its input value by 00157 *> +1 or -1). 1 <= IFST, ILST <= N. 00158 *> \endverbatim 00159 *> 00160 *> \param[out] WORK 00161 *> \verbatim 00162 *> WORK is REAL array, dimension (MAX(1,LWORK)) 00163 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00164 *> \endverbatim 00165 *> 00166 *> \param[in] LWORK 00167 *> \verbatim 00168 *> LWORK is INTEGER 00169 *> The dimension of the array WORK. 00170 *> LWORK >= 1 when N <= 1, otherwise LWORK >= 4*N + 16. 00171 *> 00172 *> If LWORK = -1, then a workspace query is assumed; the routine 00173 *> only calculates the optimal size of the WORK array, returns 00174 *> this value as the first entry of the WORK array, and no error 00175 *> message related to LWORK is issued by XERBLA. 00176 *> \endverbatim 00177 *> 00178 *> \param[out] INFO 00179 *> \verbatim 00180 *> INFO is INTEGER 00181 *> =0: successful exit. 00182 *> <0: if INFO = -i, the i-th argument had an illegal value. 00183 *> =1: The transformed matrix pair (A, B) would be too far 00184 *> from generalized Schur form; the problem is ill- 00185 *> conditioned. (A, B) may have been partially reordered, 00186 *> and ILST points to the first row of the current 00187 *> position of the block being moved. 00188 *> \endverbatim 00189 * 00190 * Authors: 00191 * ======== 00192 * 00193 *> \author Univ. of Tennessee 00194 *> \author Univ. of California Berkeley 00195 *> \author Univ. of Colorado Denver 00196 *> \author NAG Ltd. 00197 * 00198 *> \date November 2011 00199 * 00200 *> \ingroup realGEcomputational 00201 * 00202 *> \par Contributors: 00203 * ================== 00204 *> 00205 *> Bo Kagstrom and Peter Poromaa, Department of Computing Science, 00206 *> Umea University, S-901 87 Umea, Sweden. 00207 * 00208 *> \par References: 00209 * ================ 00210 *> 00211 *> \verbatim 00212 *> 00213 *> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the 00214 *> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in 00215 *> M.S. Moonen et al (eds), Linear Algebra for Large Scale and 00216 *> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218. 00217 *> \endverbatim 00218 *> 00219 * ===================================================================== 00220 SUBROUTINE STGEXC( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 00221 $ LDZ, IFST, ILST, WORK, LWORK, INFO ) 00222 * 00223 * -- LAPACK computational routine (version 3.4.0) -- 00224 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00225 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00226 * November 2011 00227 * 00228 * .. Scalar Arguments .. 00229 LOGICAL WANTQ, WANTZ 00230 INTEGER IFST, ILST, INFO, LDA, LDB, LDQ, LDZ, LWORK, N 00231 * .. 00232 * .. Array Arguments .. 00233 REAL A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 00234 $ WORK( * ), Z( LDZ, * ) 00235 * .. 00236 * 00237 * ===================================================================== 00238 * 00239 * .. Parameters .. 00240 REAL ZERO 00241 PARAMETER ( ZERO = 0.0E+0 ) 00242 * .. 00243 * .. Local Scalars .. 00244 LOGICAL LQUERY 00245 INTEGER HERE, LWMIN, NBF, NBL, NBNEXT 00246 * .. 00247 * .. External Subroutines .. 00248 EXTERNAL STGEX2, XERBLA 00249 * .. 00250 * .. Intrinsic Functions .. 00251 INTRINSIC MAX 00252 * .. 00253 * .. Executable Statements .. 00254 * 00255 * Decode and test input arguments. 00256 * 00257 INFO = 0 00258 LQUERY = ( LWORK.EQ.-1 ) 00259 IF( N.LT.0 ) THEN 00260 INFO = -3 00261 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00262 INFO = -5 00263 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00264 INFO = -7 00265 ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. ( LDQ.LT.MAX( 1, N ) ) ) THEN 00266 INFO = -9 00267 ELSE IF( LDZ.LT.1 .OR. WANTZ .AND. ( LDZ.LT.MAX( 1, N ) ) ) THEN 00268 INFO = -11 00269 ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN 00270 INFO = -12 00271 ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN 00272 INFO = -13 00273 END IF 00274 * 00275 IF( INFO.EQ.0 ) THEN 00276 IF( N.LE.1 ) THEN 00277 LWMIN = 1 00278 ELSE 00279 LWMIN = 4*N + 16 00280 END IF 00281 WORK(1) = LWMIN 00282 * 00283 IF (LWORK.LT.LWMIN .AND. .NOT.LQUERY) THEN 00284 INFO = -15 00285 END IF 00286 END IF 00287 * 00288 IF( INFO.NE.0 ) THEN 00289 CALL XERBLA( 'STGEXC', -INFO ) 00290 RETURN 00291 ELSE IF( LQUERY ) THEN 00292 RETURN 00293 END IF 00294 * 00295 * Quick return if possible 00296 * 00297 IF( N.LE.1 ) 00298 $ RETURN 00299 * 00300 * Determine the first row of the specified block and find out 00301 * if it is 1-by-1 or 2-by-2. 00302 * 00303 IF( IFST.GT.1 ) THEN 00304 IF( A( IFST, IFST-1 ).NE.ZERO ) 00305 $ IFST = IFST - 1 00306 END IF 00307 NBF = 1 00308 IF( IFST.LT.N ) THEN 00309 IF( A( IFST+1, IFST ).NE.ZERO ) 00310 $ NBF = 2 00311 END IF 00312 * 00313 * Determine the first row of the final block 00314 * and find out if it is 1-by-1 or 2-by-2. 00315 * 00316 IF( ILST.GT.1 ) THEN 00317 IF( A( ILST, ILST-1 ).NE.ZERO ) 00318 $ ILST = ILST - 1 00319 END IF 00320 NBL = 1 00321 IF( ILST.LT.N ) THEN 00322 IF( A( ILST+1, ILST ).NE.ZERO ) 00323 $ NBL = 2 00324 END IF 00325 IF( IFST.EQ.ILST ) 00326 $ RETURN 00327 * 00328 IF( IFST.LT.ILST ) THEN 00329 * 00330 * Update ILST. 00331 * 00332 IF( NBF.EQ.2 .AND. NBL.EQ.1 ) 00333 $ ILST = ILST - 1 00334 IF( NBF.EQ.1 .AND. NBL.EQ.2 ) 00335 $ ILST = ILST + 1 00336 * 00337 HERE = IFST 00338 * 00339 10 CONTINUE 00340 * 00341 * Swap with next one below. 00342 * 00343 IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN 00344 * 00345 * Current block either 1-by-1 or 2-by-2. 00346 * 00347 NBNEXT = 1 00348 IF( HERE+NBF+1.LE.N ) THEN 00349 IF( A( HERE+NBF+1, HERE+NBF ).NE.ZERO ) 00350 $ NBNEXT = 2 00351 END IF 00352 CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 00353 $ LDZ, HERE, NBF, NBNEXT, WORK, LWORK, INFO ) 00354 IF( INFO.NE.0 ) THEN 00355 ILST = HERE 00356 RETURN 00357 END IF 00358 HERE = HERE + NBNEXT 00359 * 00360 * Test if 2-by-2 block breaks into two 1-by-1 blocks. 00361 * 00362 IF( NBF.EQ.2 ) THEN 00363 IF( A( HERE+1, HERE ).EQ.ZERO ) 00364 $ NBF = 3 00365 END IF 00366 * 00367 ELSE 00368 * 00369 * Current block consists of two 1-by-1 blocks, each of which 00370 * must be swapped individually. 00371 * 00372 NBNEXT = 1 00373 IF( HERE+3.LE.N ) THEN 00374 IF( A( HERE+3, HERE+2 ).NE.ZERO ) 00375 $ NBNEXT = 2 00376 END IF 00377 CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 00378 $ LDZ, HERE+1, 1, NBNEXT, WORK, LWORK, 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, 1, 1, WORK, LWORK, INFO ) 00389 IF( INFO.NE.0 ) THEN 00390 ILST = HERE 00391 RETURN 00392 END IF 00393 HERE = HERE + 1 00394 * 00395 ELSE 00396 * 00397 * Recompute NBNEXT in case of 2-by-2 split. 00398 * 00399 IF( A( HERE+2, HERE+1 ).EQ.ZERO ) 00400 $ NBNEXT = 1 00401 IF( NBNEXT.EQ.2 ) THEN 00402 * 00403 * 2-by-2 block did not split. 00404 * 00405 CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, 00406 $ Z, LDZ, HERE, 1, NBNEXT, WORK, LWORK, 00407 $ INFO ) 00408 IF( INFO.NE.0 ) THEN 00409 ILST = HERE 00410 RETURN 00411 END IF 00412 HERE = HERE + 2 00413 ELSE 00414 * 00415 * 2-by-2 block did split. 00416 * 00417 CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, 00418 $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO ) 00419 IF( INFO.NE.0 ) THEN 00420 ILST = HERE 00421 RETURN 00422 END IF 00423 HERE = HERE + 1 00424 CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, 00425 $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO ) 00426 IF( INFO.NE.0 ) THEN 00427 ILST = HERE 00428 RETURN 00429 END IF 00430 HERE = HERE + 1 00431 END IF 00432 * 00433 END IF 00434 END IF 00435 IF( HERE.LT.ILST ) 00436 $ GO TO 10 00437 ELSE 00438 HERE = IFST 00439 * 00440 20 CONTINUE 00441 * 00442 * Swap with next one below. 00443 * 00444 IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN 00445 * 00446 * Current block either 1-by-1 or 2-by-2. 00447 * 00448 NBNEXT = 1 00449 IF( HERE.GE.3 ) THEN 00450 IF( A( HERE-1, HERE-2 ).NE.ZERO ) 00451 $ NBNEXT = 2 00452 END IF 00453 CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 00454 $ LDZ, HERE-NBNEXT, NBNEXT, NBF, WORK, LWORK, 00455 $ INFO ) 00456 IF( INFO.NE.0 ) THEN 00457 ILST = HERE 00458 RETURN 00459 END IF 00460 HERE = HERE - NBNEXT 00461 * 00462 * Test if 2-by-2 block breaks into two 1-by-1 blocks. 00463 * 00464 IF( NBF.EQ.2 ) THEN 00465 IF( A( HERE+1, HERE ).EQ.ZERO ) 00466 $ NBF = 3 00467 END IF 00468 * 00469 ELSE 00470 * 00471 * Current block consists of two 1-by-1 blocks, each of which 00472 * must be swapped individually. 00473 * 00474 NBNEXT = 1 00475 IF( HERE.GE.3 ) THEN 00476 IF( A( HERE-1, HERE-2 ).NE.ZERO ) 00477 $ NBNEXT = 2 00478 END IF 00479 CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 00480 $ LDZ, HERE-NBNEXT, NBNEXT, 1, WORK, LWORK, 00481 $ INFO ) 00482 IF( INFO.NE.0 ) THEN 00483 ILST = HERE 00484 RETURN 00485 END IF 00486 IF( NBNEXT.EQ.1 ) THEN 00487 * 00488 * Swap two 1-by-1 blocks. 00489 * 00490 CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 00491 $ LDZ, HERE, NBNEXT, 1, WORK, LWORK, INFO ) 00492 IF( INFO.NE.0 ) THEN 00493 ILST = HERE 00494 RETURN 00495 END IF 00496 HERE = HERE - 1 00497 ELSE 00498 * 00499 * Recompute NBNEXT in case of 2-by-2 split. 00500 * 00501 IF( A( HERE, HERE-1 ).EQ.ZERO ) 00502 $ NBNEXT = 1 00503 IF( NBNEXT.EQ.2 ) THEN 00504 * 00505 * 2-by-2 block did not split. 00506 * 00507 CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, 00508 $ Z, LDZ, HERE-1, 2, 1, WORK, LWORK, INFO ) 00509 IF( INFO.NE.0 ) THEN 00510 ILST = HERE 00511 RETURN 00512 END IF 00513 HERE = HERE - 2 00514 ELSE 00515 * 00516 * 2-by-2 block did split. 00517 * 00518 CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, 00519 $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO ) 00520 IF( INFO.NE.0 ) THEN 00521 ILST = HERE 00522 RETURN 00523 END IF 00524 HERE = HERE - 1 00525 CALL STGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, 00526 $ Z, LDZ, HERE, 1, 1, WORK, LWORK, INFO ) 00527 IF( INFO.NE.0 ) THEN 00528 ILST = HERE 00529 RETURN 00530 END IF 00531 HERE = HERE - 1 00532 END IF 00533 END IF 00534 END IF 00535 IF( HERE.GT.ILST ) 00536 $ GO TO 20 00537 END IF 00538 ILST = HERE 00539 WORK( 1 ) = LWMIN 00540 RETURN 00541 * 00542 * End of STGEXC 00543 * 00544 END