![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CGET22 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 CGET22( TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, W, 00012 * WORK, RWORK, RESULT ) 00013 * 00014 * .. Scalar Arguments .. 00015 * CHARACTER TRANSA, TRANSE, TRANSW 00016 * INTEGER LDA, LDE, N 00017 * .. 00018 * .. Array Arguments .. 00019 * REAL RESULT( 2 ), RWORK( * ) 00020 * COMPLEX A( LDA, * ), E( LDE, * ), W( * ), WORK( * ) 00021 * .. 00022 * 00023 * 00024 *> \par Purpose: 00025 * ============= 00026 *> 00027 *> \verbatim 00028 *> 00029 *> CGET22 does an eigenvector check. 00030 *> 00031 *> The basic test is: 00032 *> 00033 *> RESULT(1) = | A E - E W | / ( |A| |E| ulp ) 00034 *> 00035 *> using the 1-norm. It also tests the normalization of E: 00036 *> 00037 *> RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp ) 00038 *> j 00039 *> 00040 *> where E(j) is the j-th eigenvector, and m-norm is the max-norm of a 00041 *> vector. The max-norm of a complex n-vector x in this case is the 00042 *> maximum of |re(x(i)| + |im(x(i)| over i = 1, ..., n. 00043 *> \endverbatim 00044 * 00045 * Arguments: 00046 * ========== 00047 * 00048 *> \param[in] TRANSA 00049 *> \verbatim 00050 *> TRANSA is CHARACTER*1 00051 *> Specifies whether or not A is transposed. 00052 *> = 'N': No transpose 00053 *> = 'T': Transpose 00054 *> = 'C': Conjugate transpose 00055 *> \endverbatim 00056 *> 00057 *> \param[in] TRANSE 00058 *> \verbatim 00059 *> TRANSE is CHARACTER*1 00060 *> Specifies whether or not E is transposed. 00061 *> = 'N': No transpose, eigenvectors are in columns of E 00062 *> = 'T': Transpose, eigenvectors are in rows of E 00063 *> = 'C': Conjugate transpose, eigenvectors are in rows of E 00064 *> \endverbatim 00065 *> 00066 *> \param[in] TRANSW 00067 *> \verbatim 00068 *> TRANSW is CHARACTER*1 00069 *> Specifies whether or not W is transposed. 00070 *> = 'N': No transpose 00071 *> = 'T': Transpose, same as TRANSW = 'N' 00072 *> = 'C': Conjugate transpose, use -WI(j) instead of WI(j) 00073 *> \endverbatim 00074 *> 00075 *> \param[in] N 00076 *> \verbatim 00077 *> N is INTEGER 00078 *> The order of the matrix A. N >= 0. 00079 *> \endverbatim 00080 *> 00081 *> \param[in] A 00082 *> \verbatim 00083 *> A is COMPLEX array, dimension (LDA,N) 00084 *> The matrix whose eigenvectors are in E. 00085 *> \endverbatim 00086 *> 00087 *> \param[in] LDA 00088 *> \verbatim 00089 *> LDA is INTEGER 00090 *> The leading dimension of the array A. LDA >= max(1,N). 00091 *> \endverbatim 00092 *> 00093 *> \param[in] E 00094 *> \verbatim 00095 *> E is COMPLEX array, dimension (LDE,N) 00096 *> The matrix of eigenvectors. If TRANSE = 'N', the eigenvectors 00097 *> are stored in the columns of E, if TRANSE = 'T' or 'C', the 00098 *> eigenvectors are stored in the rows of E. 00099 *> \endverbatim 00100 *> 00101 *> \param[in] LDE 00102 *> \verbatim 00103 *> LDE is INTEGER 00104 *> The leading dimension of the array E. LDE >= max(1,N). 00105 *> \endverbatim 00106 *> 00107 *> \param[in] W 00108 *> \verbatim 00109 *> W is COMPLEX array, dimension (N) 00110 *> The eigenvalues of A. 00111 *> \endverbatim 00112 *> 00113 *> \param[out] WORK 00114 *> \verbatim 00115 *> WORK is COMPLEX array, dimension (N*N) 00116 *> \endverbatim 00117 *> 00118 *> \param[out] RWORK 00119 *> \verbatim 00120 *> RWORK is REAL array, dimension (N) 00121 *> \endverbatim 00122 *> 00123 *> \param[out] RESULT 00124 *> \verbatim 00125 *> RESULT is REAL array, dimension (2) 00126 *> RESULT(1) = | A E - E W | / ( |A| |E| ulp ) 00127 *> RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp ) 00128 *> \endverbatim 00129 * 00130 * Authors: 00131 * ======== 00132 * 00133 *> \author Univ. of Tennessee 00134 *> \author Univ. of California Berkeley 00135 *> \author Univ. of Colorado Denver 00136 *> \author NAG Ltd. 00137 * 00138 *> \date November 2011 00139 * 00140 *> \ingroup complex_eig 00141 * 00142 * ===================================================================== 00143 SUBROUTINE CGET22( TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, W, 00144 $ WORK, RWORK, RESULT ) 00145 * 00146 * -- LAPACK test routine (version 3.4.0) -- 00147 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00148 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00149 * November 2011 00150 * 00151 * .. Scalar Arguments .. 00152 CHARACTER TRANSA, TRANSE, TRANSW 00153 INTEGER LDA, LDE, N 00154 * .. 00155 * .. Array Arguments .. 00156 REAL RESULT( 2 ), RWORK( * ) 00157 COMPLEX A( LDA, * ), E( LDE, * ), W( * ), WORK( * ) 00158 * .. 00159 * 00160 * ===================================================================== 00161 * 00162 * .. Parameters .. 00163 REAL ZERO, ONE 00164 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00165 COMPLEX CZERO, CONE 00166 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ), 00167 $ CONE = ( 1.0E+0, 0.0E+0 ) ) 00168 * .. 00169 * .. Local Scalars .. 00170 CHARACTER NORMA, NORME 00171 INTEGER ITRNSE, ITRNSW, J, JCOL, JOFF, JROW, JVEC 00172 REAL ANORM, ENORM, ENRMAX, ENRMIN, ERRNRM, TEMP1, 00173 $ ULP, UNFL 00174 COMPLEX WTEMP 00175 * .. 00176 * .. External Functions .. 00177 LOGICAL LSAME 00178 REAL CLANGE, SLAMCH 00179 EXTERNAL LSAME, CLANGE, SLAMCH 00180 * .. 00181 * .. External Subroutines .. 00182 EXTERNAL CGEMM, CLASET 00183 * .. 00184 * .. Intrinsic Functions .. 00185 INTRINSIC ABS, AIMAG, CONJG, MAX, MIN, REAL 00186 * .. 00187 * .. Executable Statements .. 00188 * 00189 * Initialize RESULT (in case N=0) 00190 * 00191 RESULT( 1 ) = ZERO 00192 RESULT( 2 ) = ZERO 00193 IF( N.LE.0 ) 00194 $ RETURN 00195 * 00196 UNFL = SLAMCH( 'Safe minimum' ) 00197 ULP = SLAMCH( 'Precision' ) 00198 * 00199 ITRNSE = 0 00200 ITRNSW = 0 00201 NORMA = 'O' 00202 NORME = 'O' 00203 * 00204 IF( LSAME( TRANSA, 'T' ) .OR. LSAME( TRANSA, 'C' ) ) THEN 00205 NORMA = 'I' 00206 END IF 00207 * 00208 IF( LSAME( TRANSE, 'T' ) ) THEN 00209 ITRNSE = 1 00210 NORME = 'I' 00211 ELSE IF( LSAME( TRANSE, 'C' ) ) THEN 00212 ITRNSE = 2 00213 NORME = 'I' 00214 END IF 00215 * 00216 IF( LSAME( TRANSW, 'C' ) ) THEN 00217 ITRNSW = 1 00218 END IF 00219 * 00220 * Normalization of E: 00221 * 00222 ENRMIN = ONE / ULP 00223 ENRMAX = ZERO 00224 IF( ITRNSE.EQ.0 ) THEN 00225 DO 20 JVEC = 1, N 00226 TEMP1 = ZERO 00227 DO 10 J = 1, N 00228 TEMP1 = MAX( TEMP1, ABS( REAL( E( J, JVEC ) ) )+ 00229 $ ABS( AIMAG( E( J, JVEC ) ) ) ) 00230 10 CONTINUE 00231 ENRMIN = MIN( ENRMIN, TEMP1 ) 00232 ENRMAX = MAX( ENRMAX, TEMP1 ) 00233 20 CONTINUE 00234 ELSE 00235 DO 30 JVEC = 1, N 00236 RWORK( JVEC ) = ZERO 00237 30 CONTINUE 00238 * 00239 DO 50 J = 1, N 00240 DO 40 JVEC = 1, N 00241 RWORK( JVEC ) = MAX( RWORK( JVEC ), 00242 $ ABS( REAL( E( JVEC, J ) ) )+ 00243 $ ABS( AIMAG( E( JVEC, J ) ) ) ) 00244 40 CONTINUE 00245 50 CONTINUE 00246 * 00247 DO 60 JVEC = 1, N 00248 ENRMIN = MIN( ENRMIN, RWORK( JVEC ) ) 00249 ENRMAX = MAX( ENRMAX, RWORK( JVEC ) ) 00250 60 CONTINUE 00251 END IF 00252 * 00253 * Norm of A: 00254 * 00255 ANORM = MAX( CLANGE( NORMA, N, N, A, LDA, RWORK ), UNFL ) 00256 * 00257 * Norm of E: 00258 * 00259 ENORM = MAX( CLANGE( NORME, N, N, E, LDE, RWORK ), ULP ) 00260 * 00261 * Norm of error: 00262 * 00263 * Error = AE - EW 00264 * 00265 CALL CLASET( 'Full', N, N, CZERO, CZERO, WORK, N ) 00266 * 00267 JOFF = 0 00268 DO 100 JCOL = 1, N 00269 IF( ITRNSW.EQ.0 ) THEN 00270 WTEMP = W( JCOL ) 00271 ELSE 00272 WTEMP = CONJG( W( JCOL ) ) 00273 END IF 00274 * 00275 IF( ITRNSE.EQ.0 ) THEN 00276 DO 70 JROW = 1, N 00277 WORK( JOFF+JROW ) = E( JROW, JCOL )*WTEMP 00278 70 CONTINUE 00279 ELSE IF( ITRNSE.EQ.1 ) THEN 00280 DO 80 JROW = 1, N 00281 WORK( JOFF+JROW ) = E( JCOL, JROW )*WTEMP 00282 80 CONTINUE 00283 ELSE 00284 DO 90 JROW = 1, N 00285 WORK( JOFF+JROW ) = CONJG( E( JCOL, JROW ) )*WTEMP 00286 90 CONTINUE 00287 END IF 00288 JOFF = JOFF + N 00289 100 CONTINUE 00290 * 00291 CALL CGEMM( TRANSA, TRANSE, N, N, N, CONE, A, LDA, E, LDE, -CONE, 00292 $ WORK, N ) 00293 * 00294 ERRNRM = CLANGE( 'One', N, N, WORK, N, RWORK ) / ENORM 00295 * 00296 * Compute RESULT(1) (avoiding under/overflow) 00297 * 00298 IF( ANORM.GT.ERRNRM ) THEN 00299 RESULT( 1 ) = ( ERRNRM / ANORM ) / ULP 00300 ELSE 00301 IF( ANORM.LT.ONE ) THEN 00302 RESULT( 1 ) = ( MIN( ERRNRM, ANORM ) / ANORM ) / ULP 00303 ELSE 00304 RESULT( 1 ) = MIN( ERRNRM / ANORM, ONE ) / ULP 00305 END IF 00306 END IF 00307 * 00308 * Compute RESULT(2) : the normalization error in E. 00309 * 00310 RESULT( 2 ) = MAX( ABS( ENRMAX-ONE ), ABS( ENRMIN-ONE ) ) / 00311 $ ( REAL( N )*ULP ) 00312 * 00313 RETURN 00314 * 00315 * End of CGET22 00316 * 00317 END