 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DLAGV2( A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL, 00002 $ CSR, SNR ) 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 INTEGER LDA, LDB 00011 DOUBLE PRECISION CSL, CSR, SNL, SNR 00012 * .. 00013 * .. Array Arguments .. 00014 DOUBLE PRECISION A( LDA, * ), ALPHAI( 2 ), ALPHAR( 2 ), 00015 $ B( LDB, * ), BETA( 2 ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * DLAGV2 computes the Generalized Schur factorization of a real 2-by-2 00022 * matrix pencil (A,B) where B is upper triangular. This routine 00023 * computes orthogonal (rotation) matrices given by CSL, SNL and CSR, 00024 * SNR such that 00025 * 00026 * 1) if the pencil (A,B) has two real eigenvalues (include 0/0 or 1/0 00027 * types), then 00028 * 00029 * [ a11 a12 ] := [ CSL SNL ] [ a11 a12 ] [ CSR -SNR ] 00030 * [ 0 a22 ] [ -SNL CSL ] [ a21 a22 ] [ SNR CSR ] 00031 * 00032 * [ b11 b12 ] := [ CSL SNL ] [ b11 b12 ] [ CSR -SNR ] 00033 * [ 0 b22 ] [ -SNL CSL ] [ 0 b22 ] [ SNR CSR ], 00034 * 00035 * 2) if the pencil (A,B) has a pair of complex conjugate eigenvalues, 00036 * then 00037 * 00038 * [ a11 a12 ] := [ CSL SNL ] [ a11 a12 ] [ CSR -SNR ] 00039 * [ a21 a22 ] [ -SNL CSL ] [ a21 a22 ] [ SNR CSR ] 00040 * 00041 * [ b11 0 ] := [ CSL SNL ] [ b11 b12 ] [ CSR -SNR ] 00042 * [ 0 b22 ] [ -SNL CSL ] [ 0 b22 ] [ SNR CSR ] 00043 * 00044 * where b11 >= b22 > 0. 00045 * 00046 * 00047 * Arguments 00048 * ========= 00049 * 00050 * A (input/output) DOUBLE PRECISION array, dimension (LDA, 2) 00051 * On entry, the 2 x 2 matrix A. 00052 * On exit, A is overwritten by the ``A-part'' of the 00053 * generalized Schur form. 00054 * 00055 * LDA (input) INTEGER 00056 * THe leading dimension of the array A. LDA >= 2. 00057 * 00058 * B (input/output) DOUBLE PRECISION array, dimension (LDB, 2) 00059 * On entry, the upper triangular 2 x 2 matrix B. 00060 * On exit, B is overwritten by the ``B-part'' of the 00061 * generalized Schur form. 00062 * 00063 * LDB (input) INTEGER 00064 * THe leading dimension of the array B. LDB >= 2. 00065 * 00066 * ALPHAR (output) DOUBLE PRECISION array, dimension (2) 00067 * ALPHAI (output) DOUBLE PRECISION array, dimension (2) 00068 * BETA (output) DOUBLE PRECISION array, dimension (2) 00069 * (ALPHAR(k)+i*ALPHAI(k))/BETA(k) are the eigenvalues of the 00070 * pencil (A,B), k=1,2, i = sqrt(-1). Note that BETA(k) may 00071 * be zero. 00072 * 00073 * CSL (output) DOUBLE PRECISION 00074 * The cosine of the left rotation matrix. 00075 * 00076 * SNL (output) DOUBLE PRECISION 00077 * The sine of the left rotation matrix. 00078 * 00079 * CSR (output) DOUBLE PRECISION 00080 * The cosine of the right rotation matrix. 00081 * 00082 * SNR (output) DOUBLE PRECISION 00083 * The sine of the right rotation matrix. 00084 * 00085 * Further Details 00086 * =============== 00087 * 00088 * Based on contributions by 00089 * Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA 00090 * 00091 * ===================================================================== 00092 * 00093 * .. Parameters .. 00094 DOUBLE PRECISION ZERO, ONE 00095 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00096 * .. 00097 * .. Local Scalars .. 00098 DOUBLE PRECISION ANORM, ASCALE, BNORM, BSCALE, H1, H2, H3, QQ, 00099 $ R, RR, SAFMIN, SCALE1, SCALE2, T, ULP, WI, WR1, 00100 $ WR2 00101 * .. 00102 * .. External Subroutines .. 00103 EXTERNAL DLAG2, DLARTG, DLASV2, DROT 00104 * .. 00105 * .. External Functions .. 00106 DOUBLE PRECISION DLAMCH, DLAPY2 00107 EXTERNAL DLAMCH, DLAPY2 00108 * .. 00109 * .. Intrinsic Functions .. 00110 INTRINSIC ABS, MAX 00111 * .. 00112 * .. Executable Statements .. 00113 * 00114 SAFMIN = DLAMCH( 'S' ) 00115 ULP = DLAMCH( 'P' ) 00116 * 00117 * Scale A 00118 * 00119 ANORM = MAX( ABS( A( 1, 1 ) )+ABS( A( 2, 1 ) ), 00120 $ ABS( A( 1, 2 ) )+ABS( A( 2, 2 ) ), SAFMIN ) 00121 ASCALE = ONE / ANORM 00122 A( 1, 1 ) = ASCALE*A( 1, 1 ) 00123 A( 1, 2 ) = ASCALE*A( 1, 2 ) 00124 A( 2, 1 ) = ASCALE*A( 2, 1 ) 00125 A( 2, 2 ) = ASCALE*A( 2, 2 ) 00126 * 00127 * Scale B 00128 * 00129 BNORM = MAX( ABS( B( 1, 1 ) ), ABS( B( 1, 2 ) )+ABS( B( 2, 2 ) ), 00130 $ SAFMIN ) 00131 BSCALE = ONE / BNORM 00132 B( 1, 1 ) = BSCALE*B( 1, 1 ) 00133 B( 1, 2 ) = BSCALE*B( 1, 2 ) 00134 B( 2, 2 ) = BSCALE*B( 2, 2 ) 00135 * 00136 * Check if A can be deflated 00137 * 00138 IF( ABS( A( 2, 1 ) ).LE.ULP ) THEN 00139 CSL = ONE 00140 SNL = ZERO 00141 CSR = ONE 00142 SNR = ZERO 00143 A( 2, 1 ) = ZERO 00144 B( 2, 1 ) = ZERO 00145 WI = ZERO 00146 * 00147 * Check if B is singular 00148 * 00149 ELSE IF( ABS( B( 1, 1 ) ).LE.ULP ) THEN 00150 CALL DLARTG( A( 1, 1 ), A( 2, 1 ), CSL, SNL, R ) 00151 CSR = ONE 00152 SNR = ZERO 00153 CALL DROT( 2, A( 1, 1 ), LDA, A( 2, 1 ), LDA, CSL, SNL ) 00154 CALL DROT( 2, B( 1, 1 ), LDB, B( 2, 1 ), LDB, CSL, SNL ) 00155 A( 2, 1 ) = ZERO 00156 B( 1, 1 ) = ZERO 00157 B( 2, 1 ) = ZERO 00158 WI = ZERO 00159 * 00160 ELSE IF( ABS( B( 2, 2 ) ).LE.ULP ) THEN 00161 CALL DLARTG( A( 2, 2 ), A( 2, 1 ), CSR, SNR, T ) 00162 SNR = -SNR 00163 CALL DROT( 2, A( 1, 1 ), 1, A( 1, 2 ), 1, CSR, SNR ) 00164 CALL DROT( 2, B( 1, 1 ), 1, B( 1, 2 ), 1, CSR, SNR ) 00165 CSL = ONE 00166 SNL = ZERO 00167 A( 2, 1 ) = ZERO 00168 B( 2, 1 ) = ZERO 00169 B( 2, 2 ) = ZERO 00170 WI = ZERO 00171 * 00172 ELSE 00173 * 00174 * B is nonsingular, first compute the eigenvalues of (A,B) 00175 * 00176 CALL DLAG2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, WR2, 00177 $ WI ) 00178 * 00179 IF( WI.EQ.ZERO ) THEN 00180 * 00181 * two real eigenvalues, compute s*A-w*B 00182 * 00183 H1 = SCALE1*A( 1, 1 ) - WR1*B( 1, 1 ) 00184 H2 = SCALE1*A( 1, 2 ) - WR1*B( 1, 2 ) 00185 H3 = SCALE1*A( 2, 2 ) - WR1*B( 2, 2 ) 00186 * 00187 RR = DLAPY2( H1, H2 ) 00188 QQ = DLAPY2( SCALE1*A( 2, 1 ), H3 ) 00189 * 00190 IF( RR.GT.QQ ) THEN 00191 * 00192 * find right rotation matrix to zero 1,1 element of 00193 * (sA - wB) 00194 * 00195 CALL DLARTG( H2, H1, CSR, SNR, T ) 00196 * 00197 ELSE 00198 * 00199 * find right rotation matrix to zero 2,1 element of 00200 * (sA - wB) 00201 * 00202 CALL DLARTG( H3, SCALE1*A( 2, 1 ), CSR, SNR, T ) 00203 * 00204 END IF 00205 * 00206 SNR = -SNR 00207 CALL DROT( 2, A( 1, 1 ), 1, A( 1, 2 ), 1, CSR, SNR ) 00208 CALL DROT( 2, B( 1, 1 ), 1, B( 1, 2 ), 1, CSR, SNR ) 00209 * 00210 * compute inf norms of A and B 00211 * 00212 H1 = MAX( ABS( A( 1, 1 ) )+ABS( A( 1, 2 ) ), 00213 $ ABS( A( 2, 1 ) )+ABS( A( 2, 2 ) ) ) 00214 H2 = MAX( ABS( B( 1, 1 ) )+ABS( B( 1, 2 ) ), 00215 $ ABS( B( 2, 1 ) )+ABS( B( 2, 2 ) ) ) 00216 * 00217 IF( ( SCALE1*H1 ).GE.ABS( WR1 )*H2 ) THEN 00218 * 00219 * find left rotation matrix Q to zero out B(2,1) 00220 * 00221 CALL DLARTG( B( 1, 1 ), B( 2, 1 ), CSL, SNL, R ) 00222 * 00223 ELSE 00224 * 00225 * find left rotation matrix Q to zero out A(2,1) 00226 * 00227 CALL DLARTG( A( 1, 1 ), A( 2, 1 ), CSL, SNL, R ) 00228 * 00229 END IF 00230 * 00231 CALL DROT( 2, A( 1, 1 ), LDA, A( 2, 1 ), LDA, CSL, SNL ) 00232 CALL DROT( 2, B( 1, 1 ), LDB, B( 2, 1 ), LDB, CSL, SNL ) 00233 * 00234 A( 2, 1 ) = ZERO 00235 B( 2, 1 ) = ZERO 00236 * 00237 ELSE 00238 * 00239 * a pair of complex conjugate eigenvalues 00240 * first compute the SVD of the matrix B 00241 * 00242 CALL DLASV2( B( 1, 1 ), B( 1, 2 ), B( 2, 2 ), R, T, SNR, 00243 $ CSR, SNL, CSL ) 00244 * 00245 * Form (A,B) := Q(A,B)Z**T where Q is left rotation matrix and 00246 * Z is right rotation matrix computed from DLASV2 00247 * 00248 CALL DROT( 2, A( 1, 1 ), LDA, A( 2, 1 ), LDA, CSL, SNL ) 00249 CALL DROT( 2, B( 1, 1 ), LDB, B( 2, 1 ), LDB, CSL, SNL ) 00250 CALL DROT( 2, A( 1, 1 ), 1, A( 1, 2 ), 1, CSR, SNR ) 00251 CALL DROT( 2, B( 1, 1 ), 1, B( 1, 2 ), 1, CSR, SNR ) 00252 * 00253 B( 2, 1 ) = ZERO 00254 B( 1, 2 ) = ZERO 00255 * 00256 END IF 00257 * 00258 END IF 00259 * 00260 * Unscaling 00261 * 00262 A( 1, 1 ) = ANORM*A( 1, 1 ) 00263 A( 2, 1 ) = ANORM*A( 2, 1 ) 00264 A( 1, 2 ) = ANORM*A( 1, 2 ) 00265 A( 2, 2 ) = ANORM*A( 2, 2 ) 00266 B( 1, 1 ) = BNORM*B( 1, 1 ) 00267 B( 2, 1 ) = BNORM*B( 2, 1 ) 00268 B( 1, 2 ) = BNORM*B( 1, 2 ) 00269 B( 2, 2 ) = BNORM*B( 2, 2 ) 00270 * 00271 IF( WI.EQ.ZERO ) THEN 00272 ALPHAR( 1 ) = A( 1, 1 ) 00273 ALPHAR( 2 ) = A( 2, 2 ) 00274 ALPHAI( 1 ) = ZERO 00275 ALPHAI( 2 ) = ZERO 00276 BETA( 1 ) = B( 1, 1 ) 00277 BETA( 2 ) = B( 2, 2 ) 00278 ELSE 00279 ALPHAR( 1 ) = ANORM*WR1 / SCALE1 / BNORM 00280 ALPHAI( 1 ) = ANORM*WI / SCALE1 / BNORM 00281 ALPHAR( 2 ) = ALPHAR( 1 ) 00282 ALPHAI( 2 ) = -ALPHAI( 1 ) 00283 BETA( 1 ) = ONE 00284 BETA( 2 ) = ONE 00285 END IF 00286 * 00287 RETURN 00288 * 00289 * End of DLAGV2 00290 * 00291 END
1.7.3 
