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