OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
pdlltdriver.f
Go to the documentation of this file.
1 PROGRAM pdlltdriver
2*
3* -- ScaLAPACK testing driver (version 1.7) --
4* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
5* and University of California, Berkeley.
6* May 1, 1997
7*
8* Purpose
9* =======
10*
11* PDLLTDRIVER is the main test program for the DOUBLE PRECISION
12* ScaLAPACK Cholesky routines. This test driver performs an
13* A = L*L**T or A = U**T*U factorization and solve, and optionally
14* performs condition estimation and iterative refinement.
15*
16* The program must be driven by a short data file. An annotated
17* example of a data file can be obtained by deleting the first 3
18* characters from the following 18 lines:
19* 'ScaLAPACK LLt factorization input file'
20* 'Intel iPSC/860 hypercube, gamma model.'
21* 'LLT.out' output file name (if any)
22* 6 device out
23* 'U' define Lower or Upper
24* 1 number of problems sizes
25* 31 100 200 values of N
26* 1 number of NB's
27* 2 10 24 values of NB
28* 1 number of NRHS's
29* 1 values of NRHS
30* 1 Number of NBRHS's
31* 1 values of NBRHS
32* 1 number of process grids (ordered pairs of P & Q)
33* 2 values of P
34* 2 values of Q
35* 1.0 threshold
36* T (T or F) Test Cond. Est. and Iter. Ref. Routines
37*
38*
39* Internal Parameters
40* ===================
41*
42* TOTMEM INTEGER, default = 2000000
43* TOTMEM is a machine-specific parameter indicating the
44* maximum amount of available memory in bytes.
45* The user should customize TOTMEM to his platform. Remember
46* to leave room in memory for the operating system, the BLACS
47* buffer, etc. For example, on a system with 8 MB of memory
48* per process (e.g., one processor on an Intel iPSC/860), the
49* parameters we use are TOTMEM=6200000 (leaving 1.8 MB for OS,
50* code, BLACS buffer, etc). However, for PVM, we usually set
51* TOTMEM = 2000000. Some experimenting with the maximum value
52* of TOTMEM may be required.
53*
54* INTGSZ INTEGER, default = 4 bytes.
55* DBLESZ INTEGER, default = 8 bytes.
56* INTGSZ and DBLESZ indicate the length in bytes on the
57* given platform for an integer and a double precision real.
58* MEM DOUBLE PRECISION array, dimension ( TOTMEM / DBLESZ )
59*
60* All arrays used by SCALAPACK routines are allocated from
61* this array and referenced by pointers. The integer IPA,
62* for example, is a pointer to the starting element of MEM for
63* the matrix A.
64*
65* =====================================================================
66*
67* .. Parameters ..
68 INTEGER block_cyclic_2d, csrc_, ctxt_, dlen_, dtype_,
69 $ lld_, mb_, m_, nb_, n_, rsrc_
70 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
71 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
72 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
73 INTEGER dblesz, intgsz, memsiz, ntests, totmem
74 DOUBLE PRECISION padval, zero
75 PARAMETER ( dblesz = 8, intgsz = 4, totmem = 2000000,
76 $ memsiz = totmem / dblesz, ntests = 20,
77 $ padval = -9923.0d+0, zero = 0.0d+0 )
78* ..
79* .. Local Scalars ..
80 LOGICAL check, est
81 CHARACTER uplo
82 CHARACTER*6 passed
83 CHARACTER*80 outfile
84 INTEGER hh, i, iam, iaseed, ibseed, ictxt, imidpad,
85 $ info, ipa, ipa0, ipb, ipb0, ipberr, ipferr,
86 $ iprepad, ipostpad, ipw, ipw2, itemp, j, k,
87 $ kfail, KK, kpass, kskip, ktests, lcm, lcmq,
88 $ liwork, lwork, lw2, mycol, myrhs, myrow, n, nb,
89 $ nbrhs, ngrids, nmat, nnb, nnbr, nnr, nout, np,
90 $ npcol, nprocs, nprow, nq, nrhs, worksiz
91 REAL thresh
92 DOUBLE PRECISION anorm, anorm1, fresid, nops, rcond,
93 $ sresid, sresid2, tmflops
94* ..
95* .. Local Arrays ..
96 INTEGER desca( dlen_ ), descb( dlen_ ), ierr( 1 ),
97 $ nbrval( ntests ), nbval( NTESTS ),
98 $ nrval( NTESTS ), nval( ntests ),
99 $ pval( ntests ), qval( ntests )
100 DOUBLE PRECISION ctime( 2 ), mem( memsiz ), wtime( 2 )
101* ..
102* .. External Subroutines ..
103 EXTERNAL blacs_barrier, blacs_exit, blacs_gridexit,
104 $ blacs_gridinfo, blacs_gridinit, descinit,
105 $ igsum2d, blacs_pinfo, pdchekpad, pdfillpad,
106 $ pdlafchk, pdlaschk, pdlltinfo,
107 $ pdmatgen, pdpocon, pdporfs,
108 $ pdpotrf, pdpotrrv, pdpotrs, slboot,
109 $ slcombine, sltimer
110* ..
111* .. External Functions ..
112 LOGICAL lsame
113 INTEGER iceil, ilcm, numroc
114 DOUBLE PRECISION pdlansy
115 EXTERNAL iceil, ilcm, lsame, numroc, pdlansy
116* ..
117* .. Intrinsic Functions ..
118 INTRINSIC dble, max, min
119* ..
120* .. Data Statements ..
121 DATA kfail, kpass, kskip, ktests / 4*0 /
122* ..
123* .. Executable Statements ..
124*
125* Get starting information
126*
127 CALL blacs_pinfo( iam, nprocs )
128 iaseed = 100
129 ibseed = 200
130 CALL pdlltinfo( outfile, nout, uplo, nmat, nval, ntests, nnb,
131 $ nbval, ntests, nnr, nrval, ntests, nnbr, nbrval,
132 $ ntests, ngrids, pval, ntests, qval, ntests,
133 $ thresh, est, mem, iam, nprocs )
134 check = ( thresh.GE.0.0e+0 )
135*
136* Print headings
137*
138 IF( iam.EQ.0 ) THEN
139 WRITE( nout, fmt = * )
140 WRITE( nout, fmt = 9995 )
141 WRITE( nout, fmt = 9994 )
142 WRITE( nout, fmt = * )
143 END IF
144*
145* Loop over different process grids
146*
147 DO 50 i = 1, ngrids
148*
149 nprow = pval( i )
150 npcol = qval( i )
151*
152* Make sure grid information is correct
153*
154 ierr( 1 ) = 0
155 IF( nprow.LT.1 ) THEN
156 IF( iam.EQ.0 )
157 $ WRITE( nout, fmt = 9999 ) 'GRID', 'nprow', nprow
158 ierr( 1 ) = 1
159 ELSE IF( npcol.LT.1 ) THEN
160 IF( iam.EQ.0 )
161 $ WRITE( nout, fmt = 9999 ) 'GRID', 'npcol', npcol
162 ierr( 1 ) = 1
163 ELSE IF( nprow*npcol.GT.nprocs ) THEN
164 IF( iam.EQ.0 )
165 $ WRITE( nout, fmt = 9998 ) nprow*npcol, nprocs
166 ierr( 1 ) = 1
167 END IF
168*
169 IF( ierr( 1 ).GT.0 ) THEN
170 IF( iam.EQ.0 )
171 $ WRITE( nout, fmt = 9997 ) 'grid'
172 kskip = kskip + 1
173 GO TO 50
174 END IF
175*
176* Define process grid
177*
178 CALL blacs_get( -1, 0, ictxt )
179 CALL blacs_gridinit( ictxt, 'Row-major', nprow, npcol )
180 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
181*
182* Go to bottom of process grid loop if this case doesn't use my
183* process
184*
185 IF( myrow.GE.nprow .OR. mycol.GE.npcol )
186 $ GO TO 50
187*
188 DO 40 j = 1, nmat
189*
190 n = nval( j )
191*
192* Make sure matrix information is correct
193*
194 ierr( 1 ) = 0
195 IF( n.LT.1 ) THEN
196 IF( iam.EQ.0 )
197 $ WRITE( nout, fmt = 9999 ) 'MATRIX', 'N', n
198 ierr( 1 ) = 1
199 ELSE IF( n.LT.1 ) THEN
200 IF( iam.EQ.0 )
201 $ WRITE( nout, fmt = 9999 ) 'MATRIX', 'N', n
202 ierr( 1 ) = 1
203 END IF
204*
205* Check all processes for an error
206*
207 CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1, 0 )
208*
209 IF( ierr( 1 ).GT.0 ) THEN
210 IF( iam.EQ.0 )
211 $ WRITE( nout, fmt = 9997 ) 'matrix'
212 kskip = kskip + 1
213 GO TO 40
214 END IF
215*
216 DO 30 k = 1, nnb
217*
218 nb = nbval( k )
219*
220* Make sure nb is legal
221*
222 ierr( 1 ) = 0
223 IF( nb.LT.1 ) THEN
224 ierr( 1 ) = 1
225 IF( iam.EQ.0 )
226 $ WRITE( nout, fmt = 9999 ) 'nb', 'nb', NB
227 END IF
228*
229* Check all processes for an error
230*
231 CALL IGSUM2D( ICTXT, 'all', ' ', 1, 1, IERR, 1, -1, 0 )
232*
233.GT. IF( IERR( 1 )0 ) THEN
234.EQ. IF( IAM0 )
235 $ WRITE( NOUT, FMT = 9997 ) 'nb'
236 KSKIP = KSKIP + 1
237 GO TO 30
238 END IF
239*
240* Padding constants
241*
242 NP = NUMROC( N, NB, MYROW, 0, NPROW )
243 NQ = NUMROC( N, NB, MYCOL, 0, NPCOL )
244 IF( CHECK ) THEN
245 IPREPAD = MAX( NB, NP )
246 IMIDPAD = NB
247 IPOSTPAD = MAX( NB, NQ )
248 ELSE
249 IPREPAD = 0
250 IMIDPAD = 0
251 IPOSTPAD = 0
252 END IF
253*
254* Initialize the array descriptor for the matrix A
255*
256 CALL DESCINIT( DESCA, N, N, NB, NB, 0, 0, ICTXT,
257 $ MAX( 1, NP )+IMIDPAD, IERR( 1 ) )
258*
259* Check all processes for an error
260*
261 CALL IGSUM2D( ICTXT, 'all', ' ', 1, 1, IERR, 1, -1, 0 )
262*
263.LT. IF( IERR( 1 )0 ) THEN
264.EQ. IF( IAM0 )
265 $ WRITE( NOUT, FMT = 9997 ) 'descriptor'
266 KSKIP = KSKIP + 1
267 GO TO 30
268 END IF
269*
270* Assign pointers into MEM for SCALAPACK arrays, A is
271* allocated starting at position MEM( IPREPAD+1 )
272*
273 IPA = IPREPAD+1
274 IF( EST ) THEN
275 IPA0 = IPA + DESCA( LLD_ )*NQ + IPOSTPAD + IPREPAD
276 IPW = IPA0 + DESCA( LLD_ )*NQ + IPOSTPAD + IPREPAD
277 ELSE
278 IPW = IPA + DESCA( LLD_ )*NQ + IPOSTPAD + IPREPAD
279 END IF
280*
281*
282 IF( CHECK ) THEN
283*
284* Calculate the amount of workspace required by
285* the checking routines PDLAFCHK, PDPOTRRV, and
286* PDLANSY
287*
288 WORKSIZ = NP * DESCA( NB_ )
289*
290 WORKSIZ = MAX( WORKSIZ, DESCA( MB_ ) * DESCA( NB_ ) )
291*
292 LCM = ILCM( NPROW, NPCOL )
293 ITEMP = MAX( 2, 2 * NQ ) + NP
294.NE. IF( NPROWNPCOL ) THEN
295 ITEMP = ITEMP +
296 $ NB * ICEIL( ICEIL( NP, NB ), LCM / NPROW )
297 END IF
298 WORKSIZ = MAX( WORKSIZ, ITEMP )
299 WORKSIZ = WORKSIZ + IPOSTPAD
300*
301 ELSE
302*
303 WORKSIZ = IPOSTPAD
304*
305 END IF
306*
307* Check for adequate memory for problem size
308*
309 IERR( 1 ) = 0
310.GT. IF( IPW+WORKSIZMEMSIZ ) THEN
311.EQ. IF( IAM0 )
312 $ WRITE( NOUT, FMT = 9996 ) 'factorization',
313 $ ( IPW+WORKSIZ )*DBLESZ
314 IERR( 1 ) = 1
315 END IF
316*
317* Check all processes for an error
318*
319 CALL IGSUM2D( ICTXT, 'all', ' ', 1, 1, IERR, 1, -1, 0 )
320*
321.GT. IF( IERR( 1 )0 ) THEN
322.EQ. IF( IAM0 )
323 $ WRITE( NOUT, FMT = 9997 ) 'memory'
324 KSKIP = KSKIP + 1
325 GO TO 30
326 END IF
327*
328* Generate a symmetric positive definite matrix A
329*
330 CALL PDMATGEN( ICTXT, 'symm', 'diag', DESCA( M_ ),
331 $ DESCA( N_ ), DESCA( MB_ ), DESCA( NB_ ),
332 $ MEM( IPA ), DESCA( LLD_ ), DESCA( RSRC_ ),
333 $ DESCA( CSRC_ ), IASEED, 0, NP, 0, NQ,
334 $ MYROW, MYCOL, NPROW, NPCOL )
335*
336* Calculate inf-norm of A for residual error-checking
337*
338 IF( CHECK ) THEN
339 CALL PDFILLPAD( ICTXT, NP, NQ, MEM( IPA-IPREPAD ),
340 $ DESCA( LLD_ ), IPREPAD, IPOSTPAD,
341 $ PADVAL )
342 CALL PDFILLPAD( ICTXT, WORKSIZ-IPOSTPAD, 1,
343 $ MEM( IPW-IPREPAD ), WORKSIZ-IPOSTPAD,
344 $ IPREPAD, IPOSTPAD, PADVAL )
345 ANORM = PDLANSY( 'i', UPLO, N, MEM( IPA ), 1, 1,
346 $ DESCA, MEM( IPW ) )
347 ANORM1 = PDLANSY( '1', UPLO, N, MEM( IPA ), 1, 1,
348 $ DESCA, MEM( IPW ) )
349 CALL PDCHEKPAD( ICTXT, 'pdlansy', NP, NQ,
350 $ MEM( IPA-IPREPAD ), DESCA( LLD_ ),
351 $ IPREPAD, IPOSTPAD, PADVAL )
352 CALL PDCHEKPAD( ICTXT, 'pdlansy',
353 $ WORKSIZ-IPOSTPAD, 1,
354 $ MEM( IPW-IPREPAD ), WORKSIZ-IPOSTPAD,
355 $ IPREPAD, IPOSTPAD, PADVAL )
356 END IF
357*
358 IF( EST ) THEN
359 CALL PDMATGEN( ICTXT, 'symm', 'diag', DESCA( M_ ),
360 $ DESCA( N_ ), DESCA( MB_ ),
361 $ DESCA( NB_ ), MEM( IPA0 ),
362 $ DESCA( LLD_ ), DESCA( RSRC_ ),
363 $ DESCA( CSRC_ ), IASEED, 0, NP, 0, NQ,
364 $ MYROW, MYCOL, NPROW, NPCOL )
365 IF( CHECK )
366 $ CALL PDFILLPAD( ICTXT, NP, NQ,
367 $ MEM( IPA0-IPREPAD ), DESCA( LLD_ ),
368 $ IPREPAD, IPOSTPAD, PADVAL )
369 END IF
370*
371 CALL SLBOOT()
372 CALL BLACS_BARRIER( ICTXT, 'all' )
373*
374* Perform LLt factorization
375*
376 CALL SLTIMER( 1 )
377*
378 CALL PDPOTRF( UPLO, N, MEM( IPA ), 1, 1, DESCA, INFO )
379*
380 CALL SLTIMER( 1 )
381*
382.NE. IF( INFO0 ) THEN
383.EQ. IF( IAM0 )
384 $ WRITE( NOUT, FMT = * ) 'pdpotrf info=', INFO
385 KFAIL = KFAIL + 1
386 RCOND = ZERO
387 GO TO 60
388 END IF
389*
390 IF( CHECK ) THEN
391*
392* Check for memory overwrite in LLt factorization
393*
394 CALL PDCHEKPAD( ICTXT, 'pdpotrf', NP, NQ,
395 $ MEM( IPA-IPREPAD ), DESCA( LLD_ ),
396 $ IPREPAD, IPOSTPAD, PADVAL )
397 END IF
398*
399 IF( EST ) THEN
400*
401* Calculate workspace required for PDPOCON
402*
403 LWORK = MAX( 1, 2*NP ) + MAX( 1, 2*NQ ) +
404 $ MAX( 2, DESCA( NB_ )*
405 $ MAX( 1, ICEIL( NPROW-1, NPCOL ) ),
406 $ NQ + DESCA( NB_ )*
407 $ MAX( 1, ICEIL( NPCOL-1, NPROW ) ) )
408 IPW2 = IPW + LWORK + IPOSTPAD + IPREPAD
409 LIWORK = MAX( 1, NP )
410 LW2 = ICEIL( LIWORK*INTGSZ, DBLESZ ) + IPOSTPAD
411*
412 IERR( 1 ) = 0
413.GT. IF( IPW2+LW2MEMSIZ ) THEN
414.EQ. IF( IAM0 )
415 $ WRITE( NOUT, FMT = 9996 )'cond est',
416 $ ( IPW2+LW2 )*DBLESZ
417 IERR( 1 ) = 1
418 END IF
419*
420* Check all processes for an error
421*
422 CALL IGSUM2D( ICTXT, 'all', ' ', 1, 1, IERR, 1,
423 $ -1, 0 )
424*
425.GT. IF( IERR( 1 )0 ) THEN
426.EQ. IF( IAM0 )
427 $ WRITE( NOUT, FMT = 9997 ) 'memory'
428 KSKIP = KSKIP + 1
429 GO TO 60
430 END IF
431*
432 IF( CHECK ) THEN
433 CALL PDFILLPAD( ICTXT, LWORK, 1,
434 $ MEM( IPW-IPREPAD ), LWORK,
435 $ IPREPAD, IPOSTPAD, PADVAL )
436 CALL PDFILLPAD( ICTXT, LW2-IPOSTPAD, 1,
437 $ MEM( IPW2-IPREPAD ),
438 $ LW2-IPOSTPAD, IPREPAD,
439 $ IPOSTPAD, PADVAL )
440 END IF
441*
442* Compute condition number of the matrix
443*
444 CALL PDPOCON( UPLO, N, MEM( IPA ), 1, 1, DESCA,
445 $ ANORM1, RCOND, MEM( IPW ), LWORK,
446 $ MEM( IPW2 ), LIWORK, INFO )
447*
448 IF( CHECK ) THEN
449 CALL PDCHEKPAD( ICTXT, 'pdpocon', NP, NQ,
450 $ MEM( IPA-IPREPAD ), DESCA( LLD_ ),
451 $ IPREPAD, IPOSTPAD, PADVAL )
452 CALL PDCHEKPAD( ICTXT, 'pdpocon',
453 $ LWORK, 1, MEM( IPW-IPREPAD ),
454 $ LWORK, IPREPAD, IPOSTPAD,
455 $ PADVAL )
456 CALL PDCHEKPAD( ICTXT, 'pdpocon',
457 $ LW2-IPOSTPAD, 1,
458 $ MEM( IPW2-IPREPAD ), LW2-IPOSTPAD,
459 $ IPREPAD, IPOSTPAD, PADVAL )
460 END IF
461 END IF
462*
463* Loop over the different values for NRHS
464*
465 DO 20 HH = 1, NNR
466*
467 NRHS = NRVAL( HH )
468*
469 DO 10 KK = 1, NNBR
470*
471 NBRHS = NBRVAL( KK )
472*
473* Initialize Array Descriptor for rhs
474*
475 CALL DESCINIT( DESCB, N, NRHS, NB, NBRHS, 0, 0,
476 $ ICTXT, MAX( 1, NP )+IMIDPAD,
477 $ IERR( 1 ) )
478*
479* move IPW to allow room for RHS
480*
481 MYRHS = NUMROC( DESCB( N_ ), DESCB( NB_ ), MYCOL,
482 $ DESCB( CSRC_ ), NPCOL )
483 IPB = IPW
484*
485 IF( EST ) THEN
486 IPB0 = IPB + DESCB( LLD_ )*MYRHS + IPOSTPAD +
487 $ IPREPAD
488 IPFERR = IPB0 + DESCB( LLD_ )*MYRHS + IPOSTPAD
489 $ + IPREPAD
490 IPBERR = MYRHS + IPFERR + IPOSTPAD + IPREPAD
491 IPW = MYRHS + IPBERR + IPOSTPAD + IPREPAD
492 ELSE
493 IPW = IPB + DESCB( LLD_ )*MYRHS + IPOSTPAD +
494 $ IPREPAD
495 END IF
496*
497 IF( CHECK ) THEN
498*
499* Calculate the amount of workspace required by
500* the checking routines PDLASCHK
501*
502 LCMQ = LCM / NPCOL
503 WORKSIZ = MAX( WORKSIZ-IPOSTPAD,
504 $ NQ * NBRHS + NP * NBRHS +
505 $ MAX( MAX( NQ*NB, 2*NBRHS ),
506 $ NBRHS * NUMROC( NUMROC(N,NB,0,0,NPCOL),NB,
507 $ 0,0,LCMQ ) ) )
508 WORKSIZ = IPOSTPAD + WORKSIZ
509 ELSE
510 WORKSIZ = IPOSTPAD
511 END IF
512*
513 IERR( 1 ) = 0
514.GT. IF( IPW+WORKSIZMEMSIZ ) THEN
515.EQ. IF( IAM0 )
516 $ WRITE( NOUT, FMT = 9996 )'solve',
517 $ ( IPW+WORKSIZ )*DBLESZ
518 IERR( 1 ) = 1
519 END IF
520*
521* Check all processes for an error
522*
523 CALL IGSUM2D( ICTXT, 'all', ' ', 1, 1, IERR, 1,
524 $ -1, 0 )
525*
526.GT. IF( IERR( 1 )0 ) THEN
527.EQ. IF( IAM0 )
528 $ WRITE( NOUT, FMT = 9997 ) 'memory'
529 KSKIP = KSKIP + 1
530 GO TO 10
531 END IF
532*
533* Generate RHS
534*
535 CALL PDMATGEN( ICTXT, 'no', 'no', DESCB( M_ ),
536 $ DESCB( N_ ), DESCB( MB_ ),
537 $ DESCB( NB_ ), MEM( IPB ),
538 $ DESCB( LLD_ ), DESCB( RSRC_ ),
539 $ DESCB( CSRC_ ), IBSEED, 0, NP, 0,
540 $ MYRHS, MYROW, MYCOL, NPROW, NPCOL )
541*
542 IF( CHECK )
543 $ CALL PDFILLPAD( ICTXT, NP, MYRHS,
544 $ MEM( IPB-IPREPAD ),
545 $ DESCB( LLD_ ),
546 $ IPREPAD, IPOSTPAD, PADVAL )
547*
548 IF( EST ) THEN
549 CALL PDMATGEN( ICTXT, 'no', 'no', DESCB( M_ ),
550 $ DESCB( N_ ), DESCB( MB_ ),
551 $ DESCB( NB_ ), MEM( IPB0 ),
552 $ DESCB( LLD_ ), DESCB( RSRC_ ),
553 $ DESCB( CSRC_ ), IBSEED, 0, NP, 0,
554 $ MYRHS, MYROW, MYCOL, NPROW,
555 $ NPCOL )
556*
557 IF( CHECK ) THEN
558 CALL PDFILLPAD( ICTXT, NP, MYRHS,
559 $ MEM( IPB0-IPREPAD ),
560 $ DESCB( LLD_ ), IPREPAD,
561 $ IPOSTPAD, PADVAL )
562 CALL PDFILLPAD( ICTXT, 1, MYRHS,
563 $ MEM( IPFERR-IPREPAD ), 1,
564 $ IPREPAD, IPOSTPAD,
565 $ PADVAL )
566 CALL PDFILLPAD( ICTXT, 1, MYRHS,
567 $ MEM( IPBERR-IPREPAD ), 1,
568 $ IPREPAD, IPOSTPAD,
569 $ PADVAL )
570 END IF
571 END IF
572*
573 CALL BLACS_BARRIER( ICTXT, 'all' )
574 CALL SLTIMER( 2 )
575*
576* Solve linear system via Cholesky factorization
577*
578 CALL PDPOTRS( UPLO, N, NRHS, MEM( IPA ), 1, 1,
579 $ DESCA, MEM( IPB ), 1, 1, DESCB,
580 $ INFO )
581*
582 CALL SLTIMER( 2 )
583*
584 IF( CHECK ) THEN
585*
586* check for memory overwrite
587*
588 CALL PDCHEKPAD( ICTXT, 'pdpotrs', NP, NQ,
589 $ MEM( IPA-IPREPAD ),
590 $ DESCA( LLD_ ),
591 $ IPREPAD, IPOSTPAD, PADVAL )
592 CALL PDCHEKPAD( ICTXT, 'pdpotrs', NP,
593 $ MYRHS, MEM( IPB-IPREPAD ),
594 $ DESCB( LLD_ ), IPREPAD,
595 $ IPOSTPAD, PADVAL )
596*
597 CALL PDFILLPAD( ICTXT, WORKSIZ-IPOSTPAD, 1,
598 $ MEM( IPW-IPREPAD ),
599 $ WORKSIZ-IPOSTPAD, IPREPAD,
600 $ IPOSTPAD, PADVAL )
601*
602* check the solution to rhs
603*
604 CALL PDLASCHK( 'symm', 'diag', N, NRHS,
605 $ MEM( IPB ), 1, 1, DESCB,
606 $ IASEED, 1, 1, DESCA, IBSEED,
607 $ ANORM, SRESID, MEM( IPW ) )
608*
609.EQ..AND..GT. IF( IAM0 SRESIDTHRESH )
610 $ WRITE( NOUT, FMT = 9985 ) SRESID
611*
612* check for memory overwrite
613*
614 CALL PDCHEKPAD( ICTXT, 'pdlaschk', NP,
615 $ MYRHS, MEM( IPB-IPREPAD ),
616 $ DESCB( LLD_ ), IPREPAD,
617 $ IPOSTPAD, PADVAL )
618 CALL PDCHEKPAD( ICTXT, 'pdlaschk',
619 $ WORKSIZ-IPOSTPAD, 1,
620 $ MEM( IPW-IPREPAD ),
621 $ WORKSIZ-IPOSTPAD, IPREPAD,
622 $ IPOSTPAD, PADVAL )
623*
624* The second test is a NaN trap
625*
626.LE..AND. IF( ( SRESIDTHRESH )
627.EQ. $ ( (SRESID-SRESID)0.0D+0 ) ) THEN
628 KPASS = KPASS + 1
629 PASSED = 'passed'
630 ELSE
631 KFAIL = KFAIL + 1
632 PASSED = 'failed'
633 END IF
634 ELSE
635 KPASS = KPASS + 1
636 SRESID = SRESID - SRESID
637 PASSED = 'bypass'
638 END IF
639*
640 IF( EST ) THEN
641*
642* Calculate workspace required for PDPORFS
643*
644 LWORK = MAX( 1, 3*NP )
645 IPW2 = IPW + LWORK + IPOSTPAD + IPREPAD
646 LIWORK = MAX( 1, NP )
647 LW2 = ICEIL( LIWORK*INTGSZ, DBLESZ ) +
648 $ IPOSTPAD
649*
650 IERR( 1 ) = 0
651.GT. IF( IPW2+LW2MEMSIZ ) THEN
652.EQ. IF( IAM0 )
653 $ WRITE( NOUT, FMT = 9996 )
654 $ 'iter ref', ( IPW2+LW2 )*DBLESZ
655 IERR( 1 ) = 1
656 END IF
657*
658* Check all processes for an error
659*
660 CALL IGSUM2D( ICTXT, 'all', ' ', 1, 1, IERR,
661 $ 1, -1, 0 )
662*
663.GT. IF( IERR( 1 )0 ) THEN
664.EQ. IF( IAM0 )
665 $ WRITE( NOUT, FMT = 9997 )
666 $ 'memory'
667 KSKIP = KSKIP + 1
668 GO TO 10
669 END IF
670*
671 IF( CHECK ) THEN
672 CALL PDFILLPAD( ICTXT, LWORK, 1,
673 $ MEM( IPW-IPREPAD ),
674 $ LWORK, IPREPAD, IPOSTPAD,
675 $ PADVAL )
676 CALL PDFILLPAD( ICTXT, LW2-IPOSTPAD,
677 $ 1, MEM( IPW2-IPREPAD ),
678 $ LW2-IPOSTPAD,
679 $ IPREPAD, IPOSTPAD,
680 $ PADVAL )
681 END IF
682*
683* Use iterative refinement to improve the
684* computed solution
685*
686 CALL PDPORFS( UPLO, N, NRHS, MEM( IPA0 ),
687 $ 1, 1, DESCA, MEM( IPA ), 1, 1,
688 $ DESCA, MEM( IPB0 ), 1, 1,
689 $ DESCB, MEM( IPB ), 1, 1, DESCB,
690 $ MEM( IPFERR ), MEM( IPBERR ),
691 $ MEM( IPW ), LWORK, MEM( IPW2 ),
692 $ LIWORK, INFO )
693*
694* check for memory overwrite
695*
696 IF( CHECK ) THEN
697 CALL PDCHEKPAD( ICTXT, 'pdporfs', NP,
698 $ NQ, MEM( IPA0-IPREPAD ),
699 $ DESCA( LLD_ ), IPREPAD,
700 $ IPOSTPAD, PADVAL )
701 CALL PDCHEKPAD( ICTXT, 'pdporfs', NP,
702 $ NQ, MEM( IPA-IPREPAD ),
703 $ DESCA( LLD_ ), IPREPAD,
704 $ IPOSTPAD, PADVAL )
705 CALL PDCHEKPAD( ICTXT, 'pdporfs', NP,
706 $ MYRHS, MEM( IPB-IPREPAD ),
707 $ DESCB( LLD_ ), IPREPAD,
708 $ IPOSTPAD, PADVAL )
709 CALL PDCHEKPAD( ICTXT, 'pdporfs', NP,
710 $ MYRHS,
711 $ MEM( IPB0-IPREPAD ),
712 $ DESCB( LLD_ ), IPREPAD,
713 $ IPOSTPAD, PADVAL )
714 CALL PDCHEKPAD( ICTXT, 'pdporfs', 1,
715 $ MYRHS,
716 $ MEM( IPFERR-IPREPAD ), 1,
717 $ IPREPAD, IPOSTPAD,
718 $ PADVAL )
719 CALL PDCHEKPAD( ICTXT, 'pdporfs', 1,
720 $ MYRHS,
721 $ MEM( IPBERR-IPREPAD ), 1,
722 $ IPREPAD, IPOSTPAD,
723 $ PADVAL )
724 CALL PDCHEKPAD( ICTXT, 'pdporfs', LWORK,
725 $ 1, MEM( IPW-IPREPAD ),
726 $ LWORK, IPREPAD, IPOSTPAD,
727 $ PADVAL )
728 CALL PDCHEKPAD( ICTXT, 'pdporfs',
729 $ LW2-IPOSTPAD, 1,
730 $ MEM( IPW2-IPREPAD ),
731 $ LW2-IPOSTPAD,
732 $ IPREPAD, IPOSTPAD,
733 $ PADVAL )
734*
735 CALL PDFILLPAD( ICTXT, WORKSIZ-IPOSTPAD,
736 $ 1, MEM( IPW-IPREPAD ),
737 $ WORKSIZ-IPOSTPAD, IPREPAD,
738 $ IPOSTPAD, PADVAL )
739*
740* check the solution to rhs
741*
742 CALL PDLASCHK( 'symm', 'diag', N, NRHS,
743 $ MEM( IPB ), 1, 1, DESCB,
744 $ IASEED, 1, 1, DESCA,
745 $ IBSEED, ANORM, SRESID2,
746 $ MEM( IPW ) )
747*
748.EQ..AND..GT. IF( IAM0 SRESID2THRESH )
749 $ WRITE( NOUT, FMT = 9985 ) SRESID2
750*
751* check for memory overwrite
752*
753 CALL PDCHEKPAD( ICTXT, 'pdlaschk', NP,
754 $ MYRHS, MEM( IPB-IPREPAD ),
755 $ DESCB( LLD_ ), IPREPAD,
756 $ IPOSTPAD, PADVAL )
757 CALL PDCHEKPAD( ICTXT, 'pdlaschk',
758 $ WORKSIZ-IPOSTPAD, 1,
759 $ MEM( IPW-IPREPAD ),
760 $ WORKSIZ-IPOSTPAD,
761 $ IPREPAD, IPOSTPAD,
762 $ PADVAL )
763 END IF
764 END IF
765*
766* Gather maximum of all CPU and WALL clock timings
767*
768 CALL SLCOMBINE( ICTXT, 'all', '>', 'w', 2, 1,
769 $ WTIME )
770 CALL SLCOMBINE( ICTXT, 'all', '>', 'c', 2, 1,
771 $ CTIME )
772*
773* Print results
774*
775.EQ..AND..EQ. IF( MYROW0 MYCOL0 ) THEN
776*
777* 1/3 N^3 + 1/2 N^2 flops for LLt factorization
778*
779 NOPS = (DBLE(N)**3)/3.0D+0 +
780 $ (DBLE(N)**2)/2.0D+0
781*
782* nrhs * 2 N^2 flops for LLt solve.
783*
784 NOPS = NOPS + 2.0D+0*(DBLE(N)**2)*DBLE(NRHS)
785*
786* Calculate total megaflops -- factorization and
787* solve -- for WALL and CPU time, and print output
788*
789* Print WALL time if machine supports it
790*
791.GT. IF( WTIME( 1 ) + WTIME( 2 ) 0.0D+0 ) THEN
792 TMFLOPS = NOPS /
793 $ ( ( WTIME( 1 )+WTIME( 2 ) ) * 1.0D+6 )
794 ELSE
795 TMFLOPS = 0.0D+0
796 END IF
797*
798.GE. IF( WTIME( 2 )0.0D+0 )
799 $ WRITE( NOUT, FMT = 9993 ) 'wall', UPLO, N,
800 $ NB, NRHS, NBRHS, NPROW, NPCOL,
801 $ WTIME( 1 ), WTIME( 2 ), TMFLOPS,
802 $ PASSED
803*
804* Print CPU time if machine supports it
805*
806.GT. IF( CTIME( 1 )+CTIME( 2 )0.0D+0 ) THEN
807 TMFLOPS = NOPS /
808 $ ( ( CTIME( 1 )+CTIME( 2 ) ) * 1.0D+6 )
809 ELSE
810 TMFLOPS = 0.0D+0
811 END IF
812*
813.GE. IF( CTIME( 2 )0.0D+0 )
814 $ WRITE( NOUT, FMT = 9993 ) 'cpu ', UPLO, N,
815 $ NB, NRHS, NBRHS, NPROW, NPCOL,
816 $ CTIME( 1 ), CTIME( 2 ), TMFLOPS,
817 $ PASSED
818*
819 END IF
820 10 CONTINUE
821 20 CONTINUE
822*
823.AND..GT. IF( CHECK SRESIDTHRESH ) THEN
824*
825* Compute FRESID = ||A - LL'|| / (||A|| * N * eps)
826*
827 CALL PDPOTRRV( UPLO, N, MEM( IPA ), 1, 1, DESCA,
828 $ MEM( IPW ) )
829 CALL PDLAFCHK( 'symm', 'diag', N, N, MEM( IPA ), 1, 1,
830 $ DESCA, IASEED, ANORM, FRESID,
831 $ MEM( IPW ) )
832*
833* Check for memory overwrite
834*
835 CALL PDCHEKPAD( ICTXT, 'pdpotrrv', NP, NQ,
836 $ MEM( IPA-IPREPAD ), DESCA( LLD_ ),
837 $ IPREPAD, IPOSTPAD, PADVAL )
838 CALL PDCHEKPAD( ICTXT, 'pdgetrrv',
839 $ WORKSIZ-IPOSTPAD, 1,
840 $ MEM( IPW-IPREPAD ), WORKSIZ-IPOSTPAD,
841 $ IPREPAD, IPOSTPAD, PADVAL )
842*
843.EQ. IF( IAM0 ) THEN
844 IF( LSAME( UPLO, 'l' ) ) THEN
845 WRITE( NOUT, FMT = 9986 ) 'l*l''', FRESID
846 ELSE
847 WRITE( NOUT, FMT = 9986 ) 'u''*u', FRESID
848 END IF
849 END IF
850 END IF
851*
852 30 CONTINUE
853 40 CONTINUE
854 CALL BLACS_GRIDEXIT( ICTXT )
855 50 CONTINUE
856*
857* Print ending messages and close output file
858*
859 60 CONTINUE
860.EQ. IF( IAM0 ) THEN
861 KTESTS = KPASS + KFAIL + KSKIP
862 WRITE( NOUT, FMT = * )
863 WRITE( NOUT, FMT = 9992 ) KTESTS
864 IF( CHECK ) THEN
865 WRITE( NOUT, FMT = 9991 ) KPASS
866 WRITE( NOUT, FMT = 9989 ) KFAIL
867 ELSE
868 WRITE( NOUT, FMT = 9990 ) KPASS
869 END IF
870 WRITE( NOUT, FMT = 9988 ) KSKIP
871 WRITE( NOUT, FMT = * )
872 WRITE( NOUT, FMT = * )
873 WRITE( NOUT, FMT = 9987 )
874.NE..AND..NE. IF( NOUT6 NOUT0 )
875 $ CLOSE ( NOUT )
876 END IF
877*
878 CALL BLACS_EXIT( 0 )
879*
880 9999 FORMAT( 'illegal ', A6, ': ', A5, ' = ', I3,
881 $ '; it should be at least 1' )
882 9998 FORMAT( 'illegal grid: nprow*npcol = ', I4, '. it can be at most',
883 $ I4 )
884 9997 FORMAT( 'bad ', A6, ' parameters: going on to next test case.' )
885 9996 FORMAT( 'unable to perform ', A, ': need totmem of at least',
886 $ i11 )
887 9995 FORMAT( 'TIME UPLO N NB NRHS NBRHS P Q LLt Time ',
888 $ 'Slv Time MFLOPS CHECK' )
889 9994 FORMAT( '---- ---- ----- --- ---- ----- ---- ---- -------- ',
890 $ '-------- -------- ------' )
891 9993 FORMAT( a4, 4x, a1, 1x, i5, 1x, i3, 1x, i4, 1x, i5, 1x, i4, 1x,
892 $ i4, 1x, f8.2, 1x, f8.2, 1x, f8.2, 1x, a6 )
893 9992 FORMAT( 'Finished ', i6, ' tests, with the following results:' )
894 9991 FORMAT( i5, ' tests completed and passed residual checks.' )
895 9990 FORMAT( i5, ' tests completed without checking.' )
896 9989 FORMAT( i5, ' tests completed and failed residual checks.' )
897 9988 FORMAT( i5, ' tests skipped because of illegal input values.' )
898 9987 FORMAT( 'END OF TESTS.' )
899 9986 FORMAT( '||A - ', a4, '|| / (||A|| * N * eps) = ', g25.7 )
900 9985 FORMAT( '||Ax-b||/(||x||*||A||*eps*N) ', f25.7 )
901*
902 stop
903*
904* End of PDLLTDRIVER
905*
906 END
pdlafchk
subroutine pdlafchk(aform, diag, m, n, a, ia, ja, desca, iaseed, anorm, fresid, work)
Definition pdlafchk.f:3
pdmatgen
subroutine pdmatgen(ictxt, aform, diag, m, n, mb, nb, a, lda, iarow, iacol, iseed, iroff, irnum, icoff, icnum, myrow, mycol, nprow, npcol)
Definition pdmatgen.f:4
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
iceil
integer function iceil(inum, idenom)
Definition iceil.f:2
ilcm
integer function ilcm(m, n)
Definition ilcm.f:2
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
blacs_gridinit
subroutine blacs_gridinit(cntxt, c, nprow, npcol)
Definition mpi.f:745
pdpotrs
subroutine pdpotrs(uplo, n, nrhs, a, ia, ja, desca, b, ib, jb, descb, info)
Definition mpi.f:1220
pdpotrf
subroutine pdpotrf(uplo, n, a, ia, ja, desca, info)
Definition mpi.f:903
descinit
subroutine descinit(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld, info)
Definition mpi.f:777
blacs_gridexit
subroutine blacs_gridexit(cntxt)
Definition mpi.f:762
blacs_gridinfo
subroutine blacs_gridinfo(cntxt, nprow, npcol, myrow, mycol)
Definition mpi.f:754
numroc
integer function numroc(n, nb, iproc, isrcproc, nprocs)
Definition mpi.f:786
pdchekpad
subroutine pdchekpad(ictxt, mess, m, n, a, lda, ipre, ipost, chkval)
Definition pdchekpad.f:3
pdfillpad
subroutine pdfillpad(ictxt, m, n, a, lda, ipre, ipost, chkval)
Definition pdfillpad.f:2
pdgetrrv
subroutine pdgetrrv(m, n, a, ia, ja, desca, ipiv, work)
Definition pdgetrrv.f:2
pdlansy
double precision function pdlansy(norm, uplo, n, a, ia, ja, desca, work)
Definition pdlansy.f:3
pdlaschk
subroutine pdlaschk(symm, diag, n, nrhs, x, ix, jx, descx, iaseed, ia, ja, desca, ibseed, anorm, resid, work)
Definition pdlaschk.f:4
pdlltdriver
program pdlltdriver
Definition pdlltdriver.f:1
pdlltinfo
subroutine pdlltinfo(summry, nout, uplo, nmat, nval, ldnval, nnb, nbval, ldnbval, nnr, nrval, ldnrval, nnbr, nbrval, ldnbrval, ngrids, pval, ldpval, qval, ldqval, thresh, est, work, iam, nprocs)
Definition pdlltinfo.f:6
pdpocon
subroutine pdpocon(uplo, n, a, ia, ja, desca, anorm, rcond, work, lwork, iwork, liwork, info)
Definition pdpocon.f:3
pdporfs
subroutine pdporfs(uplo, n, nrhs, a, ia, ja, desca, af, iaf, jaf, descaf, b, ib, jb, descb, x, ix, jx, descx, ferr, berr, work, lwork, iwork, liwork, info)
Definition pdporfs.f:4
pdpotrrv
subroutine pdpotrrv(uplo, n, a, ia, ja, desca, work)
Definition pdpotrrv.f:2
slboot
subroutine slboot()
Definition sltimer.f:2
sltimer
subroutine sltimer(i)
Definition sltimer.f:47
slcombine
subroutine slcombine(ictxt, scope, op, timetype, n, ibeg, times)
Definition sltimer.f:267