![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SGET22 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 SGET22( TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, WR, 00012 * WI, WORK, RESULT ) 00013 * 00014 * .. Scalar Arguments .. 00015 * CHARACTER TRANSA, TRANSE, TRANSW 00016 * INTEGER LDA, LDE, N 00017 * .. 00018 * .. Array Arguments .. 00019 * REAL A( LDA, * ), E( LDE, * ), RESULT( 2 ), WI( * ), 00020 * $ WORK( * ), WR( * ) 00021 * .. 00022 * 00023 * 00024 *> \par Purpose: 00025 * ============= 00026 *> 00027 *> \verbatim 00028 *> 00029 *> SGET22 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. If an eigenvector is complex, as determined from WI(j) 00042 *> nonzero, then the max-norm of the vector ( er + i*ei ) is the maximum 00043 *> of 00044 *> |er(1)| + |ei(1)|, ... , |er(n)| + |ei(n)| 00045 *> 00046 *> W is a block diagonal matrix, with a 1 by 1 block for each real 00047 *> eigenvalue and a 2 by 2 block for each complex conjugate pair. 00048 *> If eigenvalues j and j+1 are a complex conjugate pair, so that 00049 *> WR(j) = WR(j+1) = wr and WI(j) = - WI(j+1) = wi, then the 2 by 2 00050 *> block corresponding to the pair will be: 00051 *> 00052 *> ( wr wi ) 00053 *> ( -wi wr ) 00054 *> 00055 *> Such a block multiplying an n by 2 matrix ( ur ui ) on the right 00056 *> will be the same as multiplying ur + i*ui by wr + i*wi. 00057 *> 00058 *> To handle various schemes for storage of left eigenvectors, there are 00059 *> options to use A-transpose instead of A, E-transpose instead of E, 00060 *> and/or W-transpose instead of W. 00061 *> \endverbatim 00062 * 00063 * Arguments: 00064 * ========== 00065 * 00066 *> \param[in] TRANSA 00067 *> \verbatim 00068 *> TRANSA is CHARACTER*1 00069 *> Specifies whether or not A is transposed. 00070 *> = 'N': No transpose 00071 *> = 'T': Transpose 00072 *> = 'C': Conjugate transpose (= Transpose) 00073 *> \endverbatim 00074 *> 00075 *> \param[in] TRANSE 00076 *> \verbatim 00077 *> TRANSE is CHARACTER*1 00078 *> Specifies whether or not E is transposed. 00079 *> = 'N': No transpose, eigenvectors are in columns of E 00080 *> = 'T': Transpose, eigenvectors are in rows of E 00081 *> = 'C': Conjugate transpose (= Transpose) 00082 *> \endverbatim 00083 *> 00084 *> \param[in] TRANSW 00085 *> \verbatim 00086 *> TRANSW is CHARACTER*1 00087 *> Specifies whether or not W is transposed. 00088 *> = 'N': No transpose 00089 *> = 'T': Transpose, use -WI(j) instead of WI(j) 00090 *> = 'C': Conjugate transpose, use -WI(j) instead of WI(j) 00091 *> \endverbatim 00092 *> 00093 *> \param[in] N 00094 *> \verbatim 00095 *> N is INTEGER 00096 *> The order of the matrix A. N >= 0. 00097 *> \endverbatim 00098 *> 00099 *> \param[in] A 00100 *> \verbatim 00101 *> A is REAL array, dimension (LDA,N) 00102 *> The matrix whose eigenvectors are in E. 00103 *> \endverbatim 00104 *> 00105 *> \param[in] LDA 00106 *> \verbatim 00107 *> LDA is INTEGER 00108 *> The leading dimension of the array A. LDA >= max(1,N). 00109 *> \endverbatim 00110 *> 00111 *> \param[in] E 00112 *> \verbatim 00113 *> E is REAL array, dimension (LDE,N) 00114 *> The matrix of eigenvectors. If TRANSE = 'N', the eigenvectors 00115 *> are stored in the columns of E, if TRANSE = 'T' or 'C', the 00116 *> eigenvectors are stored in the rows of E. 00117 *> \endverbatim 00118 *> 00119 *> \param[in] LDE 00120 *> \verbatim 00121 *> LDE is INTEGER 00122 *> The leading dimension of the array E. LDE >= max(1,N). 00123 *> \endverbatim 00124 *> 00125 *> \param[in] WR 00126 *> \verbatim 00127 *> WR is REAL array, dimension (N) 00128 *> \endverbatim 00129 *> 00130 *> \param[in] WI 00131 *> \verbatim 00132 *> WI is REAL array, dimension (N) 00133 *> 00134 *> The real and imaginary parts of the eigenvalues of A. 00135 *> Purely real eigenvalues are indicated by WI(j) = 0. 00136 *> Complex conjugate pairs are indicated by WR(j)=WR(j+1) and 00137 *> WI(j) = - WI(j+1) non-zero; the real part is assumed to be 00138 *> stored in the j-th row/column and the imaginary part in 00139 *> the (j+1)-th row/column. 00140 *> \endverbatim 00141 *> 00142 *> \param[out] WORK 00143 *> \verbatim 00144 *> WORK is REAL array, dimension (N*(N+1)) 00145 *> \endverbatim 00146 *> 00147 *> \param[out] RESULT 00148 *> \verbatim 00149 *> RESULT is REAL array, dimension (2) 00150 *> RESULT(1) = | A E - E W | / ( |A| |E| ulp ) 00151 *> RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp ) 00152 *> \endverbatim 00153 * 00154 * Authors: 00155 * ======== 00156 * 00157 *> \author Univ. of Tennessee 00158 *> \author Univ. of California Berkeley 00159 *> \author Univ. of Colorado Denver 00160 *> \author NAG Ltd. 00161 * 00162 *> \date November 2011 00163 * 00164 *> \ingroup single_eig 00165 * 00166 * ===================================================================== 00167 SUBROUTINE SGET22( TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, WR, 00168 $ WI, WORK, RESULT ) 00169 * 00170 * -- LAPACK test routine (version 3.4.0) -- 00171 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00172 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00173 * November 2011 00174 * 00175 * .. Scalar Arguments .. 00176 CHARACTER TRANSA, TRANSE, TRANSW 00177 INTEGER LDA, LDE, N 00178 * .. 00179 * .. Array Arguments .. 00180 REAL A( LDA, * ), E( LDE, * ), RESULT( 2 ), WI( * ), 00181 $ WORK( * ), WR( * ) 00182 * .. 00183 * 00184 * ===================================================================== 00185 * 00186 * .. Parameters .. 00187 REAL ZERO, ONE 00188 PARAMETER ( ZERO = 0.0, ONE = 1.0 ) 00189 * .. 00190 * .. Local Scalars .. 00191 CHARACTER NORMA, NORME 00192 INTEGER IECOL, IEROW, INCE, IPAIR, ITRNSE, J, JCOL, 00193 $ JVEC 00194 REAL ANORM, ENORM, ENRMAX, ENRMIN, ERRNRM, TEMP1, 00195 $ ULP, UNFL 00196 * .. 00197 * .. Local Arrays .. 00198 REAL WMAT( 2, 2 ) 00199 * .. 00200 * .. External Functions .. 00201 LOGICAL LSAME 00202 REAL SLAMCH, SLANGE 00203 EXTERNAL LSAME, SLAMCH, SLANGE 00204 * .. 00205 * .. External Subroutines .. 00206 EXTERNAL SAXPY, SGEMM, SLASET 00207 * .. 00208 * .. Intrinsic Functions .. 00209 INTRINSIC ABS, MAX, MIN, REAL 00210 * .. 00211 * .. Executable Statements .. 00212 * 00213 * Initialize RESULT (in case N=0) 00214 * 00215 RESULT( 1 ) = ZERO 00216 RESULT( 2 ) = ZERO 00217 IF( N.LE.0 ) 00218 $ RETURN 00219 * 00220 UNFL = SLAMCH( 'Safe minimum' ) 00221 ULP = SLAMCH( 'Precision' ) 00222 * 00223 ITRNSE = 0 00224 INCE = 1 00225 NORMA = 'O' 00226 NORME = 'O' 00227 * 00228 IF( LSAME( TRANSA, 'T' ) .OR. LSAME( TRANSA, 'C' ) ) THEN 00229 NORMA = 'I' 00230 END IF 00231 IF( LSAME( TRANSE, 'T' ) .OR. LSAME( TRANSE, 'C' ) ) THEN 00232 NORME = 'I' 00233 ITRNSE = 1 00234 INCE = LDE 00235 END IF 00236 * 00237 * Check normalization of E 00238 * 00239 ENRMIN = ONE / ULP 00240 ENRMAX = ZERO 00241 IF( ITRNSE.EQ.0 ) THEN 00242 * 00243 * Eigenvectors are column vectors. 00244 * 00245 IPAIR = 0 00246 DO 30 JVEC = 1, N 00247 TEMP1 = ZERO 00248 IF( IPAIR.EQ.0 .AND. JVEC.LT.N .AND. WI( JVEC ).NE.ZERO ) 00249 $ IPAIR = 1 00250 IF( IPAIR.EQ.1 ) THEN 00251 * 00252 * Complex eigenvector 00253 * 00254 DO 10 J = 1, N 00255 TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) )+ 00256 $ ABS( E( J, JVEC+1 ) ) ) 00257 10 CONTINUE 00258 ENRMIN = MIN( ENRMIN, TEMP1 ) 00259 ENRMAX = MAX( ENRMAX, TEMP1 ) 00260 IPAIR = 2 00261 ELSE IF( IPAIR.EQ.2 ) THEN 00262 IPAIR = 0 00263 ELSE 00264 * 00265 * Real eigenvector 00266 * 00267 DO 20 J = 1, N 00268 TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) ) ) 00269 20 CONTINUE 00270 ENRMIN = MIN( ENRMIN, TEMP1 ) 00271 ENRMAX = MAX( ENRMAX, TEMP1 ) 00272 IPAIR = 0 00273 END IF 00274 30 CONTINUE 00275 * 00276 ELSE 00277 * 00278 * Eigenvectors are row vectors. 00279 * 00280 DO 40 JVEC = 1, N 00281 WORK( JVEC ) = ZERO 00282 40 CONTINUE 00283 * 00284 DO 60 J = 1, N 00285 IPAIR = 0 00286 DO 50 JVEC = 1, N 00287 IF( IPAIR.EQ.0 .AND. JVEC.LT.N .AND. WI( JVEC ).NE.ZERO ) 00288 $ IPAIR = 1 00289 IF( IPAIR.EQ.1 ) THEN 00290 WORK( JVEC ) = MAX( WORK( JVEC ), 00291 $ ABS( E( J, JVEC ) )+ABS( E( J, 00292 $ JVEC+1 ) ) ) 00293 WORK( JVEC+1 ) = WORK( JVEC ) 00294 ELSE IF( IPAIR.EQ.2 ) THEN 00295 IPAIR = 0 00296 ELSE 00297 WORK( JVEC ) = MAX( WORK( JVEC ), 00298 $ ABS( E( J, JVEC ) ) ) 00299 IPAIR = 0 00300 END IF 00301 50 CONTINUE 00302 60 CONTINUE 00303 * 00304 DO 70 JVEC = 1, N 00305 ENRMIN = MIN( ENRMIN, WORK( JVEC ) ) 00306 ENRMAX = MAX( ENRMAX, WORK( JVEC ) ) 00307 70 CONTINUE 00308 END IF 00309 * 00310 * Norm of A: 00311 * 00312 ANORM = MAX( SLANGE( NORMA, N, N, A, LDA, WORK ), UNFL ) 00313 * 00314 * Norm of E: 00315 * 00316 ENORM = MAX( SLANGE( NORME, N, N, E, LDE, WORK ), ULP ) 00317 * 00318 * Norm of error: 00319 * 00320 * Error = AE - EW 00321 * 00322 CALL SLASET( 'Full', N, N, ZERO, ZERO, WORK, N ) 00323 * 00324 IPAIR = 0 00325 IEROW = 1 00326 IECOL = 1 00327 * 00328 DO 80 JCOL = 1, N 00329 IF( ITRNSE.EQ.1 ) THEN 00330 IEROW = JCOL 00331 ELSE 00332 IECOL = JCOL 00333 END IF 00334 * 00335 IF( IPAIR.EQ.0 .AND. WI( JCOL ).NE.ZERO ) 00336 $ IPAIR = 1 00337 * 00338 IF( IPAIR.EQ.1 ) THEN 00339 WMAT( 1, 1 ) = WR( JCOL ) 00340 WMAT( 2, 1 ) = -WI( JCOL ) 00341 WMAT( 1, 2 ) = WI( JCOL ) 00342 WMAT( 2, 2 ) = WR( JCOL ) 00343 CALL SGEMM( TRANSE, TRANSW, N, 2, 2, ONE, E( IEROW, IECOL ), 00344 $ LDE, WMAT, 2, ZERO, WORK( N*( JCOL-1 )+1 ), N ) 00345 IPAIR = 2 00346 ELSE IF( IPAIR.EQ.2 ) THEN 00347 IPAIR = 0 00348 * 00349 ELSE 00350 * 00351 CALL SAXPY( N, WR( JCOL ), E( IEROW, IECOL ), INCE, 00352 $ WORK( N*( JCOL-1 )+1 ), 1 ) 00353 IPAIR = 0 00354 END IF 00355 * 00356 80 CONTINUE 00357 * 00358 CALL SGEMM( TRANSA, TRANSE, N, N, N, ONE, A, LDA, E, LDE, -ONE, 00359 $ WORK, N ) 00360 * 00361 ERRNRM = SLANGE( 'One', N, N, WORK, N, WORK( N*N+1 ) ) / ENORM 00362 * 00363 * Compute RESULT(1) (avoiding under/overflow) 00364 * 00365 IF( ANORM.GT.ERRNRM ) THEN 00366 RESULT( 1 ) = ( ERRNRM / ANORM ) / ULP 00367 ELSE 00368 IF( ANORM.LT.ONE ) THEN 00369 RESULT( 1 ) = ( MIN( ERRNRM, ANORM ) / ANORM ) / ULP 00370 ELSE 00371 RESULT( 1 ) = MIN( ERRNRM / ANORM, ONE ) / ULP 00372 END IF 00373 END IF 00374 * 00375 * Compute RESULT(2) : the normalization error in E. 00376 * 00377 RESULT( 2 ) = MAX( ABS( ENRMAX-ONE ), ABS( ENRMIN-ONE ) ) / 00378 $ ( REAL( N )*ULP ) 00379 * 00380 RETURN 00381 * 00382 * End of SGET22 00383 * 00384 END