1*> \brief \b SLAGV2 computes the Generalized Schur factorization of a real 2-by-2 matrix pencil (A,B) where B is upper triangular. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download SLAGV2 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slagv2.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slagv2.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slagv2.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE SLAGV2( A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL, 22* CSR, SNR ) 23* 24* .. Scalar Arguments .. 25* INTEGER LDA, LDB 26* REAL CSL, CSR, SNL, SNR 27* .. 28* .. Array Arguments .. 29* REAL A( LDA, * ), ALPHAI( 2 ), ALPHAR( 2 ), 30* $ B( LDB, * ), BETA( 2 ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 39*> SLAGV2 computes the Generalized Schur factorization of a real 2-by-2 40*> matrix pencil (A,B) where B is upper triangular. This routine 41*> computes orthogonal (rotation) matrices given by CSL, SNL and CSR, 42*> SNR such that 43*> 44*> 1) if the pencil (A,B) has two real eigenvalues (include 0/0 or 1/0 45*> types), then 46*> 47*> [ a11 a12 ] := [ CSL SNL ] [ a11 a12 ] [ CSR -SNR ] 48*> [ 0 a22 ] [ -SNL CSL ] [ a21 a22 ] [ SNR CSR ] 49*> 50*> [ b11 b12 ] := [ CSL SNL ] [ b11 b12 ] [ CSR -SNR ] 51*> [ 0 b22 ] [ -SNL CSL ] [ 0 b22 ] [ SNR CSR ], 52*> 53*> 2) if the pencil (A,B) has a pair of complex conjugate eigenvalues, 54*> then 55*> 56*> [ a11 a12 ] := [ CSL SNL ] [ a11 a12 ] [ CSR -SNR ] 57*> [ a21 a22 ] [ -SNL CSL ] [ a21 a22 ] [ SNR CSR ] 58*> 59*> [ b11 0 ] := [ CSL SNL ] [ b11 b12 ] [ CSR -SNR ] 60*> [ 0 b22 ] [ -SNL CSL ] [ 0 b22 ] [ SNR CSR ] 61*> 62*> where b11 >= b22 > 0. 63*> 64*> \endverbatim 65* 66* Arguments: 67* ========== 68* 69*> \param[in,out] A 70*> \verbatim 71*> A is REAL array, dimension (LDA, 2) 72*> On entry, the 2 x 2 matrix A. 73*> On exit, A is overwritten by the ``A-part'' of the 74*> generalized Schur form. 75*> \endverbatim 76*> 77*> \param[in] LDA 78*> \verbatim 79*> LDA is INTEGER 80*> THe leading dimension of the array A. LDA >= 2. 81*> \endverbatim 82*> 83*> \param[in,out] B 84*> \verbatim 85*> B is REAL array, dimension (LDB, 2) 86*> On entry, the upper triangular 2 x 2 matrix B. 87*> On exit, B is overwritten by the ``B-part'' of the 88*> generalized Schur form. 89*> \endverbatim 90*> 91*> \param[in] LDB 92*> \verbatim 93*> LDB is INTEGER 94*> THe leading dimension of the array B. LDB >= 2. 95*> \endverbatim 96*> 97*> \param[out] ALPHAR 98*> \verbatim 99*> ALPHAR is REAL array, dimension (2) 100*> \endverbatim 101*> 102*> \param[out] ALPHAI 103*> \verbatim 104*> ALPHAI is REAL array, dimension (2) 105*> \endverbatim 106*> 107*> \param[out] BETA 108*> \verbatim 109*> BETA is REAL array, dimension (2) 110*> (ALPHAR(k)+i*ALPHAI(k))/BETA(k) are the eigenvalues of the 111*> pencil (A,B), k=1,2, i = sqrt(-1). Note that BETA(k) may 112*> be zero. 113*> \endverbatim 114*> 115*> \param[out] CSL 116*> \verbatim 117*> CSL is REAL 118*> The cosine of the left rotation matrix. 119*> \endverbatim 120*> 121*> \param[out] SNL 122*> \verbatim 123*> SNL is REAL 124*> The sine of the left rotation matrix. 125*> \endverbatim 126*> 127*> \param[out] CSR 128*> \verbatim 129*> CSR is REAL 130*> The cosine of the right rotation matrix. 131*> \endverbatim 132*> 133*> \param[out] SNR 134*> \verbatim 135*> SNR is REAL 136*> The sine of the right rotation matrix. 137*> \endverbatim 138* 139* Authors: 140* ======== 141* 142*> \author Univ. of Tennessee 143*> \author Univ. of California Berkeley 144*> \author Univ. of Colorado Denver 145*> \author NAG Ltd. 146* 147*> \date September 2012 148* 149*> \ingroup realOTHERauxiliary 150* 151*> \par Contributors: 152* ================== 153*> 154*> Mark Fahey, Department of Mathematics, Univ. of Kentucky, USA 155* 156* ===================================================================== 157 SUBROUTINE SLAGV2( A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL, 158 $ CSR, SNR ) 159* 160* -- LAPACK auxiliary routine (version 3.4.2) -- 161* -- LAPACK is a software package provided by Univ. of Tennessee, -- 162* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 163* September 2012 164* 165* .. Scalar Arguments .. 166 INTEGER LDA, LDB 167 REAL CSL, CSR, SNL, SNR 168* .. 169* .. Array Arguments .. 170 REAL A( LDA, * ), ALPHAI( 2 ), ALPHAR( 2 ), 171 $ B( LDB, * ), BETA( 2 ) 172* .. 173* 174* ===================================================================== 175* 176* .. Parameters .. 177 REAL ZERO, ONE 178 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 179* .. 180* .. Local Scalars .. 181 REAL ANORM, ASCALE, BNORM, BSCALE, H1, H2, H3, QQ, 182 $ R, RR, SAFMIN, SCALE1, SCALE2, T, ULP, WI, WR1, 183 $ WR2 184* .. 185* .. External Subroutines .. 186 EXTERNAL SLAG2, SLARTG, SLASV2, SROT 187* .. 188* .. External Functions .. 189 REAL SLAMCH, SLAPY2 190 EXTERNAL SLAMCH, SLAPY2 191* .. 192* .. Intrinsic Functions .. 193 INTRINSIC ABS, MAX 194* .. 195* .. Executable Statements .. 196* 197 SAFMIN = SLAMCH( 'S' ) 198 ULP = SLAMCH( 'P' ) 199* 200* Scale A 201* 202 ANORM = MAX( ABS( A( 1, 1 ) )+ABS( A( 2, 1 ) ), 203 $ ABS( A( 1, 2 ) )+ABS( A( 2, 2 ) ), SAFMIN ) 204 ASCALE = ONE / ANORM 205 A( 1, 1 ) = ASCALE*A( 1, 1 ) 206 A( 1, 2 ) = ASCALE*A( 1, 2 ) 207 A( 2, 1 ) = ASCALE*A( 2, 1 ) 208 A( 2, 2 ) = ASCALE*A( 2, 2 ) 209* 210* Scale B 211* 212 BNORM = MAX( ABS( B( 1, 1 ) ), ABS( B( 1, 2 ) )+ABS( B( 2, 2 ) ), 213 $ SAFMIN ) 214 BSCALE = ONE / BNORM 215 B( 1, 1 ) = BSCALE*B( 1, 1 ) 216 B( 1, 2 ) = BSCALE*B( 1, 2 ) 217 B( 2, 2 ) = BSCALE*B( 2, 2 ) 218* 219* Check if A can be deflated 220* 221 IF( ABS( A( 2, 1 ) ).LE.ULP ) THEN 222 CSL = ONE 223 SNL = ZERO 224 CSR = ONE 225 SNR = ZERO 226 A( 2, 1 ) = ZERO 227 B( 2, 1 ) = ZERO 228 WI = ZERO 229* 230* Check if B is singular 231* 232 ELSE IF( ABS( B( 1, 1 ) ).LE.ULP ) THEN 233 CALL SLARTG( A( 1, 1 ), A( 2, 1 ), CSL, SNL, R ) 234 CSR = ONE 235 SNR = ZERO 236 CALL SROT( 2, A( 1, 1 ), LDA, A( 2, 1 ), LDA, CSL, SNL ) 237 CALL SROT( 2, B( 1, 1 ), LDB, B( 2, 1 ), LDB, CSL, SNL ) 238 A( 2, 1 ) = ZERO 239 B( 1, 1 ) = ZERO 240 B( 2, 1 ) = ZERO 241 WI = ZERO 242* 243 ELSE IF( ABS( B( 2, 2 ) ).LE.ULP ) THEN 244 CALL SLARTG( A( 2, 2 ), A( 2, 1 ), CSR, SNR, T ) 245 SNR = -SNR 246 CALL SROT( 2, A( 1, 1 ), 1, A( 1, 2 ), 1, CSR, SNR ) 247 CALL SROT( 2, B( 1, 1 ), 1, B( 1, 2 ), 1, CSR, SNR ) 248 CSL = ONE 249 SNL = ZERO 250 A( 2, 1 ) = ZERO 251 B( 2, 1 ) = ZERO 252 B( 2, 2 ) = ZERO 253 WI = ZERO 254* 255 ELSE 256* 257* B is nonsingular, first compute the eigenvalues of (A,B) 258* 259 CALL SLAG2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, WR2, 260 $ WI ) 261* 262 IF( WI.EQ.ZERO ) THEN 263* 264* two real eigenvalues, compute s*A-w*B 265* 266 H1 = SCALE1*A( 1, 1 ) - WR1*B( 1, 1 ) 267 H2 = SCALE1*A( 1, 2 ) - WR1*B( 1, 2 ) 268 H3 = SCALE1*A( 2, 2 ) - WR1*B( 2, 2 ) 269* 270 RR = SLAPY2( H1, H2 ) 271 QQ = SLAPY2( SCALE1*A( 2, 1 ), H3 ) 272* 273 IF( RR.GT.QQ ) THEN 274* 275* find right rotation matrix to zero 1,1 element of 276* (sA - wB) 277* 278 CALL SLARTG( H2, H1, CSR, SNR, T ) 279* 280 ELSE 281* 282* find right rotation matrix to zero 2,1 element of 283* (sA - wB) 284* 285 CALL SLARTG( H3, SCALE1*A( 2, 1 ), CSR, SNR, T ) 286* 287 END IF 288* 289 SNR = -SNR 290 CALL SROT( 2, A( 1, 1 ), 1, A( 1, 2 ), 1, CSR, SNR ) 291 CALL SROT( 2, B( 1, 1 ), 1, B( 1, 2 ), 1, CSR, SNR ) 292* 293* compute inf norms of A and B 294* 295 H1 = MAX( ABS( A( 1, 1 ) )+ABS( A( 1, 2 ) ), 296 $ ABS( A( 2, 1 ) )+ABS( A( 2, 2 ) ) ) 297 H2 = MAX( ABS( B( 1, 1 ) )+ABS( B( 1, 2 ) ), 298 $ ABS( B( 2, 1 ) )+ABS( B( 2, 2 ) ) ) 299* 300 IF( ( SCALE1*H1 ).GE.ABS( WR1 )*H2 ) THEN 301* 302* find left rotation matrix Q to zero out B(2,1) 303* 304 CALL SLARTG( B( 1, 1 ), B( 2, 1 ), CSL, SNL, R ) 305* 306 ELSE 307* 308* find left rotation matrix Q to zero out A(2,1) 309* 310 CALL SLARTG( A( 1, 1 ), A( 2, 1 ), CSL, SNL, R ) 311* 312 END IF 313* 314 CALL SROT( 2, A( 1, 1 ), LDA, A( 2, 1 ), LDA, CSL, SNL ) 315 CALL SROT( 2, B( 1, 1 ), LDB, B( 2, 1 ), LDB, CSL, SNL ) 316* 317 A( 2, 1 ) = ZERO 318 B( 2, 1 ) = ZERO 319* 320 ELSE 321* 322* a pair of complex conjugate eigenvalues 323* first compute the SVD of the matrix B 324* 325 CALL SLASV2( B( 1, 1 ), B( 1, 2 ), B( 2, 2 ), R, T, SNR, 326 $ CSR, SNL, CSL ) 327* 328* Form (A,B) := Q(A,B)Z**T where Q is left rotation matrix and 329* Z is right rotation matrix computed from SLASV2 330* 331 CALL SROT( 2, A( 1, 1 ), LDA, A( 2, 1 ), LDA, CSL, SNL ) 332 CALL SROT( 2, B( 1, 1 ), LDB, B( 2, 1 ), LDB, CSL, SNL ) 333 CALL SROT( 2, A( 1, 1 ), 1, A( 1, 2 ), 1, CSR, SNR ) 334 CALL SROT( 2, B( 1, 1 ), 1, B( 1, 2 ), 1, CSR, SNR ) 335* 336 B( 2, 1 ) = ZERO 337 B( 1, 2 ) = ZERO 338* 339 END IF 340* 341 END IF 342* 343* Unscaling 344* 345 A( 1, 1 ) = ANORM*A( 1, 1 ) 346 A( 2, 1 ) = ANORM*A( 2, 1 ) 347 A( 1, 2 ) = ANORM*A( 1, 2 ) 348 A( 2, 2 ) = ANORM*A( 2, 2 ) 349 B( 1, 1 ) = BNORM*B( 1, 1 ) 350 B( 2, 1 ) = BNORM*B( 2, 1 ) 351 B( 1, 2 ) = BNORM*B( 1, 2 ) 352 B( 2, 2 ) = BNORM*B( 2, 2 ) 353* 354 IF( WI.EQ.ZERO ) THEN 355 ALPHAR( 1 ) = A( 1, 1 ) 356 ALPHAR( 2 ) = A( 2, 2 ) 357 ALPHAI( 1 ) = ZERO 358 ALPHAI( 2 ) = ZERO 359 BETA( 1 ) = B( 1, 1 ) 360 BETA( 2 ) = B( 2, 2 ) 361 ELSE 362 ALPHAR( 1 ) = ANORM*WR1 / SCALE1 / BNORM 363 ALPHAI( 1 ) = ANORM*WI / SCALE1 / BNORM 364 ALPHAR( 2 ) = ALPHAR( 1 ) 365 ALPHAI( 2 ) = -ALPHAI( 1 ) 366 BETA( 1 ) = ONE 367 BETA( 2 ) = ONE 368 END IF 369* 370 RETURN 371* 372* End of SLAGV2 373* 374 END 375