00001 SUBROUTINE CTGEX2( WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, 00002 $ LDZ, J1, INFO ) 00003 * 00004 * -- LAPACK auxiliary routine (version 3.2.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 * June 2010 00008 * 00009 * .. Scalar Arguments .. 00010 LOGICAL WANTQ, WANTZ 00011 INTEGER INFO, J1, LDA, LDB, LDQ, LDZ, N 00012 * .. 00013 * .. Array Arguments .. 00014 COMPLEX A( LDA, * ), B( LDB, * ), Q( LDQ, * ), 00015 $ Z( LDZ, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * CTGEX2 swaps adjacent diagonal 1 by 1 blocks (A11,B11) and (A22,B22) 00022 * in an upper triangular matrix pair (A, B) by an unitary equivalence 00023 * transformation. 00024 * 00025 * (A, B) must be in generalized Schur canonical form, that is, A and 00026 * B are both upper triangular. 00027 * 00028 * Optionally, the matrices Q and Z of generalized Schur vectors are 00029 * updated. 00030 * 00031 * Q(in) * A(in) * Z(in)' = Q(out) * A(out) * Z(out)' 00032 * Q(in) * B(in) * Z(in)' = Q(out) * B(out) * Z(out)' 00033 * 00034 * 00035 * Arguments 00036 * ========= 00037 * 00038 * WANTQ (input) LOGICAL 00039 * .TRUE. : update the left transformation matrix Q; 00040 * .FALSE.: do not update Q. 00041 * 00042 * WANTZ (input) LOGICAL 00043 * .TRUE. : update the right transformation matrix Z; 00044 * .FALSE.: do not update Z. 00045 * 00046 * N (input) INTEGER 00047 * The order of the matrices A and B. N >= 0. 00048 * 00049 * A (input/output) COMPLEX arrays, dimensions (LDA,N) 00050 * On entry, the matrix A in the pair (A, B). 00051 * On exit, the updated matrix A. 00052 * 00053 * LDA (input) INTEGER 00054 * The leading dimension of the array A. LDA >= max(1,N). 00055 * 00056 * B (input/output) COMPLEX arrays, dimensions (LDB,N) 00057 * On entry, the matrix B in the pair (A, B). 00058 * On exit, the updated matrix B. 00059 * 00060 * LDB (input) INTEGER 00061 * The leading dimension of the array B. LDB >= max(1,N). 00062 * 00063 * Q (input/output) COMPLEX array, dimension (LDZ,N) 00064 * If WANTQ = .TRUE, on entry, the unitary matrix Q. On exit, 00065 * the updated matrix Q. 00066 * Not referenced if WANTQ = .FALSE.. 00067 * 00068 * LDQ (input) INTEGER 00069 * The leading dimension of the array Q. LDQ >= 1; 00070 * If WANTQ = .TRUE., LDQ >= N. 00071 * 00072 * Z (input/output) COMPLEX array, dimension (LDZ,N) 00073 * If WANTZ = .TRUE, on entry, the unitary matrix Z. On exit, 00074 * the updated matrix Z. 00075 * Not referenced if WANTZ = .FALSE.. 00076 * 00077 * LDZ (input) INTEGER 00078 * The leading dimension of the array Z. LDZ >= 1; 00079 * If WANTZ = .TRUE., LDZ >= N. 00080 * 00081 * J1 (input) INTEGER 00082 * The index to the first block (A11, B11). 00083 * 00084 * INFO (output) INTEGER 00085 * =0: Successful exit. 00086 * =1: The transformed matrix pair (A, B) would be too far 00087 * from generalized Schur form; the problem is ill- 00088 * conditioned. 00089 * 00090 * 00091 * Further Details 00092 * =============== 00093 * 00094 * Based on contributions by 00095 * Bo Kagstrom and Peter Poromaa, Department of Computing Science, 00096 * Umea University, S-901 87 Umea, Sweden. 00097 * 00098 * In the current code both weak and strong stability tests are 00099 * performed. The user can omit the strong stability test by changing 00100 * the internal logical parameter WANDS to .FALSE.. See ref. [2] for 00101 * details. 00102 * 00103 * [1] B. Kagstrom; A Direct Method for Reordering Eigenvalues in the 00104 * Generalized Real Schur Form of a Regular Matrix Pair (A, B), in 00105 * M.S. Moonen et al (eds), Linear Algebra for Large Scale and 00106 * Real-Time Applications, Kluwer Academic Publ. 1993, pp 195-218. 00107 * 00108 * [2] B. Kagstrom and P. Poromaa; Computing Eigenspaces with Specified 00109 * Eigenvalues of a Regular Matrix Pair (A, B) and Condition 00110 * Estimation: Theory, Algorithms and Software, Report UMINF-94.04, 00111 * Department of Computing Science, Umea University, S-901 87 Umea, 00112 * Sweden, 1994. Also as LAPACK Working Note 87. To appear in 00113 * Numerical Algorithms, 1996. 00114 * 00115 * ===================================================================== 00116 * 00117 * .. Parameters .. 00118 COMPLEX CZERO, CONE 00119 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00120 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00121 REAL TWENTY 00122 PARAMETER ( TWENTY = 2.0E+1 ) 00123 INTEGER LDST 00124 PARAMETER ( LDST = 2 ) 00125 LOGICAL WANDS 00126 PARAMETER ( WANDS = .TRUE. ) 00127 * .. 00128 * .. Local Scalars .. 00129 LOGICAL STRONG, WEAK 00130 INTEGER I, M 00131 REAL CQ, CZ, EPS, SA, SB, SCALE, SMLNUM, SS, SUM, 00132 $ THRESH, WS 00133 COMPLEX CDUM, F, G, SQ, SZ 00134 * .. 00135 * .. Local Arrays .. 00136 COMPLEX S( LDST, LDST ), T( LDST, LDST ), WORK( 8 ) 00137 * .. 00138 * .. External Functions .. 00139 REAL SLAMCH 00140 EXTERNAL SLAMCH 00141 * .. 00142 * .. External Subroutines .. 00143 EXTERNAL CLACPY, CLARTG, CLASSQ, CROT 00144 * .. 00145 * .. Intrinsic Functions .. 00146 INTRINSIC ABS, CONJG, MAX, REAL, SQRT 00147 * .. 00148 * .. Executable Statements .. 00149 * 00150 INFO = 0 00151 * 00152 * Quick return if possible 00153 * 00154 IF( N.LE.1 ) 00155 $ RETURN 00156 * 00157 M = LDST 00158 WEAK = .FALSE. 00159 STRONG = .FALSE. 00160 * 00161 * Make a local copy of selected block in (A, B) 00162 * 00163 CALL CLACPY( 'Full', M, M, A( J1, J1 ), LDA, S, LDST ) 00164 CALL CLACPY( 'Full', M, M, B( J1, J1 ), LDB, T, LDST ) 00165 * 00166 * Compute the threshold for testing the acceptance of swapping. 00167 * 00168 EPS = SLAMCH( 'P' ) 00169 SMLNUM = SLAMCH( 'S' ) / EPS 00170 SCALE = REAL( CZERO ) 00171 SUM = REAL( CONE ) 00172 CALL CLACPY( 'Full', M, M, S, LDST, WORK, M ) 00173 CALL CLACPY( 'Full', M, M, T, LDST, WORK( M*M+1 ), M ) 00174 CALL CLASSQ( 2*M*M, WORK, 1, SCALE, SUM ) 00175 SA = SCALE*SQRT( SUM ) 00176 * 00177 * THRES has been changed from 00178 * THRESH = MAX( TEN*EPS*SA, SMLNUM ) 00179 * to 00180 * THRESH = MAX( TWENTY*EPS*SA, SMLNUM ) 00181 * on 04/01/10. 00182 * "Bug" reported by Ondra Kamenik, confirmed by Julie Langou, fixed by 00183 * Jim Demmel and Guillaume Revy. See forum post 1783. 00184 * 00185 THRESH = MAX( TWENTY*EPS*SA, SMLNUM ) 00186 * 00187 * Compute unitary QL and RQ that swap 1-by-1 and 1-by-1 blocks 00188 * using Givens rotations and perform the swap tentatively. 00189 * 00190 F = S( 2, 2 )*T( 1, 1 ) - T( 2, 2 )*S( 1, 1 ) 00191 G = S( 2, 2 )*T( 1, 2 ) - T( 2, 2 )*S( 1, 2 ) 00192 SA = ABS( S( 2, 2 ) ) 00193 SB = ABS( T( 2, 2 ) ) 00194 CALL CLARTG( G, F, CZ, SZ, CDUM ) 00195 SZ = -SZ 00196 CALL CROT( 2, S( 1, 1 ), 1, S( 1, 2 ), 1, CZ, CONJG( SZ ) ) 00197 CALL CROT( 2, T( 1, 1 ), 1, T( 1, 2 ), 1, CZ, CONJG( SZ ) ) 00198 IF( SA.GE.SB ) THEN 00199 CALL CLARTG( S( 1, 1 ), S( 2, 1 ), CQ, SQ, CDUM ) 00200 ELSE 00201 CALL CLARTG( T( 1, 1 ), T( 2, 1 ), CQ, SQ, CDUM ) 00202 END IF 00203 CALL CROT( 2, S( 1, 1 ), LDST, S( 2, 1 ), LDST, CQ, SQ ) 00204 CALL CROT( 2, T( 1, 1 ), LDST, T( 2, 1 ), LDST, CQ, SQ ) 00205 * 00206 * Weak stability test: |S21| + |T21| <= O(EPS F-norm((S, T))) 00207 * 00208 WS = ABS( S( 2, 1 ) ) + ABS( T( 2, 1 ) ) 00209 WEAK = WS.LE.THRESH 00210 IF( .NOT.WEAK ) 00211 $ GO TO 20 00212 * 00213 IF( WANDS ) THEN 00214 * 00215 * Strong stability test: 00216 * F-norm((A-QL'*S*QR, B-QL'*T*QR)) <= O(EPS*F-norm((A, B))) 00217 * 00218 CALL CLACPY( 'Full', M, M, S, LDST, WORK, M ) 00219 CALL CLACPY( 'Full', M, M, T, LDST, WORK( M*M+1 ), M ) 00220 CALL CROT( 2, WORK, 1, WORK( 3 ), 1, CZ, -CONJG( SZ ) ) 00221 CALL CROT( 2, WORK( 5 ), 1, WORK( 7 ), 1, CZ, -CONJG( SZ ) ) 00222 CALL CROT( 2, WORK, 2, WORK( 2 ), 2, CQ, -SQ ) 00223 CALL CROT( 2, WORK( 5 ), 2, WORK( 6 ), 2, CQ, -SQ ) 00224 DO 10 I = 1, 2 00225 WORK( I ) = WORK( I ) - A( J1+I-1, J1 ) 00226 WORK( I+2 ) = WORK( I+2 ) - A( J1+I-1, J1+1 ) 00227 WORK( I+4 ) = WORK( I+4 ) - B( J1+I-1, J1 ) 00228 WORK( I+6 ) = WORK( I+6 ) - B( J1+I-1, J1+1 ) 00229 10 CONTINUE 00230 SCALE = REAL( CZERO ) 00231 SUM = REAL( CONE ) 00232 CALL CLASSQ( 2*M*M, WORK, 1, SCALE, SUM ) 00233 SS = SCALE*SQRT( SUM ) 00234 STRONG = SS.LE.THRESH 00235 IF( .NOT.STRONG ) 00236 $ GO TO 20 00237 END IF 00238 * 00239 * If the swap is accepted ("weakly" and "strongly"), apply the 00240 * equivalence transformations to the original matrix pair (A,B) 00241 * 00242 CALL CROT( J1+1, A( 1, J1 ), 1, A( 1, J1+1 ), 1, CZ, CONJG( SZ ) ) 00243 CALL CROT( J1+1, B( 1, J1 ), 1, B( 1, J1+1 ), 1, CZ, CONJG( SZ ) ) 00244 CALL CROT( N-J1+1, A( J1, J1 ), LDA, A( J1+1, J1 ), LDA, CQ, SQ ) 00245 CALL CROT( N-J1+1, B( J1, J1 ), LDB, B( J1+1, J1 ), LDB, CQ, SQ ) 00246 * 00247 * Set N1 by N2 (2,1) blocks to 0 00248 * 00249 A( J1+1, J1 ) = CZERO 00250 B( J1+1, J1 ) = CZERO 00251 * 00252 * Accumulate transformations into Q and Z if requested. 00253 * 00254 IF( WANTZ ) 00255 $ CALL CROT( N, Z( 1, J1 ), 1, Z( 1, J1+1 ), 1, CZ, CONJG( SZ ) ) 00256 IF( WANTQ ) 00257 $ CALL CROT( N, Q( 1, J1 ), 1, Q( 1, J1+1 ), 1, CQ, CONJG( SQ ) ) 00258 * 00259 * Exit with INFO = 0 if swap was successfully performed. 00260 * 00261 RETURN 00262 * 00263 * Exit with INFO = 1 if swap was rejected. 00264 * 00265 20 CONTINUE 00266 INFO = 1 00267 RETURN 00268 * 00269 * End of CTGEX2 00270 * 00271 END
1.7.2