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