|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
00001 SUBROUTINE PZSDPSUBTST( WKNOWN, UPLO, N, THRESH, ABSTOL, A, COPYA, 00002 $ Z, IA, JA, DESCA, WIN, WNEW, IPREPAD, 00003 $ IPOSTPAD, WORK, LWORK, RWORK, LRWORK, 00004 $ LWORK1, IWORK, LIWORK, RESULT, TSTNRM, 00005 $ QTQNRM, NOUT ) 00006 * 00007 * -- ScaLAPACK testing routine (version 1.7) -- 00008 * University of Tennessee, Knoxville, Oak Ridge National Laboratory, 00009 * and University of California, Berkeley. 00010 * February 28, 2000 00011 * 00012 * .. Scalar Arguments .. 00013 LOGICAL WKNOWN 00014 CHARACTER UPLO 00015 INTEGER IA, IPOSTPAD, IPREPAD, JA, LIWORK, LRWORK, 00016 $ LWORK, LWORK1, N, NOUT, RESULT 00017 DOUBLE PRECISION ABSTOL, QTQNRM, THRESH, TSTNRM 00018 * .. 00019 * .. Array Arguments .. 00020 INTEGER DESCA( * ), IWORK( * ) 00021 DOUBLE PRECISION RWORK( * ), WIN( * ), WNEW( * ) 00022 COMPLEX*16 A( * ), COPYA( * ), WORK( * ), Z( * ) 00023 * .. 00024 * 00025 * Purpose 00026 * ======= 00027 * 00028 * PZSDPSUBTST calls PZHEEVD and then tests the output of 00029 * PZHEEVD 00030 * If JOBZ = 'V' then the following two tests are performed: 00031 * |AQ -QL| / (abstol + eps * norm(A) ) < N*THRESH 00032 * |QT * Q - I| / eps * norm(A) < N*THRESH 00033 * If WKNOWN then 00034 * we check to make sure that the eigenvalues match expectations 00035 * i.e. |WIN - WNEW(1+IPREPAD)| / (eps * |WIN|) < THRESH 00036 * where WIN is the array of eigenvalues as computed by 00037 * PZHEEVD when eigenvectors are requested 00038 * 00039 * Arguments 00040 * ========= 00041 * 00042 * NP = the number of rows local to a given process. 00043 * NQ = the number of columns local to a given process. 00044 * 00045 * WKNOWN (global input) INTEGER 00046 * .FALSE.: WIN does not contain the eigenvalues 00047 * .TRUE.: WIN does contain the eigenvalues 00048 * 00049 * UPLO (global input) CHARACTER*1 00050 * Specifies whether the upper or lower triangular part of the 00051 * Hermitian matrix A is stored: 00052 * = 'U': Upper triangular 00053 * = 'L': Lower triangular 00054 * 00055 * N (global input) INTEGER 00056 * Size of the matrix to be tested. (global size) 00057 * 00058 * THRESH (global input) DOUBLE PRECISION 00059 * A test will count as "failed" if the "error", computed as 00060 * described below, exceeds THRESH. Note that the error 00061 * is scaled to be O(1), so THRESH should be a reasonably 00062 * small multiple of 1, e.g., 10 or 100. In particular, 00063 * it should not depend on the precision (single vs. double) 00064 * or the size of the matrix. It must be at least zero. 00065 * 00066 * ABSTOL (global input) DOUBLE PRECISION 00067 * The absolute tolerance for the eigenvalues. An 00068 * eigenvalue is considered to be located if it has 00069 * been determined to lie in an interval whose width 00070 * is "abstol" or less. If "abstol" is less than or equal 00071 * to zero, then ulp*|T| will be used, where |T| is 00072 * the 1-norm of the matrix. If eigenvectors are 00073 * desired later by inverse iteration ("PZSTEIN"), 00074 * "abstol" MUST NOT be bigger than ulp*|T|. 00075 * 00076 * A (local workspace) COMPLEX*16 array 00077 * global dimension (N, N), local dimension (DESCA(DLEN_), NQ) 00078 * A is distributed in a block cyclic manner over both rows 00079 * and columns. 00080 * See PZHEEVD for a description of block cyclic layout. 00081 * The test matrix, which is then modified by PZHEEVD 00082 * A has already been padded front and back, use A(1+IPREPAD) 00083 * 00084 * COPYA (local input) COMPLEX*16 array, dimension(N*N) 00085 * COPYA holds a copy of the original matrix A 00086 * identical in both form and content to A 00087 * 00088 * Z (local workspace) COMPLEX*16 array, dim (N*N) 00089 * Z is distributed in the same manner as A 00090 * Z contains the eigenvector matrix 00091 * Z is used as workspace by the test routines 00092 * PZSEPCHK and PZSEPQTQ. 00093 * Z has already been padded front and back, use Z(1+IPREPAD) 00094 * 00095 * IA (global input) INTEGER 00096 * On entry, IA specifies the global row index of the submatrix 00097 * of the global matrix A, COPYA and Z to operate on. 00098 * 00099 * JA (global input) INTEGER 00100 * On entry, IA specifies the global column index of the submat 00101 * of the global matrix A, COPYA and Z to operate on. 00102 * 00103 * DESCA (global/local input) INTEGER array of dimension 8 00104 * The array descriptor for the matrix A, COPYA and Z. 00105 * 00106 * WIN (global input) DOUBLE PRECISION array, dimension (N) 00107 * If .not. WKNOWN, WIN is ignored on input 00108 * Otherwise, WIN() is taken as the standard by which the 00109 * eigenvalues are to be compared against. 00110 * 00111 * WNEW (global workspace) DOUBLE PRECISION array, dimension (N) 00112 * The eigenvalues as copmuted by this call to PZHEEVD 00113 * If JOBZ <> 'V' or RANGE <> 'A' these eigenvalues are 00114 * compared against those in WIN(). 00115 * WNEW has already been padded front and back, 00116 * use WNEW(1+IPREPAD) 00117 * 00118 * WORK (local workspace) COMPLEX*16 array, dimension (LWORK) 00119 * WORK has already been padded front and back, 00120 * use WORK(1+IPREPAD) 00121 * 00122 * LWORK (local input) INTEGER 00123 * The actual length of the array WORK after padding. 00124 * 00125 * RWORK (local workspace) DOUBLE PRECISION array, dimension (LRWORK) 00126 * RWORK has already been padded front and back, 00127 * use RWORK(1+IPREPAD) 00128 * 00129 * LRWORK (local input) INTEGER 00130 * The actual length of the array RWORK after padding. 00131 * 00132 * LWORK1 (local input) INTEGER 00133 * The amount of DOUBLE PRECISION workspace to pass to PZHEEVD 00134 * 00135 * IWORK (local workspace) INTEGER array, dimension (LIWORK) 00136 * IWORK has already been padded front and back, 00137 * use IWORK(1+IPREPAD) 00138 * 00139 * LIWORK (local input) INTEGER 00140 * The length of the array IWORK after padding. 00141 * 00142 * RESULT (global output) INTEGER 00143 * The result of this call to PZHEEVD 00144 * RESULT = -3 => This process did not participate 00145 * RESULT = 0 => All tests passed 00146 * RESULT = 1 => ONe or more tests failed 00147 * 00148 * TSTNRM (global output) DOUBLE PRECISION 00149 * |AQ- QL| / (ABSTOL+EPS*|A|)*N 00150 * 00151 * QTQNRM (global output) DOUBLE PRECISION 00152 * |QTQ -I| / N*EPS 00153 * 00154 * .. Parameters .. 00155 * 00156 INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_, 00157 $ MB_, NB_, RSRC_, CSRC_, LLD_ 00158 PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, 00159 $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 00160 $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) 00161 DOUBLE PRECISION PADVAL, FIVE, NEGONE 00162 PARAMETER ( PADVAL = 13.5285D+0, FIVE = 5.0D+0, 00163 $ NEGONE = -1.0D+0 ) 00164 COMPLEX*16 CPADVAL 00165 PARAMETER ( CPADVAL = ( 13.989D+0, 1.93D+0 ) ) 00166 INTEGER IPADVAL 00167 PARAMETER ( IPADVAL = 927 ) 00168 COMPLEX*16 CZERO, CONE, CNEGONE 00169 PARAMETER ( CZERO = 0.0D+0, CONE = 1.0D+0, 00170 $ CNEGONE = -1.0D+0 ) 00171 * .. 00172 * .. Local Scalars .. 00173 INTEGER I, IAM, INFO, ISIZEHEEVD, ISIZEHEEVX, 00174 $ ISIZESUBTST, ISIZETST, MYCOL, MYROW, NP, NPCOL, 00175 $ NPROW, NQ, RES, RSIZECHK, RSIZEHEEVD, 00176 $ RSIZEHEEVX, RSIZEQTQ, RSIZESUBTST, RSIZETST, 00177 $ SIZEHEEVD, SIZEHEEVX, SIZEMQRLEFT, 00178 $ SIZEMQRRIGHT, SIZEQRF, SIZESUBTST, SIZETMS, 00179 $ SIZETST 00180 DOUBLE PRECISION EPS, EPSNORMA, ERROR, MAXERROR, MINERROR, NORM, 00181 $ NORMWIN, SAFMIN, ULP 00182 * .. 00183 * .. Local Arrays .. 00184 INTEGER ITMP( 2 ) 00185 * .. 00186 * .. External Functions .. 00187 * 00188 INTEGER NUMROC 00189 DOUBLE PRECISION PZLANGE, PZLANHE, PDLAMCH 00190 EXTERNAL NUMROC, PZLANGE, PZLANHE, PDLAMCH 00191 * .. 00192 * .. External Subroutines .. 00193 EXTERNAL BLACS_GRIDINFO, ZLACPY, IGAMN2D, IGAMX2D, 00194 $ PZCHEKPAD, PZFILLPAD, PZGEMM, PZHEEVD, PZLASET, 00195 $ PZLASIZESEP, PZSEPCHK, PICHEKPAD, PIFILLPAD, 00196 $ PDCHEKPAD, PDFILLPAD, SLBOOT, SLTIMER 00197 * .. 00198 * .. Intrinsic Functions .. 00199 INTRINSIC ABS, MAX, MIN, DBLE 00200 * .. 00201 * .. Executable Statements .. 00202 * This is just to keep ftnchek happy 00203 IF( BLOCK_CYCLIC_2D*CSRC_*CTXT_*DLEN_*DTYPE_*LLD_*MB_*M_*NB_*N_* 00204 $ RSRC_.LT.0 )RETURN 00205 * 00206 CALL PZLASIZESEP( DESCA, IPREPAD, IPOSTPAD, SIZEMQRLEFT, 00207 $ SIZEMQRRIGHT, SIZEQRF, SIZETMS, RSIZEQTQ, 00208 $ RSIZECHK, SIZEHEEVX, RSIZEHEEVX, ISIZEHEEVX, 00209 $ SIZEHEEVD, RSIZEHEEVD, ISIZEHEEVD, SIZESUBTST, 00210 $ RSIZESUBTST, ISIZESUBTST, SIZETST, RSIZETST, 00211 $ ISIZETST ) 00212 * 00213 TSTNRM = NEGONE 00214 QTQNRM = NEGONE 00215 EPS = PDLAMCH( DESCA( CTXT_ ), 'Eps' ) 00216 SAFMIN = PDLAMCH( DESCA( CTXT_ ), 'Safe min' ) 00217 * 00218 NORMWIN = SAFMIN / EPS 00219 IF( N.GE.1 ) 00220 $ NORMWIN = MAX( ABS( WIN( 1+IPREPAD ) ), 00221 $ ABS( WIN( N+IPREPAD ) ), NORMWIN ) 00222 * 00223 DO 10 I = 1, LWORK1, 1 00224 RWORK( I+IPREPAD ) = 14.3D+0 00225 10 CONTINUE 00226 DO 20 I = 1, LIWORK, 1 00227 IWORK( I ) = 14 00228 20 CONTINUE 00229 DO 30 I = 1, LWORK, 1 00230 WORK( I+IPREPAD ) = ( 15.63D+0, 1.1D+0 ) 00231 30 CONTINUE 00232 * 00233 DO 40 I = 1, N 00234 WNEW( I+IPREPAD ) = 3.14159D+0 00235 40 CONTINUE 00236 * 00237 CALL BLACS_GRIDINFO( DESCA( CTXT_ ), NPROW, NPCOL, MYROW, MYCOL ) 00238 * 00239 IAM = 1 00240 IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) 00241 $ IAM = 0 00242 * 00243 * If this process is not involved in this test, bail out now 00244 * 00245 RESULT = -3 00246 IF( MYROW.GE.NPROW .OR. MYROW.LT.0 ) 00247 $ GO TO 60 00248 RESULT = 0 00249 * 00250 NP = NUMROC( N, DESCA( MB_ ), MYROW, 0, NPROW ) 00251 NQ = NUMROC( N, DESCA( NB_ ), MYCOL, 0, NPCOL ) 00252 * 00253 CALL ZLACPY( 'A', NP, NQ, COPYA, DESCA( LLD_ ), A( 1+IPREPAD ), 00254 $ DESCA( LLD_ ) ) 00255 * 00256 CALL PZFILLPAD( DESCA( CTXT_ ), NP, NQ, A, DESCA( LLD_ ), IPREPAD, 00257 $ IPOSTPAD, CPADVAL ) 00258 * 00259 CALL PZFILLPAD( DESCA( CTXT_ ), NP, NQ, Z, DESCA( LLD_ ), IPREPAD, 00260 $ IPOSTPAD, CPADVAL+1.0D+0 ) 00261 * 00262 CALL PDFILLPAD( DESCA( CTXT_ ), N, 1, WNEW, N, IPREPAD, IPOSTPAD, 00263 $ PADVAL+2.0D+0 ) 00264 * 00265 CALL PDFILLPAD( DESCA( CTXT_ ), LWORK1, 1, RWORK, LWORK1, IPREPAD, 00266 $ IPOSTPAD, PADVAL+4.0D+0 ) 00267 * 00268 CALL PIFILLPAD( DESCA( CTXT_ ), LIWORK, 1, IWORK, LIWORK, IPREPAD, 00269 $ IPOSTPAD, IPADVAL ) 00270 * 00271 CALL PZFILLPAD( DESCA( CTXT_ ), LWORK, 1, WORK, LWORK, IPREPAD, 00272 $ IPOSTPAD, CPADVAL+4.1D+0 ) 00273 * 00274 CALL SLBOOT 00275 CALL SLTIMER( 1 ) 00276 CALL SLTIMER( 6 ) 00277 * 00278 CALL PZHEEVD( 'V', UPLO, N, A( 1+IPREPAD ), IA, JA, DESCA, 00279 $ WNEW( 1+IPREPAD ), Z( 1+IPREPAD ), IA, JA, DESCA, 00280 $ WORK( 1+IPREPAD ), SIZEHEEVD, RWORK( 1+IPREPAD ), 00281 $ LWORK1, IWORK( 1+IPREPAD ), LIWORK, INFO ) 00282 CALL SLTIMER( 6 ) 00283 CALL SLTIMER( 1 ) 00284 * 00285 IF( THRESH.LE.0 ) THEN 00286 RESULT = 0 00287 ELSE 00288 CALL PZCHEKPAD( DESCA( CTXT_ ), 'PZHEEVD-A', NP, NQ, A, 00289 $ DESCA( LLD_ ), IPREPAD, IPOSTPAD, CPADVAL ) 00290 * 00291 CALL PZCHEKPAD( DESCA( CTXT_ ), 'PZHEEVD-Z', NP, NQ, Z, 00292 $ DESCA( LLD_ ), IPREPAD, IPOSTPAD, 00293 $ CPADVAL+1.0D+0 ) 00294 * 00295 CALL PDCHEKPAD( DESCA( CTXT_ ), 'PZHEEVD-WNEW', N, 1, WNEW, N, 00296 $ IPREPAD, IPOSTPAD, PADVAL+2.0D+0 ) 00297 * 00298 CALL PDCHEKPAD( DESCA( CTXT_ ), 'PZHEEVD-rWORK', LWORK1, 1, 00299 $ RWORK, LWORK1, IPREPAD, IPOSTPAD, 00300 $ PADVAL+4.0D+0 ) 00301 * 00302 CALL PZCHEKPAD( DESCA( CTXT_ ), 'PZHEEVD-WORK', LWORK, 1, WORK, 00303 $ LWORK, IPREPAD, IPOSTPAD, CPADVAL+4.1D+0 ) 00304 * 00305 CALL PICHEKPAD( DESCA( CTXT_ ), 'PZHEEVD-IWORK', LIWORK, 1, 00306 $ IWORK, LIWORK, IPREPAD, IPOSTPAD, IPADVAL ) 00307 * 00308 * Check INFO 00309 * 00310 * Make sure that all processes return the same value of INFO 00311 * 00312 ITMP( 1 ) = INFO 00313 ITMP( 2 ) = INFO 00314 * 00315 CALL IGAMN2D( DESCA( CTXT_ ), 'a', ' ', 1, 1, ITMP, 1, 1, 1, 00316 $ -1, -1, 0 ) 00317 CALL IGAMX2D( DESCA( CTXT_ ), 'a', ' ', 1, 1, ITMP( 2 ), 1, 1, 00318 $ 1, -1, -1, 0 ) 00319 * 00320 * 00321 IF( ITMP( 1 ).NE.ITMP( 2 ) ) THEN 00322 IF( IAM.EQ.0 ) 00323 $ WRITE( NOUT, FMT = * ) 00324 $ 'Different processes return different INFO' 00325 RESULT = 1 00326 ELSE IF( INFO.NE.0 ) THEN 00327 IF( IAM.EQ.0 ) 00328 $ WRITE( NOUT, FMT = 9996 )INFO 00329 RESULT = 1 00330 END IF 00331 * 00332 * Compute eps * norm(A) 00333 * 00334 IF( N.EQ.0 ) THEN 00335 EPSNORMA = EPS 00336 ELSE 00337 EPSNORMA = PZLANHE( 'I', UPLO, N, COPYA, IA, JA, DESCA, 00338 $ RWORK )*EPS 00339 END IF 00340 * 00341 * Note that a couple key variables get redefined in PZSEPCHK 00342 * as described by this table: 00343 * 00344 * PZSEPTST name PZSEPCHK name 00345 * ------------- ------------- 00346 * COPYA A 00347 * Z Q 00348 * A C 00349 * 00350 * Perform the |AQ - QE| test 00351 * 00352 CALL PDFILLPAD( DESCA( CTXT_ ), RSIZECHK, 1, RWORK, RSIZECHK, 00353 $ IPREPAD, IPOSTPAD, 4.3D+0 ) 00354 * 00355 CALL PZSEPCHK( N, N, COPYA, IA, JA, DESCA, 00356 $ MAX( ABSTOL+EPSNORMA, SAFMIN ), THRESH, 00357 $ Z( 1+IPREPAD ), IA, JA, DESCA, A( 1+IPREPAD ), 00358 $ IA, JA, DESCA, WNEW( 1+IPREPAD ), 00359 $ RWORK( 1+IPREPAD ), RSIZECHK, TSTNRM, RES ) 00360 * 00361 CALL PDCHEKPAD( DESCA( CTXT_ ), 'PZSDPCHK-rWORK', RSIZECHK, 1, 00362 $ RWORK, RSIZECHK, IPREPAD, IPOSTPAD, 4.3D+0 ) 00363 * 00364 IF( RES.NE.0 ) THEN 00365 RESULT = 1 00366 WRITE( NOUT, FMT = 9995 ) 00367 END IF 00368 * 00369 * Perform the |QTQ - I| test 00370 * 00371 CALL PDFILLPAD( DESCA( CTXT_ ), RSIZEQTQ, 1, RWORK, RSIZEQTQ, 00372 $ IPREPAD, IPOSTPAD, 4.3D+0 ) 00373 * 00374 * 00375 RES = 0 00376 ULP = PDLAMCH( DESCA( CTXT_ ), 'P' ) 00377 CALL PZLASET( 'A', N, N, CZERO, CONE, A( 1+IPREPAD ), IA, JA, 00378 $ DESCA ) 00379 CALL PZGEMM( 'Conjugate transpose', 'N', N, N, N, CNEGONE, 00380 $ Z( 1+IPREPAD ), IA, JA, DESCA, Z( 1+IPREPAD ), IA, 00381 $ JA, DESCA, CONE, A( 1+IPREPAD ), IA, JA, DESCA ) 00382 NORM = PZLANGE( '1', N, N, A( 1+IPREPAD ), IA, JA, DESCA, 00383 $ WORK( 1+IPREPAD ) ) 00384 QTQNRM = NORM / ( DBLE( MAX( N, 1 ) )*ULP ) 00385 IF( QTQNRM.GT.THRESH ) THEN 00386 RES = 1 00387 END IF 00388 CALL PDCHEKPAD( DESCA( CTXT_ ), 'PZSEPQTQ-rWORK', RSIZEQTQ, 1, 00389 $ RWORK, RSIZEQTQ, IPREPAD, IPOSTPAD, 4.3D+0 ) 00390 * 00391 IF( RES.NE.0 ) THEN 00392 RESULT = 1 00393 WRITE( NOUT, FMT = 9994 ) 00394 END IF 00395 * 00396 IF( INFO.NE.0 ) THEN 00397 IF( IAM.EQ.0 ) 00398 $ WRITE( NOUT, FMT = 9998 )INFO 00399 RESULT = 1 00400 END IF 00401 * 00402 IF( INFO.NE.0 ) THEN 00403 IF( IAM.EQ.0 ) 00404 $ WRITE( NOUT, FMT = 9998 )INFO 00405 RESULT = 1 00406 END IF 00407 END IF 00408 * 00409 * Check to make sure that we have the right eigenvalues 00410 * 00411 IF( WKNOWN .AND. N.GT.0 ) THEN 00412 * 00413 * Find the largest difference between the computed 00414 * and expected eigenvalues 00415 * 00416 MINERROR = NORMWIN 00417 MAXERROR = 0.0D+00 00418 * 00419 DO 50 I = 1, N 00420 ERROR = ABS( WIN( I+IPREPAD )-WNEW( I+IPREPAD ) ) 00421 MAXERROR = MAX( MAXERROR, ERROR ) 00422 50 CONTINUE 00423 MINERROR = MIN( MAXERROR, MINERROR ) 00424 * 00425 IF( MINERROR.GT.NORMWIN*FIVE*THRESH*EPS ) THEN 00426 IF( IAM.EQ.0 ) 00427 $ WRITE( NOUT, FMT = 9997 )MINERROR, NORMWIN 00428 RESULT = 1 00429 END IF 00430 END IF 00431 * 00432 * 00433 * All processes should report the same result 00434 * 00435 CALL IGAMX2D( DESCA( CTXT_ ), 'a', ' ', 1, 1, RESULT, 1, 1, 1, -1, 00436 $ -1, 0 ) 00437 * 00438 60 CONTINUE 00439 * 00440 RETURN 00441 * 00442 9999 FORMAT( 'PZHEEVD returned INFO=', I7 ) 00443 9998 FORMAT( 'PZSEPQTQ returned INFO=', I7 ) 00444 9997 FORMAT( 'PZSDPSUBTST minerror =', D11.2, ' normwin=', D11.2 ) 00445 9996 FORMAT( 'PZHEEVD returned INFO=', I7, 00446 $ ' despite adequate workspace' ) 00447 9995 FORMAT( 'PZHEEVD failed the |AQ -QE| test' ) 00448 9994 FORMAT( 'PZHEEVD failed the |QTQ -I| test' ) 00449 * 00450 * End of PZSDPSUBTST 00451 * 00452 END
