1*> \brief \b SLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary to avoid over-/underflow. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download SLAG2 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slag2.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slag2.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slag2.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE SLAG2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, 22* WR2, WI ) 23* 24* .. Scalar Arguments .. 25* INTEGER LDA, LDB 26* REAL SAFMIN, SCALE1, SCALE2, WI, WR1, WR2 27* .. 28* .. Array Arguments .. 29* REAL A( LDA, * ), B( LDB, * ) 30* .. 31* 32* 33*> \par Purpose: 34* ============= 35*> 36*> \verbatim 37*> 38*> SLAG2 computes the eigenvalues of a 2 x 2 generalized eigenvalue 39*> problem A - w B, with scaling as necessary to avoid over-/underflow. 40*> 41*> The scaling factor "s" results in a modified eigenvalue equation 42*> 43*> s A - w B 44*> 45*> where s is a non-negative scaling factor chosen so that w, w B, 46*> and s A do not overflow and, if possible, do not underflow, either. 47*> \endverbatim 48* 49* Arguments: 50* ========== 51* 52*> \param[in] A 53*> \verbatim 54*> A is REAL array, dimension (LDA, 2) 55*> On entry, the 2 x 2 matrix A. It is assumed that its 1-norm 56*> is less than 1/SAFMIN. Entries less than 57*> sqrt(SAFMIN)*norm(A) are subject to being treated as zero. 58*> \endverbatim 59*> 60*> \param[in] LDA 61*> \verbatim 62*> LDA is INTEGER 63*> The leading dimension of the array A. LDA >= 2. 64*> \endverbatim 65*> 66*> \param[in] B 67*> \verbatim 68*> B is REAL array, dimension (LDB, 2) 69*> On entry, the 2 x 2 upper triangular matrix B. It is 70*> assumed that the one-norm of B is less than 1/SAFMIN. The 71*> diagonals should be at least sqrt(SAFMIN) times the largest 72*> element of B (in absolute value); if a diagonal is smaller 73*> than that, then +/- sqrt(SAFMIN) will be used instead of 74*> that diagonal. 75*> \endverbatim 76*> 77*> \param[in] LDB 78*> \verbatim 79*> LDB is INTEGER 80*> The leading dimension of the array B. LDB >= 2. 81*> \endverbatim 82*> 83*> \param[in] SAFMIN 84*> \verbatim 85*> SAFMIN is REAL 86*> The smallest positive number s.t. 1/SAFMIN does not 87*> overflow. (This should always be SLAMCH('S') -- it is an 88*> argument in order to avoid having to call SLAMCH frequently.) 89*> \endverbatim 90*> 91*> \param[out] SCALE1 92*> \verbatim 93*> SCALE1 is REAL 94*> A scaling factor used to avoid over-/underflow in the 95*> eigenvalue equation which defines the first eigenvalue. If 96*> the eigenvalues are complex, then the eigenvalues are 97*> ( WR1 +/- WI i ) / SCALE1 (which may lie outside the 98*> exponent range of the machine), SCALE1=SCALE2, and SCALE1 99*> will always be positive. If the eigenvalues are real, then 100*> the first (real) eigenvalue is WR1 / SCALE1 , but this may 101*> overflow or underflow, and in fact, SCALE1 may be zero or 102*> less than the underflow threshold if the exact eigenvalue 103*> is sufficiently large. 104*> \endverbatim 105*> 106*> \param[out] SCALE2 107*> \verbatim 108*> SCALE2 is REAL 109*> A scaling factor used to avoid over-/underflow in the 110*> eigenvalue equation which defines the second eigenvalue. If 111*> the eigenvalues are complex, then SCALE2=SCALE1. If the 112*> eigenvalues are real, then the second (real) eigenvalue is 113*> WR2 / SCALE2 , but this may overflow or underflow, and in 114*> fact, SCALE2 may be zero or less than the underflow 115*> threshold if the exact eigenvalue is sufficiently large. 116*> \endverbatim 117*> 118*> \param[out] WR1 119*> \verbatim 120*> WR1 is REAL 121*> If the eigenvalue is real, then WR1 is SCALE1 times the 122*> eigenvalue closest to the (2,2) element of A B**(-1). If the 123*> eigenvalue is complex, then WR1=WR2 is SCALE1 times the real 124*> part of the eigenvalues. 125*> \endverbatim 126*> 127*> \param[out] WR2 128*> \verbatim 129*> WR2 is REAL 130*> If the eigenvalue is real, then WR2 is SCALE2 times the 131*> other eigenvalue. If the eigenvalue is complex, then 132*> WR1=WR2 is SCALE1 times the real part of the eigenvalues. 133*> \endverbatim 134*> 135*> \param[out] WI 136*> \verbatim 137*> WI is REAL 138*> If the eigenvalue is real, then WI is zero. If the 139*> eigenvalue is complex, then WI is SCALE1 times the imaginary 140*> part of the eigenvalues. WI will always be non-negative. 141*> \endverbatim 142* 143* Authors: 144* ======== 145* 146*> \author Univ. of Tennessee 147*> \author Univ. of California Berkeley 148*> \author Univ. of Colorado Denver 149*> \author NAG Ltd. 150* 151*> \ingroup realOTHERauxiliary 152* 153* ===================================================================== 154 SUBROUTINE SLAG2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, 155 $ WR2, WI ) 156* 157* -- LAPACK auxiliary routine -- 158* -- LAPACK is a software package provided by Univ. of Tennessee, -- 159* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 160* 161* .. Scalar Arguments .. 162 INTEGER LDA, LDB 163 REAL SAFMIN, SCALE1, SCALE2, WI, WR1, WR2 164* .. 165* .. Array Arguments .. 166 REAL A( LDA, * ), B( LDB, * ) 167* .. 168* 169* ===================================================================== 170* 171* .. Parameters .. 172 REAL ZERO, ONE, TWO 173 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 ) 174 REAL HALF 175 PARAMETER ( HALF = ONE / TWO ) 176 REAL FUZZY1 177 PARAMETER ( FUZZY1 = ONE+1.0E-5 ) 178* .. 179* .. Local Scalars .. 180 REAL A11, A12, A21, A22, ABI22, ANORM, AS11, AS12, 181 $ AS22, ASCALE, B11, B12, B22, BINV11, BINV22, 182 $ BMIN, BNORM, BSCALE, BSIZE, C1, C2, C3, C4, C5, 183 $ DIFF, DISCR, PP, QQ, R, RTMAX, RTMIN, S1, S2, 184 $ SAFMAX, SHIFT, SS, SUM, WABS, WBIG, WDET, 185 $ WSCALE, WSIZE, WSMALL 186* .. 187* .. Intrinsic Functions .. 188 INTRINSIC ABS, MAX, MIN, SIGN, SQRT 189* .. 190* .. Executable Statements .. 191* 192 RTMIN = SQRT( SAFMIN ) 193 RTMAX = ONE / RTMIN 194 SAFMAX = ONE / SAFMIN 195* 196* Scale A 197* 198 ANORM = MAX( ABS( A( 1, 1 ) )+ABS( A( 2, 1 ) ), 199 $ ABS( A( 1, 2 ) )+ABS( A( 2, 2 ) ), SAFMIN ) 200 ASCALE = ONE / ANORM 201 A11 = ASCALE*A( 1, 1 ) 202 A21 = ASCALE*A( 2, 1 ) 203 A12 = ASCALE*A( 1, 2 ) 204 A22 = ASCALE*A( 2, 2 ) 205* 206* Perturb B if necessary to insure non-singularity 207* 208 B11 = B( 1, 1 ) 209 B12 = B( 1, 2 ) 210 B22 = B( 2, 2 ) 211 BMIN = RTMIN*MAX( ABS( B11 ), ABS( B12 ), ABS( B22 ), RTMIN ) 212 IF( ABS( B11 ).LT.BMIN ) 213 $ B11 = SIGN( BMIN, B11 ) 214 IF( ABS( B22 ).LT.BMIN ) 215 $ B22 = SIGN( BMIN, B22 ) 216* 217* Scale B 218* 219 BNORM = MAX( ABS( B11 ), ABS( B12 )+ABS( B22 ), SAFMIN ) 220 BSIZE = MAX( ABS( B11 ), ABS( B22 ) ) 221 BSCALE = ONE / BSIZE 222 B11 = B11*BSCALE 223 B12 = B12*BSCALE 224 B22 = B22*BSCALE 225* 226* Compute larger eigenvalue by method described by C. van Loan 227* 228* ( AS is A shifted by -SHIFT*B ) 229* 230 BINV11 = ONE / B11 231 BINV22 = ONE / B22 232 S1 = A11*BINV11 233 S2 = A22*BINV22 234 IF( ABS( S1 ).LE.ABS( S2 ) ) THEN 235 AS12 = A12 - S1*B12 236 AS22 = A22 - S1*B22 237 SS = A21*( BINV11*BINV22 ) 238 ABI22 = AS22*BINV22 - SS*B12 239 PP = HALF*ABI22 240 SHIFT = S1 241 ELSE 242 AS12 = A12 - S2*B12 243 AS11 = A11 - S2*B11 244 SS = A21*( BINV11*BINV22 ) 245 ABI22 = -SS*B12 246 PP = HALF*( AS11*BINV11+ABI22 ) 247 SHIFT = S2 248 END IF 249 QQ = SS*AS12 250 IF( ABS( PP*RTMIN ).GE.ONE ) THEN 251 DISCR = ( RTMIN*PP )**2 + QQ*SAFMIN 252 R = SQRT( ABS( DISCR ) )*RTMAX 253 ELSE 254 IF( PP**2+ABS( QQ ).LE.SAFMIN ) THEN 255 DISCR = ( RTMAX*PP )**2 + QQ*SAFMAX 256 R = SQRT( ABS( DISCR ) )*RTMIN 257 ELSE 258 DISCR = PP**2 + QQ 259 R = SQRT( ABS( DISCR ) ) 260 END IF 261 END IF 262* 263* Note: the test of R in the following IF is to cover the case when 264* DISCR is small and negative and is flushed to zero during 265* the calculation of R. On machines which have a consistent 266* flush-to-zero threshold and handle numbers above that 267* threshold correctly, it would not be necessary. 268* 269 IF( DISCR.GE.ZERO .OR. R.EQ.ZERO ) THEN 270 SUM = PP + SIGN( R, PP ) 271 DIFF = PP - SIGN( R, PP ) 272 WBIG = SHIFT + SUM 273* 274* Compute smaller eigenvalue 275* 276 WSMALL = SHIFT + DIFF 277 IF( HALF*ABS( WBIG ).GT.MAX( ABS( WSMALL ), SAFMIN ) ) THEN 278 WDET = ( A11*A22-A12*A21 )*( BINV11*BINV22 ) 279 WSMALL = WDET / WBIG 280 END IF 281* 282* Choose (real) eigenvalue closest to 2,2 element of A*B**(-1) 283* for WR1. 284* 285 IF( PP.GT.ABI22 ) THEN 286 WR1 = MIN( WBIG, WSMALL ) 287 WR2 = MAX( WBIG, WSMALL ) 288 ELSE 289 WR1 = MAX( WBIG, WSMALL ) 290 WR2 = MIN( WBIG, WSMALL ) 291 END IF 292 WI = ZERO 293 ELSE 294* 295* Complex eigenvalues 296* 297 WR1 = SHIFT + PP 298 WR2 = WR1 299 WI = R 300 END IF 301* 302* Further scaling to avoid underflow and overflow in computing 303* SCALE1 and overflow in computing w*B. 304* 305* This scale factor (WSCALE) is bounded from above using C1 and C2, 306* and from below using C3 and C4. 307* C1 implements the condition s A must never overflow. 308* C2 implements the condition w B must never overflow. 309* C3, with C2, 310* implement the condition that s A - w B must never overflow. 311* C4 implements the condition s should not underflow. 312* C5 implements the condition max(s,|w|) should be at least 2. 313* 314 C1 = BSIZE*( SAFMIN*MAX( ONE, ASCALE ) ) 315 C2 = SAFMIN*MAX( ONE, BNORM ) 316 C3 = BSIZE*SAFMIN 317 IF( ASCALE.LE.ONE .AND. BSIZE.LE.ONE ) THEN 318 C4 = MIN( ONE, ( ASCALE / SAFMIN )*BSIZE ) 319 ELSE 320 C4 = ONE 321 END IF 322 IF( ASCALE.LE.ONE .OR. BSIZE.LE.ONE ) THEN 323 C5 = MIN( ONE, ASCALE*BSIZE ) 324 ELSE 325 C5 = ONE 326 END IF 327* 328* Scale first eigenvalue 329* 330 WABS = ABS( WR1 ) + ABS( WI ) 331 WSIZE = MAX( SAFMIN, C1, FUZZY1*( WABS*C2+C3 ), 332 $ MIN( C4, HALF*MAX( WABS, C5 ) ) ) 333 IF( WSIZE.NE.ONE ) THEN 334 WSCALE = ONE / WSIZE 335 IF( WSIZE.GT.ONE ) THEN 336 SCALE1 = ( MAX( ASCALE, BSIZE )*WSCALE )* 337 $ MIN( ASCALE, BSIZE ) 338 ELSE 339 SCALE1 = ( MIN( ASCALE, BSIZE )*WSCALE )* 340 $ MAX( ASCALE, BSIZE ) 341 END IF 342 WR1 = WR1*WSCALE 343 IF( WI.NE.ZERO ) THEN 344 WI = WI*WSCALE 345 WR2 = WR1 346 SCALE2 = SCALE1 347 END IF 348 ELSE 349 SCALE1 = ASCALE*BSIZE 350 SCALE2 = SCALE1 351 END IF 352* 353* Scale second eigenvalue (if real) 354* 355 IF( WI.EQ.ZERO ) THEN 356 WSIZE = MAX( SAFMIN, C1, FUZZY1*( ABS( WR2 )*C2+C3 ), 357 $ MIN( C4, HALF*MAX( ABS( WR2 ), C5 ) ) ) 358 IF( WSIZE.NE.ONE ) THEN 359 WSCALE = ONE / WSIZE 360 IF( WSIZE.GT.ONE ) THEN 361 SCALE2 = ( MAX( ASCALE, BSIZE )*WSCALE )* 362 $ MIN( ASCALE, BSIZE ) 363 ELSE 364 SCALE2 = ( MIN( ASCALE, BSIZE )*WSCALE )* 365 $ MAX( ASCALE, BSIZE ) 366 END IF 367 WR2 = WR2*WSCALE 368 ELSE 369 SCALE2 = ASCALE*BSIZE 370 END IF 371 END IF 372* 373* End of SLAG2 374* 375 RETURN 376 END 377