1 SUBROUTINE PZGET22( TRANSA, TRANSE, TRANSW, N, A, DESCA, E, DESCE, 2 $ W, WORK, DESCW, RWORK, RESULT ) 3* 4* -- ScaLAPACK testing routine (version 1.7) -- 5* University of Tennessee, Knoxville, Oak Ridge National Laboratory, 6* and University of California, Berkeley. 7* August 14, 2001 8* 9* .. Scalar Arguments .. 10 CHARACTER TRANSA, TRANSE, TRANSW 11 INTEGER N 12* .. 13* .. Array Arguments .. 14 INTEGER DESCA( * ), DESCE( * ), DESCW( * ) 15 DOUBLE PRECISION RESULT( 2 ), RWORK( * ) 16 COMPLEX*16 A( * ), E( * ), W( * ), WORK( * ) 17* .. 18* 19* Purpose 20* ======= 21* 22* PZGET22 does an eigenvector check. 23* 24* The basic test is: 25* 26* RESULT(1) = | A E - E W | / ( |A| |E| ulp ) 27* 28* using the 1-norm. It also tests the normalization of E: 29* 30* RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp ) 31* j 32* 33* where E(j) is the j-th eigenvector, and m-norm is the max-norm of a 34* vector. The max-norm of a complex n-vector x in this case is the 35* maximum of |re(x(i)| + |im(x(i)| over i = 1, ..., n. 36* 37* Arguments 38* ========== 39* 40* TRANSA (input) CHARACTER*1 41* Specifies whether or not A is transposed. 42* = 'N': No transpose 43* = 'T': Transpose 44* = 'C': Conjugate transpose 45* 46* TRANSE (input) CHARACTER*1 47* Specifies whether or not E is transposed. 48* = 'N': No transpose, eigenvectors are in columns of E 49* = 'T': Transpose, eigenvectors are in rows of E 50* = 'C': Conjugate transpose, eigenvectors are in rows of E 51* 52* TRANSW (input) CHARACTER*1 53* Specifies whether or not W is transposed. 54* = 'N': No transpose 55* = 'T': Transpose, same as TRANSW = 'N' 56* = 'C': Conjugate transpose, use -WI(j) instead of WI(j) 57* 58* N (input) INTEGER 59* The order of the matrix A. N >= 0. 60* 61* A (input) COMPLEX*16 array, dimension (*) 62* The matrix whose eigenvectors are in E. 63* 64* DESCA (input) INTEGER array, dimension(*) 65* 66* E (input) COMPLEX*16 array, dimension (*) 67* The matrix of eigenvectors. If TRANSE = 'N', the eigenvectors 68* are stored in the columns of E, if TRANSE = 'T' or 'C', the 69* eigenvectors are stored in the rows of E. 70* 71* DESCE (input) INTEGER array, dimension(*) 72* 73* W (input) COMPLEX*16 array, dimension (N) 74* The eigenvalues of A. 75* 76* WORK (workspace) COMPLEX*16 array, dimension (*) 77* DESCW (input) INTEGER array, dimension(*) 78* 79* RWORK (workspace) DOUBLE PRECISION array, dimension (N) 80* 81* RESULT (output) DOUBLE PRECISION array, dimension (2) 82* RESULT(1) = | A E - E W | / ( |A| |E| ulp ) 83* RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp ) 84* j 85* Further Details 86* =============== 87* 88* Contributed by Mark Fahey, June, 2000 89* 90* ===================================================================== 91* 92* .. Parameters .. 93 INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_, 94 $ MB_, NB_, RSRC_, CSRC_, LLD_ 95 PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, 96 $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 97 $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) 98 DOUBLE PRECISION ZERO, ONE 99 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 100 COMPLEX*16 CZERO, CONE 101 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), 102 $ CONE = ( 1.0D+0, 0.0D+0 ) ) 103* .. 104* .. Local Scalars .. 105 CHARACTER NORMA, NORME 106 INTEGER ICOL, II, IROW, ITRNSE, ITRNSW, J, JCOL, JJ, 107 $ JROW, JVEC, LDA, LDE, LDW, MB, MYCOL, MYROW, 108 $ NB, NPCOL, NPROW, CONTXT, CA, CSRC, RA, RSRC 109 DOUBLE PRECISION ANORM, ENORM, ENRMAX, ENRMIN, ERRNRM, TEMP1, 110 $ ULP, UNFL 111 COMPLEX*16 CDUM, WTEMP 112* .. 113* .. External Functions .. 114 LOGICAL LSAME 115 DOUBLE PRECISION PDLAMCH, PZLANGE 116 EXTERNAL LSAME, PDLAMCH, PZLANGE 117* .. 118* .. External Subroutines .. 119 EXTERNAL BLACS_GRIDINFO, DGAMN2D, DGAMX2D, INFOG2L, 120 $ PZAXPY, PZGEMM, PZLASET 121* .. 122* .. Intrinsic Functions .. 123 INTRINSIC ABS, DBLE, DCONJG, DIMAG, MAX, MIN 124* .. 125* .. Statement Functions .. 126 DOUBLE PRECISION CABS1 127* .. 128* .. Statement Function definitions .. 129 CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) ) 130* .. 131* .. Executable Statements .. 132* 133* Initialize RESULT (in case N=0) 134* 135 RESULT( 1 ) = ZERO 136 RESULT( 2 ) = ZERO 137 IF( N.LE.0 ) 138 $ RETURN 139* 140 CONTXT = DESCA( CTXT_ ) 141 RSRC = DESCA( RSRC_ ) 142 CSRC = DESCA( CSRC_ ) 143 NB = DESCA( NB_ ) 144 MB = DESCA( MB_ ) 145 LDA = DESCA( LLD_ ) 146 LDE = DESCE( LLD_ ) 147 LDW = DESCW( LLD_ ) 148 CALL BLACS_GRIDINFO( CONTXT, NPROW, NPCOL, MYROW, MYCOL ) 149* 150 UNFL = PDLAMCH( CONTXT, 'Safe minimum' ) 151 ULP = PDLAMCH( CONTXT, 'Precision' ) 152* 153 ITRNSE = 0 154 ITRNSW = 0 155 NORMA = 'O' 156 NORME = 'O' 157* 158 IF( LSAME( TRANSA, 'T' ) .OR. LSAME( TRANSA, 'C' ) ) THEN 159 NORMA = 'I' 160 END IF 161* 162 IF( LSAME( TRANSE, 'T' ) ) THEN 163 ITRNSE = 1 164 NORME = 'I' 165 ELSE IF( LSAME( TRANSE, 'C' ) ) THEN 166 ITRNSE = 2 167 NORME = 'I' 168 END IF 169* 170 IF( LSAME( TRANSW, 'C' ) ) THEN 171 ITRNSW = 1 172 END IF 173* 174* Normalization of E: 175* 176 ENRMIN = ONE / ULP 177 ENRMAX = ZERO 178 IF( ITRNSE.EQ.0 ) THEN 179 DO 20 JVEC = 1, N 180 TEMP1 = ZERO 181 DO 10 J = 1, N 182 CALL INFOG2L( J, JVEC, DESCE, NPROW, NPCOL, MYROW, MYCOL, 183 $ IROW, ICOL, II, JJ ) 184 IF( ( MYROW.EQ.II ) .AND. ( MYCOL.EQ.JJ ) ) THEN 185 TEMP1 = MAX( TEMP1, CABS1( E( ( ICOL-1 )*LDE+ 186 $ IROW ) ) ) 187 END IF 188 10 CONTINUE 189 IF( MYCOL.EQ.JJ ) THEN 190 CALL DGAMX2D( CONTXT, 'Col', ' ', 1, 1, TEMP1, 1, RA, CA, 191 $ -1, -1, -1 ) 192 ENRMIN = MIN( ENRMIN, TEMP1 ) 193 ENRMAX = MAX( ENRMAX, TEMP1 ) 194 END IF 195 20 CONTINUE 196 CALL DGAMX2D( CONTXT, 'Row', ' ', 1, 1, ENRMAX, 1, RA, CA, -1, 197 $ -1, -1 ) 198 CALL DGAMN2D( CONTXT, 'Row', ' ', 1, 1, ENRMIN, 1, RA, CA, -1, 199 $ -1, -1 ) 200 ELSE 201 DO 40 J = 1, N 202 TEMP1 = ZERO 203 DO 30 JVEC = 1, N 204 CALL INFOG2L( J, JVEC, DESCE, NPROW, NPCOL, MYROW, MYCOL, 205 $ IROW, ICOL, II, JJ ) 206 IF( ( MYROW.EQ.II ) .AND. ( MYCOL.EQ.JJ ) ) THEN 207 TEMP1 = MAX( TEMP1, CABS1( E( ( ICOL-1 )*LDE+ 208 $ IROW ) ) ) 209 END IF 210 30 CONTINUE 211 IF( MYROW.EQ.II ) THEN 212 CALL DGAMX2D( CONTXT, 'Row', ' ', 1, 1, TEMP1, 1, RA, CA, 213 $ -1, -1, -1 ) 214 ENRMIN = MIN( ENRMIN, TEMP1 ) 215 ENRMAX = MAX( ENRMAX, TEMP1 ) 216 END IF 217 40 CONTINUE 218 CALL DGAMX2D( CONTXT, 'Row', ' ', 1, 1, ENRMAX, 1, RA, CA, -1, 219 $ -1, -1 ) 220 CALL DGAMN2D( CONTXT, 'Row', ' ', 1, 1, ENRMIN, 1, RA, CA, -1, 221 $ -1, -1 ) 222 END IF 223* 224* Norm of A: 225* 226 ANORM = MAX( PZLANGE( NORMA, N, N, A, 1, 1, DESCA, RWORK ), UNFL ) 227* 228* Norm of E: 229* 230 ENORM = MAX( PZLANGE( NORME, N, N, E, 1, 1, DESCE, RWORK ), ULP ) 231* 232* Norm of error: 233* 234* Error = AE - EW 235* 236 CALL PZLASET( 'Full', N, N, CZERO, CZERO, WORK, 1, 1, DESCW ) 237* 238 DO 60 JCOL = 1, N 239 IF( ITRNSW.EQ.0 ) THEN 240 WTEMP = W( JCOL ) 241 ELSE 242 WTEMP = DCONJG( W( JCOL ) ) 243 END IF 244* 245 IF( ITRNSE.EQ.0 ) THEN 246 CALL PZAXPY( N, WTEMP, E, 1, JCOL, DESCE, 1, WORK, 1, JCOL, 247 $ DESCW, 1 ) 248 ELSE IF( ITRNSE.EQ.1 ) THEN 249 CALL PZAXPY( N, WTEMP, E, JCOL, 1, DESCE, N, WORK, 1, JCOL, 250 $ DESCW, 1 ) 251 ELSE 252 CALL PZAXPY( N, DCONJG( WTEMP ), E, JCOL, 1, DESCE, N, WORK, 253 $ 1, JCOL, DESCW, 1 ) 254 DO 50 JROW = 1, N 255 CALL INFOG2L( JROW, JCOL, DESCW, NPROW, NPCOL, MYROW, 256 $ MYCOL, IROW, ICOL, II, JJ ) 257 IF( ( MYROW.EQ.II ) .AND. ( MYCOL.EQ.JJ ) ) THEN 258 WORK( ( JCOL-1 )*LDW+JROW ) 259 $ = DCONJG( WORK( ( JCOL-1 )*LDW+JROW ) ) 260 END IF 261 50 CONTINUE 262 END IF 263 60 CONTINUE 264* 265 CALL PZGEMM( TRANSA, TRANSE, N, N, N, CONE, A, 1, 1, DESCA, E, 1, 266 $ 1, DESCE, -CONE, WORK, 1, 1, DESCW ) 267* 268 ERRNRM = PZLANGE( 'One', N, N, WORK, 1, 1, DESCW, RWORK ) / ENORM 269* 270* Compute RESULT(1) (avoiding under/overflow) 271* 272 IF( ANORM.GT.ERRNRM ) THEN 273 RESULT( 1 ) = ( ERRNRM / ANORM ) / ULP 274 ELSE 275 IF( ANORM.LT.ONE ) THEN 276 RESULT( 1 ) = ( MIN( ERRNRM, ANORM ) / ANORM ) / ULP 277 ELSE 278 RESULT( 1 ) = MIN( ERRNRM / ANORM, ONE ) / ULP 279 END IF 280 END IF 281* 282* Compute RESULT(2) : the normalization error in E. 283* 284 RESULT( 2 ) = MAX( ABS( ENRMAX-ONE ), ABS( ENRMIN-ONE ) ) / 285 $ ( DBLE( N )*ULP ) 286* 287 RETURN 288* 289* End of PZGET22 290* 291 END 292