 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DLAG2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, 00002 $ WR2, WI ) 00003 * 00004 * -- LAPACK auxiliary routine (version 3.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 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 INTEGER LDA, LDB 00011 DOUBLE PRECISION SAFMIN, SCALE1, SCALE2, WI, WR1, WR2 00012 * .. 00013 * .. Array Arguments .. 00014 DOUBLE PRECISION A( LDA, * ), B( LDB, * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * DLAG2 computes the eigenvalues of a 2 x 2 generalized eigenvalue 00021 * problem A - w B, with scaling as necessary to avoid over-/underflow. 00022 * 00023 * The scaling factor "s" results in a modified eigenvalue equation 00024 * 00025 * s A - w B 00026 * 00027 * where s is a non-negative scaling factor chosen so that w, w B, 00028 * and s A do not overflow and, if possible, do not underflow, either. 00029 * 00030 * Arguments 00031 * ========= 00032 * 00033 * A (input) DOUBLE PRECISION array, dimension (LDA, 2) 00034 * On entry, the 2 x 2 matrix A. It is assumed that its 1-norm 00035 * is less than 1/SAFMIN. Entries less than 00036 * sqrt(SAFMIN)*norm(A) are subject to being treated as zero. 00037 * 00038 * LDA (input) INTEGER 00039 * The leading dimension of the array A. LDA >= 2. 00040 * 00041 * B (input) DOUBLE PRECISION array, dimension (LDB, 2) 00042 * On entry, the 2 x 2 upper triangular matrix B. It is 00043 * assumed that the one-norm of B is less than 1/SAFMIN. The 00044 * diagonals should be at least sqrt(SAFMIN) times the largest 00045 * element of B (in absolute value); if a diagonal is smaller 00046 * than that, then +/- sqrt(SAFMIN) will be used instead of 00047 * that diagonal. 00048 * 00049 * LDB (input) INTEGER 00050 * The leading dimension of the array B. LDB >= 2. 00051 * 00052 * SAFMIN (input) DOUBLE PRECISION 00053 * The smallest positive number s.t. 1/SAFMIN does not 00054 * overflow. (This should always be DLAMCH('S') -- it is an 00055 * argument in order to avoid having to call DLAMCH frequently.) 00056 * 00057 * SCALE1 (output) DOUBLE PRECISION 00058 * A scaling factor used to avoid over-/underflow in the 00059 * eigenvalue equation which defines the first eigenvalue. If 00060 * the eigenvalues are complex, then the eigenvalues are 00061 * ( WR1 +/- WI i ) / SCALE1 (which may lie outside the 00062 * exponent range of the machine), SCALE1=SCALE2, and SCALE1 00063 * will always be positive. If the eigenvalues are real, then 00064 * the first (real) eigenvalue is WR1 / SCALE1 , but this may 00065 * overflow or underflow, and in fact, SCALE1 may be zero or 00066 * less than the underflow threshhold if the exact eigenvalue 00067 * is sufficiently large. 00068 * 00069 * SCALE2 (output) DOUBLE PRECISION 00070 * A scaling factor used to avoid over-/underflow in the 00071 * eigenvalue equation which defines the second eigenvalue. If 00072 * the eigenvalues are complex, then SCALE2=SCALE1. If the 00073 * eigenvalues are real, then the second (real) eigenvalue is 00074 * WR2 / SCALE2 , but this may overflow or underflow, and in 00075 * fact, SCALE2 may be zero or less than the underflow 00076 * threshhold if the exact eigenvalue is sufficiently large. 00077 * 00078 * WR1 (output) DOUBLE PRECISION 00079 * If the eigenvalue is real, then WR1 is SCALE1 times the 00080 * eigenvalue closest to the (2,2) element of A B**(-1). If the 00081 * eigenvalue is complex, then WR1=WR2 is SCALE1 times the real 00082 * part of the eigenvalues. 00083 * 00084 * WR2 (output) DOUBLE PRECISION 00085 * If the eigenvalue is real, then WR2 is SCALE2 times the 00086 * other eigenvalue. If the eigenvalue is complex, then 00087 * WR1=WR2 is SCALE1 times the real part of the eigenvalues. 00088 * 00089 * WI (output) DOUBLE PRECISION 00090 * If the eigenvalue is real, then WI is zero. If the 00091 * eigenvalue is complex, then WI is SCALE1 times the imaginary 00092 * part of the eigenvalues. WI will always be non-negative. 00093 * 00094 * ===================================================================== 00095 * 00096 * .. Parameters .. 00097 DOUBLE PRECISION ZERO, ONE, TWO 00098 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 ) 00099 DOUBLE PRECISION HALF 00100 PARAMETER ( HALF = ONE / TWO ) 00101 DOUBLE PRECISION FUZZY1 00102 PARAMETER ( FUZZY1 = ONE+1.0D-5 ) 00103 * .. 00104 * .. Local Scalars .. 00105 DOUBLE PRECISION A11, A12, A21, A22, ABI22, ANORM, AS11, AS12, 00106 $ AS22, ASCALE, B11, B12, B22, BINV11, BINV22, 00107 $ BMIN, BNORM, BSCALE, BSIZE, C1, C2, C3, C4, C5, 00108 $ DIFF, DISCR, PP, QQ, R, RTMAX, RTMIN, S1, S2, 00109 $ SAFMAX, SHIFT, SS, SUM, WABS, WBIG, WDET, 00110 $ WSCALE, WSIZE, WSMALL 00111 * .. 00112 * .. Intrinsic Functions .. 00113 INTRINSIC ABS, MAX, MIN, SIGN, SQRT 00114 * .. 00115 * .. Executable Statements .. 00116 * 00117 RTMIN = SQRT( SAFMIN ) 00118 RTMAX = ONE / RTMIN 00119 SAFMAX = ONE / SAFMIN 00120 * 00121 * Scale A 00122 * 00123 ANORM = MAX( ABS( A( 1, 1 ) )+ABS( A( 2, 1 ) ), 00124 $ ABS( A( 1, 2 ) )+ABS( A( 2, 2 ) ), SAFMIN ) 00125 ASCALE = ONE / ANORM 00126 A11 = ASCALE*A( 1, 1 ) 00127 A21 = ASCALE*A( 2, 1 ) 00128 A12 = ASCALE*A( 1, 2 ) 00129 A22 = ASCALE*A( 2, 2 ) 00130 * 00131 * Perturb B if necessary to insure non-singularity 00132 * 00133 B11 = B( 1, 1 ) 00134 B12 = B( 1, 2 ) 00135 B22 = B( 2, 2 ) 00136 BMIN = RTMIN*MAX( ABS( B11 ), ABS( B12 ), ABS( B22 ), RTMIN ) 00137 IF( ABS( B11 ).LT.BMIN ) 00138 $ B11 = SIGN( BMIN, B11 ) 00139 IF( ABS( B22 ).LT.BMIN ) 00140 $ B22 = SIGN( BMIN, B22 ) 00141 * 00142 * Scale B 00143 * 00144 BNORM = MAX( ABS( B11 ), ABS( B12 )+ABS( B22 ), SAFMIN ) 00145 BSIZE = MAX( ABS( B11 ), ABS( B22 ) ) 00146 BSCALE = ONE / BSIZE 00147 B11 = B11*BSCALE 00148 B12 = B12*BSCALE 00149 B22 = B22*BSCALE 00150 * 00151 * Compute larger eigenvalue by method described by C. van Loan 00152 * 00153 * ( AS is A shifted by -SHIFT*B ) 00154 * 00155 BINV11 = ONE / B11 00156 BINV22 = ONE / B22 00157 S1 = A11*BINV11 00158 S2 = A22*BINV22 00159 IF( ABS( S1 ).LE.ABS( S2 ) ) THEN 00160 AS12 = A12 - S1*B12 00161 AS22 = A22 - S1*B22 00162 SS = A21*( BINV11*BINV22 ) 00163 ABI22 = AS22*BINV22 - SS*B12 00164 PP = HALF*ABI22 00165 SHIFT = S1 00166 ELSE 00167 AS12 = A12 - S2*B12 00168 AS11 = A11 - S2*B11 00169 SS = A21*( BINV11*BINV22 ) 00170 ABI22 = -SS*B12 00171 PP = HALF*( AS11*BINV11+ABI22 ) 00172 SHIFT = S2 00173 END IF 00174 QQ = SS*AS12 00175 IF( ABS( PP*RTMIN ).GE.ONE ) THEN 00176 DISCR = ( RTMIN*PP )**2 + QQ*SAFMIN 00177 R = SQRT( ABS( DISCR ) )*RTMAX 00178 ELSE 00179 IF( PP**2+ABS( QQ ).LE.SAFMIN ) THEN 00180 DISCR = ( RTMAX*PP )**2 + QQ*SAFMAX 00181 R = SQRT( ABS( DISCR ) )*RTMIN 00182 ELSE 00183 DISCR = PP**2 + QQ 00184 R = SQRT( ABS( DISCR ) ) 00185 END IF 00186 END IF 00187 * 00188 * Note: the test of R in the following IF is to cover the case when 00189 * DISCR is small and negative and is flushed to zero during 00190 * the calculation of R. On machines which have a consistent 00191 * flush-to-zero threshhold and handle numbers above that 00192 * threshhold correctly, it would not be necessary. 00193 * 00194 IF( DISCR.GE.ZERO .OR. R.EQ.ZERO ) THEN 00195 SUM = PP + SIGN( R, PP ) 00196 DIFF = PP - SIGN( R, PP ) 00197 WBIG = SHIFT + SUM 00198 * 00199 * Compute smaller eigenvalue 00200 * 00201 WSMALL = SHIFT + DIFF 00202 IF( HALF*ABS( WBIG ).GT.MAX( ABS( WSMALL ), SAFMIN ) ) THEN 00203 WDET = ( A11*A22-A12*A21 )*( BINV11*BINV22 ) 00204 WSMALL = WDET / WBIG 00205 END IF 00206 * 00207 * Choose (real) eigenvalue closest to 2,2 element of A*B**(-1) 00208 * for WR1. 00209 * 00210 IF( PP.GT.ABI22 ) THEN 00211 WR1 = MIN( WBIG, WSMALL ) 00212 WR2 = MAX( WBIG, WSMALL ) 00213 ELSE 00214 WR1 = MAX( WBIG, WSMALL ) 00215 WR2 = MIN( WBIG, WSMALL ) 00216 END IF 00217 WI = ZERO 00218 ELSE 00219 * 00220 * Complex eigenvalues 00221 * 00222 WR1 = SHIFT + PP 00223 WR2 = WR1 00224 WI = R 00225 END IF 00226 * 00227 * Further scaling to avoid underflow and overflow in computing 00228 * SCALE1 and overflow in computing w*B. 00229 * 00230 * This scale factor (WSCALE) is bounded from above using C1 and C2, 00231 * and from below using C3 and C4. 00232 * C1 implements the condition s A must never overflow. 00233 * C2 implements the condition w B must never overflow. 00234 * C3, with C2, 00235 * implement the condition that s A - w B must never overflow. 00236 * C4 implements the condition s should not underflow. 00237 * C5 implements the condition max(s,|w|) should be at least 2. 00238 * 00239 C1 = BSIZE*( SAFMIN*MAX( ONE, ASCALE ) ) 00240 C2 = SAFMIN*MAX( ONE, BNORM ) 00241 C3 = BSIZE*SAFMIN 00242 IF( ASCALE.LE.ONE .AND. BSIZE.LE.ONE ) THEN 00243 C4 = MIN( ONE, ( ASCALE / SAFMIN )*BSIZE ) 00244 ELSE 00245 C4 = ONE 00246 END IF 00247 IF( ASCALE.LE.ONE .OR. BSIZE.LE.ONE ) THEN 00248 C5 = MIN( ONE, ASCALE*BSIZE ) 00249 ELSE 00250 C5 = ONE 00251 END IF 00252 * 00253 * Scale first eigenvalue 00254 * 00255 WABS = ABS( WR1 ) + ABS( WI ) 00256 WSIZE = MAX( SAFMIN, C1, FUZZY1*( WABS*C2+C3 ), 00257 $ MIN( C4, HALF*MAX( WABS, C5 ) ) ) 00258 IF( WSIZE.NE.ONE ) THEN 00259 WSCALE = ONE / WSIZE 00260 IF( WSIZE.GT.ONE ) THEN 00261 SCALE1 = ( MAX( ASCALE, BSIZE )*WSCALE )* 00262 $ MIN( ASCALE, BSIZE ) 00263 ELSE 00264 SCALE1 = ( MIN( ASCALE, BSIZE )*WSCALE )* 00265 $ MAX( ASCALE, BSIZE ) 00266 END IF 00267 WR1 = WR1*WSCALE 00268 IF( WI.NE.ZERO ) THEN 00269 WI = WI*WSCALE 00270 WR2 = WR1 00271 SCALE2 = SCALE1 00272 END IF 00273 ELSE 00274 SCALE1 = ASCALE*BSIZE 00275 SCALE2 = SCALE1 00276 END IF 00277 * 00278 * Scale second eigenvalue (if real) 00279 * 00280 IF( WI.EQ.ZERO ) THEN 00281 WSIZE = MAX( SAFMIN, C1, FUZZY1*( ABS( WR2 )*C2+C3 ), 00282 $ MIN( C4, HALF*MAX( ABS( WR2 ), C5 ) ) ) 00283 IF( WSIZE.NE.ONE ) THEN 00284 WSCALE = ONE / WSIZE 00285 IF( WSIZE.GT.ONE ) THEN 00286 SCALE2 = ( MAX( ASCALE, BSIZE )*WSCALE )* 00287 $ MIN( ASCALE, BSIZE ) 00288 ELSE 00289 SCALE2 = ( MIN( ASCALE, BSIZE )*WSCALE )* 00290 $ MAX( ASCALE, BSIZE ) 00291 END IF 00292 WR2 = WR2*WSCALE 00293 ELSE 00294 SCALE2 = ASCALE*BSIZE 00295 END IF 00296 END IF 00297 * 00298 * End of DLAG2 00299 * 00300 RETURN 00301 END
1.7.3 
