|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
00001 PROGRAM PSLLTDRIVER 00002 * 00003 * -- ScaLAPACK testing driver (version 1.7) -- 00004 * University of Tennessee, Knoxville, Oak Ridge National Laboratory, 00005 * and University of California, Berkeley. 00006 * May 1, 1997 00007 * 00008 * Purpose 00009 * ======= 00010 * 00011 * PSLLTDRIVER is the main test program for the REAL 00012 * ScaLAPACK Cholesky routines. This test driver performs an 00013 * A = L*L**T or A = U**T*U factorization and solve, and optionally 00014 * performs condition estimation and iterative refinement. 00015 * 00016 * The program must be driven by a short data file. An annotated 00017 * example of a data file can be obtained by deleting the first 3 00018 * characters from the following 18 lines: 00019 * 'ScaLAPACK LLt factorization input file' 00020 * 'Intel iPSC/860 hypercube, gamma model.' 00021 * 'LLT.out' output file name (if any) 00022 * 6 device out 00023 * 'U' define Lower or Upper 00024 * 1 number of problems sizes 00025 * 31 100 200 values of N 00026 * 1 number of NB's 00027 * 2 10 24 values of NB 00028 * 1 number of NRHS's 00029 * 1 values of NRHS 00030 * 1 Number of NBRHS's 00031 * 1 values of NBRHS 00032 * 1 number of process grids (ordered pairs of P & Q) 00033 * 2 values of P 00034 * 2 values of Q 00035 * 1.0 threshold 00036 * T (T or F) Test Cond. Est. and Iter. Ref. Routines 00037 * 00038 * 00039 * Internal Parameters 00040 * =================== 00041 * 00042 * TOTMEM INTEGER, default = 2000000 00043 * TOTMEM is a machine-specific parameter indicating the 00044 * maximum amount of available memory in bytes. 00045 * The user should customize TOTMEM to his platform. Remember 00046 * to leave room in memory for the operating system, the BLACS 00047 * buffer, etc. For example, on a system with 8 MB of memory 00048 * per process (e.g., one processor on an Intel iPSC/860), the 00049 * parameters we use are TOTMEM=6200000 (leaving 1.8 MB for OS, 00050 * code, BLACS buffer, etc). However, for PVM, we usually set 00051 * TOTMEM = 2000000. Some experimenting with the maximum value 00052 * of TOTMEM may be required. 00053 * 00054 * INTGSZ INTEGER, default = 4 bytes. 00055 * REALSZ INTEGER, default = 4 bytes. 00056 * INTGSZ and REALSZ indicate the length in bytes on the 00057 * given platform for an integer and a single precision real. 00058 * MEM REAL array, dimension ( TOTMEM / REALSZ ) 00059 * 00060 * All arrays used by SCALAPACK routines are allocated from 00061 * this array and referenced by pointers. The integer IPA, 00062 * for example, is a pointer to the starting element of MEM for 00063 * the matrix A. 00064 * 00065 * ===================================================================== 00066 * 00067 * .. Parameters .. 00068 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_, 00069 $ LLD_, MB_, M_, NB_, N_, RSRC_ 00070 PARAMETER ( BLOCK_CYCLIC_2D = 1, DLEN_ = 9, DTYPE_ = 1, 00071 $ CTXT_ = 2, M_ = 3, N_ = 4, MB_ = 5, NB_ = 6, 00072 $ RSRC_ = 7, CSRC_ = 8, LLD_ = 9 ) 00073 INTEGER INTGSZ, MEMSIZ, NTESTS, REALSZ, TOTMEM 00074 REAL PADVAL, ZERO 00075 PARAMETER ( INTGSZ = 4, REALSZ = 4, TOTMEM = 2000000, 00076 $ MEMSIZ = TOTMEM / REALSZ, NTESTS = 20, 00077 $ PADVAL = -9923.0E+0, ZERO = 0.0E+0 ) 00078 * .. 00079 * .. Local Scalars .. 00080 LOGICAL CHECK, EST 00081 CHARACTER UPLO 00082 CHARACTER*6 PASSED 00083 CHARACTER*80 OUTFILE 00084 INTEGER HH, I, IAM, IASEED, IBSEED, ICTXT, IMIDPAD, 00085 $ INFO, IPA, IPA0, IPB, IPB0, IPBERR, IPFERR, 00086 $ IPREPAD, IPOSTPAD, IPW, IPW2, ITEMP, J, K, 00087 $ KFAIL, KK, KPASS, KSKIP, KTESTS, LCM, LCMQ, 00088 $ LIWORK, LWORK, LW2, MYCOL, MYRHS, MYROW, N, NB, 00089 $ NBRHS, NGRIDS, NMAT, NNB, NNBR, NNR, NOUT, NP, 00090 $ NPCOL, NPROCS, NPROW, NQ, NRHS, WORKSIZ 00091 REAL ANORM, ANORM1, FRESID, RCOND, SRESID, SRESID2, 00092 $ THRESH 00093 DOUBLE PRECISION NOPS, TMFLOPS 00094 * .. 00095 * .. Local Arrays .. 00096 INTEGER DESCA( DLEN_ ), DESCB( DLEN_ ), IERR( 1 ), 00097 $ NBRVAL( NTESTS ), NBVAL( NTESTS ), 00098 $ NRVAL( NTESTS ), NVAL( NTESTS ), 00099 $ PVAL( NTESTS ), QVAL( NTESTS ) 00100 REAL MEM( MEMSIZ ) 00101 DOUBLE PRECISION CTIME( 2 ), WTIME( 2 ) 00102 * .. 00103 * .. External Subroutines .. 00104 EXTERNAL BLACS_BARRIER, BLACS_EXIT, BLACS_GRIDEXIT, 00105 $ BLACS_GRIDINFO, BLACS_GRIDINIT, DESCINIT, 00106 $ IGSUM2D, BLACS_PINFO, PSCHEKPAD, PSFILLPAD, 00107 $ PSLAFCHK, PSLASCHK, PSLLTINFO, 00108 $ PSMATGEN, PSPOCON, PSPORFS, 00109 $ PSPOTRF, PSPOTRRV, PSPOTRS, SLBOOT, 00110 $ SLCOMBINE, SLTIMER 00111 * .. 00112 * .. External Functions .. 00113 LOGICAL LSAME 00114 INTEGER ICEIL, ILCM, NUMROC 00115 REAL PSLANSY 00116 EXTERNAL ICEIL, ILCM, LSAME, NUMROC, PSLANSY 00117 * .. 00118 * .. Intrinsic Functions .. 00119 INTRINSIC DBLE, MAX, MIN 00120 * .. 00121 * .. Data Statements .. 00122 DATA KFAIL, KPASS, KSKIP, KTESTS / 4*0 / 00123 * .. 00124 * .. Executable Statements .. 00125 * 00126 * Get starting information 00127 * 00128 CALL BLACS_PINFO( IAM, NPROCS ) 00129 IASEED = 100 00130 IBSEED = 200 00131 CALL PSLLTINFO( OUTFILE, NOUT, UPLO, NMAT, NVAL, NTESTS, NNB, 00132 $ NBVAL, NTESTS, NNR, NRVAL, NTESTS, NNBR, NBRVAL, 00133 $ NTESTS, NGRIDS, PVAL, NTESTS, QVAL, NTESTS, 00134 $ THRESH, EST, MEM, IAM, NPROCS ) 00135 CHECK = ( THRESH.GE.0.0E+0 ) 00136 * 00137 * Print headings 00138 * 00139 IF( IAM.EQ.0 ) THEN 00140 WRITE( NOUT, FMT = * ) 00141 WRITE( NOUT, FMT = 9995 ) 00142 WRITE( NOUT, FMT = 9994 ) 00143 WRITE( NOUT, FMT = * ) 00144 END IF 00145 * 00146 * Loop over different process grids 00147 * 00148 DO 50 I = 1, NGRIDS 00149 * 00150 NPROW = PVAL( I ) 00151 NPCOL = QVAL( I ) 00152 * 00153 * Make sure grid information is correct 00154 * 00155 IERR( 1 ) = 0 00156 IF( NPROW.LT.1 ) THEN 00157 IF( IAM.EQ.0 ) 00158 $ WRITE( NOUT, FMT = 9999 ) 'GRID', 'nprow', NPROW 00159 IERR( 1 ) = 1 00160 ELSE IF( NPCOL.LT.1 ) THEN 00161 IF( IAM.EQ.0 ) 00162 $ WRITE( NOUT, FMT = 9999 ) 'GRID', 'npcol', NPCOL 00163 IERR( 1 ) = 1 00164 ELSE IF( NPROW*NPCOL.GT.NPROCS ) THEN 00165 IF( IAM.EQ.0 ) 00166 $ WRITE( NOUT, FMT = 9998 ) NPROW*NPCOL, NPROCS 00167 IERR( 1 ) = 1 00168 END IF 00169 * 00170 IF( IERR( 1 ).GT.0 ) THEN 00171 IF( IAM.EQ.0 ) 00172 $ WRITE( NOUT, FMT = 9997 ) 'grid' 00173 KSKIP = KSKIP + 1 00174 GO TO 50 00175 END IF 00176 * 00177 * Define process grid 00178 * 00179 CALL BLACS_GET( -1, 0, ICTXT ) 00180 CALL BLACS_GRIDINIT( ICTXT, 'Row-major', NPROW, NPCOL ) 00181 CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL ) 00182 * 00183 * Go to bottom of process grid loop if this case doesn't use my 00184 * process 00185 * 00186 IF( MYROW.GE.NPROW .OR. MYCOL.GE.NPCOL ) 00187 $ GO TO 50 00188 * 00189 DO 40 J = 1, NMAT 00190 * 00191 N = NVAL( J ) 00192 * 00193 * Make sure matrix information is correct 00194 * 00195 IERR( 1 ) = 0 00196 IF( N.LT.1 ) THEN 00197 IF( IAM.EQ.0 ) 00198 $ WRITE( NOUT, FMT = 9999 ) 'MATRIX', 'N', N 00199 IERR( 1 ) = 1 00200 ELSE IF( N.LT.1 ) THEN 00201 IF( IAM.EQ.0 ) 00202 $ WRITE( NOUT, FMT = 9999 ) 'MATRIX', 'N', N 00203 IERR( 1 ) = 1 00204 END IF 00205 * 00206 * Check all processes for an error 00207 * 00208 CALL IGSUM2D( ICTXT, 'All', ' ', 1, 1, IERR, 1, -1, 0 ) 00209 * 00210 IF( IERR( 1 ).GT.0 ) THEN 00211 IF( IAM.EQ.0 ) 00212 $ WRITE( NOUT, FMT = 9997 ) 'matrix' 00213 KSKIP = KSKIP + 1 00214 GO TO 40 00215 END IF 00216 * 00217 DO 30 K = 1, NNB 00218 * 00219 NB = NBVAL( K ) 00220 * 00221 * Make sure nb is legal 00222 * 00223 IERR( 1 ) = 0 00224 IF( NB.LT.1 ) THEN 00225 IERR( 1 ) = 1 00226 IF( IAM.EQ.0 ) 00227 $ WRITE( NOUT, FMT = 9999 ) 'NB', 'NB', NB 00228 END IF 00229 * 00230 * Check all processes for an error 00231 * 00232 CALL IGSUM2D( ICTXT, 'All', ' ', 1, 1, IERR, 1, -1, 0 ) 00233 * 00234 IF( IERR( 1 ).GT.0 ) THEN 00235 IF( IAM.EQ.0 ) 00236 $ WRITE( NOUT, FMT = 9997 ) 'NB' 00237 KSKIP = KSKIP + 1 00238 GO TO 30 00239 END IF 00240 * 00241 * Padding constants 00242 * 00243 NP = NUMROC( N, NB, MYROW, 0, NPROW ) 00244 NQ = NUMROC( N, NB, MYCOL, 0, NPCOL ) 00245 IF( CHECK ) THEN 00246 IPREPAD = MAX( NB, NP ) 00247 IMIDPAD = NB 00248 IPOSTPAD = MAX( NB, NQ ) 00249 ELSE 00250 IPREPAD = 0 00251 IMIDPAD = 0 00252 IPOSTPAD = 0 00253 END IF 00254 * 00255 * Initialize the array descriptor for the matrix A 00256 * 00257 CALL DESCINIT( DESCA, N, N, NB, NB, 0, 0, ICTXT, 00258 $ MAX( 1, NP )+IMIDPAD, IERR( 1 ) ) 00259 * 00260 * Check all processes for an error 00261 * 00262 CALL IGSUM2D( ICTXT, 'All', ' ', 1, 1, IERR, 1, -1, 0 ) 00263 * 00264 IF( IERR( 1 ).LT.0 ) THEN 00265 IF( IAM.EQ.0 ) 00266 $ WRITE( NOUT, FMT = 9997 ) 'descriptor' 00267 KSKIP = KSKIP + 1 00268 GO TO 30 00269 END IF 00270 * 00271 * Assign pointers into MEM for SCALAPACK arrays, A is 00272 * allocated starting at position MEM( IPREPAD+1 ) 00273 * 00274 IPA = IPREPAD+1 00275 IF( EST ) THEN 00276 IPA0 = IPA + DESCA( LLD_ )*NQ + IPOSTPAD + IPREPAD 00277 IPW = IPA0 + DESCA( LLD_ )*NQ + IPOSTPAD + IPREPAD 00278 ELSE 00279 IPW = IPA + DESCA( LLD_ )*NQ + IPOSTPAD + IPREPAD 00280 END IF 00281 * 00282 * 00283 IF( CHECK ) THEN 00284 * 00285 * Calculate the amount of workspace required by 00286 * the checking routines PSLAFCHK, PSPOTRRV, and 00287 * PSLANSY 00288 * 00289 WORKSIZ = NP * DESCA( NB_ ) 00290 * 00291 WORKSIZ = MAX( WORKSIZ, DESCA( MB_ ) * DESCA( NB_ ) ) 00292 * 00293 LCM = ILCM( NPROW, NPCOL ) 00294 ITEMP = MAX( 2, 2 * NQ ) + NP 00295 IF( NPROW.NE.NPCOL ) THEN 00296 ITEMP = ITEMP + 00297 $ NB * ICEIL( ICEIL( NP, NB ), LCM / NPROW ) 00298 END IF 00299 WORKSIZ = MAX( WORKSIZ, ITEMP ) 00300 WORKSIZ = WORKSIZ + IPOSTPAD 00301 * 00302 ELSE 00303 * 00304 WORKSIZ = IPOSTPAD 00305 * 00306 END IF 00307 * 00308 * Check for adequate memory for problem size 00309 * 00310 IERR( 1 ) = 0 00311 IF( IPW+WORKSIZ.GT.MEMSIZ ) THEN 00312 IF( IAM.EQ.0 ) 00313 $ WRITE( NOUT, FMT = 9996 ) 'factorization', 00314 $ ( IPW+WORKSIZ )*REALSZ 00315 IERR( 1 ) = 1 00316 END IF 00317 * 00318 * Check all processes for an error 00319 * 00320 CALL IGSUM2D( ICTXT, 'All', ' ', 1, 1, IERR, 1, -1, 0 ) 00321 * 00322 IF( IERR( 1 ).GT.0 ) THEN 00323 IF( IAM.EQ.0 ) 00324 $ WRITE( NOUT, FMT = 9997 ) 'MEMORY' 00325 KSKIP = KSKIP + 1 00326 GO TO 30 00327 END IF 00328 * 00329 * Generate a symmetric positive definite matrix A 00330 * 00331 CALL PSMATGEN( ICTXT, 'Symm', 'Diag', DESCA( M_ ), 00332 $ DESCA( N_ ), DESCA( MB_ ), DESCA( NB_ ), 00333 $ MEM( IPA ), DESCA( LLD_ ), DESCA( RSRC_ ), 00334 $ DESCA( CSRC_ ), IASEED, 0, NP, 0, NQ, 00335 $ MYROW, MYCOL, NPROW, NPCOL ) 00336 * 00337 * Calculate inf-norm of A for residual error-checking 00338 * 00339 IF( CHECK ) THEN 00340 CALL PSFILLPAD( ICTXT, NP, NQ, MEM( IPA-IPREPAD ), 00341 $ DESCA( LLD_ ), IPREPAD, IPOSTPAD, 00342 $ PADVAL ) 00343 CALL PSFILLPAD( ICTXT, WORKSIZ-IPOSTPAD, 1, 00344 $ MEM( IPW-IPREPAD ), WORKSIZ-IPOSTPAD, 00345 $ IPREPAD, IPOSTPAD, PADVAL ) 00346 ANORM = PSLANSY( 'I', UPLO, N, MEM( IPA ), 1, 1, 00347 $ DESCA, MEM( IPW ) ) 00348 ANORM1 = PSLANSY( '1', UPLO, N, MEM( IPA ), 1, 1, 00349 $ DESCA, MEM( IPW ) ) 00350 CALL PSCHEKPAD( ICTXT, 'PSLANSY', NP, NQ, 00351 $ MEM( IPA-IPREPAD ), DESCA( LLD_ ), 00352 $ IPREPAD, IPOSTPAD, PADVAL ) 00353 CALL PSCHEKPAD( ICTXT, 'PSLANSY', 00354 $ WORKSIZ-IPOSTPAD, 1, 00355 $ MEM( IPW-IPREPAD ), WORKSIZ-IPOSTPAD, 00356 $ IPREPAD, IPOSTPAD, PADVAL ) 00357 END IF 00358 * 00359 IF( EST ) THEN 00360 CALL PSMATGEN( ICTXT, 'Symm', 'Diag', DESCA( M_ ), 00361 $ DESCA( N_ ), DESCA( MB_ ), 00362 $ DESCA( NB_ ), MEM( IPA0 ), 00363 $ DESCA( LLD_ ), DESCA( RSRC_ ), 00364 $ DESCA( CSRC_ ), IASEED, 0, NP, 0, NQ, 00365 $ MYROW, MYCOL, NPROW, NPCOL ) 00366 IF( CHECK ) 00367 $ CALL PSFILLPAD( ICTXT, NP, NQ, 00368 $ MEM( IPA0-IPREPAD ), DESCA( LLD_ ), 00369 $ IPREPAD, IPOSTPAD, PADVAL ) 00370 END IF 00371 * 00372 CALL SLBOOT() 00373 CALL BLACS_BARRIER( ICTXT, 'All' ) 00374 * 00375 * Perform LLt factorization 00376 * 00377 CALL SLTIMER( 1 ) 00378 * 00379 CALL PSPOTRF( UPLO, N, MEM( IPA ), 1, 1, DESCA, INFO ) 00380 * 00381 CALL SLTIMER( 1 ) 00382 * 00383 IF( INFO.NE.0 ) THEN 00384 IF( IAM.EQ.0 ) 00385 $ WRITE( NOUT, FMT = * ) 'PSPOTRF INFO=', INFO 00386 KFAIL = KFAIL + 1 00387 RCOND = ZERO 00388 GO TO 60 00389 END IF 00390 * 00391 IF( CHECK ) THEN 00392 * 00393 * Check for memory overwrite in LLt factorization 00394 * 00395 CALL PSCHEKPAD( ICTXT, 'PSPOTRF', NP, NQ, 00396 $ MEM( IPA-IPREPAD ), DESCA( LLD_ ), 00397 $ IPREPAD, IPOSTPAD, PADVAL ) 00398 END IF 00399 * 00400 IF( EST ) THEN 00401 * 00402 * Calculate workspace required for PSPOCON 00403 * 00404 LWORK = MAX( 1, 2*NP ) + MAX( 1, 2*NQ ) + 00405 $ MAX( 2, DESCA( NB_ )* 00406 $ MAX( 1, ICEIL( NPROW-1, NPCOL ) ), 00407 $ NQ + DESCA( NB_ )* 00408 $ MAX( 1, ICEIL( NPCOL-1, NPROW ) ) ) 00409 IPW2 = IPW + LWORK + IPOSTPAD + IPREPAD 00410 LIWORK = MAX( 1, NP ) 00411 LW2 = ICEIL( LIWORK*INTGSZ, REALSZ ) + IPOSTPAD 00412 * 00413 IERR( 1 ) = 0 00414 IF( IPW2+LW2.GT.MEMSIZ ) THEN 00415 IF( IAM.EQ.0 ) 00416 $ WRITE( NOUT, FMT = 9996 )'cond est', 00417 $ ( IPW2+LW2 )*REALSZ 00418 IERR( 1 ) = 1 00419 END IF 00420 * 00421 * Check all processes for an error 00422 * 00423 CALL IGSUM2D( ICTXT, 'All', ' ', 1, 1, IERR, 1, 00424 $ -1, 0 ) 00425 * 00426 IF( IERR( 1 ).GT.0 ) THEN 00427 IF( IAM.EQ.0 ) 00428 $ WRITE( NOUT, FMT = 9997 ) 'MEMORY' 00429 KSKIP = KSKIP + 1 00430 GO TO 60 00431 END IF 00432 * 00433 IF( CHECK ) THEN 00434 CALL PSFILLPAD( ICTXT, LWORK, 1, 00435 $ MEM( IPW-IPREPAD ), LWORK, 00436 $ IPREPAD, IPOSTPAD, PADVAL ) 00437 CALL PSFILLPAD( ICTXT, LW2-IPOSTPAD, 1, 00438 $ MEM( IPW2-IPREPAD ), 00439 $ LW2-IPOSTPAD, IPREPAD, 00440 $ IPOSTPAD, PADVAL ) 00441 END IF 00442 * 00443 * Compute condition number of the matrix 00444 * 00445 CALL PSPOCON( UPLO, N, MEM( IPA ), 1, 1, DESCA, 00446 $ ANORM1, RCOND, MEM( IPW ), LWORK, 00447 $ MEM( IPW2 ), LIWORK, INFO ) 00448 * 00449 IF( CHECK ) THEN 00450 CALL PSCHEKPAD( ICTXT, 'PSPOCON', NP, NQ, 00451 $ MEM( IPA-IPREPAD ), DESCA( LLD_ ), 00452 $ IPREPAD, IPOSTPAD, PADVAL ) 00453 CALL PSCHEKPAD( ICTXT, 'PSPOCON', 00454 $ LWORK, 1, MEM( IPW-IPREPAD ), 00455 $ LWORK, IPREPAD, IPOSTPAD, 00456 $ PADVAL ) 00457 CALL PSCHEKPAD( ICTXT, 'PSPOCON', 00458 $ LW2-IPOSTPAD, 1, 00459 $ MEM( IPW2-IPREPAD ), LW2-IPOSTPAD, 00460 $ IPREPAD, IPOSTPAD, PADVAL ) 00461 END IF 00462 END IF 00463 * 00464 * Loop over the different values for NRHS 00465 * 00466 DO 20 HH = 1, NNR 00467 * 00468 NRHS = NRVAL( HH ) 00469 * 00470 DO 10 KK = 1, NNBR 00471 * 00472 NBRHS = NBRVAL( KK ) 00473 * 00474 * Initialize Array Descriptor for rhs 00475 * 00476 CALL DESCINIT( DESCB, N, NRHS, NB, NBRHS, 0, 0, 00477 $ ICTXT, MAX( 1, NP )+IMIDPAD, 00478 $ IERR( 1 ) ) 00479 * 00480 * move IPW to allow room for RHS 00481 * 00482 MYRHS = NUMROC( DESCB( N_ ), DESCB( NB_ ), MYCOL, 00483 $ DESCB( CSRC_ ), NPCOL ) 00484 IPB = IPW 00485 * 00486 IF( EST ) THEN 00487 IPB0 = IPB + DESCB( LLD_ )*MYRHS + IPOSTPAD + 00488 $ IPREPAD 00489 IPFERR = IPB0 + DESCB( LLD_ )*MYRHS + IPOSTPAD 00490 $ + IPREPAD 00491 IPBERR = MYRHS + IPFERR + IPOSTPAD + IPREPAD 00492 IPW = MYRHS + IPBERR + IPOSTPAD + IPREPAD 00493 ELSE 00494 IPW = IPB + DESCB( LLD_ )*MYRHS + IPOSTPAD + 00495 $ IPREPAD 00496 END IF 00497 * 00498 IF( CHECK ) THEN 00499 * 00500 * Calculate the amount of workspace required by 00501 * the checking routines PSLASCHK 00502 * 00503 LCMQ = LCM / NPCOL 00504 WORKSIZ = MAX( WORKSIZ-IPOSTPAD, 00505 $ NQ * NBRHS + NP * NBRHS + 00506 $ MAX( MAX( NQ*NB, 2*NBRHS ), 00507 $ NBRHS * NUMROC( NUMROC(N,NB,0,0,NPCOL),NB, 00508 $ 0,0,LCMQ ) ) ) 00509 WORKSIZ = IPOSTPAD + WORKSIZ 00510 ELSE 00511 WORKSIZ = IPOSTPAD 00512 END IF 00513 * 00514 IERR( 1 ) = 0 00515 IF( IPW+WORKSIZ.GT.MEMSIZ ) THEN 00516 IF( IAM.EQ.0 ) 00517 $ WRITE( NOUT, FMT = 9996 )'solve', 00518 $ ( IPW+WORKSIZ )*REALSZ 00519 IERR( 1 ) = 1 00520 END IF 00521 * 00522 * Check all processes for an error 00523 * 00524 CALL IGSUM2D( ICTXT, 'All', ' ', 1, 1, IERR, 1, 00525 $ -1, 0 ) 00526 * 00527 IF( IERR( 1 ).GT.0 ) THEN 00528 IF( IAM.EQ.0 ) 00529 $ WRITE( NOUT, FMT = 9997 ) 'MEMORY' 00530 KSKIP = KSKIP + 1 00531 GO TO 10 00532 END IF 00533 * 00534 * Generate RHS 00535 * 00536 CALL PSMATGEN( ICTXT, 'No', 'No', DESCB( M_ ), 00537 $ DESCB( N_ ), DESCB( MB_ ), 00538 $ DESCB( NB_ ), MEM( IPB ), 00539 $ DESCB( LLD_ ), DESCB( RSRC_ ), 00540 $ DESCB( CSRC_ ), IBSEED, 0, NP, 0, 00541 $ MYRHS, MYROW, MYCOL, NPROW, NPCOL ) 00542 * 00543 IF( CHECK ) 00544 $ CALL PSFILLPAD( ICTXT, NP, MYRHS, 00545 $ MEM( IPB-IPREPAD ), 00546 $ DESCB( LLD_ ), 00547 $ IPREPAD, IPOSTPAD, PADVAL ) 00548 * 00549 IF( EST ) THEN 00550 CALL PSMATGEN( ICTXT, 'No', 'No', DESCB( M_ ), 00551 $ DESCB( N_ ), DESCB( MB_ ), 00552 $ DESCB( NB_ ), MEM( IPB0 ), 00553 $ DESCB( LLD_ ), DESCB( RSRC_ ), 00554 $ DESCB( CSRC_ ), IBSEED, 0, NP, 0, 00555 $ MYRHS, MYROW, MYCOL, NPROW, 00556 $ NPCOL ) 00557 * 00558 IF( CHECK ) THEN 00559 CALL PSFILLPAD( ICTXT, NP, MYRHS, 00560 $ MEM( IPB0-IPREPAD ), 00561 $ DESCB( LLD_ ), IPREPAD, 00562 $ IPOSTPAD, PADVAL ) 00563 CALL PSFILLPAD( ICTXT, 1, MYRHS, 00564 $ MEM( IPFERR-IPREPAD ), 1, 00565 $ IPREPAD, IPOSTPAD, 00566 $ PADVAL ) 00567 CALL PSFILLPAD( ICTXT, 1, MYRHS, 00568 $ MEM( IPBERR-IPREPAD ), 1, 00569 $ IPREPAD, IPOSTPAD, 00570 $ PADVAL ) 00571 END IF 00572 END IF 00573 * 00574 CALL BLACS_BARRIER( ICTXT, 'All' ) 00575 CALL SLTIMER( 2 ) 00576 * 00577 * Solve linear system via Cholesky factorization 00578 * 00579 CALL PSPOTRS( UPLO, N, NRHS, MEM( IPA ), 1, 1, 00580 $ DESCA, MEM( IPB ), 1, 1, DESCB, 00581 $ INFO ) 00582 * 00583 CALL SLTIMER( 2 ) 00584 * 00585 IF( CHECK ) THEN 00586 * 00587 * check for memory overwrite 00588 * 00589 CALL PSCHEKPAD( ICTXT, 'PSPOTRS', NP, NQ, 00590 $ MEM( IPA-IPREPAD ), 00591 $ DESCA( LLD_ ), 00592 $ IPREPAD, IPOSTPAD, PADVAL ) 00593 CALL PSCHEKPAD( ICTXT, 'PSPOTRS', NP, 00594 $ MYRHS, MEM( IPB-IPREPAD ), 00595 $ DESCB( LLD_ ), IPREPAD, 00596 $ IPOSTPAD, PADVAL ) 00597 * 00598 CALL PSFILLPAD( ICTXT, WORKSIZ-IPOSTPAD, 1, 00599 $ MEM( IPW-IPREPAD ), 00600 $ WORKSIZ-IPOSTPAD, IPREPAD, 00601 $ IPOSTPAD, PADVAL ) 00602 * 00603 * check the solution to rhs 00604 * 00605 CALL PSLASCHK( 'Symm', 'Diag', N, NRHS, 00606 $ MEM( IPB ), 1, 1, DESCB, 00607 $ IASEED, 1, 1, DESCA, IBSEED, 00608 $ ANORM, SRESID, MEM( IPW ) ) 00609 * 00610 IF( IAM.EQ.0 .AND. SRESID.GT.THRESH ) 00611 $ WRITE( NOUT, FMT = 9985 ) SRESID 00612 * 00613 * check for memory overwrite 00614 * 00615 CALL PSCHEKPAD( ICTXT, 'PSLASCHK', NP, 00616 $ MYRHS, MEM( IPB-IPREPAD ), 00617 $ DESCB( LLD_ ), IPREPAD, 00618 $ IPOSTPAD, PADVAL ) 00619 CALL PSCHEKPAD( ICTXT, 'PSLASCHK', 00620 $ WORKSIZ-IPOSTPAD, 1, 00621 $ MEM( IPW-IPREPAD ), 00622 $ WORKSIZ-IPOSTPAD, IPREPAD, 00623 $ IPOSTPAD, PADVAL ) 00624 * 00625 * The second test is a NaN trap 00626 * 00627 IF( ( SRESID.LE.THRESH ).AND. 00628 $ ( (SRESID-SRESID).EQ.0.0E+0 ) ) THEN 00629 KPASS = KPASS + 1 00630 PASSED = 'PASSED' 00631 ELSE 00632 KFAIL = KFAIL + 1 00633 PASSED = 'FAILED' 00634 END IF 00635 ELSE 00636 KPASS = KPASS + 1 00637 SRESID = SRESID - SRESID 00638 PASSED = 'BYPASS' 00639 END IF 00640 * 00641 IF( EST ) THEN 00642 * 00643 * Calculate workspace required for PSPORFS 00644 * 00645 LWORK = MAX( 1, 3*NP ) 00646 IPW2 = IPW + LWORK + IPOSTPAD + IPREPAD 00647 LIWORK = MAX( 1, NP ) 00648 LW2 = ICEIL( LIWORK*INTGSZ, REALSZ ) + 00649 $ IPOSTPAD 00650 * 00651 IERR( 1 ) = 0 00652 IF( IPW2+LW2.GT.MEMSIZ ) THEN 00653 IF( IAM.EQ.0 ) 00654 $ WRITE( NOUT, FMT = 9996 ) 00655 $ 'iter ref', ( IPW2+LW2 )*REALSZ 00656 IERR( 1 ) = 1 00657 END IF 00658 * 00659 * Check all processes for an error 00660 * 00661 CALL IGSUM2D( ICTXT, 'All', ' ', 1, 1, IERR, 00662 $ 1, -1, 0 ) 00663 * 00664 IF( IERR( 1 ).GT.0 ) THEN 00665 IF( IAM.EQ.0 ) 00666 $ WRITE( NOUT, FMT = 9997 ) 00667 $ 'MEMORY' 00668 KSKIP = KSKIP + 1 00669 GO TO 10 00670 END IF 00671 * 00672 IF( CHECK ) THEN 00673 CALL PSFILLPAD( ICTXT, LWORK, 1, 00674 $ MEM( IPW-IPREPAD ), 00675 $ LWORK, IPREPAD, IPOSTPAD, 00676 $ PADVAL ) 00677 CALL PSFILLPAD( ICTXT, LW2-IPOSTPAD, 00678 $ 1, MEM( IPW2-IPREPAD ), 00679 $ LW2-IPOSTPAD, 00680 $ IPREPAD, IPOSTPAD, 00681 $ PADVAL ) 00682 END IF 00683 * 00684 * Use iterative refinement to improve the 00685 * computed solution 00686 * 00687 CALL PSPORFS( UPLO, N, NRHS, MEM( IPA0 ), 00688 $ 1, 1, DESCA, MEM( IPA ), 1, 1, 00689 $ DESCA, MEM( IPB0 ), 1, 1, 00690 $ DESCB, MEM( IPB ), 1, 1, DESCB, 00691 $ MEM( IPFERR ), MEM( IPBERR ), 00692 $ MEM( IPW ), LWORK, MEM( IPW2 ), 00693 $ LIWORK, INFO ) 00694 * 00695 * check for memory overwrite 00696 * 00697 IF( CHECK ) THEN 00698 CALL PSCHEKPAD( ICTXT, 'PSPORFS', NP, 00699 $ NQ, MEM( IPA0-IPREPAD ), 00700 $ DESCA( LLD_ ), IPREPAD, 00701 $ IPOSTPAD, PADVAL ) 00702 CALL PSCHEKPAD( ICTXT, 'PSPORFS', NP, 00703 $ NQ, MEM( IPA-IPREPAD ), 00704 $ DESCA( LLD_ ), IPREPAD, 00705 $ IPOSTPAD, PADVAL ) 00706 CALL PSCHEKPAD( ICTXT, 'PSPORFS', NP, 00707 $ MYRHS, MEM( IPB-IPREPAD ), 00708 $ DESCB( LLD_ ), IPREPAD, 00709 $ IPOSTPAD, PADVAL ) 00710 CALL PSCHEKPAD( ICTXT, 'PSPORFS', NP, 00711 $ MYRHS, 00712 $ MEM( IPB0-IPREPAD ), 00713 $ DESCB( LLD_ ), IPREPAD, 00714 $ IPOSTPAD, PADVAL ) 00715 CALL PSCHEKPAD( ICTXT, 'PSPORFS', 1, 00716 $ MYRHS, 00717 $ MEM( IPFERR-IPREPAD ), 1, 00718 $ IPREPAD, IPOSTPAD, 00719 $ PADVAL ) 00720 CALL PSCHEKPAD( ICTXT, 'PSPORFS', 1, 00721 $ MYRHS, 00722 $ MEM( IPBERR-IPREPAD ), 1, 00723 $ IPREPAD, IPOSTPAD, 00724 $ PADVAL ) 00725 CALL PSCHEKPAD( ICTXT, 'PSPORFS', LWORK, 00726 $ 1, MEM( IPW-IPREPAD ), 00727 $ LWORK, IPREPAD, IPOSTPAD, 00728 $ PADVAL ) 00729 CALL PSCHEKPAD( ICTXT, 'PSPORFS', 00730 $ LW2-IPOSTPAD, 1, 00731 $ MEM( IPW2-IPREPAD ), 00732 $ LW2-IPOSTPAD, 00733 $ IPREPAD, IPOSTPAD, 00734 $ PADVAL ) 00735 * 00736 CALL PSFILLPAD( ICTXT, WORKSIZ-IPOSTPAD, 00737 $ 1, MEM( IPW-IPREPAD ), 00738 $ WORKSIZ-IPOSTPAD, IPREPAD, 00739 $ IPOSTPAD, PADVAL ) 00740 * 00741 * check the solution to rhs 00742 * 00743 CALL PSLASCHK( 'Symm', 'Diag', N, NRHS, 00744 $ MEM( IPB ), 1, 1, DESCB, 00745 $ IASEED, 1, 1, DESCA, 00746 $ IBSEED, ANORM, SRESID2, 00747 $ MEM( IPW ) ) 00748 * 00749 IF( IAM.EQ.0 .AND. SRESID2.GT.THRESH ) 00750 $ WRITE( NOUT, FMT = 9985 ) SRESID2 00751 * 00752 * check for memory overwrite 00753 * 00754 CALL PSCHEKPAD( ICTXT, 'PSLASCHK', NP, 00755 $ MYRHS, MEM( IPB-IPREPAD ), 00756 $ DESCB( LLD_ ), IPREPAD, 00757 $ IPOSTPAD, PADVAL ) 00758 CALL PSCHEKPAD( ICTXT, 'PSLASCHK', 00759 $ WORKSIZ-IPOSTPAD, 1, 00760 $ MEM( IPW-IPREPAD ), 00761 $ WORKSIZ-IPOSTPAD, 00762 $ IPREPAD, IPOSTPAD, 00763 $ PADVAL ) 00764 END IF 00765 END IF 00766 * 00767 * Gather maximum of all CPU and WALL clock timings 00768 * 00769 CALL SLCOMBINE( ICTXT, 'All', '>', 'W', 2, 1, 00770 $ WTIME ) 00771 CALL SLCOMBINE( ICTXT, 'All', '>', 'C', 2, 1, 00772 $ CTIME ) 00773 * 00774 * Print results 00775 * 00776 IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN 00777 * 00778 * 1/3 N^3 + 1/2 N^2 flops for LLt factorization 00779 * 00780 NOPS = (DBLE(N)**3)/3.0D+0 + 00781 $ (DBLE(N)**2)/2.0D+0 00782 * 00783 * nrhs * 2 N^2 flops for LLt solve. 00784 * 00785 NOPS = NOPS + 2.0D+0*(DBLE(N)**2)*DBLE(NRHS) 00786 * 00787 * Calculate total megaflops -- factorization and 00788 * solve -- for WALL and CPU time, and print output 00789 * 00790 * Print WALL time if machine supports it 00791 * 00792 IF( WTIME( 1 ) + WTIME( 2 ) .GT. 0.0D+0 ) THEN 00793 TMFLOPS = NOPS / 00794 $ ( ( WTIME( 1 )+WTIME( 2 ) ) * 1.0D+6 ) 00795 ELSE 00796 TMFLOPS = 0.0D+0 00797 END IF 00798 * 00799 IF( WTIME( 2 ).GE.0.0D+0 ) 00800 $ WRITE( NOUT, FMT = 9993 ) 'WALL', UPLO, N, 00801 $ NB, NRHS, NBRHS, NPROW, NPCOL, 00802 $ WTIME( 1 ), WTIME( 2 ), TMFLOPS, 00803 $ PASSED 00804 * 00805 * Print CPU time if machine supports it 00806 * 00807 IF( CTIME( 1 )+CTIME( 2 ).GT.0.0D+0 ) THEN 00808 TMFLOPS = NOPS / 00809 $ ( ( CTIME( 1 )+CTIME( 2 ) ) * 1.0D+6 ) 00810 ELSE 00811 TMFLOPS = 0.0D+0 00812 END IF 00813 * 00814 IF( CTIME( 2 ).GE.0.0D+0 ) 00815 $ WRITE( NOUT, FMT = 9993 ) 'CPU ', UPLO, N, 00816 $ NB, NRHS, NBRHS, NPROW, NPCOL, 00817 $ CTIME( 1 ), CTIME( 2 ), TMFLOPS, 00818 $ PASSED 00819 * 00820 END IF 00821 10 CONTINUE 00822 20 CONTINUE 00823 * 00824 IF( CHECK .AND. SRESID.GT.THRESH ) THEN 00825 * 00826 * Compute FRESID = ||A - LL'|| / (||A|| * N * eps) 00827 * 00828 CALL PSPOTRRV( UPLO, N, MEM( IPA ), 1, 1, DESCA, 00829 $ MEM( IPW ) ) 00830 CALL PSLAFCHK( 'Symm', 'Diag', N, N, MEM( IPA ), 1, 1, 00831 $ DESCA, IASEED, ANORM, FRESID, 00832 $ MEM( IPW ) ) 00833 * 00834 * Check for memory overwrite 00835 * 00836 CALL PSCHEKPAD( ICTXT, 'PSPOTRRV', NP, NQ, 00837 $ MEM( IPA-IPREPAD ), DESCA( LLD_ ), 00838 $ IPREPAD, IPOSTPAD, PADVAL ) 00839 CALL PSCHEKPAD( ICTXT, 'PSGETRRV', 00840 $ WORKSIZ-IPOSTPAD, 1, 00841 $ MEM( IPW-IPREPAD ), WORKSIZ-IPOSTPAD, 00842 $ IPREPAD, IPOSTPAD, PADVAL ) 00843 * 00844 IF( IAM.EQ.0 ) THEN 00845 IF( LSAME( UPLO, 'L' ) ) THEN 00846 WRITE( NOUT, FMT = 9986 ) 'L*L''', FRESID 00847 ELSE 00848 WRITE( NOUT, FMT = 9986 ) 'U''*U', FRESID 00849 END IF 00850 END IF 00851 END IF 00852 * 00853 30 CONTINUE 00854 40 CONTINUE 00855 CALL BLACS_GRIDEXIT( ICTXT ) 00856 50 CONTINUE 00857 * 00858 * Print ending messages and close output file 00859 * 00860 60 CONTINUE 00861 IF( IAM.EQ.0 ) THEN 00862 KTESTS = KPASS + KFAIL + KSKIP 00863 WRITE( NOUT, FMT = * ) 00864 WRITE( NOUT, FMT = 9992 ) KTESTS 00865 IF( CHECK ) THEN 00866 WRITE( NOUT, FMT = 9991 ) KPASS 00867 WRITE( NOUT, FMT = 9989 ) KFAIL 00868 ELSE 00869 WRITE( NOUT, FMT = 9990 ) KPASS 00870 END IF 00871 WRITE( NOUT, FMT = 9988 ) KSKIP 00872 WRITE( NOUT, FMT = * ) 00873 WRITE( NOUT, FMT = * ) 00874 WRITE( NOUT, FMT = 9987 ) 00875 IF( NOUT.NE.6 .AND. NOUT.NE.0 ) 00876 $ CLOSE ( NOUT ) 00877 END IF 00878 * 00879 CALL BLACS_EXIT( 0 ) 00880 * 00881 9999 FORMAT( 'ILLEGAL ', A6, ': ', A5, ' = ', I3, 00882 $ '; It should be at least 1' ) 00883 9998 FORMAT( 'ILLEGAL GRID: nprow*npcol = ', I4, '. It can be at most', 00884 $ I4 ) 00885 9997 FORMAT( 'Bad ', A6, ' parameters: going on to next test case.' ) 00886 9996 FORMAT( 'Unable to perform ', A, ': need TOTMEM of at least', 00887 $ I11 ) 00888 9995 FORMAT( 'TIME UPLO N NB NRHS NBRHS P Q LLt Time ', 00889 $ 'Slv Time MFLOPS CHECK' ) 00890 9994 FORMAT( '---- ---- ----- --- ---- ----- ---- ---- -------- ', 00891 $ '-------- -------- ------' ) 00892 9993 FORMAT( A4, 4X, A1, 1X, I5, 1X, I3, 1X, I4, 1X, I5, 1X, I4, 1X, 00893 $ I4, 1X, F8.2, 1X, F8.2, 1X, F8.2, 1X, A6 ) 00894 9992 FORMAT( 'Finished ', I6, ' tests, with the following results:' ) 00895 9991 FORMAT( I5, ' tests completed and passed residual checks.' ) 00896 9990 FORMAT( I5, ' tests completed without checking.' ) 00897 9989 FORMAT( I5, ' tests completed and failed residual checks.' ) 00898 9988 FORMAT( I5, ' tests skipped because of illegal input values.' ) 00899 9987 FORMAT( 'END OF TESTS.' ) 00900 9986 FORMAT( '||A - ', A4, '|| / (||A|| * N * eps) = ', G25.7 ) 00901 9985 FORMAT( '||Ax-b||/(||x||*||A||*eps*N) ', F25.7 ) 00902 * 00903 STOP 00904 * 00905 * End of PSLLTDRIVER 00906 * 00907 END
