SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
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 IF( ierr( 1 ).GT.0 ) THEN
234 IF( iam.EQ.0 )
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 IF( ierr( 1 ).LT.0 ) THEN
264 IF( iam.EQ.0 )
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 IF( nprow.NE.npcol ) 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 IF( ipw+worksiz.GT.memsiz ) THEN
311 IF( iam.EQ.0 )
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 IF( ierr( 1 ).GT.0 ) THEN
322 IF( iam.EQ.0 )
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 IF( info.NE.0 ) THEN
383 IF( iam.EQ.0 )
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 IF( ipw2+lw2.GT.memsiz ) THEN
414 IF( iam.EQ.0 )
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 IF( ierr( 1 ).GT.0 ) THEN
426 IF( iam.EQ.0 )
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 IF( ipw+worksiz.GT.memsiz ) THEN
515 IF( iam.EQ.0 )
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 IF( ierr( 1 ).GT.0 ) THEN
527 IF( iam.EQ.0 )
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 IF( iam.EQ.0 .AND. sresid.GT.thresh )
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 IF( ( sresid.LE.thresh ).AND.
627 $ ( (sresid-sresid).EQ.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 IF( ipw2+lw2.GT.memsiz ) THEN
652 IF( iam.EQ.0 )
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 IF( ierr( 1 ).GT.0 ) THEN
664 IF( iam.EQ.0 )
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 IF( iam.EQ.0 .AND. sresid2.GT.thresh )
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 IF( myrow.EQ.0 .AND. mycol.EQ.0 ) 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 IF( wtime( 1 ) + wtime( 2 ) .GT. 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 IF( wtime( 2 ).GE.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 IF( ctime( 1 )+ctime( 2 ).GT.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 IF( ctime( 2 ).GE.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 IF( check .AND. sresid.GT.thresh ) 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 IF( iam.EQ.0 ) 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 IF( iam.EQ.0 ) 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 IF( nout.NE.6 .AND. nout.NE.0 )
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
descinit
subroutine descinit(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld, info)
Definition descinit.f:3
iceil
integer function iceil(inum, idenom)
Definition iceil.f:2
ilcm
integer function ilcm(m, n)
Definition ilcm.f:2
numroc
integer function numroc(n, nb, iproc, isrcproc, nprocs)
Definition numroc.f:2
max
#define max(A, B)
Definition pcgemr.c:180
min
#define min(A, B)
Definition pcgemr.c:181
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
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
pdpotrf
subroutine pdpotrf(uplo, n, a, ia, ja, desca, info)
Definition pdpotrf.f:2
pdpotrrv
subroutine pdpotrrv(uplo, n, a, ia, ja, desca, work)
Definition pdpotrrv.f:2
pdpotrs
subroutine pdpotrs(uplo, n, nrhs, a, ia, ja, desca, b, ib, jb, descb, info)
Definition pdpotrs.f:3
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
lsame
logical function lsame(ca, cb)
Definition tools.f:1724