 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DGET22( TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, WR, 00002 $ WI, 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 CHARACTER TRANSA, TRANSE, TRANSW 00010 INTEGER LDA, LDE, N 00011 * .. 00012 * .. Array Arguments .. 00013 DOUBLE PRECISION A( LDA, * ), E( LDE, * ), RESULT( 2 ), WI( * ), 00014 $ WORK( * ), WR( * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * DGET22 does an eigenvector check. 00021 * 00022 * The basic test is: 00023 * 00024 * RESULT(1) = | A E - E W | / ( |A| |E| ulp ) 00025 * 00026 * using the 1-norm. It also tests the normalization of E: 00027 * 00028 * RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp ) 00029 * j 00030 * 00031 * where E(j) is the j-th eigenvector, and m-norm is the max-norm of a 00032 * vector. If an eigenvector is complex, as determined from WI(j) 00033 * nonzero, then the max-norm of the vector ( er + i*ei ) is the maximum 00034 * of 00035 * |er(1)| + |ei(1)|, ... , |er(n)| + |ei(n)| 00036 * 00037 * W is a block diagonal matrix, with a 1 by 1 block for each real 00038 * eigenvalue and a 2 by 2 block for each complex conjugate pair. 00039 * If eigenvalues j and j+1 are a complex conjugate pair, so that 00040 * WR(j) = WR(j+1) = wr and WI(j) = - WI(j+1) = wi, then the 2 by 2 00041 * block corresponding to the pair will be: 00042 * 00043 * ( wr wi ) 00044 * ( -wi wr ) 00045 * 00046 * Such a block multiplying an n by 2 matrix ( ur ui ) on the right 00047 * will be the same as multiplying ur + i*ui by wr + i*wi. 00048 * 00049 * To handle various schemes for storage of left eigenvectors, there are 00050 * options to use A-transpose instead of A, E-transpose instead of E, 00051 * and/or W-transpose instead of W. 00052 * 00053 * Arguments 00054 * ========== 00055 * 00056 * TRANSA (input) CHARACTER*1 00057 * Specifies whether or not A is transposed. 00058 * = 'N': No transpose 00059 * = 'T': Transpose 00060 * = 'C': Conjugate transpose (= Transpose) 00061 * 00062 * TRANSE (input) CHARACTER*1 00063 * Specifies whether or not E is transposed. 00064 * = 'N': No transpose, eigenvectors are in columns of E 00065 * = 'T': Transpose, eigenvectors are in rows of E 00066 * = 'C': Conjugate transpose (= Transpose) 00067 * 00068 * TRANSW (input) CHARACTER*1 00069 * Specifies whether or not W is transposed. 00070 * = 'N': No transpose 00071 * = 'T': Transpose, use -WI(j) instead of WI(j) 00072 * = 'C': Conjugate transpose, use -WI(j) instead of WI(j) 00073 * 00074 * N (input) INTEGER 00075 * The order of the matrix A. N >= 0. 00076 * 00077 * A (input) DOUBLE PRECISION array, dimension (LDA,N) 00078 * The matrix whose eigenvectors are in E. 00079 * 00080 * LDA (input) INTEGER 00081 * The leading dimension of the array A. LDA >= max(1,N). 00082 * 00083 * E (input) DOUBLE PRECISION array, dimension (LDE,N) 00084 * The matrix of eigenvectors. If TRANSE = 'N', the eigenvectors 00085 * are stored in the columns of E, if TRANSE = 'T' or 'C', the 00086 * eigenvectors are stored in the rows of E. 00087 * 00088 * LDE (input) INTEGER 00089 * The leading dimension of the array E. LDE >= max(1,N). 00090 * 00091 * WR (input) DOUBLE PRECISION array, dimension (N) 00092 * WI (input) DOUBLE PRECISION array, dimension (N) 00093 * The real and imaginary parts of the eigenvalues of A. 00094 * Purely real eigenvalues are indicated by WI(j) = 0. 00095 * Complex conjugate pairs are indicated by WR(j)=WR(j+1) and 00096 * WI(j) = - WI(j+1) non-zero; the real part is assumed to be 00097 * stored in the j-th row/column and the imaginary part in 00098 * the (j+1)-th row/column. 00099 * 00100 * WORK (workspace) DOUBLE PRECISION array, dimension (N*(N+1)) 00101 * 00102 * RESULT (output) DOUBLE PRECISION array, dimension (2) 00103 * RESULT(1) = | A E - E W | / ( |A| |E| ulp ) 00104 * RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp ) 00105 * 00106 * ===================================================================== 00107 * 00108 * .. Parameters .. 00109 DOUBLE PRECISION ZERO, ONE 00110 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00111 * .. 00112 * .. Local Scalars .. 00113 CHARACTER NORMA, NORME 00114 INTEGER IECOL, IEROW, INCE, IPAIR, ITRNSE, J, JCOL, 00115 $ JVEC 00116 DOUBLE PRECISION ANORM, ENORM, ENRMAX, ENRMIN, ERRNRM, TEMP1, 00117 $ ULP, UNFL 00118 * .. 00119 * .. Local Arrays .. 00120 DOUBLE PRECISION WMAT( 2, 2 ) 00121 * .. 00122 * .. External Functions .. 00123 LOGICAL LSAME 00124 DOUBLE PRECISION DLAMCH, DLANGE 00125 EXTERNAL LSAME, DLAMCH, DLANGE 00126 * .. 00127 * .. External Subroutines .. 00128 EXTERNAL DAXPY, DGEMM, DLASET 00129 * .. 00130 * .. Intrinsic Functions .. 00131 INTRINSIC ABS, DBLE, MAX, MIN 00132 * .. 00133 * .. Executable Statements .. 00134 * 00135 * Initialize RESULT (in case N=0) 00136 * 00137 RESULT( 1 ) = ZERO 00138 RESULT( 2 ) = ZERO 00139 IF( N.LE.0 ) 00140 $ RETURN 00141 * 00142 UNFL = DLAMCH( 'Safe minimum' ) 00143 ULP = DLAMCH( 'Precision' ) 00144 * 00145 ITRNSE = 0 00146 INCE = 1 00147 NORMA = 'O' 00148 NORME = 'O' 00149 * 00150 IF( LSAME( TRANSA, 'T' ) .OR. LSAME( TRANSA, 'C' ) ) THEN 00151 NORMA = 'I' 00152 END IF 00153 IF( LSAME( TRANSE, 'T' ) .OR. LSAME( TRANSE, 'C' ) ) THEN 00154 NORME = 'I' 00155 ITRNSE = 1 00156 INCE = LDE 00157 END IF 00158 * 00159 * Check normalization of E 00160 * 00161 ENRMIN = ONE / ULP 00162 ENRMAX = ZERO 00163 IF( ITRNSE.EQ.0 ) THEN 00164 * 00165 * Eigenvectors are column vectors. 00166 * 00167 IPAIR = 0 00168 DO 30 JVEC = 1, N 00169 TEMP1 = ZERO 00170 IF( IPAIR.EQ.0 .AND. JVEC.LT.N .AND. WI( JVEC ).NE.ZERO ) 00171 $ IPAIR = 1 00172 IF( IPAIR.EQ.1 ) THEN 00173 * 00174 * Complex eigenvector 00175 * 00176 DO 10 J = 1, N 00177 TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) )+ 00178 $ ABS( E( J, JVEC+1 ) ) ) 00179 10 CONTINUE 00180 ENRMIN = MIN( ENRMIN, TEMP1 ) 00181 ENRMAX = MAX( ENRMAX, TEMP1 ) 00182 IPAIR = 2 00183 ELSE IF( IPAIR.EQ.2 ) THEN 00184 IPAIR = 0 00185 ELSE 00186 * 00187 * Real eigenvector 00188 * 00189 DO 20 J = 1, N 00190 TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) ) ) 00191 20 CONTINUE 00192 ENRMIN = MIN( ENRMIN, TEMP1 ) 00193 ENRMAX = MAX( ENRMAX, TEMP1 ) 00194 IPAIR = 0 00195 END IF 00196 30 CONTINUE 00197 * 00198 ELSE 00199 * 00200 * Eigenvectors are row vectors. 00201 * 00202 DO 40 JVEC = 1, N 00203 WORK( JVEC ) = ZERO 00204 40 CONTINUE 00205 * 00206 DO 60 J = 1, N 00207 IPAIR = 0 00208 DO 50 JVEC = 1, N 00209 IF( IPAIR.EQ.0 .AND. JVEC.LT.N .AND. WI( JVEC ).NE.ZERO ) 00210 $ IPAIR = 1 00211 IF( IPAIR.EQ.1 ) THEN 00212 WORK( JVEC ) = MAX( WORK( JVEC ), 00213 $ ABS( E( J, JVEC ) )+ABS( E( J, 00214 $ JVEC+1 ) ) ) 00215 WORK( JVEC+1 ) = WORK( JVEC ) 00216 ELSE IF( IPAIR.EQ.2 ) THEN 00217 IPAIR = 0 00218 ELSE 00219 WORK( JVEC ) = MAX( WORK( JVEC ), 00220 $ ABS( E( J, JVEC ) ) ) 00221 IPAIR = 0 00222 END IF 00223 50 CONTINUE 00224 60 CONTINUE 00225 * 00226 DO 70 JVEC = 1, N 00227 ENRMIN = MIN( ENRMIN, WORK( JVEC ) ) 00228 ENRMAX = MAX( ENRMAX, WORK( JVEC ) ) 00229 70 CONTINUE 00230 END IF 00231 * 00232 * Norm of A: 00233 * 00234 ANORM = MAX( DLANGE( NORMA, N, N, A, LDA, WORK ), UNFL ) 00235 * 00236 * Norm of E: 00237 * 00238 ENORM = MAX( DLANGE( NORME, N, N, E, LDE, WORK ), ULP ) 00239 * 00240 * Norm of error: 00241 * 00242 * Error = AE - EW 00243 * 00244 CALL DLASET( 'Full', N, N, ZERO, ZERO, WORK, N ) 00245 * 00246 IPAIR = 0 00247 IEROW = 1 00248 IECOL = 1 00249 * 00250 DO 80 JCOL = 1, N 00251 IF( ITRNSE.EQ.1 ) THEN 00252 IEROW = JCOL 00253 ELSE 00254 IECOL = JCOL 00255 END IF 00256 * 00257 IF( IPAIR.EQ.0 .AND. WI( JCOL ).NE.ZERO ) 00258 $ IPAIR = 1 00259 * 00260 IF( IPAIR.EQ.1 ) THEN 00261 WMAT( 1, 1 ) = WR( JCOL ) 00262 WMAT( 2, 1 ) = -WI( JCOL ) 00263 WMAT( 1, 2 ) = WI( JCOL ) 00264 WMAT( 2, 2 ) = WR( JCOL ) 00265 CALL DGEMM( TRANSE, TRANSW, N, 2, 2, ONE, E( IEROW, IECOL ), 00266 $ LDE, WMAT, 2, ZERO, WORK( N*( JCOL-1 )+1 ), N ) 00267 IPAIR = 2 00268 ELSE IF( IPAIR.EQ.2 ) THEN 00269 IPAIR = 0 00270 * 00271 ELSE 00272 * 00273 CALL DAXPY( N, WR( JCOL ), E( IEROW, IECOL ), INCE, 00274 $ WORK( N*( JCOL-1 )+1 ), 1 ) 00275 IPAIR = 0 00276 END IF 00277 * 00278 80 CONTINUE 00279 * 00280 CALL DGEMM( TRANSA, TRANSE, N, N, N, ONE, A, LDA, E, LDE, -ONE, 00281 $ WORK, N ) 00282 * 00283 ERRNRM = DLANGE( 'One', N, N, WORK, N, WORK( N*N+1 ) ) / ENORM 00284 * 00285 * Compute RESULT(1) (avoiding under/overflow) 00286 * 00287 IF( ANORM.GT.ERRNRM ) THEN 00288 RESULT( 1 ) = ( ERRNRM / ANORM ) / ULP 00289 ELSE 00290 IF( ANORM.LT.ONE ) THEN 00291 RESULT( 1 ) = ( MIN( ERRNRM, ANORM ) / ANORM ) / ULP 00292 ELSE 00293 RESULT( 1 ) = MIN( ERRNRM / ANORM, ONE ) / ULP 00294 END IF 00295 END IF 00296 * 00297 * Compute RESULT(2) : the normalization error in E. 00298 * 00299 RESULT( 2 ) = MAX( ABS( ENRMAX-ONE ), ABS( ENRMIN-ONE ) ) / 00300 $ ( DBLE( N )*ULP ) 00301 * 00302 RETURN 00303 * 00304 * End of DGET22 00305 * 00306 END
1.7.3 
