![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SLAG2 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SLAG2 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slag2.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slag2.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slag2.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SLAG2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, 00022 * WR2, WI ) 00023 * 00024 * .. Scalar Arguments .. 00025 * INTEGER LDA, LDB 00026 * REAL SAFMIN, SCALE1, SCALE2, WI, WR1, WR2 00027 * .. 00028 * .. Array Arguments .. 00029 * REAL A( LDA, * ), B( LDB, * ) 00030 * .. 00031 * 00032 * 00033 *> \par Purpose: 00034 * ============= 00035 *> 00036 *> \verbatim 00037 *> 00038 *> SLAG2 computes the eigenvalues of a 2 x 2 generalized eigenvalue 00039 *> problem A - w B, with scaling as necessary to avoid over-/underflow. 00040 *> 00041 *> The scaling factor "s" results in a modified eigenvalue equation 00042 *> 00043 *> s A - w B 00044 *> 00045 *> where s is a non-negative scaling factor chosen so that w, w B, 00046 *> and s A do not overflow and, if possible, do not underflow, either. 00047 *> \endverbatim 00048 * 00049 * Arguments: 00050 * ========== 00051 * 00052 *> \param[in] A 00053 *> \verbatim 00054 *> A is REAL array, dimension (LDA, 2) 00055 *> On entry, the 2 x 2 matrix A. It is assumed that its 1-norm 00056 *> is less than 1/SAFMIN. Entries less than 00057 *> sqrt(SAFMIN)*norm(A) are subject to being treated as zero. 00058 *> \endverbatim 00059 *> 00060 *> \param[in] LDA 00061 *> \verbatim 00062 *> LDA is INTEGER 00063 *> The leading dimension of the array A. LDA >= 2. 00064 *> \endverbatim 00065 *> 00066 *> \param[in] B 00067 *> \verbatim 00068 *> B is REAL array, dimension (LDB, 2) 00069 *> On entry, the 2 x 2 upper triangular matrix B. It is 00070 *> assumed that the one-norm of B is less than 1/SAFMIN. The 00071 *> diagonals should be at least sqrt(SAFMIN) times the largest 00072 *> element of B (in absolute value); if a diagonal is smaller 00073 *> than that, then +/- sqrt(SAFMIN) will be used instead of 00074 *> that diagonal. 00075 *> \endverbatim 00076 *> 00077 *> \param[in] LDB 00078 *> \verbatim 00079 *> LDB is INTEGER 00080 *> The leading dimension of the array B. LDB >= 2. 00081 *> \endverbatim 00082 *> 00083 *> \param[in] SAFMIN 00084 *> \verbatim 00085 *> SAFMIN is REAL 00086 *> The smallest positive number s.t. 1/SAFMIN does not 00087 *> overflow. (This should always be SLAMCH('S') -- it is an 00088 *> argument in order to avoid having to call SLAMCH frequently.) 00089 *> \endverbatim 00090 *> 00091 *> \param[out] SCALE1 00092 *> \verbatim 00093 *> SCALE1 is REAL 00094 *> A scaling factor used to avoid over-/underflow in the 00095 *> eigenvalue equation which defines the first eigenvalue. If 00096 *> the eigenvalues are complex, then the eigenvalues are 00097 *> ( WR1 +/- WI i ) / SCALE1 (which may lie outside the 00098 *> exponent range of the machine), SCALE1=SCALE2, and SCALE1 00099 *> will always be positive. If the eigenvalues are real, then 00100 *> the first (real) eigenvalue is WR1 / SCALE1 , but this may 00101 *> overflow or underflow, and in fact, SCALE1 may be zero or 00102 *> less than the underflow threshhold if the exact eigenvalue 00103 *> is sufficiently large. 00104 *> \endverbatim 00105 *> 00106 *> \param[out] SCALE2 00107 *> \verbatim 00108 *> SCALE2 is REAL 00109 *> A scaling factor used to avoid over-/underflow in the 00110 *> eigenvalue equation which defines the second eigenvalue. If 00111 *> the eigenvalues are complex, then SCALE2=SCALE1. If the 00112 *> eigenvalues are real, then the second (real) eigenvalue is 00113 *> WR2 / SCALE2 , but this may overflow or underflow, and in 00114 *> fact, SCALE2 may be zero or less than the underflow 00115 *> threshhold if the exact eigenvalue is sufficiently large. 00116 *> \endverbatim 00117 *> 00118 *> \param[out] WR1 00119 *> \verbatim 00120 *> WR1 is REAL 00121 *> If the eigenvalue is real, then WR1 is SCALE1 times the 00122 *> eigenvalue closest to the (2,2) element of A B**(-1). If the 00123 *> eigenvalue is complex, then WR1=WR2 is SCALE1 times the real 00124 *> part of the eigenvalues. 00125 *> \endverbatim 00126 *> 00127 *> \param[out] WR2 00128 *> \verbatim 00129 *> WR2 is REAL 00130 *> If the eigenvalue is real, then WR2 is SCALE2 times the 00131 *> other eigenvalue. If the eigenvalue is complex, then 00132 *> WR1=WR2 is SCALE1 times the real part of the eigenvalues. 00133 *> \endverbatim 00134 *> 00135 *> \param[out] WI 00136 *> \verbatim 00137 *> WI is REAL 00138 *> If the eigenvalue is real, then WI is zero. If the 00139 *> eigenvalue is complex, then WI is SCALE1 times the imaginary 00140 *> part of the eigenvalues. WI will always be non-negative. 00141 *> \endverbatim 00142 * 00143 * Authors: 00144 * ======== 00145 * 00146 *> \author Univ. of Tennessee 00147 *> \author Univ. of California Berkeley 00148 *> \author Univ. of Colorado Denver 00149 *> \author NAG Ltd. 00150 * 00151 *> \date November 2011 00152 * 00153 *> \ingroup realOTHERauxiliary 00154 * 00155 * ===================================================================== 00156 SUBROUTINE SLAG2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, 00157 $ WR2, WI ) 00158 * 00159 * -- LAPACK auxiliary routine (version 3.4.0) -- 00160 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00161 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00162 * November 2011 00163 * 00164 * .. Scalar Arguments .. 00165 INTEGER LDA, LDB 00166 REAL SAFMIN, SCALE1, SCALE2, WI, WR1, WR2 00167 * .. 00168 * .. Array Arguments .. 00169 REAL A( LDA, * ), B( LDB, * ) 00170 * .. 00171 * 00172 * ===================================================================== 00173 * 00174 * .. Parameters .. 00175 REAL ZERO, ONE, TWO 00176 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TWO = 2.0E+0 ) 00177 REAL HALF 00178 PARAMETER ( HALF = ONE / TWO ) 00179 REAL FUZZY1 00180 PARAMETER ( FUZZY1 = ONE+1.0E-5 ) 00181 * .. 00182 * .. Local Scalars .. 00183 REAL A11, A12, A21, A22, ABI22, ANORM, AS11, AS12, 00184 $ AS22, ASCALE, B11, B12, B22, BINV11, BINV22, 00185 $ BMIN, BNORM, BSCALE, BSIZE, C1, C2, C3, C4, C5, 00186 $ DIFF, DISCR, PP, QQ, R, RTMAX, RTMIN, S1, S2, 00187 $ SAFMAX, SHIFT, SS, SUM, WABS, WBIG, WDET, 00188 $ WSCALE, WSIZE, WSMALL 00189 * .. 00190 * .. Intrinsic Functions .. 00191 INTRINSIC ABS, MAX, MIN, SIGN, SQRT 00192 * .. 00193 * .. Executable Statements .. 00194 * 00195 RTMIN = SQRT( SAFMIN ) 00196 RTMAX = ONE / RTMIN 00197 SAFMAX = ONE / SAFMIN 00198 * 00199 * Scale A 00200 * 00201 ANORM = MAX( ABS( A( 1, 1 ) )+ABS( A( 2, 1 ) ), 00202 $ ABS( A( 1, 2 ) )+ABS( A( 2, 2 ) ), SAFMIN ) 00203 ASCALE = ONE / ANORM 00204 A11 = ASCALE*A( 1, 1 ) 00205 A21 = ASCALE*A( 2, 1 ) 00206 A12 = ASCALE*A( 1, 2 ) 00207 A22 = ASCALE*A( 2, 2 ) 00208 * 00209 * Perturb B if necessary to insure non-singularity 00210 * 00211 B11 = B( 1, 1 ) 00212 B12 = B( 1, 2 ) 00213 B22 = B( 2, 2 ) 00214 BMIN = RTMIN*MAX( ABS( B11 ), ABS( B12 ), ABS( B22 ), RTMIN ) 00215 IF( ABS( B11 ).LT.BMIN ) 00216 $ B11 = SIGN( BMIN, B11 ) 00217 IF( ABS( B22 ).LT.BMIN ) 00218 $ B22 = SIGN( BMIN, B22 ) 00219 * 00220 * Scale B 00221 * 00222 BNORM = MAX( ABS( B11 ), ABS( B12 )+ABS( B22 ), SAFMIN ) 00223 BSIZE = MAX( ABS( B11 ), ABS( B22 ) ) 00224 BSCALE = ONE / BSIZE 00225 B11 = B11*BSCALE 00226 B12 = B12*BSCALE 00227 B22 = B22*BSCALE 00228 * 00229 * Compute larger eigenvalue by method described by C. van Loan 00230 * 00231 * ( AS is A shifted by -SHIFT*B ) 00232 * 00233 BINV11 = ONE / B11 00234 BINV22 = ONE / B22 00235 S1 = A11*BINV11 00236 S2 = A22*BINV22 00237 IF( ABS( S1 ).LE.ABS( S2 ) ) THEN 00238 AS12 = A12 - S1*B12 00239 AS22 = A22 - S1*B22 00240 SS = A21*( BINV11*BINV22 ) 00241 ABI22 = AS22*BINV22 - SS*B12 00242 PP = HALF*ABI22 00243 SHIFT = S1 00244 ELSE 00245 AS12 = A12 - S2*B12 00246 AS11 = A11 - S2*B11 00247 SS = A21*( BINV11*BINV22 ) 00248 ABI22 = -SS*B12 00249 PP = HALF*( AS11*BINV11+ABI22 ) 00250 SHIFT = S2 00251 END IF 00252 QQ = SS*AS12 00253 IF( ABS( PP*RTMIN ).GE.ONE ) THEN 00254 DISCR = ( RTMIN*PP )**2 + QQ*SAFMIN 00255 R = SQRT( ABS( DISCR ) )*RTMAX 00256 ELSE 00257 IF( PP**2+ABS( QQ ).LE.SAFMIN ) THEN 00258 DISCR = ( RTMAX*PP )**2 + QQ*SAFMAX 00259 R = SQRT( ABS( DISCR ) )*RTMIN 00260 ELSE 00261 DISCR = PP**2 + QQ 00262 R = SQRT( ABS( DISCR ) ) 00263 END IF 00264 END IF 00265 * 00266 * Note: the test of R in the following IF is to cover the case when 00267 * DISCR is small and negative and is flushed to zero during 00268 * the calculation of R. On machines which have a consistent 00269 * flush-to-zero threshhold and handle numbers above that 00270 * threshhold correctly, it would not be necessary. 00271 * 00272 IF( DISCR.GE.ZERO .OR. R.EQ.ZERO ) THEN 00273 SUM = PP + SIGN( R, PP ) 00274 DIFF = PP - SIGN( R, PP ) 00275 WBIG = SHIFT + SUM 00276 * 00277 * Compute smaller eigenvalue 00278 * 00279 WSMALL = SHIFT + DIFF 00280 IF( HALF*ABS( WBIG ).GT.MAX( ABS( WSMALL ), SAFMIN ) ) THEN 00281 WDET = ( A11*A22-A12*A21 )*( BINV11*BINV22 ) 00282 WSMALL = WDET / WBIG 00283 END IF 00284 * 00285 * Choose (real) eigenvalue closest to 2,2 element of A*B**(-1) 00286 * for WR1. 00287 * 00288 IF( PP.GT.ABI22 ) THEN 00289 WR1 = MIN( WBIG, WSMALL ) 00290 WR2 = MAX( WBIG, WSMALL ) 00291 ELSE 00292 WR1 = MAX( WBIG, WSMALL ) 00293 WR2 = MIN( WBIG, WSMALL ) 00294 END IF 00295 WI = ZERO 00296 ELSE 00297 * 00298 * Complex eigenvalues 00299 * 00300 WR1 = SHIFT + PP 00301 WR2 = WR1 00302 WI = R 00303 END IF 00304 * 00305 * Further scaling to avoid underflow and overflow in computing 00306 * SCALE1 and overflow in computing w*B. 00307 * 00308 * This scale factor (WSCALE) is bounded from above using C1 and C2, 00309 * and from below using C3 and C4. 00310 * C1 implements the condition s A must never overflow. 00311 * C2 implements the condition w B must never overflow. 00312 * C3, with C2, 00313 * implement the condition that s A - w B must never overflow. 00314 * C4 implements the condition s should not underflow. 00315 * C5 implements the condition max(s,|w|) should be at least 2. 00316 * 00317 C1 = BSIZE*( SAFMIN*MAX( ONE, ASCALE ) ) 00318 C2 = SAFMIN*MAX( ONE, BNORM ) 00319 C3 = BSIZE*SAFMIN 00320 IF( ASCALE.LE.ONE .AND. BSIZE.LE.ONE ) THEN 00321 C4 = MIN( ONE, ( ASCALE / SAFMIN )*BSIZE ) 00322 ELSE 00323 C4 = ONE 00324 END IF 00325 IF( ASCALE.LE.ONE .OR. BSIZE.LE.ONE ) THEN 00326 C5 = MIN( ONE, ASCALE*BSIZE ) 00327 ELSE 00328 C5 = ONE 00329 END IF 00330 * 00331 * Scale first eigenvalue 00332 * 00333 WABS = ABS( WR1 ) + ABS( WI ) 00334 WSIZE = MAX( SAFMIN, C1, FUZZY1*( WABS*C2+C3 ), 00335 $ MIN( C4, HALF*MAX( WABS, C5 ) ) ) 00336 IF( WSIZE.NE.ONE ) THEN 00337 WSCALE = ONE / WSIZE 00338 IF( WSIZE.GT.ONE ) THEN 00339 SCALE1 = ( MAX( ASCALE, BSIZE )*WSCALE )* 00340 $ MIN( ASCALE, BSIZE ) 00341 ELSE 00342 SCALE1 = ( MIN( ASCALE, BSIZE )*WSCALE )* 00343 $ MAX( ASCALE, BSIZE ) 00344 END IF 00345 WR1 = WR1*WSCALE 00346 IF( WI.NE.ZERO ) THEN 00347 WI = WI*WSCALE 00348 WR2 = WR1 00349 SCALE2 = SCALE1 00350 END IF 00351 ELSE 00352 SCALE1 = ASCALE*BSIZE 00353 SCALE2 = SCALE1 00354 END IF 00355 * 00356 * Scale second eigenvalue (if real) 00357 * 00358 IF( WI.EQ.ZERO ) THEN 00359 WSIZE = MAX( SAFMIN, C1, FUZZY1*( ABS( WR2 )*C2+C3 ), 00360 $ MIN( C4, HALF*MAX( ABS( WR2 ), C5 ) ) ) 00361 IF( WSIZE.NE.ONE ) THEN 00362 WSCALE = ONE / WSIZE 00363 IF( WSIZE.GT.ONE ) THEN 00364 SCALE2 = ( MAX( ASCALE, BSIZE )*WSCALE )* 00365 $ MIN( ASCALE, BSIZE ) 00366 ELSE 00367 SCALE2 = ( MIN( ASCALE, BSIZE )*WSCALE )* 00368 $ MAX( ASCALE, BSIZE ) 00369 END IF 00370 WR2 = WR2*WSCALE 00371 ELSE 00372 SCALE2 = ASCALE*BSIZE 00373 END IF 00374 END IF 00375 * 00376 * End of SLAG2 00377 * 00378 RETURN 00379 END