|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
00001 * 00002 * 00003 SUBROUTINE PCGSEPSUBTST( WKNOWN, IBTYPE, JOBZ, RANGE, UPLO, N, VL, 00004 $ VU, IL, IU, THRESH, ABSTOL, A, COPYA, B, 00005 $ COPYB, Z, IA, JA, DESCA, WIN, WNEW, 00006 $ IFAIL, ICLUSTR, GAP, IPREPAD, IPOSTPAD, 00007 $ WORK, LWORK, RWORK, LRWORK, LWORK1, 00008 $ IWORK, LIWORK, RESULT, TSTNRM, QTQNRM, 00009 $ NOUT ) 00010 * 00011 * -- ScaLAPACK routine (version 1.7) -- 00012 * University of Tennessee, Knoxville, Oak Ridge National Laboratory, 00013 * and University of California, Berkeley. 00014 * May 1, 1997 00015 * 00016 * .. Scalar Arguments .. 00017 LOGICAL WKNOWN 00018 CHARACTER JOBZ, RANGE, UPLO 00019 INTEGER IA, IBTYPE, IL, IPOSTPAD, IPREPAD, IU, JA, 00020 $ LIWORK, LRWORK, LWORK, LWORK1, N, NOUT, RESULT 00021 REAL ABSTOL, QTQNRM, THRESH, TSTNRM, VL, VU 00022 * .. 00023 * .. Array Arguments .. 00024 INTEGER DESCA( * ), ICLUSTR( * ), IFAIL( * ), 00025 $ IWORK( * ) 00026 REAL GAP( * ), RWORK( * ), WIN( * ), WNEW( * ) 00027 COMPLEX A( * ), B( * ), COPYA( * ), COPYB( * ), 00028 $ WORK( * ), Z( * ) 00029 * .. 00030 * 00031 * Purpose 00032 * ======= 00033 * 00034 * PCGSEPSUBTST calls PCHEGVX and then tests the output of 00035 * PCHEGVX 00036 * If JOBZ = 'V' then the following two tests are performed: 00037 * |AQ -QL| / (abstol + eps * norm(A) ) < THRESH 00038 * |QT * Q - I| / eps * norm(A) < THRESH 00039 * If WKNOWN then 00040 * we check to make sure that the eigenvalues match expectations 00041 * i.e. |WIN - WNEW(1+IPREPAD)| / (eps * |WIN|) < THRESH 00042 * where WIN is the array of eigenvalues as computed by 00043 * PCHEGVX when eigenvectors are requested 00044 * 00045 * Arguments 00046 * ========= 00047 * 00048 * NP = the number of rows local to a given process. 00049 * NQ = the number of columns local to a given process. 00050 * 00051 * WKNOWN (global input) INTEGER 00052 * .FALSE.: WIN does not contain the eigenvalues 00053 * .TRUE.: WIN does contain the eigenvalues 00054 * 00055 * IBTYPE (global input) INTEGER 00056 * Specifies the problem type to be solved: 00057 * = 1: sub( A )*x = (lambda)*sub( B )*x 00058 * = 2: sub( A )*sub( B )*x = (lambda)*x 00059 * = 3: sub( B )*sub( A )*x = (lambda)*x 00060 * 00061 * 00062 * JOBZ (global input) CHARACTER*1 00063 * Specifies whether or not to compute the eigenvectors: 00064 * = 'N': Compute eigenvalues only. 00065 * = 'V': Compute eigenvalues and eigenvectors. 00066 * Must be 'V' on first call to PCGSEPSUBTST 00067 * 00068 * RANGE (global input) CHARACTER*1 00069 * = 'A': all eigenvalues will be found. 00070 * = 'V': all eigenvalues in the interval [VL,VU] 00071 * will be found. 00072 * = 'I': the IL-th through IU-th eigenvalues will be found. 00073 * Must be 'A' on first call to PCGSEPSUBTST 00074 * 00075 * UPLO (global input) CHARACTER*1 00076 * Specifies whether the upper or lower triangular part of the 00077 * Hermitian matrix A is stored: 00078 * = 'U': Upper triangular 00079 * = 'L': Lower triangular 00080 * 00081 * N (global input) INTEGER 00082 * Size of the matrix to be tested. (global size) 00083 * 00084 * VL (global input) REAL 00085 * If RANGE='V', the lower bound of the interval to be searched 00086 * for eigenvalues. Not referenced if RANGE = 'A' or 'I'. 00087 * 00088 * VU (global input) REAL 00089 * If RANGE='V', the upper bound of the interval to be searched 00090 * for eigenvalues. Not referenced if RANGE = 'A' or 'I'. 00091 * 00092 * IL (global input) INTEGER 00093 * If RANGE='I', the index (from smallest to largest) of the 00094 * smallest eigenvalue to be returned. IL >= 1. 00095 * Not referenced if RANGE = 'A' or 'V'. 00096 * 00097 * IU (global input) INTEGER 00098 * If RANGE='I', the index (from smallest to largest) of the 00099 * largest eigenvalue to be returned. min(IL,N) <= IU <= N. 00100 * Not referenced if RANGE = 'A' or 'V'. 00101 * 00102 * THRESH (global input) REAL 00103 * A test will count as "failed" if the "error", computed as 00104 * described below, exceeds THRESH. Note that the error 00105 * is scaled to be O(1), so THRESH should be a reasonably 00106 * small multiple of 1, e.g., 10 or 100. In particular, 00107 * it should not depend on the precision (single vs. double) 00108 * or the size of the matrix. It must be at least zero. 00109 * 00110 * ABSTOL (global input) REAL 00111 * The absolute tolerance for the eigenvalues. An 00112 * eigenvalue is considered to be located if it has 00113 * been determined to lie in an interval whose width 00114 * is "abstol" or less. If "abstol" is less than or equal 00115 * to zero, then ulp*|T| will be used, where |T| is 00116 * the 1-norm of the matrix. If eigenvectors are 00117 * desired later by inverse iteration ("PCSTEIN"), 00118 * "abstol" MUST NOT be bigger than ulp*|T|. 00119 * 00120 * A (local workspace) COMPLEX array 00121 * global dimension (N, N), local dimension (DESCA(DLEN_), NQ) 00122 * A is distributed in a block cyclic manner over both rows 00123 * and columns. 00124 * See PCHEGVX for a description of block cyclic layout. 00125 * The test matrix, which is then modified by PCHEGVX 00126 * A has already been padded front and back, use A(1+IPREPAD) 00127 * 00128 * COPYA (local input) COMPLEX array, dimension(N*N) 00129 * COPYA holds a copy of the original matrix A 00130 * identical in both form and content to A 00131 * 00132 * B (local workspace) COMPLEX array, dim (N*N) 00133 * global dimension (N, N), local dimension (LDA, NQ) 00134 * A is distributed in a block cyclic manner over both rows 00135 * and columns. 00136 * The B test matrix, which is then modified by PCHEGVX 00137 * 00138 * COPYB (local input) COMPLEX array, dim (N, N) 00139 * COPYB is used to hold an identical copy of the array B 00140 * identical in both form and content to B 00141 * 00142 * Z (local workspace) COMPLEX array, dim (N*N) 00143 * Z is distributed in the same manner as A 00144 * Z contains the eigenvector matrix 00145 * Z is used as workspace by the test routines 00146 * PCGSEPCHK and PCSEPQTQ. 00147 * Z has already been padded front and back, use Z(1+IPREPAD) 00148 * 00149 * IA (global input) INTEGER 00150 * On entry, IA specifies the global row index of the submatrix 00151 * of the global matrix A, COPYA and Z to operate on. 00152 * 00153 * JA (global input) INTEGER 00154 * On entry, IA specifies the global column index of the submat 00155 * of the global matrix A, COPYA and Z to operate on. 00156 * 00157 * DESCA (global/local input) INTEGER array of dimension 8 00158 * The array descriptor for the matrix A, COPYA and Z. 00159 * 00160 * WIN (global input) REAL array, dimension (N) 00161 * If .not. WKNOWN, WIN is ignored on input 00162 * Otherwise, WIN() is taken as the standard by which the 00163 * eigenvalues are to be compared against. 00164 * 00165 * WNEW (global workspace) REAL array, dimension (N) 00166 * The eigenvalues as copmuted by this call to PCHEGVX 00167 * If JOBZ <> 'V' or RANGE <> 'A' these eigenvalues are 00168 * compared against those in WIN(). 00169 * WNEW has already been padded front and back, 00170 * use WNEW(1+IPREPAD) 00171 * 00172 * IFAIL (global output) INTEGER array, dimension (N) 00173 * If JOBZ = 'V', then on normal exit, the first M elements of 00174 * IFAIL are zero. If INFO > 0 on exit, then IFAIL contains the 00175 * indices of the eigenvectors that failed to converge. 00176 * If JOBZ = 'N', then IFAIL is not referenced. 00177 * IFAIL has already been padded front and back, 00178 * use IFAIL(1+IPREPAD) 00179 * 00180 * ICLUSTR (global workspace) integer array, dimension (2*NPROW*NPCOL) 00181 * 00182 * GAP (global workspace) REAL array, 00183 * dimension (NPROW*NPCOL) 00184 * 00185 * WORK (local workspace) COMPLEX array, dimension (LWORK) 00186 * WORK has already been padded front and back, 00187 * use WORK(1+IPREPAD) 00188 * 00189 * LWORK (local input) INTEGER 00190 * The actual length of the array WORK after padding. 00191 * 00192 * RWORK (local workspace) REAL array, dimension (LRWORK) 00193 * RWORK has already been padded front and back, 00194 * use RWORK(1+IPREPAD) 00195 * 00196 * LRWORK (local input) INTEGER 00197 * The actual length of the array RWORK after padding. 00198 * 00199 * LWORK1 (local input) INTEGER 00200 * The amount of real workspace to pass to PCHEGVX 00201 * 00202 * IWORK (local workspace) INTEGER array, dimension (LIWORK) 00203 * IWORK has already been padded front and back, 00204 * use IWORK(1+IPREPAD) 00205 * 00206 * LIWORK (local input) INTEGER 00207 * The length of the array IWORK after padding. 00208 * 00209 * RESULT (global output) INTEGER 00210 * The result of this call to PCHEGVX 00211 * RESULT = -3 => This process did not participate 00212 * RESULT = 0 => All tests passed 00213 * RESULT = 1 => ONe or more tests failed 00214 * 00215 * TSTNRM (global output) REAL 00216 * |AQ- QL| / |A|*N*EPS 00217 * 00218 * QTQNRM (global output) REAL 00219 * |QTQ -I| / N*EPS 00220 * 00221 * .. Parameters .. 00222 * 00223 INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_, 00224 $ MB_, NB_, RSRC_, CSRC_, LLD_ 00225 PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, 00226 $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 00227 $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) 00228 REAL PADVAL, FIVE, NEGONE 00229 PARAMETER ( PADVAL = 13.5285E+0, FIVE = 5.0E+0, 00230 $ NEGONE = -1.0E+0 ) 00231 COMPLEX CPADVAL 00232 PARAMETER ( CPADVAL = ( 13.989E+0, 1.93E+0 ) ) 00233 INTEGER IPADVAL 00234 PARAMETER ( IPADVAL = 927 ) 00235 * .. 00236 * .. Local Scalars .. 00237 LOGICAL MISSLARGEST, MISSSMALLEST 00238 INTEGER I, IAM, INDIWRK, INFO, ISIZEHEEVX, ISIZESUBTST, 00239 $ ISIZETST, J, M, MAXEIGS, MAXIL, MAXIU, MAXSIZE, 00240 $ MINIL, MQ, MYCOL, MYIL, MYROW, NCLUSTERS, NP, 00241 $ NPCOL, NPROW, NQ, NZ, OLDIL, OLDIU, OLDNZ, RES, 00242 $ RSIZECHK, RSIZEHEEVX, RSIZEQTQ, RSIZESUBTST, 00243 $ RSIZETST, SIZEHEEVX, SIZEMQRLEFT, SIZEMQRRIGHT, 00244 $ SIZEQRF, SIZESUBTST, SIZETMS, SIZETST, VALSIZE, 00245 $ VECSIZE 00246 REAL EPS, ERROR, MAXERROR, MAXVU, MINERROR, MINVL, 00247 $ NORMWIN, OLDVL, OLDVU, ORFAC, SAFMIN 00248 * .. 00249 * .. Local Arrays .. 00250 INTEGER DESCZ( DLEN_ ), DSEED( 4 ), ITMP( 2 ) 00251 * .. 00252 * .. External Functions .. 00253 * 00254 LOGICAL LSAME 00255 INTEGER NUMROC 00256 REAL PSLAMCH 00257 EXTERNAL LSAME, NUMROC, PSLAMCH 00258 * .. 00259 * .. External Subroutines .. 00260 EXTERNAL BLACS_GRIDINFO, CLACPY, DESCINIT, IGAMN2D, 00261 $ IGAMX2D, PCCHEKPAD, PCELSET, PCFILLPAD, 00262 $ PCGSEPCHK, PCHEGVX, PCLASIZEGSEP, 00263 $ PCLASIZEHEEVX, PICHEKPAD, PIFILLPAD, PSCHEKPAD, 00264 $ PSFILLPAD, SGAMN2D, SGAMX2D, SLBOOT, SLTIMER 00265 * .. 00266 * .. Intrinsic Functions .. 00267 INTRINSIC ABS, MAX, MIN, MOD 00268 * .. 00269 * .. Executable Statements .. 00270 * This is just to keep ftnchek happy 00271 IF( BLOCK_CYCLIC_2D*CSRC_*CTXT_*DLEN_*DTYPE_*LLD_*MB_*M_*NB_*N_* 00272 $ RSRC_.LT.0 )RETURN 00273 CALL PCLASIZEGSEP( DESCA, IPREPAD, IPOSTPAD, SIZEMQRLEFT, 00274 $ SIZEMQRRIGHT, SIZEQRF, SIZETMS, RSIZEQTQ, 00275 $ RSIZECHK, SIZEHEEVX, RSIZEHEEVX, ISIZEHEEVX, 00276 $ SIZESUBTST, RSIZESUBTST, ISIZESUBTST, SIZETST, 00277 $ RSIZETST, ISIZETST ) 00278 * 00279 TSTNRM = NEGONE 00280 QTQNRM = NEGONE 00281 EPS = PSLAMCH( DESCA( CTXT_ ), 'Eps' ) 00282 SAFMIN = PSLAMCH( DESCA( CTXT_ ), 'Safe min' ) 00283 * 00284 NORMWIN = SAFMIN / EPS 00285 IF( N.GE.1 ) 00286 $ NORMWIN = MAX( ABS( WIN( 1 ) ), ABS( WIN( N ) ), NORMWIN ) 00287 * 00288 * Make sure that we aren't using information from previous calls 00289 * 00290 NZ = -13 00291 OLDNZ = NZ 00292 OLDIL = IL 00293 OLDIU = IU 00294 OLDVL = VL 00295 OLDVU = VU 00296 * 00297 DO 10 I = 1, LWORK1, 1 00298 RWORK( I+IPREPAD ) = 14.3E+0 00299 10 CONTINUE 00300 DO 20 I = 1, LIWORK, 1 00301 IWORK( I+IPREPAD ) = 14 00302 20 CONTINUE 00303 DO 30 I = 1, LWORK, 1 00304 WORK( I+IPREPAD ) = ( 15.63E+0, 1.1E+0 ) 00305 30 CONTINUE 00306 * 00307 DO 40 I = 1, N 00308 WNEW( I+IPREPAD ) = 3.14159E+0 00309 40 CONTINUE 00310 * 00311 ICLUSTR( 1+IPREPAD ) = 139 00312 * 00313 IF( LSAME( JOBZ, 'N' ) ) THEN 00314 MAXEIGS = 0 00315 ELSE 00316 IF( LSAME( RANGE, 'A' ) ) THEN 00317 MAXEIGS = N 00318 ELSE IF( LSAME( RANGE, 'I' ) ) THEN 00319 MAXEIGS = IU - IL + 1 00320 ELSE 00321 MINVL = VL - NORMWIN*FIVE*EPS - ABSTOL 00322 MAXVU = VU + NORMWIN*FIVE*EPS + ABSTOL 00323 MINIL = 1 00324 MAXIU = 0 00325 DO 50 I = 1, N 00326 IF( WIN( I ).LT.MINVL ) 00327 $ MINIL = MINIL + 1 00328 IF( WIN( I ).LE.MAXVU ) 00329 $ MAXIU = MAXIU + 1 00330 50 CONTINUE 00331 * 00332 MAXEIGS = MAXIU - MINIL + 1 00333 END IF 00334 END IF 00335 * 00336 * 00337 CALL DESCINIT( DESCZ, DESCA( M_ ), DESCA( N_ ), DESCA( MB_ ), 00338 $ DESCA( NB_ ), DESCA( RSRC_ ), DESCA( CSRC_ ), 00339 $ DESCA( CTXT_ ), DESCA( LLD_ ), INFO ) 00340 * 00341 CALL BLACS_GRIDINFO( DESCA( CTXT_ ), NPROW, NPCOL, MYROW, MYCOL ) 00342 INDIWRK = 1 + IPREPAD + NPROW*NPCOL + 1 00343 * 00344 IAM = 1 00345 IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) 00346 $ IAM = 0 00347 * 00348 * If this process is not involved in this test, bail out now 00349 * 00350 RESULT = -3 00351 IF( MYROW.GE.NPROW .OR. MYROW.LT.0 ) 00352 $ GO TO 160 00353 RESULT = 0 00354 * 00355 * 00356 * DSEED is not used in this call to PCLASIZEHEEVX, the 00357 * following line just makes ftnchek happy. 00358 * 00359 DSEED( 1 ) = 1 00360 * 00361 CALL PCLASIZEHEEVX( WKNOWN, RANGE, N, DESCA, VL, VU, IL, IU, 00362 $ DSEED, WIN, MAXSIZE, VECSIZE, VALSIZE ) 00363 * 00364 NP = NUMROC( N, DESCA( MB_ ), MYROW, 0, NPROW ) 00365 NQ = NUMROC( N, DESCA( NB_ ), MYCOL, 0, NPCOL ) 00366 MQ = NUMROC( MAXEIGS, DESCA( NB_ ), MYCOL, 0, NPCOL ) 00367 * 00368 CALL CLACPY( 'A', NP, NQ, COPYA, DESCA( LLD_ ), A( 1+IPREPAD ), 00369 $ DESCA( LLD_ ) ) 00370 * 00371 CALL CLACPY( 'A', NP, NQ, COPYB, DESCA( LLD_ ), B( 1+IPREPAD ), 00372 $ DESCA( LLD_ ) ) 00373 * 00374 CALL PCFILLPAD( DESCA( CTXT_ ), NP, NQ, B, DESCA( LLD_ ), IPREPAD, 00375 $ IPOSTPAD, CPADVAL+1.0E+2 ) 00376 * 00377 CALL PCFILLPAD( DESCA( CTXT_ ), NP, NQ, A, DESCA( LLD_ ), IPREPAD, 00378 $ IPOSTPAD, CPADVAL ) 00379 * 00380 CALL PCFILLPAD( DESCZ( CTXT_ ), NP, MQ, Z, DESCZ( LLD_ ), IPREPAD, 00381 $ IPOSTPAD, CPADVAL+1.0E+0 ) 00382 * 00383 CALL PSFILLPAD( DESCA( CTXT_ ), N, 1, WNEW, N, IPREPAD, IPOSTPAD, 00384 $ PADVAL+2.0E+0 ) 00385 * 00386 CALL PSFILLPAD( DESCA( CTXT_ ), NPROW*NPCOL, 1, GAP, NPROW*NPCOL, 00387 $ IPREPAD, IPOSTPAD, PADVAL+3.0E+0 ) 00388 * 00389 CALL PSFILLPAD( DESCA( CTXT_ ), LWORK1, 1, RWORK, LWORK1, IPREPAD, 00390 $ IPOSTPAD, PADVAL+4.0E+0 ) 00391 * 00392 CALL PIFILLPAD( DESCA( CTXT_ ), LIWORK, 1, IWORK, LIWORK, IPREPAD, 00393 $ IPOSTPAD, IPADVAL ) 00394 * 00395 CALL PIFILLPAD( DESCA( CTXT_ ), N, 1, IFAIL, N, IPREPAD, IPOSTPAD, 00396 $ IPADVAL ) 00397 * 00398 CALL PIFILLPAD( DESCA( CTXT_ ), 2*NPROW*NPCOL, 1, ICLUSTR, 00399 $ 2*NPROW*NPCOL, IPREPAD, IPOSTPAD, IPADVAL ) 00400 * 00401 CALL PCFILLPAD( DESCA( CTXT_ ), LWORK, 1, WORK, LWORK, IPREPAD, 00402 $ IPOSTPAD, CPADVAL+4.1E+0 ) 00403 * 00404 * Make sure that PCHEGVX does not cheat (i.e. use answers 00405 * already computed.) 00406 * 00407 DO 70 I = 1, N, 1 00408 DO 60 J = 1, MAXEIGS, 1 00409 CALL PCELSET( Z( 1+IPREPAD ), I, J, DESCA, 00410 $ ( 13.0E+0, 1.34E+0 ) ) 00411 60 CONTINUE 00412 70 CONTINUE 00413 * 00414 ORFAC = -1.0E+0 00415 * 00416 CALL SLBOOT 00417 CALL SLTIMER( 1 ) 00418 CALL SLTIMER( 6 ) 00419 CALL PCHEGVX( IBTYPE, JOBZ, RANGE, UPLO, N, A( 1+IPREPAD ), IA, 00420 $ JA, DESCA, B( 1+IPREPAD ), IA, JA, DESCA, VL, VU, 00421 $ IL, IU, ABSTOL, M, NZ, WNEW( 1+IPREPAD ), ORFAC, 00422 $ Z( 1+IPREPAD ), IA, JA, DESCA, WORK( 1+IPREPAD ), 00423 $ SIZEHEEVX, RWORK( 1+IPREPAD ), LWORK1, 00424 $ IWORK( 1+IPREPAD ), LIWORK, IFAIL( 1+IPREPAD ), 00425 $ ICLUSTR( 1+IPREPAD ), GAP( 1+IPREPAD ), INFO ) 00426 CALL SLTIMER( 6 ) 00427 CALL SLTIMER( 1 ) 00428 * 00429 IF( THRESH.LE.0 ) THEN 00430 RESULT = 0 00431 ELSE 00432 CALL PCCHEKPAD( DESCA( CTXT_ ), 'PCHEGVX-B', NP, NQ, B, 00433 $ DESCA( LLD_ ), IPREPAD, IPOSTPAD, 00434 $ CPADVAL+1.0E+2 ) 00435 * 00436 CALL PCCHEKPAD( DESCA( CTXT_ ), 'PCHEGVX-A', NP, NQ, A, 00437 $ DESCA( LLD_ ), IPREPAD, IPOSTPAD, CPADVAL ) 00438 * 00439 CALL PCCHEKPAD( DESCZ( CTXT_ ), 'PCHEGVX-Z', NP, MQ, Z, 00440 $ DESCZ( LLD_ ), IPREPAD, IPOSTPAD, 00441 $ CPADVAL+1.0E+0 ) 00442 * 00443 CALL PSCHEKPAD( DESCA( CTXT_ ), 'PCHEGVX-WNEW', N, 1, WNEW, N, 00444 $ IPREPAD, IPOSTPAD, PADVAL+2.0E+0 ) 00445 * 00446 CALL PSCHEKPAD( DESCA( CTXT_ ), 'PCHEGVX-GAP', NPROW*NPCOL, 1, 00447 $ GAP, NPROW*NPCOL, IPREPAD, IPOSTPAD, 00448 $ PADVAL+3.0E+0 ) 00449 * 00450 CALL PSCHEKPAD( DESCA( CTXT_ ), 'PCHEGVX-rWORK', LWORK1, 1, 00451 $ RWORK, LWORK1, IPREPAD, IPOSTPAD, 00452 $ PADVAL+4.0E+0 ) 00453 * 00454 CALL PCCHEKPAD( DESCA( CTXT_ ), 'PCHEGVX-WORK', LWORK, 1, WORK, 00455 $ LWORK, IPREPAD, IPOSTPAD, CPADVAL+4.1E+0 ) 00456 * 00457 CALL PICHEKPAD( DESCA( CTXT_ ), 'PCHEGVX-IWORK', LIWORK, 1, 00458 $ IWORK, LIWORK, IPREPAD, IPOSTPAD, IPADVAL ) 00459 * 00460 CALL PICHEKPAD( DESCA( CTXT_ ), 'PCHEGVX-IFAIL', N, 1, IFAIL, 00461 $ N, IPREPAD, IPOSTPAD, IPADVAL ) 00462 * 00463 CALL PICHEKPAD( DESCA( CTXT_ ), 'PCHEGVX-ICLUSTR', 00464 $ 2*NPROW*NPCOL, 1, ICLUSTR, 2*NPROW*NPCOL, 00465 $ IPREPAD, IPOSTPAD, IPADVAL ) 00466 * 00467 * 00468 * Since we now know the spectrum, we can potentially reduce MAXSIZE. 00469 * 00470 IF( LSAME( RANGE, 'A' ) ) THEN 00471 CALL PCLASIZEHEEVX( .TRUE., RANGE, N, DESCA, VL, VU, IL, IU, 00472 $ DSEED, WNEW( 1+IPREPAD ), MAXSIZE, 00473 $ VECSIZE, VALSIZE ) 00474 END IF 00475 * 00476 * 00477 * Check INFO 00478 * 00479 * 00480 * Make sure that all processes return the same value of INFO 00481 * 00482 ITMP( 1 ) = INFO 00483 ITMP( 2 ) = INFO 00484 * 00485 CALL IGAMN2D( DESCA( CTXT_ ), 'a', ' ', 1, 1, ITMP, 1, 1, 1, 00486 $ -1, -1, 0 ) 00487 CALL IGAMX2D( DESCA( CTXT_ ), 'a', ' ', 1, 1, ITMP( 2 ), 1, 1, 00488 $ 1, -1, -1, 0 ) 00489 * 00490 * 00491 IF( ITMP( 1 ).NE.ITMP( 2 ) ) THEN 00492 IF( IAM.EQ.0 ) 00493 $ WRITE( NOUT, FMT = * ) 00494 $ 'Different processes return different INFO' 00495 RESULT = 1 00496 ELSE IF( MOD( INFO, 2 ).EQ.1 .OR. INFO.GT.7 .OR. INFO.LT.0 ) 00497 $ THEN 00498 IF( IAM.EQ.0 ) 00499 $ WRITE( NOUT, FMT = 9999 )INFO 00500 RESULT = 1 00501 ELSE IF( MOD( INFO / 2, 2 ).EQ.1 .AND. LWORK1.GE.MAXSIZE ) THEN 00502 IF( IAM.EQ.0 ) 00503 $ WRITE( NOUT, FMT = 9996 )INFO 00504 RESULT = 1 00505 ELSE IF( MOD( INFO / 4, 2 ).EQ.1 .AND. LWORK1.GE.VECSIZE ) THEN 00506 IF( IAM.EQ.0 ) 00507 $ WRITE( NOUT, FMT = 9996 )INFO 00508 RESULT = 1 00509 END IF 00510 * 00511 * 00512 IF( LSAME( JOBZ, 'V' ) .AND. ( ICLUSTR( 1+IPREPAD ).NE. 00513 $ 0 ) .AND. ( MOD( INFO / 2, 2 ).NE.1 ) ) THEN 00514 IF( IAM.EQ.0 ) 00515 $ WRITE( NOUT, FMT = 9995 ) 00516 RESULT = 1 00517 END IF 00518 * 00519 * Check M 00520 * 00521 IF( ( M.LT.0 ) .OR. ( M.GT.N ) ) THEN 00522 IF( IAM.EQ.0 ) 00523 $ WRITE( NOUT, FMT = 9994 ) 00524 RESULT = 1 00525 ELSE IF( LSAME( RANGE, 'A' ) .AND. ( M.NE.N ) ) THEN 00526 IF( IAM.EQ.0 ) 00527 $ WRITE( NOUT, FMT = 9993 ) 00528 RESULT = 1 00529 ELSE IF( LSAME( RANGE, 'I' ) .AND. ( M.NE.IU-IL+1 ) ) THEN 00530 IF( IAM.EQ.0 ) 00531 $ WRITE( NOUT, FMT = 9992 ) 00532 RESULT = 1 00533 ELSE IF( LSAME( JOBZ, 'V' ) .AND. 00534 $ ( .NOT.( LSAME( RANGE, 'V' ) ) ) .AND. ( M.NE.NZ ) ) 00535 $ THEN 00536 IF( IAM.EQ.0 ) 00537 $ WRITE( NOUT, FMT = 9991 ) 00538 RESULT = 1 00539 END IF 00540 * 00541 * Check NZ 00542 * 00543 IF( LSAME( JOBZ, 'V' ) ) THEN 00544 IF( LSAME( RANGE, 'V' ) ) THEN 00545 IF( NZ.GT.M ) THEN 00546 IF( IAM.EQ.0 ) 00547 $ WRITE( NOUT, FMT = 9990 ) 00548 RESULT = 1 00549 END IF 00550 IF( NZ.LT.M .AND. MOD( INFO / 4, 2 ).NE.1 ) THEN 00551 IF( IAM.EQ.0 ) 00552 $ WRITE( NOUT, FMT = 9989 ) 00553 RESULT = 1 00554 END IF 00555 ELSE 00556 IF( NZ.NE.M ) THEN 00557 IF( IAM.EQ.0 ) 00558 $ WRITE( NOUT, FMT = 9988 ) 00559 RESULT = 1 00560 END IF 00561 END IF 00562 END IF 00563 IF( RESULT.EQ.0 ) THEN 00564 * 00565 * Make sure that all processes return the same # of eigenvalues 00566 * 00567 ITMP( 1 ) = M 00568 ITMP( 2 ) = M 00569 * 00570 CALL IGAMN2D( DESCA( CTXT_ ), 'a', ' ', 1, 1, ITMP, 1, 1, 1, 00571 $ -1, -1, 0 ) 00572 CALL IGAMX2D( DESCA( CTXT_ ), 'a', ' ', 1, 1, ITMP( 2 ), 1, 00573 $ 1, 1, -1, -1, 0 ) 00574 * 00575 IF( ITMP( 1 ).NE.ITMP( 2 ) ) THEN 00576 IF( IAM.EQ.0 ) 00577 $ WRITE( NOUT, FMT = 9987 ) 00578 RESULT = 1 00579 ELSE 00580 * 00581 * Make sure that different processes return the same eigenvalues 00582 * 00583 DO 80 I = 1, M 00584 RWORK( I ) = WNEW( I+IPREPAD ) 00585 RWORK( I+M ) = WNEW( I+IPREPAD ) 00586 80 CONTINUE 00587 * 00588 CALL SGAMN2D( DESCA( CTXT_ ), 'a', ' ', M, 1, RWORK, M, 00589 $ 1, 1, -1, -1, 0 ) 00590 CALL SGAMX2D( DESCA( CTXT_ ), 'a', ' ', M, 1, 00591 $ RWORK( 1+M ), M, 1, 1, -1, -1, 0 ) 00592 * 00593 DO 90 I = 1, M 00594 * 00595 IF( RESULT.EQ.0 .AND. ( ABS( RWORK( I )-RWORK( M+ 00596 $ I ) ).GT.FIVE*EPS*ABS( RWORK( I ) ) ) ) THEN 00597 IF( IAM.EQ.0 ) 00598 $ WRITE( NOUT, FMT = 9986 ) 00599 RESULT = 1 00600 END IF 00601 90 CONTINUE 00602 END IF 00603 END IF 00604 * 00605 * Make sure that all processes return the same # of clusters 00606 * 00607 IF( LSAME( JOBZ, 'V' ) ) THEN 00608 NCLUSTERS = 0 00609 DO 100 I = 0, NPROW*NPCOL - 1 00610 IF( ICLUSTR( 1+IPREPAD+2*I ).EQ.0 ) 00611 $ GO TO 110 00612 NCLUSTERS = NCLUSTERS + 1 00613 100 CONTINUE 00614 110 CONTINUE 00615 ITMP( 1 ) = NCLUSTERS 00616 ITMP( 2 ) = NCLUSTERS 00617 * 00618 CALL IGAMN2D( DESCA( CTXT_ ), 'a', ' ', 1, 1, ITMP, 1, 1, 1, 00619 $ -1, -1, 0 ) 00620 CALL IGAMX2D( DESCA( CTXT_ ), 'a', ' ', 1, 1, ITMP( 2 ), 1, 00621 $ 1, 1, -1, -1, 0 ) 00622 * 00623 IF( ITMP( 1 ).NE.ITMP( 2 ) ) THEN 00624 IF( IAM.EQ.0 ) 00625 $ WRITE( NOUT, FMT = 9985 ) 00626 RESULT = 1 00627 ELSE 00628 * 00629 * Make sure that different processes return the same clusters 00630 * 00631 DO 120 I = 1, NCLUSTERS 00632 IWORK( INDIWRK+I ) = ICLUSTR( I+IPREPAD ) 00633 IWORK( INDIWRK+I+NCLUSTERS ) = ICLUSTR( I+IPREPAD ) 00634 120 CONTINUE 00635 CALL IGAMN2D( DESCA( CTXT_ ), 'a', ' ', NCLUSTERS*2+1, 1, 00636 $ IWORK( INDIWRK+1 ), NCLUSTERS*2+1, 1, 1, 00637 $ -1, -1, 0 ) 00638 CALL IGAMX2D( DESCA( CTXT_ ), 'a', ' ', NCLUSTERS*2+1, 1, 00639 $ IWORK( INDIWRK+1+NCLUSTERS ), 00640 $ NCLUSTERS*2+1, 1, 1, -1, -1, 0 ) 00641 * 00642 * 00643 DO 130 I = 1, NCLUSTERS 00644 IF( RESULT.EQ.0 .AND. IWORK( INDIWRK+I ).NE. 00645 $ IWORK( INDIWRK+NCLUSTERS+I ) ) THEN 00646 IF( IAM.EQ.0 ) 00647 $ WRITE( NOUT, FMT = 9984 ) 00648 RESULT = 1 00649 END IF 00650 130 CONTINUE 00651 * 00652 IF( ICLUSTR( 1+IPREPAD+NCLUSTERS*2 ).NE.0 ) THEN 00653 IF( IAM.EQ.0 ) 00654 $ WRITE( NOUT, FMT = 9983 ) 00655 RESULT = 1 00656 END IF 00657 END IF 00658 END IF 00659 * 00660 * 00661 CALL IGAMX2D( DESCA( CTXT_ ), 'a', ' ', 1, 1, RESULT, 1, 1, 1, 00662 $ -1, -1, 0 ) 00663 IF( RESULT.NE.0 ) 00664 $ GO TO 160 00665 * 00666 * Note that a couple key variables get redefined in PCGSEPCHK 00667 * as described by this table: 00668 * 00669 * PCGSEPTST name PCGSEPCHK name 00670 * ------------- ------------- 00671 * COPYA A 00672 * Z Q 00673 * B B 00674 * A C 00675 * 00676 * 00677 IF( LSAME( JOBZ, 'V' ) ) THEN 00678 * 00679 * Perform the residual check 00680 * 00681 CALL PSFILLPAD( DESCA( CTXT_ ), RSIZECHK, 1, RWORK, 00682 $ RSIZECHK, IPREPAD, IPOSTPAD, 4.3E+0 ) 00683 * 00684 CALL PCGSEPCHK( IBTYPE, N, NZ, COPYA, IA, JA, DESCA, COPYB, 00685 $ IA, JA, DESCA, THRESH, Z( 1+IPREPAD ), IA, 00686 $ JA, DESCZ, A( 1+IPREPAD ), IA, JA, DESCA, 00687 $ WNEW( 1+IPREPAD ), RWORK( 1+IPREPAD ), 00688 $ RSIZECHK, TSTNRM, RES ) 00689 * 00690 CALL PSCHEKPAD( DESCA( CTXT_ ), 'PCGSEPCHK-rWORK', RSIZECHK, 00691 $ 1, RWORK, RSIZECHK, IPREPAD, IPOSTPAD, 00692 $ 4.3E+0 ) 00693 * 00694 IF( RES.NE.0 ) 00695 $ RESULT = 1 00696 END IF 00697 * 00698 * Check to make sure that we have the right eigenvalues 00699 * 00700 IF( WKNOWN ) THEN 00701 * 00702 * Set up MYIL if necessary 00703 * 00704 MYIL = IL 00705 * 00706 IF( LSAME( RANGE, 'V' ) ) THEN 00707 MYIL = 1 00708 MINIL = 1 00709 MAXIL = N - M + 1 00710 ELSE 00711 IF( LSAME( RANGE, 'A' ) ) THEN 00712 MYIL = 1 00713 END IF 00714 MINIL = MYIL 00715 MAXIL = MYIL 00716 END IF 00717 * 00718 * Find the largest difference between the computed 00719 * and expected eigenvalues 00720 * 00721 MINERROR = NORMWIN 00722 * 00723 DO 150 MYIL = MINIL, MAXIL 00724 MAXERROR = 0 00725 * 00726 * Make sure that we aren't skipping any important eigenvalues 00727 * 00728 MISSSMALLEST = .TRUE. 00729 IF( .NOT.LSAME( RANGE, 'V' ) .OR. ( MYIL.EQ.1 ) ) 00730 $ MISSSMALLEST = .FALSE. 00731 IF( MISSSMALLEST .AND. ( WIN( MYIL-1 ).LT.VL+NORMWIN* 00732 $ FIVE*THRESH*EPS ) )MISSSMALLEST = .FALSE. 00733 MISSLARGEST = .TRUE. 00734 IF( .NOT.LSAME( RANGE, 'V' ) .OR. ( MYIL.EQ.MAXIL ) ) 00735 $ MISSLARGEST = .FALSE. 00736 IF( MISSLARGEST .AND. ( WIN( MYIL+M ).GT.VU-NORMWIN*FIVE* 00737 $ THRESH*EPS ) )MISSLARGEST = .FALSE. 00738 IF( .NOT.MISSSMALLEST ) THEN 00739 IF( .NOT.MISSLARGEST ) THEN 00740 * 00741 * Make sure that the eigenvalues that we report are OK 00742 * 00743 DO 140 I = 1, M 00744 ERROR = ABS( WIN( I+MYIL-1 )-WNEW( I+IPREPAD ) ) 00745 MAXERROR = MAX( MAXERROR, ERROR ) 00746 140 CONTINUE 00747 * 00748 MINERROR = MIN( MAXERROR, MINERROR ) 00749 END IF 00750 END IF 00751 150 CONTINUE 00752 * 00753 * 00754 * If JOBZ = 'V' and RANGE='A', we might be comparing 00755 * against our estimate of what the eigenvalues ought to 00756 * be, rather than comparing against what PxHEGVX computed 00757 * last time around, so we have to be more generous. 00758 * 00759 IF( LSAME( JOBZ, 'V' ) .AND. LSAME( RANGE, 'A' ) ) THEN 00760 IF( MINERROR.GT.NORMWIN*FIVE*FIVE*THRESH*EPS ) THEN 00761 IF( IAM.EQ.0 ) 00762 $ WRITE( NOUT, FMT = 9997 )MINERROR, NORMWIN 00763 RESULT = 1 00764 END IF 00765 ELSE 00766 IF( MINERROR.GT.NORMWIN*FIVE*THRESH*EPS ) THEN 00767 IF( IAM.EQ.0 ) 00768 $ WRITE( NOUT, FMT = 9997 )MINERROR, NORMWIN 00769 RESULT = 1 00770 END IF 00771 END IF 00772 END IF 00773 * 00774 * 00775 * Make sure that the IL, IU, VL and VU were not altered 00776 * 00777 IF( IL.NE.OLDIL .OR. IU.NE.OLDIU .OR. VL.NE.OLDVL .OR. VU.NE. 00778 $ OLDVU ) THEN 00779 IF( IAM.EQ.0 ) 00780 $ WRITE( NOUT, FMT = 9982 ) 00781 RESULT = 1 00782 END IF 00783 * 00784 IF( LSAME( JOBZ, 'N' ) .AND. ( NZ.NE.OLDNZ ) ) THEN 00785 IF( IAM.EQ.0 ) 00786 $ WRITE( NOUT, FMT = 9981 ) 00787 RESULT = 1 00788 END IF 00789 * 00790 END IF 00791 * 00792 * All processes should report the same result 00793 * 00794 CALL IGAMX2D( DESCA( CTXT_ ), 'a', ' ', 1, 1, RESULT, 1, 1, 1, -1, 00795 $ -1, 0 ) 00796 * 00797 160 CONTINUE 00798 * 00799 * 00800 RETURN 00801 * 00802 9999 FORMAT( 'PCHEGVX returned INFO=', I7 ) 00803 9998 FORMAT( 'PCSEPQTQ returned INFO=', I7 ) 00804 9997 FORMAT( 'PCGSEPSUBTST minerror =', D11.2, ' normwin=', D11.2 ) 00805 9996 FORMAT( 'PCHEGVX returned INFO=', I7, 00806 $ ' despite adequate workspace' ) 00807 9995 FORMAT( 'ICLUSTR(1).NE.0 but mod(INFO/2,2).NE.1' ) 00808 9994 FORMAT( 'M not in the range 0 to N' ) 00809 9993 FORMAT( 'M not equal to N' ) 00810 9992 FORMAT( 'M not equal to IU-IL+1' ) 00811 9991 FORMAT( 'M not equal to NZ' ) 00812 9990 FORMAT( 'NZ > M' ) 00813 9989 FORMAT( 'NZ < M' ) 00814 9988 FORMAT( 'NZ not equal to M' ) 00815 9987 FORMAT( 'Different processes return different values for M' ) 00816 9986 FORMAT( 'Different processes return different eigenvalues' ) 00817 9985 FORMAT( 'Different processes return ', 00818 $ 'different numbers of clusters' ) 00819 9984 FORMAT( 'Different processes return different clusters' ) 00820 9983 FORMAT( 'ICLUSTR not zero terminated' ) 00821 9982 FORMAT( 'IL, IU, VL or VU altered by PCHEGVX' ) 00822 9981 FORMAT( 'NZ altered by PCHEGVX with JOBZ=N' ) 00823 * 00824 * End of PCGSEPSUBTST 00825 * 00826 END