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