![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DGET52 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 DGET52( LEFT, N, A, LDA, B, LDB, E, LDE, ALPHAR, 00012 * ALPHAI, BETA, WORK, RESULT ) 00013 * 00014 * .. Scalar Arguments .. 00015 * LOGICAL LEFT 00016 * INTEGER LDA, LDB, LDE, N 00017 * .. 00018 * .. Array Arguments .. 00019 * DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ), 00020 * $ B( LDB, * ), BETA( * ), E( LDE, * ), 00021 * $ RESULT( 2 ), WORK( * ) 00022 * .. 00023 * 00024 * 00025 *> \par Purpose: 00026 * ============= 00027 *> 00028 *> \verbatim 00029 *> 00030 *> DGET52 does an eigenvector check for the generalized eigenvalue 00031 *> problem. 00032 *> 00033 *> The basic test for right eigenvectors is: 00034 *> 00035 *> | b(j) A E(j) - a(j) B E(j) | 00036 *> RESULT(1) = max ------------------------------- 00037 *> j n ulp max( |b(j) A|, |a(j) B| ) 00038 *> 00039 *> using the 1-norm. Here, a(j)/b(j) = w is the j-th generalized 00040 *> eigenvalue of A - w B, or, equivalently, b(j)/a(j) = m is the j-th 00041 *> generalized eigenvalue of m A - B. 00042 *> 00043 *> For real eigenvalues, the test is straightforward. For complex 00044 *> eigenvalues, E(j) and a(j) are complex, represented by 00045 *> Er(j) + i*Ei(j) and ar(j) + i*ai(j), resp., so the test for that 00046 *> eigenvector becomes 00047 *> 00048 *> max( |Wr|, |Wi| ) 00049 *> -------------------------------------------- 00050 *> n ulp max( |b(j) A|, (|ar(j)|+|ai(j)|) |B| ) 00051 *> 00052 *> where 00053 *> 00054 *> Wr = b(j) A Er(j) - ar(j) B Er(j) + ai(j) B Ei(j) 00055 *> 00056 *> Wi = b(j) A Ei(j) - ai(j) B Er(j) - ar(j) B Ei(j) 00057 *> 00058 *> T T _ 00059 *> For left eigenvectors, A , B , a, and b are used. 00060 *> 00061 *> DGET52 also tests the normalization of E. Each eigenvector is 00062 *> supposed to be normalized so that the maximum "absolute value" 00063 *> of its elements is 1, where in this case, "absolute value" 00064 *> of a complex value x is |Re(x)| + |Im(x)| ; let us call this 00065 *> maximum "absolute value" norm of a vector v M(v). 00066 *> if a(j)=b(j)=0, then the eigenvector is set to be the jth coordinate 00067 *> vector. The normalization test is: 00068 *> 00069 *> RESULT(2) = max | M(v(j)) - 1 | / ( n ulp ) 00070 *> eigenvectors v(j) 00071 *> \endverbatim 00072 * 00073 * Arguments: 00074 * ========== 00075 * 00076 *> \param[in] LEFT 00077 *> \verbatim 00078 *> LEFT is LOGICAL 00079 *> =.TRUE.: The eigenvectors in the columns of E are assumed 00080 *> to be *left* eigenvectors. 00081 *> =.FALSE.: The eigenvectors in the columns of E are assumed 00082 *> to be *right* eigenvectors. 00083 *> \endverbatim 00084 *> 00085 *> \param[in] N 00086 *> \verbatim 00087 *> N is INTEGER 00088 *> The size of the matrices. If it is zero, DGET52 does 00089 *> nothing. It must be at least zero. 00090 *> \endverbatim 00091 *> 00092 *> \param[in] A 00093 *> \verbatim 00094 *> A is DOUBLE PRECISION array, dimension (LDA, N) 00095 *> The matrix A. 00096 *> \endverbatim 00097 *> 00098 *> \param[in] LDA 00099 *> \verbatim 00100 *> LDA is INTEGER 00101 *> The leading dimension of A. It must be at least 1 00102 *> and at least N. 00103 *> \endverbatim 00104 *> 00105 *> \param[in] B 00106 *> \verbatim 00107 *> B is DOUBLE PRECISION array, dimension (LDB, N) 00108 *> The matrix B. 00109 *> \endverbatim 00110 *> 00111 *> \param[in] LDB 00112 *> \verbatim 00113 *> LDB is INTEGER 00114 *> The leading dimension of B. It must be at least 1 00115 *> and at least N. 00116 *> \endverbatim 00117 *> 00118 *> \param[in] E 00119 *> \verbatim 00120 *> E is DOUBLE PRECISION array, dimension (LDE, N) 00121 *> The matrix of eigenvectors. It must be O( 1 ). Complex 00122 *> eigenvalues and eigenvectors always come in pairs, the 00123 *> eigenvalue and its conjugate being stored in adjacent 00124 *> elements of ALPHAR, ALPHAI, and BETA. Thus, if a(j)/b(j) 00125 *> and a(j+1)/b(j+1) are a complex conjugate pair of 00126 *> generalized eigenvalues, then E(,j) contains the real part 00127 *> of the eigenvector and E(,j+1) contains the imaginary part. 00128 *> Note that whether E(,j) is a real eigenvector or part of a 00129 *> complex one is specified by whether ALPHAI(j) is zero or not. 00130 *> \endverbatim 00131 *> 00132 *> \param[in] LDE 00133 *> \verbatim 00134 *> LDE is INTEGER 00135 *> The leading dimension of E. It must be at least 1 and at 00136 *> least N. 00137 *> \endverbatim 00138 *> 00139 *> \param[in] ALPHAR 00140 *> \verbatim 00141 *> ALPHAR is DOUBLE PRECISION array, dimension (N) 00142 *> The real parts of the values a(j) as described above, which, 00143 *> along with b(j), define the generalized eigenvalues. 00144 *> Complex eigenvalues always come in complex conjugate pairs 00145 *> a(j)/b(j) and a(j+1)/b(j+1), which are stored in adjacent 00146 *> elements in ALPHAR, ALPHAI, and BETA. Thus, if the j-th 00147 *> and (j+1)-st eigenvalues form a pair, ALPHAR(j+1)/BETA(j+1) 00148 *> is assumed to be equal to ALPHAR(j)/BETA(j). 00149 *> \endverbatim 00150 *> 00151 *> \param[in] ALPHAI 00152 *> \verbatim 00153 *> ALPHAI is DOUBLE PRECISION array, dimension (N) 00154 *> The imaginary parts of the values a(j) as described above, 00155 *> which, along with b(j), define the generalized eigenvalues. 00156 *> If ALPHAI(j)=0, then the eigenvalue is real, otherwise it 00157 *> is part of a complex conjugate pair. Complex eigenvalues 00158 *> always come in complex conjugate pairs a(j)/b(j) and 00159 *> a(j+1)/b(j+1), which are stored in adjacent elements in 00160 *> ALPHAR, ALPHAI, and BETA. Thus, if the j-th and (j+1)-st 00161 *> eigenvalues form a pair, ALPHAI(j+1)/BETA(j+1) is assumed to 00162 *> be equal to -ALPHAI(j)/BETA(j). Also, nonzero values in 00163 *> ALPHAI are assumed to always come in adjacent pairs. 00164 *> \endverbatim 00165 *> 00166 *> \param[in] BETA 00167 *> \verbatim 00168 *> BETA is DOUBLE PRECISION array, dimension (N) 00169 *> The values b(j) as described above, which, along with a(j), 00170 *> define the generalized eigenvalues. 00171 *> \endverbatim 00172 *> 00173 *> \param[out] WORK 00174 *> \verbatim 00175 *> WORK is DOUBLE PRECISION array, dimension (N**2+N) 00176 *> \endverbatim 00177 *> 00178 *> \param[out] RESULT 00179 *> \verbatim 00180 *> RESULT is DOUBLE PRECISION array, dimension (2) 00181 *> The values computed by the test described above. If A E or 00182 *> B E is likely to overflow, then RESULT(1:2) is set to 00183 *> 10 / ulp. 00184 *> \endverbatim 00185 * 00186 * Authors: 00187 * ======== 00188 * 00189 *> \author Univ. of Tennessee 00190 *> \author Univ. of California Berkeley 00191 *> \author Univ. of Colorado Denver 00192 *> \author NAG Ltd. 00193 * 00194 *> \date November 2011 00195 * 00196 *> \ingroup double_eig 00197 * 00198 * ===================================================================== 00199 SUBROUTINE DGET52( LEFT, N, A, LDA, B, LDB, E, LDE, ALPHAR, 00200 $ ALPHAI, BETA, WORK, RESULT ) 00201 * 00202 * -- LAPACK test routine (version 3.4.0) -- 00203 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00204 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00205 * November 2011 00206 * 00207 * .. Scalar Arguments .. 00208 LOGICAL LEFT 00209 INTEGER LDA, LDB, LDE, N 00210 * .. 00211 * .. Array Arguments .. 00212 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ), 00213 $ B( LDB, * ), BETA( * ), E( LDE, * ), 00214 $ RESULT( 2 ), WORK( * ) 00215 * .. 00216 * 00217 * ===================================================================== 00218 * 00219 * .. Parameters .. 00220 DOUBLE PRECISION ZERO, ONE, TEN 00221 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TEN = 10.0D0 ) 00222 * .. 00223 * .. Local Scalars .. 00224 LOGICAL ILCPLX 00225 CHARACTER NORMAB, TRANS 00226 INTEGER J, JVEC 00227 DOUBLE PRECISION ABMAX, ACOEF, ALFMAX, ANORM, BCOEFI, BCOEFR, 00228 $ BETMAX, BNORM, ENORM, ENRMER, ERRNRM, SAFMAX, 00229 $ SAFMIN, SALFI, SALFR, SBETA, SCALE, TEMP1, ULP 00230 * .. 00231 * .. External Functions .. 00232 DOUBLE PRECISION DLAMCH, DLANGE 00233 EXTERNAL DLAMCH, DLANGE 00234 * .. 00235 * .. External Subroutines .. 00236 EXTERNAL DGEMV 00237 * .. 00238 * .. Intrinsic Functions .. 00239 INTRINSIC ABS, DBLE, MAX 00240 * .. 00241 * .. Executable Statements .. 00242 * 00243 RESULT( 1 ) = ZERO 00244 RESULT( 2 ) = ZERO 00245 IF( N.LE.0 ) 00246 $ RETURN 00247 * 00248 SAFMIN = DLAMCH( 'Safe minimum' ) 00249 SAFMAX = ONE / SAFMIN 00250 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) 00251 * 00252 IF( LEFT ) THEN 00253 TRANS = 'T' 00254 NORMAB = 'I' 00255 ELSE 00256 TRANS = 'N' 00257 NORMAB = 'O' 00258 END IF 00259 * 00260 * Norm of A, B, and E: 00261 * 00262 ANORM = MAX( DLANGE( NORMAB, N, N, A, LDA, WORK ), SAFMIN ) 00263 BNORM = MAX( DLANGE( NORMAB, N, N, B, LDB, WORK ), SAFMIN ) 00264 ENORM = MAX( DLANGE( 'O', N, N, E, LDE, WORK ), ULP ) 00265 ALFMAX = SAFMAX / MAX( ONE, BNORM ) 00266 BETMAX = SAFMAX / MAX( ONE, ANORM ) 00267 * 00268 * Compute error matrix. 00269 * Column i = ( b(i) A - a(i) B ) E(i) / max( |a(i) B| |b(i) A| ) 00270 * 00271 ILCPLX = .FALSE. 00272 DO 10 JVEC = 1, N 00273 IF( ILCPLX ) THEN 00274 * 00275 * 2nd Eigenvalue/-vector of pair -- do nothing 00276 * 00277 ILCPLX = .FALSE. 00278 ELSE 00279 SALFR = ALPHAR( JVEC ) 00280 SALFI = ALPHAI( JVEC ) 00281 SBETA = BETA( JVEC ) 00282 IF( SALFI.EQ.ZERO ) THEN 00283 * 00284 * Real eigenvalue and -vector 00285 * 00286 ABMAX = MAX( ABS( SALFR ), ABS( SBETA ) ) 00287 IF( ABS( SALFR ).GT.ALFMAX .OR. ABS( SBETA ).GT. 00288 $ BETMAX .OR. ABMAX.LT.ONE ) THEN 00289 SCALE = ONE / MAX( ABMAX, SAFMIN ) 00290 SALFR = SCALE*SALFR 00291 SBETA = SCALE*SBETA 00292 END IF 00293 SCALE = ONE / MAX( ABS( SALFR )*BNORM, 00294 $ ABS( SBETA )*ANORM, SAFMIN ) 00295 ACOEF = SCALE*SBETA 00296 BCOEFR = SCALE*SALFR 00297 CALL DGEMV( TRANS, N, N, ACOEF, A, LDA, E( 1, JVEC ), 1, 00298 $ ZERO, WORK( N*( JVEC-1 )+1 ), 1 ) 00299 CALL DGEMV( TRANS, N, N, -BCOEFR, B, LDA, E( 1, JVEC ), 00300 $ 1, ONE, WORK( N*( JVEC-1 )+1 ), 1 ) 00301 ELSE 00302 * 00303 * Complex conjugate pair 00304 * 00305 ILCPLX = .TRUE. 00306 IF( JVEC.EQ.N ) THEN 00307 RESULT( 1 ) = TEN / ULP 00308 RETURN 00309 END IF 00310 ABMAX = MAX( ABS( SALFR )+ABS( SALFI ), ABS( SBETA ) ) 00311 IF( ABS( SALFR )+ABS( SALFI ).GT.ALFMAX .OR. 00312 $ ABS( SBETA ).GT.BETMAX .OR. ABMAX.LT.ONE ) THEN 00313 SCALE = ONE / MAX( ABMAX, SAFMIN ) 00314 SALFR = SCALE*SALFR 00315 SALFI = SCALE*SALFI 00316 SBETA = SCALE*SBETA 00317 END IF 00318 SCALE = ONE / MAX( ( ABS( SALFR )+ABS( SALFI ) )*BNORM, 00319 $ ABS( SBETA )*ANORM, SAFMIN ) 00320 ACOEF = SCALE*SBETA 00321 BCOEFR = SCALE*SALFR 00322 BCOEFI = SCALE*SALFI 00323 IF( LEFT ) THEN 00324 BCOEFI = -BCOEFI 00325 END IF 00326 * 00327 CALL DGEMV( TRANS, N, N, ACOEF, A, LDA, E( 1, JVEC ), 1, 00328 $ ZERO, WORK( N*( JVEC-1 )+1 ), 1 ) 00329 CALL DGEMV( TRANS, N, N, -BCOEFR, B, LDA, E( 1, JVEC ), 00330 $ 1, ONE, WORK( N*( JVEC-1 )+1 ), 1 ) 00331 CALL DGEMV( TRANS, N, N, BCOEFI, B, LDA, E( 1, JVEC+1 ), 00332 $ 1, ONE, WORK( N*( JVEC-1 )+1 ), 1 ) 00333 * 00334 CALL DGEMV( TRANS, N, N, ACOEF, A, LDA, E( 1, JVEC+1 ), 00335 $ 1, ZERO, WORK( N*JVEC+1 ), 1 ) 00336 CALL DGEMV( TRANS, N, N, -BCOEFI, B, LDA, E( 1, JVEC ), 00337 $ 1, ONE, WORK( N*JVEC+1 ), 1 ) 00338 CALL DGEMV( TRANS, N, N, -BCOEFR, B, LDA, E( 1, JVEC+1 ), 00339 $ 1, ONE, WORK( N*JVEC+1 ), 1 ) 00340 END IF 00341 END IF 00342 10 CONTINUE 00343 * 00344 ERRNRM = DLANGE( 'One', N, N, WORK, N, WORK( N**2+1 ) ) / ENORM 00345 * 00346 * Compute RESULT(1) 00347 * 00348 RESULT( 1 ) = ERRNRM / ULP 00349 * 00350 * Normalization of E: 00351 * 00352 ENRMER = ZERO 00353 ILCPLX = .FALSE. 00354 DO 40 JVEC = 1, N 00355 IF( ILCPLX ) THEN 00356 ILCPLX = .FALSE. 00357 ELSE 00358 TEMP1 = ZERO 00359 IF( ALPHAI( JVEC ).EQ.ZERO ) THEN 00360 DO 20 J = 1, N 00361 TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) ) ) 00362 20 CONTINUE 00363 ENRMER = MAX( ENRMER, TEMP1-ONE ) 00364 ELSE 00365 ILCPLX = .TRUE. 00366 DO 30 J = 1, N 00367 TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) )+ 00368 $ ABS( E( J, JVEC+1 ) ) ) 00369 30 CONTINUE 00370 ENRMER = MAX( ENRMER, TEMP1-ONE ) 00371 END IF 00372 END IF 00373 40 CONTINUE 00374 * 00375 * Compute RESULT(2) : the normalization error in E. 00376 * 00377 RESULT( 2 ) = ENRMER / ( DBLE( N )*ULP ) 00378 * 00379 RETURN 00380 * 00381 * End of DGET52 00382 * 00383 END