![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CGET52 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 CGET52( LEFT, N, A, LDA, B, LDB, E, LDE, ALPHA, BETA, 00012 * WORK, RWORK, RESULT ) 00013 * 00014 * .. Scalar Arguments .. 00015 * LOGICAL LEFT 00016 * INTEGER LDA, LDB, LDE, N 00017 * .. 00018 * .. Array Arguments .. 00019 * REAL RESULT( 2 ), RWORK( * ) 00020 * COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ), 00021 * $ BETA( * ), E( LDE, * ), WORK( * ) 00022 * .. 00023 * 00024 * 00025 *> \par Purpose: 00026 * ============= 00027 *> 00028 *> \verbatim 00029 *> 00030 *> CGET52 does an eigenvector check for the generalized eigenvalue 00031 *> problem. 00032 *> 00033 *> The basic test for right eigenvectors is: 00034 *> 00035 *> | b(i) A E(i) - a(i) B E(i) | 00036 *> RESULT(1) = max ------------------------------- 00037 *> i n ulp max( |b(i) A|, |a(i) B| ) 00038 *> 00039 *> using the 1-norm. Here, a(i)/b(i) = w is the i-th generalized 00040 *> eigenvalue of A - w B, or, equivalently, b(i)/a(i) = m is the i-th 00041 *> generalized eigenvalue of m A - B. 00042 *> 00043 *> H H _ _ 00044 *> For left eigenvectors, A , B , a, and b are used. 00045 *> 00046 *> CGET52 also tests the normalization of E. Each eigenvector is 00047 *> supposed to be normalized so that the maximum "absolute value" 00048 *> of its elements is 1, where in this case, "absolute value" 00049 *> of a complex value x is |Re(x)| + |Im(x)| ; let us call this 00050 *> maximum "absolute value" norm of a vector v M(v). 00051 *> if a(i)=b(i)=0, then the eigenvector is set to be the jth coordinate 00052 *> vector. The normalization test is: 00053 *> 00054 *> RESULT(2) = max | M(v(i)) - 1 | / ( n ulp ) 00055 *> eigenvectors v(i) 00056 *> \endverbatim 00057 * 00058 * Arguments: 00059 * ========== 00060 * 00061 *> \param[in] LEFT 00062 *> \verbatim 00063 *> LEFT is LOGICAL 00064 *> =.TRUE.: The eigenvectors in the columns of E are assumed 00065 *> to be *left* eigenvectors. 00066 *> =.FALSE.: The eigenvectors in the columns of E are assumed 00067 *> to be *right* eigenvectors. 00068 *> \endverbatim 00069 *> 00070 *> \param[in] N 00071 *> \verbatim 00072 *> N is INTEGER 00073 *> The size of the matrices. If it is zero, CGET52 does 00074 *> nothing. It must be at least zero. 00075 *> \endverbatim 00076 *> 00077 *> \param[in] A 00078 *> \verbatim 00079 *> A is COMPLEX array, dimension (LDA, N) 00080 *> The matrix A. 00081 *> \endverbatim 00082 *> 00083 *> \param[in] LDA 00084 *> \verbatim 00085 *> LDA is INTEGER 00086 *> The leading dimension of A. It must be at least 1 00087 *> and at least N. 00088 *> \endverbatim 00089 *> 00090 *> \param[in] B 00091 *> \verbatim 00092 *> B is COMPLEX array, dimension (LDB, N) 00093 *> The matrix B. 00094 *> \endverbatim 00095 *> 00096 *> \param[in] LDB 00097 *> \verbatim 00098 *> LDB is INTEGER 00099 *> The leading dimension of B. It must be at least 1 00100 *> and at least N. 00101 *> \endverbatim 00102 *> 00103 *> \param[in] E 00104 *> \verbatim 00105 *> E is COMPLEX array, dimension (LDE, N) 00106 *> The matrix of eigenvectors. It must be O( 1 ). 00107 *> \endverbatim 00108 *> 00109 *> \param[in] LDE 00110 *> \verbatim 00111 *> LDE is INTEGER 00112 *> The leading dimension of E. It must be at least 1 and at 00113 *> least N. 00114 *> \endverbatim 00115 *> 00116 *> \param[in] ALPHA 00117 *> \verbatim 00118 *> ALPHA is COMPLEX array, dimension (N) 00119 *> The values a(i) as described above, which, along with b(i), 00120 *> define the generalized eigenvalues. 00121 *> \endverbatim 00122 *> 00123 *> \param[in] BETA 00124 *> \verbatim 00125 *> BETA is COMPLEX array, dimension (N) 00126 *> The values b(i) as described above, which, along with a(i), 00127 *> define the generalized eigenvalues. 00128 *> \endverbatim 00129 *> 00130 *> \param[out] WORK 00131 *> \verbatim 00132 *> WORK is COMPLEX array, dimension (N**2) 00133 *> \endverbatim 00134 *> 00135 *> \param[out] RWORK 00136 *> \verbatim 00137 *> RWORK is REAL array, dimension (N) 00138 *> \endverbatim 00139 *> 00140 *> \param[out] RESULT 00141 *> \verbatim 00142 *> RESULT is REAL array, dimension (2) 00143 *> The values computed by the test described above. If A E or 00144 *> B E is likely to overflow, then RESULT(1:2) is set to 00145 *> 10 / ulp. 00146 *> \endverbatim 00147 * 00148 * Authors: 00149 * ======== 00150 * 00151 *> \author Univ. of Tennessee 00152 *> \author Univ. of California Berkeley 00153 *> \author Univ. of Colorado Denver 00154 *> \author NAG Ltd. 00155 * 00156 *> \date November 2011 00157 * 00158 *> \ingroup complex_eig 00159 * 00160 * ===================================================================== 00161 SUBROUTINE CGET52( LEFT, N, A, LDA, B, LDB, E, LDE, ALPHA, BETA, 00162 $ WORK, RWORK, RESULT ) 00163 * 00164 * -- LAPACK test routine (version 3.4.0) -- 00165 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00166 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00167 * November 2011 00168 * 00169 * .. Scalar Arguments .. 00170 LOGICAL LEFT 00171 INTEGER LDA, LDB, LDE, N 00172 * .. 00173 * .. Array Arguments .. 00174 REAL RESULT( 2 ), RWORK( * ) 00175 COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ), 00176 $ BETA( * ), E( LDE, * ), WORK( * ) 00177 * .. 00178 * 00179 * ===================================================================== 00180 * 00181 * .. Parameters .. 00182 REAL ZERO, ONE 00183 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00184 COMPLEX CZERO, CONE 00185 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00186 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00187 * .. 00188 * .. Local Scalars .. 00189 CHARACTER NORMAB, TRANS 00190 INTEGER J, JVEC 00191 REAL ABMAX, ALFMAX, ANORM, BETMAX, BNORM, ENORM, 00192 $ ENRMER, ERRNRM, SAFMAX, SAFMIN, SCALE, TEMP1, 00193 $ ULP 00194 COMPLEX ACOEFF, ALPHAI, BCOEFF, BETAI, X 00195 * .. 00196 * .. External Functions .. 00197 REAL CLANGE, SLAMCH 00198 EXTERNAL CLANGE, SLAMCH 00199 * .. 00200 * .. External Subroutines .. 00201 EXTERNAL CGEMV 00202 * .. 00203 * .. Intrinsic Functions .. 00204 INTRINSIC ABS, AIMAG, CONJG, MAX, REAL 00205 * .. 00206 * .. Statement Functions .. 00207 REAL ABS1 00208 * .. 00209 * .. Statement Function definitions .. 00210 ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) ) 00211 * .. 00212 * .. Executable Statements .. 00213 * 00214 RESULT( 1 ) = ZERO 00215 RESULT( 2 ) = ZERO 00216 IF( N.LE.0 ) 00217 $ RETURN 00218 * 00219 SAFMIN = SLAMCH( 'Safe minimum' ) 00220 SAFMAX = ONE / SAFMIN 00221 ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' ) 00222 * 00223 IF( LEFT ) THEN 00224 TRANS = 'C' 00225 NORMAB = 'I' 00226 ELSE 00227 TRANS = 'N' 00228 NORMAB = 'O' 00229 END IF 00230 * 00231 * Norm of A, B, and E: 00232 * 00233 ANORM = MAX( CLANGE( NORMAB, N, N, A, LDA, RWORK ), SAFMIN ) 00234 BNORM = MAX( CLANGE( NORMAB, N, N, B, LDB, RWORK ), SAFMIN ) 00235 ENORM = MAX( CLANGE( 'O', N, N, E, LDE, RWORK ), ULP ) 00236 ALFMAX = SAFMAX / MAX( ONE, BNORM ) 00237 BETMAX = SAFMAX / MAX( ONE, ANORM ) 00238 * 00239 * Compute error matrix. 00240 * Column i = ( b(i) A - a(i) B ) E(i) / max( |a(i) B| |b(i) A| ) 00241 * 00242 DO 10 JVEC = 1, N 00243 ALPHAI = ALPHA( JVEC ) 00244 BETAI = BETA( JVEC ) 00245 ABMAX = MAX( ABS1( ALPHAI ), ABS1( BETAI ) ) 00246 IF( ABS1( ALPHAI ).GT.ALFMAX .OR. ABS1( BETAI ).GT.BETMAX .OR. 00247 $ ABMAX.LT.ONE ) THEN 00248 SCALE = ONE / MAX( ABMAX, SAFMIN ) 00249 ALPHAI = SCALE*ALPHAI 00250 BETAI = SCALE*BETAI 00251 END IF 00252 SCALE = ONE / MAX( ABS1( ALPHAI )*BNORM, ABS1( BETAI )*ANORM, 00253 $ SAFMIN ) 00254 ACOEFF = SCALE*BETAI 00255 BCOEFF = SCALE*ALPHAI 00256 IF( LEFT ) THEN 00257 ACOEFF = CONJG( ACOEFF ) 00258 BCOEFF = CONJG( BCOEFF ) 00259 END IF 00260 CALL CGEMV( TRANS, N, N, ACOEFF, A, LDA, E( 1, JVEC ), 1, 00261 $ CZERO, WORK( N*( JVEC-1 )+1 ), 1 ) 00262 CALL CGEMV( TRANS, N, N, -BCOEFF, B, LDA, E( 1, JVEC ), 1, 00263 $ CONE, WORK( N*( JVEC-1 )+1 ), 1 ) 00264 10 CONTINUE 00265 * 00266 ERRNRM = CLANGE( 'One', N, N, WORK, N, RWORK ) / ENORM 00267 * 00268 * Compute RESULT(1) 00269 * 00270 RESULT( 1 ) = ERRNRM / ULP 00271 * 00272 * Normalization of E: 00273 * 00274 ENRMER = ZERO 00275 DO 30 JVEC = 1, N 00276 TEMP1 = ZERO 00277 DO 20 J = 1, N 00278 TEMP1 = MAX( TEMP1, ABS1( E( J, JVEC ) ) ) 00279 20 CONTINUE 00280 ENRMER = MAX( ENRMER, TEMP1-ONE ) 00281 30 CONTINUE 00282 * 00283 * Compute RESULT(2) : the normalization error in E. 00284 * 00285 RESULT( 2 ) = ENRMER / ( REAL( N )*ULP ) 00286 * 00287 RETURN 00288 * 00289 * End of CGET52 00290 * 00291 END