![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZTGEX2 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZTGEX2 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztgex2.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztgex2.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztgex2.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 00022 * LDZ, J1, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * LOGICAL WANTQ, WANTZ 00026 * INTEGER INFO, J1, LDA, LDB, LDQ, LDZ, N 00027 * .. 00028 * .. Array Arguments .. 00029 * COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 00030 * $ Z( LDZ, * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> ZTGEX2 swaps adjacent diagonal 1 by 1 blocks (A11,B11) and (A22,B22) 00040 *> in an upper triangular matrix pair (A, B) by an unitary equivalence 00041 *> transformation. 00042 *> 00043 *> (A, B) must be in generalized Schur canonical form, that is, A and 00044 *> B are both upper triangular. 00045 *> 00046 *> Optionally, the matrices Q and Z of generalized Schur vectors are 00047 *> updated. 00048 *> 00049 *> Q(in) * A(in) * Z(in)**H = Q(out) * A(out) * Z(out)**H 00050 *> Q(in) * B(in) * Z(in)**H = Q(out) * B(out) * Z(out)**H 00051 *> 00052 *> \endverbatim 00053 * 00054 * Arguments: 00055 * ========== 00056 * 00057 *> \param[in] WANTQ 00058 *> \verbatim 00059 *> WANTQ is LOGICAL 00060 *> .TRUE. : update the left transformation matrix Q; 00061 *> .FALSE.: do not update Q. 00062 *> \endverbatim 00063 *> 00064 *> \param[in] WANTZ 00065 *> \verbatim 00066 *> WANTZ is LOGICAL 00067 *> .TRUE. : update the right transformation matrix Z; 00068 *> .FALSE.: do not update Z. 00069 *> \endverbatim 00070 *> 00071 *> \param[in] N 00072 *> \verbatim 00073 *> N is INTEGER 00074 *> The order of the matrices A and B. N >= 0. 00075 *> \endverbatim 00076 *> 00077 *> \param[in,out] A 00078 *> \verbatim 00079 *> A is COMPLEX*16 arrays, dimensions (LDA,N) 00080 *> On entry, the matrix A in the pair (A, B). 00081 *> On exit, the updated matrix A. 00082 *> \endverbatim 00083 *> 00084 *> \param[in] LDA 00085 *> \verbatim 00086 *> LDA is INTEGER 00087 *> The leading dimension of the array A. LDA >= max(1,N). 00088 *> \endverbatim 00089 *> 00090 *> \param[in,out] B 00091 *> \verbatim 00092 *> B is COMPLEX*16 arrays, dimensions (LDB,N) 00093 *> On entry, the matrix B in the pair (A, B). 00094 *> On exit, the updated matrix B. 00095 *> \endverbatim 00096 *> 00097 *> \param[in] LDB 00098 *> \verbatim 00099 *> LDB is INTEGER 00100 *> The leading dimension of the array B. LDB >= max(1,N). 00101 *> \endverbatim 00102 *> 00103 *> \param[in,out] Q 00104 *> \verbatim 00105 *> Q is COMPLEX*16 array, dimension (LDZ,N) 00106 *> If WANTQ = .TRUE, on entry, the unitary matrix Q. On exit, 00107 *> the updated matrix Q. 00108 *> Not referenced if WANTQ = .FALSE.. 00109 *> \endverbatim 00110 *> 00111 *> \param[in] LDQ 00112 *> \verbatim 00113 *> LDQ is INTEGER 00114 *> The leading dimension of the array Q. LDQ >= 1; 00115 *> If WANTQ = .TRUE., LDQ >= N. 00116 *> \endverbatim 00117 *> 00118 *> \param[in,out] Z 00119 *> \verbatim 00120 *> Z is COMPLEX*16 array, dimension (LDZ,N) 00121 *> If WANTZ = .TRUE, on entry, the unitary matrix Z. On exit, 00122 *> the updated matrix Z. 00123 *> Not referenced if WANTZ = .FALSE.. 00124 *> \endverbatim 00125 *> 00126 *> \param[in] LDZ 00127 *> \verbatim 00128 *> LDZ is INTEGER 00129 *> The leading dimension of the array Z. LDZ >= 1; 00130 *> If WANTZ = .TRUE., LDZ >= N. 00131 *> \endverbatim 00132 *> 00133 *> \param[in] J1 00134 *> \verbatim 00135 *> J1 is INTEGER 00136 *> The index to the first block (A11, B11). 00137 *> \endverbatim 00138 *> 00139 *> \param[out] INFO 00140 *> \verbatim 00141 *> INFO is INTEGER 00142 *> =0: Successful exit. 00143 *> =1: The transformed matrix pair (A, B) would be too far 00144 *> from generalized Schur form; the problem is ill- 00145 *> conditioned. 00146 *> \endverbatim 00147 * 00148 * Authors: 00149 * ======== 00150 * 00151 *> \author Univ. of Tennessee 00152 *> \author Univ. of California Berkeley 00153 *> \author Univ. of Colorado Denver 00154 *> \author NAG Ltd. 00155 * 00156 *> \date November 2011 00157 * 00158 *> \ingroup complex16GEauxiliary 00159 * 00160 *> \par Further Details: 00161 * ===================== 00162 *> 00163 *> In the current code both weak and strong stability tests are 00164 *> performed. The user can omit the strong stability test by changing 00165 *> the internal logical parameter WANDS to .FALSE.. See ref. [2] for 00166 *> details. 00167 * 00168 *> \par Contributors: 00169 * ================== 00170 *> 00171 *> Bo Kagstrom and Peter Poromaa, Department of Computing Science, 00172 *> Umea University, S-901 87 Umea, Sweden. 00173 * 00174 *> \par References: 00175 * ================ 00176 *> 00177 *> [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the 00178 *> Generalized Real Schur Form of a Regular Matrix Pair (A, B), in 00179 *> M.S. Moonen et al (eds), Linear Algebra for Large Scale and 00180 *> Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218. 00181 *> \n 00182 *> [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified 00183 *> Eigenvalues of a Regular Matrix Pair (A, B) and Condition 00184 *> Estimation: Theory, Algorithms and Software, Report UMINF-94.04, 00185 *> Department of Computing Science, Umea University, S-901 87 Umea, 00186 *> Sweden, 1994. Also as LAPACK Working Note 87. To appear in 00187 *> Numerical Algorithms, 1996. 00188 *> 00189 * ===================================================================== 00190 SUBROUTINE ZTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 00191 $ LDZ, J1, INFO ) 00192 * 00193 * -- LAPACK auxiliary routine (version 3.4.0) -- 00194 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00195 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00196 * November 2011 00197 * 00198 * .. Scalar Arguments .. 00199 LOGICAL WANTQ, WANTZ 00200 INTEGER INFO, J1, LDA, LDB, LDQ, LDZ, N 00201 * .. 00202 * .. Array Arguments .. 00203 COMPLEX*16 A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 00204 $ Z( LDZ, * ) 00205 * .. 00206 * 00207 * ===================================================================== 00208 * 00209 * .. Parameters .. 00210 COMPLEX*16 CZERO, CONE 00211 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 00212 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 00213 DOUBLE PRECISION TWENTY 00214 PARAMETER ( TWENTY = 2.0D+1 ) 00215 INTEGER LDST 00216 PARAMETER ( LDST = 2 ) 00217 LOGICAL WANDS 00218 PARAMETER ( WANDS = .TRUE. ) 00219 * .. 00220 * .. Local Scalars .. 00221 LOGICAL DTRONG, WEAK 00222 INTEGER I, M 00223 DOUBLE PRECISION CQ, CZ, EPS, SA, SB, SCALE, SMLNUM, SS, SUM, 00224 $ THRESH, WS 00225 COMPLEX*16 CDUM, F, G, SQ, SZ 00226 * .. 00227 * .. Local Arrays .. 00228 COMPLEX*16 S( LDST, LDST ), T( LDST, LDST ), WORK( 8 ) 00229 * .. 00230 * .. External Functions .. 00231 DOUBLE PRECISION DLAMCH 00232 EXTERNAL DLAMCH 00233 * .. 00234 * .. External Subroutines .. 00235 EXTERNAL ZLACPY, ZLARTG, ZLASSQ, ZROT 00236 * .. 00237 * .. Intrinsic Functions .. 00238 INTRINSIC ABS, DBLE, DCONJG, MAX, SQRT 00239 * .. 00240 * .. Executable Statements .. 00241 * 00242 INFO = 0 00243 * 00244 * Quick return if possible 00245 * 00246 IF( N.LE.1 ) 00247 $ RETURN 00248 * 00249 M = LDST 00250 WEAK = .FALSE. 00251 DTRONG = .FALSE. 00252 * 00253 * Make a local copy of selected block in (A, B) 00254 * 00255 CALL ZLACPY( 'Full', M, M, A( J1, J1 ), LDA, S, LDST ) 00256 CALL ZLACPY( 'Full', M, M, B( J1, J1 ), LDB, T, LDST ) 00257 * 00258 * Compute the threshold for testing the acceptance of swapping. 00259 * 00260 EPS = DLAMCH( 'P' ) 00261 SMLNUM = DLAMCH( 'S' ) / EPS 00262 SCALE = DBLE( CZERO ) 00263 SUM = DBLE( CONE ) 00264 CALL ZLACPY( 'Full', M, M, S, LDST, WORK, M ) 00265 CALL ZLACPY( 'Full', M, M, T, LDST, WORK( M*M+1 ), M ) 00266 CALL ZLASSQ( 2*M*M, WORK, 1, SCALE, SUM ) 00267 SA = SCALE*SQRT( SUM ) 00268 * 00269 * THRES has been changed from 00270 * THRESH = MAX( TEN*EPS*SA, SMLNUM ) 00271 * to 00272 * THRESH = MAX( TWENTY*EPS*SA, SMLNUM ) 00273 * on 04/01/10. 00274 * "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by 00275 * Jim Demmel and Guillaume Revy. See forum post 1783. 00276 * 00277 THRESH = MAX( TWENTY*EPS*SA, SMLNUM ) 00278 * 00279 * Compute unitary QL and RQ that swap 1-by-1 and 1-by-1 blocks 00280 * using Givens rotations and perform the swap tentatively. 00281 * 00282 F = S( 2, 2 )*T( 1, 1 ) - T( 2, 2 )*S( 1, 1 ) 00283 G = S( 2, 2 )*T( 1, 2 ) - T( 2, 2 )*S( 1, 2 ) 00284 SA = ABS( S( 2, 2 ) ) 00285 SB = ABS( T( 2, 2 ) ) 00286 CALL ZLARTG( G, F, CZ, SZ, CDUM ) 00287 SZ = -SZ 00288 CALL ZROT( 2, S( 1, 1 ), 1, S( 1, 2 ), 1, CZ, DCONJG( SZ ) ) 00289 CALL ZROT( 2, T( 1, 1 ), 1, T( 1, 2 ), 1, CZ, DCONJG( SZ ) ) 00290 IF( SA.GE.SB ) THEN 00291 CALL ZLARTG( S( 1, 1 ), S( 2, 1 ), CQ, SQ, CDUM ) 00292 ELSE 00293 CALL ZLARTG( T( 1, 1 ), T( 2, 1 ), CQ, SQ, CDUM ) 00294 END IF 00295 CALL ZROT( 2, S( 1, 1 ), LDST, S( 2, 1 ), LDST, CQ, SQ ) 00296 CALL ZROT( 2, T( 1, 1 ), LDST, T( 2, 1 ), LDST, CQ, SQ ) 00297 * 00298 * Weak stability test: |S21| + |T21| <= O(EPS F-norm((S, T))) 00299 * 00300 WS = ABS( S( 2, 1 ) ) + ABS( T( 2, 1 ) ) 00301 WEAK = WS.LE.THRESH 00302 IF( .NOT.WEAK ) 00303 $ GO TO 20 00304 * 00305 IF( WANDS ) THEN 00306 * 00307 * Strong stability test: 00308 * F-norm((A-QL**H*S*QR, B-QL**H*T*QR)) <= O(EPS*F-norm((A, B))) 00309 * 00310 CALL ZLACPY( 'Full', M, M, S, LDST, WORK, M ) 00311 CALL ZLACPY( 'Full', M, M, T, LDST, WORK( M*M+1 ), M ) 00312 CALL ZROT( 2, WORK, 1, WORK( 3 ), 1, CZ, -DCONJG( SZ ) ) 00313 CALL ZROT( 2, WORK( 5 ), 1, WORK( 7 ), 1, CZ, -DCONJG( SZ ) ) 00314 CALL ZROT( 2, WORK, 2, WORK( 2 ), 2, CQ, -SQ ) 00315 CALL ZROT( 2, WORK( 5 ), 2, WORK( 6 ), 2, CQ, -SQ ) 00316 DO 10 I = 1, 2 00317 WORK( I ) = WORK( I ) - A( J1+I-1, J1 ) 00318 WORK( I+2 ) = WORK( I+2 ) - A( J1+I-1, J1+1 ) 00319 WORK( I+4 ) = WORK( I+4 ) - B( J1+I-1, J1 ) 00320 WORK( I+6 ) = WORK( I+6 ) - B( J1+I-1, J1+1 ) 00321 10 CONTINUE 00322 SCALE = DBLE( CZERO ) 00323 SUM = DBLE( CONE ) 00324 CALL ZLASSQ( 2*M*M, WORK, 1, SCALE, SUM ) 00325 SS = SCALE*SQRT( SUM ) 00326 DTRONG = SS.LE.THRESH 00327 IF( .NOT.DTRONG ) 00328 $ GO TO 20 00329 END IF 00330 * 00331 * If the swap is accepted ("weakly" and "strongly"), apply the 00332 * equivalence transformations to the original matrix pair (A,B) 00333 * 00334 CALL ZROT( J1+1, A( 1, J1 ), 1, A( 1, J1+1 ), 1, CZ, 00335 $ DCONJG( SZ ) ) 00336 CALL ZROT( J1+1, B( 1, J1 ), 1, B( 1, J1+1 ), 1, CZ, 00337 $ DCONJG( SZ ) ) 00338 CALL ZROT( N-J1+1, A( J1, J1 ), LDA, A( J1+1, J1 ), LDA, CQ, SQ ) 00339 CALL ZROT( N-J1+1, B( J1, J1 ), LDB, B( J1+1, J1 ), LDB, CQ, SQ ) 00340 * 00341 * Set N1 by N2 (2,1) blocks to 0 00342 * 00343 A( J1+1, J1 ) = CZERO 00344 B( J1+1, J1 ) = CZERO 00345 * 00346 * Accumulate transformations into Q and Z if requested. 00347 * 00348 IF( WANTZ ) 00349 $ CALL ZROT( N, Z( 1, J1 ), 1, Z( 1, J1+1 ), 1, CZ, 00350 $ DCONJG( SZ ) ) 00351 IF( WANTQ ) 00352 $ CALL ZROT( N, Q( 1, J1 ), 1, Q( 1, J1+1 ), 1, CQ, 00353 $ DCONJG( SQ ) ) 00354 * 00355 * Exit with INFO = 0 if swap was successfully performed. 00356 * 00357 RETURN 00358 * 00359 * Exit with INFO = 1 if swap was rejected. 00360 * 00361 20 CONTINUE 00362 INFO = 1 00363 RETURN 00364 * 00365 * End of ZTGEX2 00366 * 00367 END