1*> \brief \b ZGET22 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8* Definition: 9* =========== 10* 11* SUBROUTINE ZGET22( TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, W, 12* WORK, RWORK, RESULT ) 13* 14* .. Scalar Arguments .. 15* CHARACTER TRANSA, TRANSE, TRANSW 16* INTEGER LDA, LDE, N 17* .. 18* .. Array Arguments .. 19* DOUBLE PRECISION RESULT( 2 ), RWORK( * ) 20* COMPLEX*16 A( LDA, * ), E( LDE, * ), W( * ), WORK( * ) 21* .. 22* 23* 24*> \par Purpose: 25* ============= 26*> 27*> \verbatim 28*> 29*> ZGET22 does an eigenvector check. 30*> 31*> The basic test is: 32*> 33*> RESULT(1) = | A E - E W | / ( |A| |E| ulp ) 34*> 35*> using the 1-norm. It also tests the normalization of E: 36*> 37*> RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp ) 38*> j 39*> 40*> where E(j) is the j-th eigenvector, and m-norm is the max-norm of a 41*> vector. The max-norm of a complex n-vector x in this case is the 42*> maximum of |re(x(i)| + |im(x(i)| over i = 1, ..., n. 43*> \endverbatim 44* 45* Arguments: 46* ========== 47* 48*> \param[in] TRANSA 49*> \verbatim 50*> TRANSA is CHARACTER*1 51*> Specifies whether or not A is transposed. 52*> = 'N': No transpose 53*> = 'T': Transpose 54*> = 'C': Conjugate transpose 55*> \endverbatim 56*> 57*> \param[in] TRANSE 58*> \verbatim 59*> TRANSE is CHARACTER*1 60*> Specifies whether or not E is transposed. 61*> = 'N': No transpose, eigenvectors are in columns of E 62*> = 'T': Transpose, eigenvectors are in rows of E 63*> = 'C': Conjugate transpose, eigenvectors are in rows of E 64*> \endverbatim 65*> 66*> \param[in] TRANSW 67*> \verbatim 68*> TRANSW is CHARACTER*1 69*> Specifies whether or not W is transposed. 70*> = 'N': No transpose 71*> = 'T': Transpose, same as TRANSW = 'N' 72*> = 'C': Conjugate transpose, use -WI(j) instead of WI(j) 73*> \endverbatim 74*> 75*> \param[in] N 76*> \verbatim 77*> N is INTEGER 78*> The order of the matrix A. N >= 0. 79*> \endverbatim 80*> 81*> \param[in] A 82*> \verbatim 83*> A is COMPLEX*16 array, dimension (LDA,N) 84*> The matrix whose eigenvectors are in E. 85*> \endverbatim 86*> 87*> \param[in] LDA 88*> \verbatim 89*> LDA is INTEGER 90*> The leading dimension of the array A. LDA >= max(1,N). 91*> \endverbatim 92*> 93*> \param[in] E 94*> \verbatim 95*> E is COMPLEX*16 array, dimension (LDE,N) 96*> The matrix of eigenvectors. If TRANSE = 'N', the eigenvectors 97*> are stored in the columns of E, if TRANSE = 'T' or 'C', the 98*> eigenvectors are stored in the rows of E. 99*> \endverbatim 100*> 101*> \param[in] LDE 102*> \verbatim 103*> LDE is INTEGER 104*> The leading dimension of the array E. LDE >= max(1,N). 105*> \endverbatim 106*> 107*> \param[in] W 108*> \verbatim 109*> W is COMPLEX*16 array, dimension (N) 110*> The eigenvalues of A. 111*> \endverbatim 112*> 113*> \param[out] WORK 114*> \verbatim 115*> WORK is COMPLEX*16 array, dimension (N*N) 116*> \endverbatim 117*> 118*> \param[out] RWORK 119*> \verbatim 120*> RWORK is DOUBLE PRECISION array, dimension (N) 121*> \endverbatim 122*> 123*> \param[out] RESULT 124*> \verbatim 125*> RESULT is DOUBLE PRECISION array, dimension (2) 126*> RESULT(1) = | A E - E W | / ( |A| |E| ulp ) 127*> RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp ) 128*> j 129*> \endverbatim 130* 131* Authors: 132* ======== 133* 134*> \author Univ. of Tennessee 135*> \author Univ. of California Berkeley 136*> \author Univ. of Colorado Denver 137*> \author NAG Ltd. 138* 139*> \ingroup complex16_eig 140* 141* ===================================================================== 142 SUBROUTINE ZGET22( TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, W, 143 $ WORK, RWORK, RESULT ) 144* 145* -- LAPACK test routine -- 146* -- LAPACK is a software package provided by Univ. of Tennessee, -- 147* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 148* 149* .. Scalar Arguments .. 150 CHARACTER TRANSA, TRANSE, TRANSW 151 INTEGER LDA, LDE, N 152* .. 153* .. Array Arguments .. 154 DOUBLE PRECISION RESULT( 2 ), RWORK( * ) 155 COMPLEX*16 A( LDA, * ), E( LDE, * ), W( * ), WORK( * ) 156* .. 157* 158* ===================================================================== 159* 160* .. Parameters .. 161 DOUBLE PRECISION ZERO, ONE 162 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 163 COMPLEX*16 CZERO, CONE 164 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 165 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 166* .. 167* .. Local Scalars .. 168 CHARACTER NORMA, NORME 169 INTEGER ITRNSE, ITRNSW, J, JCOL, JOFF, JROW, JVEC 170 DOUBLE PRECISION ANORM, ENORM, ENRMAX, ENRMIN, ERRNRM, TEMP1, 171 $ ULP, UNFL 172 COMPLEX*16 WTEMP 173* .. 174* .. External Functions .. 175 LOGICAL LSAME 176 DOUBLE PRECISION DLAMCH, ZLANGE 177 EXTERNAL LSAME, DLAMCH, ZLANGE 178* .. 179* .. External Subroutines .. 180 EXTERNAL ZGEMM, ZLASET 181* .. 182* .. Intrinsic Functions .. 183 INTRINSIC ABS, DBLE, DCONJG, DIMAG, MAX, MIN 184* .. 185* .. Executable Statements .. 186* 187* Initialize RESULT (in case N=0) 188* 189 RESULT( 1 ) = ZERO 190 RESULT( 2 ) = ZERO 191 IF( N.LE.0 ) 192 $ RETURN 193* 194 UNFL = DLAMCH( 'Safe minimum' ) 195 ULP = DLAMCH( 'Precision' ) 196* 197 ITRNSE = 0 198 ITRNSW = 0 199 NORMA = 'O' 200 NORME = 'O' 201* 202 IF( LSAME( TRANSA, 'T' ) .OR. LSAME( TRANSA, 'C' ) ) THEN 203 NORMA = 'I' 204 END IF 205* 206 IF( LSAME( TRANSE, 'T' ) ) THEN 207 ITRNSE = 1 208 NORME = 'I' 209 ELSE IF( LSAME( TRANSE, 'C' ) ) THEN 210 ITRNSE = 2 211 NORME = 'I' 212 END IF 213* 214 IF( LSAME( TRANSW, 'C' ) ) THEN 215 ITRNSW = 1 216 END IF 217* 218* Normalization of E: 219* 220 ENRMIN = ONE / ULP 221 ENRMAX = ZERO 222 IF( ITRNSE.EQ.0 ) THEN 223 DO 20 JVEC = 1, N 224 TEMP1 = ZERO 225 DO 10 J = 1, N 226 TEMP1 = MAX( TEMP1, ABS( DBLE( E( J, JVEC ) ) )+ 227 $ ABS( DIMAG( E( J, JVEC ) ) ) ) 228 10 CONTINUE 229 ENRMIN = MIN( ENRMIN, TEMP1 ) 230 ENRMAX = MAX( ENRMAX, TEMP1 ) 231 20 CONTINUE 232 ELSE 233 DO 30 JVEC = 1, N 234 RWORK( JVEC ) = ZERO 235 30 CONTINUE 236* 237 DO 50 J = 1, N 238 DO 40 JVEC = 1, N 239 RWORK( JVEC ) = MAX( RWORK( JVEC ), 240 $ ABS( DBLE( E( JVEC, J ) ) )+ 241 $ ABS( DIMAG( E( JVEC, J ) ) ) ) 242 40 CONTINUE 243 50 CONTINUE 244* 245 DO 60 JVEC = 1, N 246 ENRMIN = MIN( ENRMIN, RWORK( JVEC ) ) 247 ENRMAX = MAX( ENRMAX, RWORK( JVEC ) ) 248 60 CONTINUE 249 END IF 250* 251* Norm of A: 252* 253 ANORM = MAX( ZLANGE( NORMA, N, N, A, LDA, RWORK ), UNFL ) 254* 255* Norm of E: 256* 257 ENORM = MAX( ZLANGE( NORME, N, N, E, LDE, RWORK ), ULP ) 258* 259* Norm of error: 260* 261* Error = AE - EW 262* 263 CALL ZLASET( 'Full', N, N, CZERO, CZERO, WORK, N ) 264* 265 JOFF = 0 266 DO 100 JCOL = 1, N 267 IF( ITRNSW.EQ.0 ) THEN 268 WTEMP = W( JCOL ) 269 ELSE 270 WTEMP = DCONJG( W( JCOL ) ) 271 END IF 272* 273 IF( ITRNSE.EQ.0 ) THEN 274 DO 70 JROW = 1, N 275 WORK( JOFF+JROW ) = E( JROW, JCOL )*WTEMP 276 70 CONTINUE 277 ELSE IF( ITRNSE.EQ.1 ) THEN 278 DO 80 JROW = 1, N 279 WORK( JOFF+JROW ) = E( JCOL, JROW )*WTEMP 280 80 CONTINUE 281 ELSE 282 DO 90 JROW = 1, N 283 WORK( JOFF+JROW ) = DCONJG( E( JCOL, JROW ) )*WTEMP 284 90 CONTINUE 285 END IF 286 JOFF = JOFF + N 287 100 CONTINUE 288* 289 CALL ZGEMM( TRANSA, TRANSE, N, N, N, CONE, A, LDA, E, LDE, -CONE, 290 $ WORK, N ) 291* 292 ERRNRM = ZLANGE( 'One', N, N, WORK, N, RWORK ) / ENORM 293* 294* Compute RESULT(1) (avoiding under/overflow) 295* 296 IF( ANORM.GT.ERRNRM ) THEN 297 RESULT( 1 ) = ( ERRNRM / ANORM ) / ULP 298 ELSE 299 IF( ANORM.LT.ONE ) THEN 300 RESULT( 1 ) = ONE / ULP 301 ELSE 302 RESULT( 1 ) = MIN( ERRNRM / ANORM, ONE ) / ULP 303 END IF 304 END IF 305* 306* Compute RESULT(2) : the normalization error in E. 307* 308 RESULT( 2 ) = MAX( ABS( ENRMAX-ONE ), ABS( ENRMIN-ONE ) ) / 309 $ ( DBLE( N )*ULP ) 310* 311 RETURN 312* 313* End of ZGET22 314* 315 END 316