![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DTREXC 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DTREXC + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrexc.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrexc.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrexc.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DTREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK, 00022 * INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER COMPQ 00026 * INTEGER IFST, ILST, INFO, LDQ, LDT, N 00027 * .. 00028 * .. Array Arguments .. 00029 * DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * ) 00030 * .. 00031 * 00032 * 00033 *> \par Purpose: 00034 * ============= 00035 *> 00036 *> \verbatim 00037 *> 00038 *> DTREXC reorders the real Schur factorization of a real matrix 00039 *> A = Q*T*Q**T, so that the diagonal block of T with row index IFST is 00040 *> moved to row ILST. 00041 *> 00042 *> The real Schur form T is reordered by an orthogonal similarity 00043 *> transformation Z**T*T*Z, and optionally the matrix Q of Schur vectors 00044 *> is updated by postmultiplying it with Z. 00045 *> 00046 *> T must be in Schur canonical form (as returned by DHSEQR), that is, 00047 *> block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each 00048 *> 2-by-2 diagonal block has its diagonal elements equal and its 00049 *> off-diagonal elements of opposite sign. 00050 *> \endverbatim 00051 * 00052 * Arguments: 00053 * ========== 00054 * 00055 *> \param[in] COMPQ 00056 *> \verbatim 00057 *> COMPQ is CHARACTER*1 00058 *> = 'V': update the matrix Q of Schur vectors; 00059 *> = 'N': do not update Q. 00060 *> \endverbatim 00061 *> 00062 *> \param[in] N 00063 *> \verbatim 00064 *> N is INTEGER 00065 *> The order of the matrix T. N >= 0. 00066 *> \endverbatim 00067 *> 00068 *> \param[in,out] T 00069 *> \verbatim 00070 *> T is DOUBLE PRECISION array, dimension (LDT,N) 00071 *> On entry, the upper quasi-triangular matrix T, in Schur 00072 *> Schur canonical form. 00073 *> On exit, the reordered upper quasi-triangular matrix, again 00074 *> in Schur canonical form. 00075 *> \endverbatim 00076 *> 00077 *> \param[in] LDT 00078 *> \verbatim 00079 *> LDT is INTEGER 00080 *> The leading dimension of the array T. LDT >= max(1,N). 00081 *> \endverbatim 00082 *> 00083 *> \param[in,out] Q 00084 *> \verbatim 00085 *> Q is DOUBLE PRECISION array, dimension (LDQ,N) 00086 *> On entry, if COMPQ = 'V', the matrix Q of Schur vectors. 00087 *> On exit, if COMPQ = 'V', Q has been postmultiplied by the 00088 *> orthogonal transformation matrix Z which reorders T. 00089 *> If COMPQ = 'N', Q is not referenced. 00090 *> \endverbatim 00091 *> 00092 *> \param[in] LDQ 00093 *> \verbatim 00094 *> LDQ is INTEGER 00095 *> The leading dimension of the array Q. LDQ >= max(1,N). 00096 *> \endverbatim 00097 *> 00098 *> \param[in,out] IFST 00099 *> \verbatim 00100 *> IFST is INTEGER 00101 *> \endverbatim 00102 *> 00103 *> \param[in,out] ILST 00104 *> \verbatim 00105 *> ILST is INTEGER 00106 *> 00107 *> Specify the reordering of the diagonal blocks of T. 00108 *> The block with row index IFST is moved to row ILST, by a 00109 *> sequence of transpositions between adjacent blocks. 00110 *> On exit, if IFST pointed on entry to the second row of a 00111 *> 2-by-2 block, it is changed to point to the first row; ILST 00112 *> always points to the first row of the block in its final 00113 *> position (which may differ from its input value by +1 or -1). 00114 *> 1 <= IFST <= N; 1 <= ILST <= N. 00115 *> \endverbatim 00116 *> 00117 *> \param[out] WORK 00118 *> \verbatim 00119 *> WORK is DOUBLE PRECISION array, dimension (N) 00120 *> \endverbatim 00121 *> 00122 *> \param[out] INFO 00123 *> \verbatim 00124 *> INFO is INTEGER 00125 *> = 0: successful exit 00126 *> < 0: if INFO = -i, the i-th argument had an illegal value 00127 *> = 1: two adjacent blocks were too close to swap (the problem 00128 *> is very ill-conditioned); T may have been partially 00129 *> reordered, and ILST points to the first row of the 00130 *> current position of the block being moved. 00131 *> \endverbatim 00132 * 00133 * Authors: 00134 * ======== 00135 * 00136 *> \author Univ. of Tennessee 00137 *> \author Univ. of California Berkeley 00138 *> \author Univ. of Colorado Denver 00139 *> \author NAG Ltd. 00140 * 00141 *> \date November 2011 00142 * 00143 *> \ingroup doubleOTHERcomputational 00144 * 00145 * ===================================================================== 00146 SUBROUTINE DTREXC( COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, WORK, 00147 $ INFO ) 00148 * 00149 * -- LAPACK computational routine (version 3.4.0) -- 00150 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00151 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00152 * November 2011 00153 * 00154 * .. Scalar Arguments .. 00155 CHARACTER COMPQ 00156 INTEGER IFST, ILST, INFO, LDQ, LDT, N 00157 * .. 00158 * .. Array Arguments .. 00159 DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * ) 00160 * .. 00161 * 00162 * ===================================================================== 00163 * 00164 * .. Parameters .. 00165 DOUBLE PRECISION ZERO 00166 PARAMETER ( ZERO = 0.0D+0 ) 00167 * .. 00168 * .. Local Scalars .. 00169 LOGICAL WANTQ 00170 INTEGER HERE, NBF, NBL, NBNEXT 00171 * .. 00172 * .. External Functions .. 00173 LOGICAL LSAME 00174 EXTERNAL LSAME 00175 * .. 00176 * .. External Subroutines .. 00177 EXTERNAL DLAEXC, XERBLA 00178 * .. 00179 * .. Intrinsic Functions .. 00180 INTRINSIC MAX 00181 * .. 00182 * .. Executable Statements .. 00183 * 00184 * Decode and test the input arguments. 00185 * 00186 INFO = 0 00187 WANTQ = LSAME( COMPQ, 'V' ) 00188 IF( .NOT.WANTQ .AND. .NOT.LSAME( COMPQ, 'N' ) ) THEN 00189 INFO = -1 00190 ELSE IF( N.LT.0 ) THEN 00191 INFO = -2 00192 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN 00193 INFO = -4 00194 ELSE IF( LDQ.LT.1 .OR. ( WANTQ .AND. LDQ.LT.MAX( 1, N ) ) ) THEN 00195 INFO = -6 00196 ELSE IF( IFST.LT.1 .OR. IFST.GT.N ) THEN 00197 INFO = -7 00198 ELSE IF( ILST.LT.1 .OR. ILST.GT.N ) THEN 00199 INFO = -8 00200 END IF 00201 IF( INFO.NE.0 ) THEN 00202 CALL XERBLA( 'DTREXC', -INFO ) 00203 RETURN 00204 END IF 00205 * 00206 * Quick return if possible 00207 * 00208 IF( N.LE.1 ) 00209 $ RETURN 00210 * 00211 * Determine the first row of specified block 00212 * and find out it is 1 by 1 or 2 by 2. 00213 * 00214 IF( IFST.GT.1 ) THEN 00215 IF( T( IFST, IFST-1 ).NE.ZERO ) 00216 $ IFST = IFST - 1 00217 END IF 00218 NBF = 1 00219 IF( IFST.LT.N ) THEN 00220 IF( T( IFST+1, IFST ).NE.ZERO ) 00221 $ NBF = 2 00222 END IF 00223 * 00224 * Determine the first row of the final block 00225 * and find out it is 1 by 1 or 2 by 2. 00226 * 00227 IF( ILST.GT.1 ) THEN 00228 IF( T( ILST, ILST-1 ).NE.ZERO ) 00229 $ ILST = ILST - 1 00230 END IF 00231 NBL = 1 00232 IF( ILST.LT.N ) THEN 00233 IF( T( ILST+1, ILST ).NE.ZERO ) 00234 $ NBL = 2 00235 END IF 00236 * 00237 IF( IFST.EQ.ILST ) 00238 $ RETURN 00239 * 00240 IF( IFST.LT.ILST ) THEN 00241 * 00242 * Update ILST 00243 * 00244 IF( NBF.EQ.2 .AND. NBL.EQ.1 ) 00245 $ ILST = ILST - 1 00246 IF( NBF.EQ.1 .AND. NBL.EQ.2 ) 00247 $ ILST = ILST + 1 00248 * 00249 HERE = IFST 00250 * 00251 10 CONTINUE 00252 * 00253 * Swap block with next one below 00254 * 00255 IF( NBF.EQ.1 .OR. NBF.EQ.2 ) THEN 00256 * 00257 * Current block either 1 by 1 or 2 by 2 00258 * 00259 NBNEXT = 1 00260 IF( HERE+NBF+1.LE.N ) THEN 00261 IF( T( HERE+NBF+1, HERE+NBF ).NE.ZERO ) 00262 $ NBNEXT = 2 00263 END IF 00264 CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBF, NBNEXT, 00265 $ WORK, INFO ) 00266 IF( INFO.NE.0 ) THEN 00267 ILST = HERE 00268 RETURN 00269 END IF 00270 HERE = HERE + NBNEXT 00271 * 00272 * Test if 2 by 2 block breaks into two 1 by 1 blocks 00273 * 00274 IF( NBF.EQ.2 ) THEN 00275 IF( T( HERE+1, HERE ).EQ.ZERO ) 00276 $ NBF = 3 00277 END IF 00278 * 00279 ELSE 00280 * 00281 * Current block consists of two 1 by 1 blocks each of which 00282 * must be swapped individually 00283 * 00284 NBNEXT = 1 00285 IF( HERE+3.LE.N ) THEN 00286 IF( T( HERE+3, HERE+2 ).NE.ZERO ) 00287 $ NBNEXT = 2 00288 END IF 00289 CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, NBNEXT, 00290 $ WORK, INFO ) 00291 IF( INFO.NE.0 ) THEN 00292 ILST = HERE 00293 RETURN 00294 END IF 00295 IF( NBNEXT.EQ.1 ) THEN 00296 * 00297 * Swap two 1 by 1 blocks, no problems possible 00298 * 00299 CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, NBNEXT, 00300 $ WORK, INFO ) 00301 HERE = HERE + 1 00302 ELSE 00303 * 00304 * Recompute NBNEXT in case 2 by 2 split 00305 * 00306 IF( T( HERE+2, HERE+1 ).EQ.ZERO ) 00307 $ NBNEXT = 1 00308 IF( NBNEXT.EQ.2 ) THEN 00309 * 00310 * 2 by 2 Block did not split 00311 * 00312 CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 00313 $ NBNEXT, WORK, INFO ) 00314 IF( INFO.NE.0 ) THEN 00315 ILST = HERE 00316 RETURN 00317 END IF 00318 HERE = HERE + 2 00319 ELSE 00320 * 00321 * 2 by 2 Block did split 00322 * 00323 CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1, 00324 $ WORK, INFO ) 00325 CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE+1, 1, 1, 00326 $ WORK, INFO ) 00327 HERE = HERE + 2 00328 END IF 00329 END IF 00330 END IF 00331 IF( HERE.LT.ILST ) 00332 $ GO TO 10 00333 * 00334 ELSE 00335 * 00336 HERE = IFST 00337 20 CONTINUE 00338 * 00339 * Swap block with next one above 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( T( HERE-1, HERE-2 ).NE.ZERO ) 00348 $ NBNEXT = 2 00349 END IF 00350 CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT, 00351 $ NBF, WORK, INFO ) 00352 IF( INFO.NE.0 ) THEN 00353 ILST = HERE 00354 RETURN 00355 END IF 00356 HERE = HERE - NBNEXT 00357 * 00358 * Test if 2 by 2 block breaks into two 1 by 1 blocks 00359 * 00360 IF( NBF.EQ.2 ) THEN 00361 IF( T( HERE+1, HERE ).EQ.ZERO ) 00362 $ NBF = 3 00363 END IF 00364 * 00365 ELSE 00366 * 00367 * Current block consists of two 1 by 1 blocks each of which 00368 * must be swapped individually 00369 * 00370 NBNEXT = 1 00371 IF( HERE.GE.3 ) THEN 00372 IF( T( HERE-1, HERE-2 ).NE.ZERO ) 00373 $ NBNEXT = 2 00374 END IF 00375 CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-NBNEXT, NBNEXT, 00376 $ 1, WORK, INFO ) 00377 IF( INFO.NE.0 ) THEN 00378 ILST = HERE 00379 RETURN 00380 END IF 00381 IF( NBNEXT.EQ.1 ) THEN 00382 * 00383 * Swap two 1 by 1 blocks, no problems possible 00384 * 00385 CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, NBNEXT, 1, 00386 $ WORK, INFO ) 00387 HERE = HERE - 1 00388 ELSE 00389 * 00390 * Recompute NBNEXT in case 2 by 2 split 00391 * 00392 IF( T( HERE, HERE-1 ).EQ.ZERO ) 00393 $ NBNEXT = 1 00394 IF( NBNEXT.EQ.2 ) THEN 00395 * 00396 * 2 by 2 Block did not split 00397 * 00398 CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 2, 1, 00399 $ WORK, INFO ) 00400 IF( INFO.NE.0 ) THEN 00401 ILST = HERE 00402 RETURN 00403 END IF 00404 HERE = HERE - 2 00405 ELSE 00406 * 00407 * 2 by 2 Block did split 00408 * 00409 CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE, 1, 1, 00410 $ WORK, INFO ) 00411 CALL DLAEXC( WANTQ, N, T, LDT, Q, LDQ, HERE-1, 1, 1, 00412 $ WORK, INFO ) 00413 HERE = HERE - 2 00414 END IF 00415 END IF 00416 END IF 00417 IF( HERE.GT.ILST ) 00418 $ GO TO 20 00419 END IF 00420 ILST = HERE 00421 * 00422 RETURN 00423 * 00424 * End of DTREXC 00425 * 00426 END