1 SUBROUTINE DLAGV2( A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL, 2 $ CSR, SNR ) 3* 4* -- LAPACK auxiliary routine (version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* June 30, 1999 8* 9* .. Scalar Arguments .. 10 INTEGER LDA, LDB 11 DOUBLE PRECISION CSL, CSR, SNL, SNR 12* .. 13* .. Array Arguments .. 14 DOUBLE PRECISION A( LDA, * ), ALPHAI( 2 ), ALPHAR( 2 ), 15 $ B( LDB, * ), BETA( 2 ) 16* .. 17* 18* Purpose 19* ======= 20* 21* DLAGV2 computes the Generalized Schur factorization of a real 2-by-2 22* matrix pencil (A,B) where B is upper triangular. This routine 23* computes orthogonal (rotation) matrices given by CSL, SNL and CSR, 24* SNR such that 25* 26* 1) if the pencil (A,B) has two real eigenvalues (include 0/0 or 1/0 27* types), then 28* 29* [ a11 a12 ] := [ CSL SNL ] [ a11 a12 ] [ CSR -SNR ] 30* [ 0 a22 ] [ -SNL CSL ] [ a21 a22 ] [ SNR CSR ] 31* 32* [ b11 b12 ] := [ CSL SNL ] [ b11 b12 ] [ CSR -SNR ] 33* [ 0 b22 ] [ -SNL CSL ] [ 0 b22 ] [ SNR CSR ], 34* 35* 2) if the pencil (A,B) has a pair of complex conjugate eigenvalues, 36* then 37* 38* [ a11 a12 ] := [ CSL SNL ] [ a11 a12 ] [ CSR -SNR ] 39* [ a21 a22 ] [ -SNL CSL ] [ a21 a22 ] [ SNR CSR ] 40* 41* [ b11 0 ] := [ CSL SNL ] [ b11 b12 ] [ CSR -SNR ] 42* [ 0 b22 ] [ -SNL CSL ] [ 0 b22 ] [ SNR CSR ] 43* 44* where b11 >= b22 > 0. 45* 46* 47* Arguments 48* ========= 49* 50* A (input/output) DOUBLE PRECISION array, dimension (LDA, 2) 51* On entry, the 2 x 2 matrix A. 52* On exit, A is overwritten by the ``A-part'' of the 53* generalized Schur form. 54* 55* LDA (input) INTEGER 56* THe leading dimension of the array A. LDA >= 2. 57* 58* B (input/output) DOUBLE PRECISION array, dimension (LDB, 2) 59* On entry, the upper triangular 2 x 2 matrix B. 60* On exit, B is overwritten by the ``B-part'' of the 61* generalized Schur form. 62* 63* LDB (input) INTEGER 64* THe leading dimension of the array B. LDB >= 2. 65* 66* ALPHAR (output) DOUBLE PRECISION array, dimension (2) 67* ALPHAI (output) DOUBLE PRECISION array, dimension (2) 68* BETA (output) DOUBLE PRECISION array, dimension (2) 69* (ALPHAR(k)+i*ALPHAI(k))/BETA(k) are the eigenvalues of the 70* pencil (A,B), k=1,2, i = sqrt(-1). Note that BETA(k) may 71* be zero. 72* 73* CSL (output) DOUBLE PRECISION 74* The cosine of the left rotation matrix. 75* 76* SNL (output) DOUBLE PRECISION 77* The sine of the left rotation matrix. 78* 79* CSR (output) DOUBLE PRECISION 80* The cosine of the right rotation matrix. 81* 82* SNR (output) DOUBLE PRECISION 83* The sine of the right rotation matrix. 84* 85* Further Details 86* =============== 87* 88* Based on contributions by 89* Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA 90* 91* ===================================================================== 92* 93* .. Parameters .. 94 DOUBLE PRECISION ZERO, ONE 95 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 96* .. 97* .. Local Scalars .. 98 DOUBLE PRECISION ANORM, ASCALE, BNORM, BSCALE, H1, H2, H3, QQ, 99 $ R, RR, SAFMIN, SCALE1, SCALE2, T, ULP, WI, WR1, 100 $ WR2 101* .. 102* .. External Subroutines .. 103 EXTERNAL DLAG2, DLARTG, DLASV2, DROT 104* .. 105* .. External Functions .. 106 DOUBLE PRECISION DLAMCH, DLAPY2 107 EXTERNAL DLAMCH, DLAPY2 108* .. 109* .. Intrinsic Functions .. 110 INTRINSIC ABS, MAX 111* .. 112* .. Executable Statements .. 113* 114 SAFMIN = DLAMCH( 'S' ) 115 ULP = DLAMCH( 'P' ) 116* 117* Scale A 118* 119 ANORM = MAX( ABS( A( 1, 1 ) )+ABS( A( 2, 1 ) ), 120 $ ABS( A( 1, 2 ) )+ABS( A( 2, 2 ) ), SAFMIN ) 121 ASCALE = ONE / ANORM 122 A( 1, 1 ) = ASCALE*A( 1, 1 ) 123 A( 1, 2 ) = ASCALE*A( 1, 2 ) 124 A( 2, 1 ) = ASCALE*A( 2, 1 ) 125 A( 2, 2 ) = ASCALE*A( 2, 2 ) 126* 127* Scale B 128* 129 BNORM = MAX( ABS( B( 1, 1 ) ), ABS( B( 1, 2 ) )+ABS( B( 2, 2 ) ), 130 $ SAFMIN ) 131 BSCALE = ONE / BNORM 132 B( 1, 1 ) = BSCALE*B( 1, 1 ) 133 B( 1, 2 ) = BSCALE*B( 1, 2 ) 134 B( 2, 2 ) = BSCALE*B( 2, 2 ) 135* 136* Check if A can be deflated 137* 138 IF( ABS( A( 2, 1 ) ).LE.ULP ) THEN 139 CSL = ONE 140 SNL = ZERO 141 CSR = ONE 142 SNR = ZERO 143 A( 2, 1 ) = ZERO 144 B( 2, 1 ) = ZERO 145* 146* Check if B is singular 147* 148 ELSE IF( ABS( B( 1, 1 ) ).LE.ULP ) THEN 149 CALL DLARTG( A( 1, 1 ), A( 2, 1 ), CSL, SNL, R ) 150 CSR = ONE 151 SNR = ZERO 152 CALL DROT( 2, A( 1, 1 ), LDA, A( 2, 1 ), LDA, CSL, SNL ) 153 CALL DROT( 2, B( 1, 1 ), LDB, B( 2, 1 ), LDB, CSL, SNL ) 154 A( 2, 1 ) = ZERO 155 B( 1, 1 ) = ZERO 156 B( 2, 1 ) = ZERO 157* 158 ELSE IF( ABS( B( 2, 2 ) ).LE.ULP ) THEN 159 CALL DLARTG( A( 2, 2 ), A( 2, 1 ), CSR, SNR, T ) 160 SNR = -SNR 161 CALL DROT( 2, A( 1, 1 ), 1, A( 1, 2 ), 1, CSR, SNR ) 162 CALL DROT( 2, B( 1, 1 ), 1, B( 1, 2 ), 1, CSR, SNR ) 163 CSL = ONE 164 SNL = ZERO 165 A( 2, 1 ) = ZERO 166 B( 2, 1 ) = ZERO 167 B( 2, 2 ) = ZERO 168* 169 ELSE 170* 171* B is nonsingular, first compute the eigenvalues of (A,B) 172* 173 CALL DLAG2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, WR2, 174 $ WI ) 175* 176 IF( WI.EQ.ZERO ) THEN 177* 178* two real eigenvalues, compute s*A-w*B 179* 180 H1 = SCALE1*A( 1, 1 ) - WR1*B( 1, 1 ) 181 H2 = SCALE1*A( 1, 2 ) - WR1*B( 1, 2 ) 182 H3 = SCALE1*A( 2, 2 ) - WR1*B( 2, 2 ) 183* 184 RR = DLAPY2( H1, H2 ) 185 QQ = DLAPY2( SCALE1*A( 2, 1 ), H3 ) 186* 187 IF( RR.GT.QQ ) THEN 188* 189* find right rotation matrix to zero 1,1 element of 190* (sA - wB) 191* 192 CALL DLARTG( H2, H1, CSR, SNR, T ) 193* 194 ELSE 195* 196* find right rotation matrix to zero 2,1 element of 197* (sA - wB) 198* 199 CALL DLARTG( H3, SCALE1*A( 2, 1 ), CSR, SNR, T ) 200* 201 END IF 202* 203 SNR = -SNR 204 CALL DROT( 2, A( 1, 1 ), 1, A( 1, 2 ), 1, CSR, SNR ) 205 CALL DROT( 2, B( 1, 1 ), 1, B( 1, 2 ), 1, CSR, SNR ) 206* 207* compute inf norms of A and B 208* 209 H1 = MAX( ABS( A( 1, 1 ) )+ABS( A( 1, 2 ) ), 210 $ ABS( A( 2, 1 ) )+ABS( A( 2, 2 ) ) ) 211 H2 = MAX( ABS( B( 1, 1 ) )+ABS( B( 1, 2 ) ), 212 $ ABS( B( 2, 1 ) )+ABS( B( 2, 2 ) ) ) 213* 214 IF( ( SCALE1*H1 ).GE.ABS( WR1 )*H2 ) THEN 215* 216* find left rotation matrix Q to zero out B(2,1) 217* 218 CALL DLARTG( B( 1, 1 ), B( 2, 1 ), CSL, SNL, R ) 219* 220 ELSE 221* 222* find left rotation matrix Q to zero out A(2,1) 223* 224 CALL DLARTG( A( 1, 1 ), A( 2, 1 ), CSL, SNL, R ) 225* 226 END IF 227* 228 CALL DROT( 2, A( 1, 1 ), LDA, A( 2, 1 ), LDA, CSL, SNL ) 229 CALL DROT( 2, B( 1, 1 ), LDB, B( 2, 1 ), LDB, CSL, SNL ) 230* 231 A( 2, 1 ) = ZERO 232 B( 2, 1 ) = ZERO 233* 234 ELSE 235* 236* a pair of complex conjugate eigenvalues 237* first compute the SVD of the matrix B 238* 239 CALL DLASV2( B( 1, 1 ), B( 1, 2 ), B( 2, 2 ), R, T, SNR, 240 $ CSR, SNL, CSL ) 241* 242* Form (A,B) := Q(A,B)Z' where Q is left rotation matrix and 243* Z is right rotation matrix computed from DLASV2 244* 245 CALL DROT( 2, A( 1, 1 ), LDA, A( 2, 1 ), LDA, CSL, SNL ) 246 CALL DROT( 2, B( 1, 1 ), LDB, B( 2, 1 ), LDB, CSL, SNL ) 247 CALL DROT( 2, A( 1, 1 ), 1, A( 1, 2 ), 1, CSR, SNR ) 248 CALL DROT( 2, B( 1, 1 ), 1, B( 1, 2 ), 1, CSR, SNR ) 249* 250 B( 2, 1 ) = ZERO 251 B( 1, 2 ) = ZERO 252* 253 END IF 254* 255 END IF 256* 257* Unscaling 258* 259 A( 1, 1 ) = ANORM*A( 1, 1 ) 260 A( 2, 1 ) = ANORM*A( 2, 1 ) 261 A( 1, 2 ) = ANORM*A( 1, 2 ) 262 A( 2, 2 ) = ANORM*A( 2, 2 ) 263 B( 1, 1 ) = BNORM*B( 1, 1 ) 264 B( 2, 1 ) = BNORM*B( 2, 1 ) 265 B( 1, 2 ) = BNORM*B( 1, 2 ) 266 B( 2, 2 ) = BNORM*B( 2, 2 ) 267* 268 IF( WI.EQ.ZERO ) THEN 269 ALPHAR( 1 ) = A( 1, 1 ) 270 ALPHAR( 2 ) = A( 2, 2 ) 271 ALPHAI( 1 ) = ZERO 272 ALPHAI( 2 ) = ZERO 273 BETA( 1 ) = B( 1, 1 ) 274 BETA( 2 ) = B( 2, 2 ) 275 ELSE 276 ALPHAR( 1 ) = ANORM*WR1 / SCALE1 / BNORM 277 ALPHAI( 1 ) = ANORM*WI / SCALE1 / BNORM 278 ALPHAR( 2 ) = ALPHAR( 1 ) 279 ALPHAI( 2 ) = -ALPHAI( 1 ) 280 BETA( 1 ) = ONE 281 BETA( 2 ) = ONE 282 END IF 283* 284 10 CONTINUE 285* 286 RETURN 287* 288* End of DLAGV2 289* 290 END 291