![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DGET53 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 * Definition: 00009 * =========== 00010 * 00011 * SUBROUTINE DGET53( A, LDA, B, LDB, SCALE, WR, WI, RESULT, INFO ) 00012 * 00013 * .. Scalar Arguments .. 00014 * INTEGER INFO, LDA, LDB 00015 * DOUBLE PRECISION RESULT, SCALE, WI, WR 00016 * .. 00017 * .. Array Arguments .. 00018 * DOUBLE PRECISION A( LDA, * ), B( LDB, * ) 00019 * .. 00020 * 00021 * 00022 *> \par Purpose: 00023 * ============= 00024 *> 00025 *> \verbatim 00026 *> 00027 *> DGET53 checks the generalized eigenvalues computed by DLAG2. 00028 *> 00029 *> The basic test for an eigenvalue is: 00030 *> 00031 *> | det( s A - w B ) | 00032 *> RESULT = --------------------------------------------------- 00033 *> ulp max( s norm(A), |w| norm(B) )*norm( s A - w B ) 00034 *> 00035 *> Two "safety checks" are performed: 00036 *> 00037 *> (1) ulp*max( s*norm(A), |w|*norm(B) ) must be at least 00038 *> safe_minimum. This insures that the test performed is 00039 *> not essentially det(0*A + 0*B)=0. 00040 *> 00041 *> (2) s*norm(A) + |w|*norm(B) must be less than 1/safe_minimum. 00042 *> This insures that s*A - w*B will not overflow. 00043 *> 00044 *> If these tests are not passed, then s and w are scaled and 00045 *> tested anyway, if this is possible. 00046 *> \endverbatim 00047 * 00048 * Arguments: 00049 * ========== 00050 * 00051 *> \param[in] A 00052 *> \verbatim 00053 *> A is DOUBLE PRECISION array, dimension (LDA, 2) 00054 *> The 2x2 matrix A. 00055 *> \endverbatim 00056 *> 00057 *> \param[in] LDA 00058 *> \verbatim 00059 *> LDA is INTEGER 00060 *> The leading dimension of A. It must be at least 2. 00061 *> \endverbatim 00062 *> 00063 *> \param[in] B 00064 *> \verbatim 00065 *> B is DOUBLE PRECISION array, dimension (LDB, N) 00066 *> The 2x2 upper-triangular matrix B. 00067 *> \endverbatim 00068 *> 00069 *> \param[in] LDB 00070 *> \verbatim 00071 *> LDB is INTEGER 00072 *> The leading dimension of B. It must be at least 2. 00073 *> \endverbatim 00074 *> 00075 *> \param[in] SCALE 00076 *> \verbatim 00077 *> SCALE is DOUBLE PRECISION 00078 *> The "scale factor" s in the formula s A - w B . It is 00079 *> assumed to be non-negative. 00080 *> \endverbatim 00081 *> 00082 *> \param[in] WR 00083 *> \verbatim 00084 *> WR is DOUBLE PRECISION 00085 *> The real part of the eigenvalue w in the formula 00086 *> s A - w B . 00087 *> \endverbatim 00088 *> 00089 *> \param[in] WI 00090 *> \verbatim 00091 *> WI is DOUBLE PRECISION 00092 *> The imaginary part of the eigenvalue w in the formula 00093 *> s A - w B . 00094 *> \endverbatim 00095 *> 00096 *> \param[out] RESULT 00097 *> \verbatim 00098 *> RESULT is DOUBLE PRECISION 00099 *> If INFO is 2 or less, the value computed by the test 00100 *> described above. 00101 *> If INFO=3, this will just be 1/ulp. 00102 *> \endverbatim 00103 *> 00104 *> \param[out] INFO 00105 *> \verbatim 00106 *> INFO is INTEGER 00107 *> =0: The input data pass the "safety checks". 00108 *> =1: s*norm(A) + |w|*norm(B) > 1/safe_minimum. 00109 *> =2: ulp*max( s*norm(A), |w|*norm(B) ) < safe_minimum 00110 *> =3: same as INFO=2, but s and w could not be scaled so 00111 *> as to compute the test. 00112 *> \endverbatim 00113 * 00114 * Authors: 00115 * ======== 00116 * 00117 *> \author Univ. of Tennessee 00118 *> \author Univ. of California Berkeley 00119 *> \author Univ. of Colorado Denver 00120 *> \author NAG Ltd. 00121 * 00122 *> \date November 2011 00123 * 00124 *> \ingroup double_eig 00125 * 00126 * ===================================================================== 00127 SUBROUTINE DGET53( A, LDA, B, LDB, SCALE, WR, WI, RESULT, INFO ) 00128 * 00129 * -- LAPACK test routine (version 3.4.0) -- 00130 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00131 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00132 * November 2011 00133 * 00134 * .. Scalar Arguments .. 00135 INTEGER INFO, LDA, LDB 00136 DOUBLE PRECISION RESULT, SCALE, WI, WR 00137 * .. 00138 * .. Array Arguments .. 00139 DOUBLE PRECISION A( LDA, * ), B( LDB, * ) 00140 * .. 00141 * 00142 * ===================================================================== 00143 * 00144 * .. Parameters .. 00145 DOUBLE PRECISION ZERO, ONE 00146 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00147 * .. 00148 * .. Local Scalars .. 00149 DOUBLE PRECISION ABSW, ANORM, BNORM, CI11, CI12, CI22, CNORM, 00150 $ CR11, CR12, CR21, CR22, CSCALE, DETI, DETR, S1, 00151 $ SAFMIN, SCALES, SIGMIN, TEMP, ULP, WIS, WRS 00152 * .. 00153 * .. External Functions .. 00154 DOUBLE PRECISION DLAMCH 00155 EXTERNAL DLAMCH 00156 * .. 00157 * .. Intrinsic Functions .. 00158 INTRINSIC ABS, MAX, SQRT 00159 * .. 00160 * .. Executable Statements .. 00161 * 00162 * Initialize 00163 * 00164 INFO = 0 00165 RESULT = ZERO 00166 SCALES = SCALE 00167 WRS = WR 00168 WIS = WI 00169 * 00170 * Machine constants and norms 00171 * 00172 SAFMIN = DLAMCH( 'Safe minimum' ) 00173 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) 00174 ABSW = ABS( WRS ) + ABS( WIS ) 00175 ANORM = MAX( ABS( A( 1, 1 ) )+ABS( A( 2, 1 ) ), 00176 $ ABS( A( 1, 2 ) )+ABS( A( 2, 2 ) ), SAFMIN ) 00177 BNORM = MAX( ABS( B( 1, 1 ) ), ABS( B( 1, 2 ) )+ABS( B( 2, 2 ) ), 00178 $ SAFMIN ) 00179 * 00180 * Check for possible overflow. 00181 * 00182 TEMP = ( SAFMIN*BNORM )*ABSW + ( SAFMIN*ANORM )*SCALES 00183 IF( TEMP.GE.ONE ) THEN 00184 * 00185 * Scale down to avoid overflow 00186 * 00187 INFO = 1 00188 TEMP = ONE / TEMP 00189 SCALES = SCALES*TEMP 00190 WRS = WRS*TEMP 00191 WIS = WIS*TEMP 00192 ABSW = ABS( WRS ) + ABS( WIS ) 00193 END IF 00194 S1 = MAX( ULP*MAX( SCALES*ANORM, ABSW*BNORM ), 00195 $ SAFMIN*MAX( SCALES, ABSW ) ) 00196 * 00197 * Check for W and SCALE essentially zero. 00198 * 00199 IF( S1.LT.SAFMIN ) THEN 00200 INFO = 2 00201 IF( SCALES.LT.SAFMIN .AND. ABSW.LT.SAFMIN ) THEN 00202 INFO = 3 00203 RESULT = ONE / ULP 00204 RETURN 00205 END IF 00206 * 00207 * Scale up to avoid underflow 00208 * 00209 TEMP = ONE / MAX( SCALES*ANORM+ABSW*BNORM, SAFMIN ) 00210 SCALES = SCALES*TEMP 00211 WRS = WRS*TEMP 00212 WIS = WIS*TEMP 00213 ABSW = ABS( WRS ) + ABS( WIS ) 00214 S1 = MAX( ULP*MAX( SCALES*ANORM, ABSW*BNORM ), 00215 $ SAFMIN*MAX( SCALES, ABSW ) ) 00216 IF( S1.LT.SAFMIN ) THEN 00217 INFO = 3 00218 RESULT = ONE / ULP 00219 RETURN 00220 END IF 00221 END IF 00222 * 00223 * Compute C = s A - w B 00224 * 00225 CR11 = SCALES*A( 1, 1 ) - WRS*B( 1, 1 ) 00226 CI11 = -WIS*B( 1, 1 ) 00227 CR21 = SCALES*A( 2, 1 ) 00228 CR12 = SCALES*A( 1, 2 ) - WRS*B( 1, 2 ) 00229 CI12 = -WIS*B( 1, 2 ) 00230 CR22 = SCALES*A( 2, 2 ) - WRS*B( 2, 2 ) 00231 CI22 = -WIS*B( 2, 2 ) 00232 * 00233 * Compute the smallest singular value of s A - w B: 00234 * 00235 * |det( s A - w B )| 00236 * sigma_min = ------------------ 00237 * norm( s A - w B ) 00238 * 00239 CNORM = MAX( ABS( CR11 )+ABS( CI11 )+ABS( CR21 ), 00240 $ ABS( CR12 )+ABS( CI12 )+ABS( CR22 )+ABS( CI22 ), SAFMIN ) 00241 CSCALE = ONE / SQRT( CNORM ) 00242 DETR = ( CSCALE*CR11 )*( CSCALE*CR22 ) - 00243 $ ( CSCALE*CI11 )*( CSCALE*CI22 ) - 00244 $ ( CSCALE*CR12 )*( CSCALE*CR21 ) 00245 DETI = ( CSCALE*CR11 )*( CSCALE*CI22 ) + 00246 $ ( CSCALE*CI11 )*( CSCALE*CR22 ) - 00247 $ ( CSCALE*CI12 )*( CSCALE*CR21 ) 00248 SIGMIN = ABS( DETR ) + ABS( DETI ) 00249 RESULT = SIGMIN / S1 00250 RETURN 00251 * 00252 * End of DGET53 00253 * 00254 END