SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
•All Classes Files Functions Variables Typedefs Macros
pslltdriver.f
Go to the documentation of this file.
1 PROGRAM pslltdriver
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* PSLLTDRIVER is the main test program for the REAL
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* REALSZ INTEGER, default = 4 bytes.
56* INTGSZ and REALSZ indicate the length in bytes on the
57* given platform for an integer and a single precision real.
58* MEM REAL array, dimension ( TOTMEM / REALSZ )
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 intgsz, memsiz, ntests, realsz, totmem
74 REAL padval, zero
75 parameter( intgsz = 4, realsz = 4, totmem = 2000000,
76 $ memsiz = totmem / realsz, ntests = 20,
77 $ padval = -9923.0e+0, zero = 0.0e+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 anorm, anorm1, fresid, rcond, sresid, sresid2,
92 $ thresh
93 DOUBLE PRECISION nops, 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 REAL mem( memsiz )
101 DOUBLE PRECISION ctime( 2 ), wtime( 2 )
102* ..
103* .. External Subroutines ..
104 EXTERNAL blacs_barrier, blacs_exit, blacs_gridexit,
105 $ blacs_gridinfo, blacs_gridinit, descinit,
106 $ igsum2d, blacs_pinfo, pschekpad, psfillpad,
107 $ pslafchk, pslaschk, pslltinfo,
108 $ psmatgen, pspocon, psporfs,
109 $ pspotrf, pspotrrv, pspotrs, slboot,
110 $ slcombine, sltimer
111* ..
112* .. External Functions ..
113 LOGICAL lsame
114 INTEGER iceil, ilcm, numroc
115 REAL pslansy
116 EXTERNAL iceil, ilcm, lsame, numroc, pslansy
117* ..
118* .. Intrinsic Functions ..
119 INTRINSIC dble, max, min
120* ..
121* .. Data Statements ..
122 DATA kfail, kpass, kskip, ktests / 4*0 /
123* ..
124* .. Executable Statements ..
125*
126* Get starting information
127*
128 CALL blacs_pinfo( iam, nprocs )
129 iaseed = 100
130 ibseed = 200
131 CALL pslltinfo( outfile, nout, uplo, nmat, nval, ntests, nnb,
132 $ nbval, ntests, nnr, nrval, ntests, nnbr, nbrval,
133 $ ntests, ngrids, pval, ntests, qval, ntests,
134 $ thresh, est, mem, iam, nprocs )
135 check = ( thresh.GE.0.0e+0 )
136*
137* Print headings
138*
139 IF( iam.EQ.0 ) THEN
140 WRITE( nout, fmt = * )
141 WRITE( nout, fmt = 9995 )
142 WRITE( nout, fmt = 9994 )
143 WRITE( nout, fmt = * )
144 END IF
145*
146* Loop over different process grids
147*
148 DO 50 i = 1, ngrids
149*
150 nprow = pval( i )
151 npcol = qval( i )
152*
153* Make sure grid information is correct
154*
155 ierr( 1 ) = 0
156 IF( nprow.LT.1 ) THEN
157 IF( iam.EQ.0 )
158 $ WRITE( nout, fmt = 9999 ) 'GRID', 'nprow', nprow
159 ierr( 1 ) = 1
160 ELSE IF( npcol.LT.1 ) THEN
161 IF( iam.EQ.0 )
162 $ WRITE( nout, fmt = 9999 ) 'GRID', 'npcol', npcol
163 ierr( 1 ) = 1
164 ELSE IF( nprow*npcol.GT.nprocs ) THEN
165 IF( iam.EQ.0 )
166 $ WRITE( nout, fmt = 9998 ) nprow*npcol, nprocs
167 ierr( 1 ) = 1
168 END IF
169*
170 IF( ierr( 1 ).GT.0 ) THEN
171 IF( iam.EQ.0 )
172 $ WRITE( nout, fmt = 9997 ) 'grid'
173 kskip = kskip + 1
174 GO TO 50
175 END IF
176*
177* Define process grid
178*
179 CALL blacs_get( -1, 0, ictxt )
180 CALL blacs_gridinit( ictxt, 'Row-major', nprow, npcol )
181 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
182*
183* Go to bottom of process grid loop if this case doesn't use my
184* process
185*
186 IF( myrow.GE.nprow .OR. mycol.GE.npcol )
187 $ GO TO 50
188*
189 DO 40 j = 1, nmat
190*
191 n = nval( j )
192*
193* Make sure matrix information is correct
194*
195 ierr( 1 ) = 0
196 IF( n.LT.1 ) THEN
197 IF( iam.EQ.0 )
198 $ WRITE( nout, fmt = 9999 ) 'MATRIX', 'N', n
199 ierr( 1 ) = 1
200 ELSE IF( n.LT.1 ) THEN
201 IF( iam.EQ.0 )
202 $ WRITE( nout, fmt = 9999 ) 'MATRIX', 'N', n
203 ierr( 1 ) = 1
204 END IF
205*
206* Check all processes for an error
207*
208 CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1, 0 )
209*
210 IF( ierr( 1 ).GT.0 ) THEN
211 IF( iam.EQ.0 )
212 $ WRITE( nout, fmt = 9997 ) 'matrix'
213 kskip = kskip + 1
214 GO TO 40
215 END IF
216*
217 DO 30 k = 1, nnb
218*
219 nb = nbval( k )
220*
221* Make sure nb is legal
222*
223 ierr( 1 ) = 0
224 IF( nb.LT.1 ) THEN
225 ierr( 1 ) = 1
226 IF( iam.EQ.0 )
227 $ WRITE( nout, fmt = 9999 ) 'NB', 'NB', nb
228 END IF
229*
230* Check all processes for an error
231*
232 CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1, 0 )
233*
234 IF( ierr( 1 ).GT.0 ) THEN
235 IF( iam.EQ.0 )
236 $ WRITE( nout, fmt = 9997 ) 'NB'
237 kskip = kskip + 1
238 GO TO 30
239 END IF
240*
241* Padding constants
242*
243 np = numroc( n, nb, myrow, 0, nprow )
244 nq = numroc( n, nb, mycol, 0, npcol )
245 IF( check ) THEN
246 iprepad = max( nb, np )
247 imidpad = nb
248 ipostpad = max( nb, nq )
249 ELSE
250 iprepad = 0
251 imidpad = 0
252 ipostpad = 0
253 END IF
254*
255* Initialize the array descriptor for the matrix A
256*
257 CALL descinit( desca, n, n, nb, nb, 0, 0, ictxt,
258 $ max( 1, np )+imidpad, ierr( 1 ) )
259*
260* Check all processes for an error
261*
262 CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1, 0 )
263*
264 IF( ierr( 1 ).LT.0 ) THEN
265 IF( iam.EQ.0 )
266 $ WRITE( nout, fmt = 9997 ) 'descriptor'
267 kskip = kskip + 1
268 GO TO 30
269 END IF
270*
271* Assign pointers into MEM for SCALAPACK arrays, A is
272* allocated starting at position MEM( IPREPAD+1 )
273*
274 ipa = iprepad+1
275 IF( est ) THEN
276 ipa0 = ipa + desca( lld_ )*nq + ipostpad + iprepad
277 ipw = ipa0 + desca( lld_ )*nq + ipostpad + iprepad
278 ELSE
279 ipw = ipa + desca( lld_ )*nq + ipostpad + iprepad
280 END IF
281*
282*
283 IF( check ) THEN
284*
285* Calculate the amount of workspace required by
286* the checking routines PSLAFCHK, PSPOTRRV, and
287* PSLANSY
288*
289 worksiz = np * desca( nb_ )
290*
291 worksiz = max( worksiz, desca( mb_ ) * desca( nb_ ) )
292*
293 lcm = ilcm( nprow, npcol )
294 itemp = max( 2, 2 * nq ) + np
295 IF( nprow.NE.npcol ) THEN
296 itemp = itemp +
297 $ nb * iceil( iceil( np, nb ), lcm / nprow )
298 END IF
299 worksiz = max( worksiz, itemp )
300 worksiz = worksiz + ipostpad
301*
302 ELSE
303*
304 worksiz = ipostpad
305*
306 END IF
307*
308* Check for adequate memory for problem size
309*
310 ierr( 1 ) = 0
311 IF( ipw+worksiz.GT.memsiz ) THEN
312 IF( iam.EQ.0 )
313 $ WRITE( nout, fmt = 9996 ) 'factorization',
314 $ ( ipw+worksiz )*realsz
315 ierr( 1 ) = 1
316 END IF
317*
318* Check all processes for an error
319*
320 CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1, 0 )
321*
322 IF( ierr( 1 ).GT.0 ) THEN
323 IF( iam.EQ.0 )
324 $ WRITE( nout, fmt = 9997 ) 'MEMORY'
325 kskip = kskip + 1
326 GO TO 30
327 END IF
328*
329* Generate a symmetric positive definite matrix A
330*
331 CALL psmatgen( ictxt, 'Symm', 'Diag', desca( m_ ),
332 $ desca( n_ ), desca( mb_ ), desca( nb_ ),
333 $ mem( ipa ), desca( lld_ ), desca( rsrc_ ),
334 $ desca( csrc_ ), iaseed, 0, np, 0, nq,
335 $ myrow, mycol, nprow, npcol )
336*
337* Calculate inf-norm of A for residual error-checking
338*
339 IF( check ) THEN
340 CALL psfillpad( ictxt, np, nq, mem( ipa-iprepad ),
341 $ desca( lld_ ), iprepad, ipostpad,
342 $ padval )
343 CALL psfillpad( ictxt, worksiz-ipostpad, 1,
344 $ mem( ipw-iprepad ), worksiz-ipostpad,
345 $ iprepad, ipostpad, padval )
346 anorm = pslansy( 'I', uplo, n, mem( ipa ), 1, 1,
347 $ desca, mem( ipw ) )
348 anorm1 = pslansy( '1', uplo, n, mem( ipa ), 1, 1,
349 $ desca, mem( ipw ) )
350 CALL pschekpad( ictxt, 'PSLANSY', np, nq,
351 $ mem( ipa-iprepad ), desca( lld_ ),
352 $ iprepad, ipostpad, padval )
353 CALL pschekpad( ictxt, 'PSLANSY',
354 $ worksiz-ipostpad, 1,
355 $ mem( ipw-iprepad ), worksiz-ipostpad,
356 $ iprepad, ipostpad, padval )
357 END IF
358*
359 IF( est ) THEN
360 CALL psmatgen( ictxt, 'Symm', 'Diag', desca( m_ ),
361 $ desca( n_ ), desca( mb_ ),
362 $ desca( nb_ ), mem( ipa0 ),
363 $ desca( lld_ ), desca( rsrc_ ),
364 $ desca( csrc_ ), iaseed, 0, np, 0, nq,
365 $ myrow, mycol, nprow, npcol )
366 IF( check )
367 $ CALL psfillpad( ictxt, np, nq,
368 $ mem( ipa0-iprepad ), desca( lld_ ),
369 $ iprepad, ipostpad, padval )
370 END IF
371*
372 CALL slboot()
373 CALL blacs_barrier( ictxt, 'All' )
374*
375* Perform LLt factorization
376*
377 CALL sltimer( 1 )
378*
379 CALL pspotrf( uplo, n, mem( ipa ), 1, 1, desca, info )
380*
381 CALL sltimer( 1 )
382*
383 IF( info.NE.0 ) THEN
384 IF( iam.EQ.0 )
385 $ WRITE( nout, fmt = * ) 'PSPOTRF INFO=', info
386 kfail = kfail + 1
387 rcond = zero
388 GO TO 60
389 END IF
390*
391 IF( check ) THEN
392*
393* Check for memory overwrite in LLt factorization
394*
395 CALL pschekpad( ictxt, 'PSPOTRF', np, nq,
396 $ mem( ipa-iprepad ), desca( lld_ ),
397 $ iprepad, ipostpad, padval )
398 END IF
399*
400 IF( est ) THEN
401*
402* Calculate workspace required for PSPOCON
403*
404 lwork = max( 1, 2*np ) + max( 1, 2*nq ) +
405 $ max( 2, desca( nb_ )*
406 $ max( 1, iceil( nprow-1, npcol ) ),
407 $ nq + desca( nb_ )*
408 $ max( 1, iceil( npcol-1, nprow ) ) )
409 ipw2 = ipw + lwork + ipostpad + iprepad
410 liwork = max( 1, np )
411 lw2 = iceil( liwork*intgsz, realsz ) + ipostpad
412*
413 ierr( 1 ) = 0
414 IF( ipw2+lw2.GT.memsiz ) THEN
415 IF( iam.EQ.0 )
416 $ WRITE( nout, fmt = 9996 )'cond est',
417 $ ( ipw2+lw2 )*realsz
418 ierr( 1 ) = 1
419 END IF
420*
421* Check all processes for an error
422*
423 CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1,
424 $ -1, 0 )
425*
426 IF( ierr( 1 ).GT.0 ) THEN
427 IF( iam.EQ.0 )
428 $ WRITE( nout, fmt = 9997 ) 'MEMORY'
429 kskip = kskip + 1
430 GO TO 60
431 END IF
432*
433 IF( check ) THEN
434 CALL psfillpad( ictxt, lwork, 1,
435 $ mem( ipw-iprepad ), lwork,
436 $ iprepad, ipostpad, padval )
437 CALL psfillpad( ictxt, lw2-ipostpad, 1,
438 $ mem( ipw2-iprepad ),
439 $ lw2-ipostpad, iprepad,
440 $ ipostpad, padval )
441 END IF
442*
443* Compute condition number of the matrix
444*
445 CALL pspocon( uplo, n, mem( ipa ), 1, 1, desca,
446 $ anorm1, rcond, mem( ipw ), lwork,
447 $ mem( ipw2 ), liwork, info )
448*
449 IF( check ) THEN
450 CALL pschekpad( ictxt, 'PSPOCON', np, nq,
451 $ mem( ipa-iprepad ), desca( lld_ ),
452 $ iprepad, ipostpad, padval )
453 CALL pschekpad( ictxt, 'PSPOCON',
454 $ lwork, 1, mem( ipw-iprepad ),
455 $ lwork, iprepad, ipostpad,
456 $ padval )
457 CALL pschekpad( ictxt, 'PSPOCON',
458 $ lw2-ipostpad, 1,
459 $ mem( ipw2-iprepad ), lw2-ipostpad,
460 $ iprepad, ipostpad, padval )
461 END IF
462 END IF
463*
464* Loop over the different values for NRHS
465*
466 DO 20 hh = 1, nnr
467*
468 nrhs = nrval( hh )
469*
470 DO 10 kk = 1, nnbr
471*
472 nbrhs = nbrval( kk )
473*
474* Initialize Array Descriptor for rhs
475*
476 CALL descinit( descb, n, nrhs, nb, nbrhs, 0, 0,
477 $ ictxt, max( 1, np )+imidpad,
478 $ ierr( 1 ) )
479*
480* move IPW to allow room for RHS
481*
482 myrhs = numroc( descb( n_ ), descb( nb_ ), mycol,
483 $ descb( csrc_ ), npcol )
484 ipb = ipw
485*
486 IF( est ) THEN
487 ipb0 = ipb + descb( lld_ )*myrhs + ipostpad +
488 $ iprepad
489 ipferr = ipb0 + descb( lld_ )*myrhs + ipostpad
490 $ + iprepad
491 ipberr = myrhs + ipferr + ipostpad + iprepad
492 ipw = myrhs + ipberr + ipostpad + iprepad
493 ELSE
494 ipw = ipb + descb( lld_ )*myrhs + ipostpad +
495 $ iprepad
496 END IF
497*
498 IF( check ) THEN
499*
500* Calculate the amount of workspace required by
501* the checking routines PSLASCHK
502*
503 lcmq = lcm / npcol
504 worksiz = max( worksiz-ipostpad,
505 $ nq * nbrhs + np * nbrhs +
506 $ max( max( nq*nb, 2*nbrhs ),
507 $ nbrhs * numroc( numroc(n,nb,0,0,npcol),nb,
508 $ 0,0,lcmq ) ) )
509 worksiz = ipostpad + worksiz
510 ELSE
511 worksiz = ipostpad
512 END IF
513*
514 ierr( 1 ) = 0
515 IF( ipw+worksiz.GT.memsiz ) THEN
516 IF( iam.EQ.0 )
517 $ WRITE( nout, fmt = 9996 )'solve',
518 $ ( ipw+worksiz )*realsz
519 ierr( 1 ) = 1
520 END IF
521*
522* Check all processes for an error
523*
524 CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1,
525 $ -1, 0 )
526*
527 IF( ierr( 1 ).GT.0 ) THEN
528 IF( iam.EQ.0 )
529 $ WRITE( nout, fmt = 9997 ) 'MEMORY'
530 kskip = kskip + 1
531 GO TO 10
532 END IF
533*
534* Generate RHS
535*
536 CALL psmatgen( ictxt, 'No', 'No', descb( m_ ),
537 $ descb( n_ ), descb( mb_ ),
538 $ descb( nb_ ), mem( ipb ),
539 $ descb( lld_ ), descb( rsrc_ ),
540 $ descb( csrc_ ), ibseed, 0, np, 0,
541 $ myrhs, myrow, mycol, nprow, npcol )
542*
543 IF( check )
544 $ CALL psfillpad( ictxt, np, myrhs,
545 $ mem( ipb-iprepad ),
546 $ descb( lld_ ),
547 $ iprepad, ipostpad, padval )
548*
549 IF( est ) THEN
550 CALL psmatgen( ictxt, 'No', 'No', descb( m_ ),
551 $ descb( n_ ), descb( mb_ ),
552 $ descb( nb_ ), mem( ipb0 ),
553 $ descb( lld_ ), descb( rsrc_ ),
554 $ descb( csrc_ ), ibseed, 0, np, 0,
555 $ myrhs, myrow, mycol, nprow,
556 $ npcol )
557*
558 IF( check ) THEN
559 CALL psfillpad( ictxt, np, myrhs,
560 $ mem( ipb0-iprepad ),
561 $ descb( lld_ ), iprepad,
562 $ ipostpad, padval )
563 CALL psfillpad( ictxt, 1, myrhs,
564 $ mem( ipferr-iprepad ), 1,
565 $ iprepad, ipostpad,
566 $ padval )
567 CALL psfillpad( ictxt, 1, myrhs,
568 $ mem( ipberr-iprepad ), 1,
569 $ iprepad, ipostpad,
570 $ padval )
571 END IF
572 END IF
573*
574 CALL blacs_barrier( ictxt, 'All' )
575 CALL sltimer( 2 )
576*
577* Solve linear system via Cholesky factorization
578*
579 CALL pspotrs( uplo, n, nrhs, mem( ipa ), 1, 1,
580 $ desca, mem( ipb ), 1, 1, descb,
581 $ info )
582*
583 CALL sltimer( 2 )
584*
585 IF( check ) THEN
586*
587* check for memory overwrite
588*
589 CALL pschekpad( ictxt, 'PSPOTRS', np, nq,
590 $ mem( ipa-iprepad ),
591 $ desca( lld_ ),
592 $ iprepad, ipostpad, padval )
593 CALL pschekpad( ictxt, 'PSPOTRS', np,
594 $ myrhs, mem( ipb-iprepad ),
595 $ descb( lld_ ), iprepad,
596 $ ipostpad, padval )
597*
598 CALL psfillpad( ictxt, worksiz-ipostpad, 1,
599 $ mem( ipw-iprepad ),
600 $ worksiz-ipostpad, iprepad,
601 $ ipostpad, padval )
602*
603* check the solution to rhs
604*
605 CALL pslaschk( 'Symm', 'Diag', n, nrhs,
606 $ mem( ipb ), 1, 1, descb,
607 $ iaseed, 1, 1, desca, ibseed,
608 $ anorm, sresid, mem( ipw ) )
609*
610 IF( iam.EQ.0 .AND. sresid.GT.thresh )
611 $ WRITE( nout, fmt = 9985 ) sresid
612*
613* check for memory overwrite
614*
615 CALL pschekpad( ictxt, 'PSLASCHK', np,
616 $ myrhs, mem( ipb-iprepad ),
617 $ descb( lld_ ), iprepad,
618 $ ipostpad, padval )
619 CALL pschekpad( ictxt, 'PSLASCHK',
620 $ worksiz-ipostpad, 1,
621 $ mem( ipw-iprepad ),
622 $ worksiz-ipostpad, iprepad,
623 $ ipostpad, padval )
624*
625* The second test is a NaN trap
626*
627 IF( ( sresid.LE.thresh ).AND.
628 $ ( (sresid-sresid).EQ.0.0e+0 ) ) THEN
629 kpass = kpass + 1
630 passed = 'PASSED'
631 ELSE
632 kfail = kfail + 1
633 passed = 'FAILED'
634 END IF
635 ELSE
636 kpass = kpass + 1
637 sresid = sresid - sresid
638 passed = 'BYPASS'
639 END IF
640*
641 IF( est ) THEN
642*
643* Calculate workspace required for PSPORFS
644*
645 lwork = max( 1, 3*np )
646 ipw2 = ipw + lwork + ipostpad + iprepad
647 liwork = max( 1, np )
648 lw2 = iceil( liwork*intgsz, realsz ) +
649 $ ipostpad
650*
651 ierr( 1 ) = 0
652 IF( ipw2+lw2.GT.memsiz ) THEN
653 IF( iam.EQ.0 )
654 $ WRITE( nout, fmt = 9996 )
655 $ 'iter ref', ( ipw2+lw2 )*realsz
656 ierr( 1 ) = 1
657 END IF
658*
659* Check all processes for an error
660*
661 CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr,
662 $ 1, -1, 0 )
663*
664 IF( ierr( 1 ).GT.0 ) THEN
665 IF( iam.EQ.0 )
666 $ WRITE( nout, fmt = 9997 )
667 $ 'MEMORY'
668 kskip = kskip + 1
669 GO TO 10
670 END IF
671*
672 IF( check ) THEN
673 CALL psfillpad( ictxt, lwork, 1,
674 $ mem( ipw-iprepad ),
675 $ lwork, iprepad, ipostpad,
676 $ padval )
677 CALL psfillpad( ictxt, lw2-ipostpad,
678 $ 1, mem( ipw2-iprepad ),
679 $ lw2-ipostpad,
680 $ iprepad, ipostpad,
681 $ padval )
682 END IF
683*
684* Use iterative refinement to improve the
685* computed solution
686*
687 CALL psporfs( uplo, n, nrhs, mem( ipa0 ),
688 $ 1, 1, desca, mem( ipa ), 1, 1,
689 $ desca, mem( ipb0 ), 1, 1,
690 $ descb, mem( ipb ), 1, 1, descb,
691 $ mem( ipferr ), mem( ipberr ),
692 $ mem( ipw ), lwork, mem( ipw2 ),
693 $ liwork, info )
694*
695* check for memory overwrite
696*
697 IF( check ) THEN
698 CALL pschekpad( ictxt, 'PSPORFS', np,
699 $ nq, mem( ipa0-iprepad ),
700 $ desca( lld_ ), iprepad,
701 $ ipostpad, padval )
702 CALL pschekpad( ictxt, 'PSPORFS', np,
703 $ nq, mem( ipa-iprepad ),
704 $ desca( lld_ ), iprepad,
705 $ ipostpad, padval )
706 CALL pschekpad( ictxt, 'PSPORFS', np,
707 $ myrhs, mem( ipb-iprepad ),
708 $ descb( lld_ ), iprepad,
709 $ ipostpad, padval )
710 CALL pschekpad( ictxt, 'PSPORFS', np,
711 $ myrhs,
712 $ mem( ipb0-iprepad ),
713 $ descb( lld_ ), iprepad,
714 $ ipostpad, padval )
715 CALL pschekpad( ictxt, 'PSPORFS', 1,
716 $ myrhs,
717 $ mem( ipferr-iprepad ), 1,
718 $ iprepad, ipostpad,
719 $ padval )
720 CALL pschekpad( ictxt, 'PSPORFS', 1,
721 $ myrhs,
722 $ mem( ipberr-iprepad ), 1,
723 $ iprepad, ipostpad,
724 $ padval )
725 CALL pschekpad( ictxt, 'PSPORFS', lwork,
726 $ 1, mem( ipw-iprepad ),
727 $ lwork, iprepad, ipostpad,
728 $ padval )
729 CALL pschekpad( ictxt, 'PSPORFS',
730 $ lw2-ipostpad, 1,
731 $ mem( ipw2-iprepad ),
732 $ lw2-ipostpad,
733 $ iprepad, ipostpad,
734 $ padval )
735*
736 CALL psfillpad( ictxt, worksiz-ipostpad,
737 $ 1, mem( ipw-iprepad ),
738 $ worksiz-ipostpad, iprepad,
739 $ ipostpad, padval )
740*
741* check the solution to rhs
742*
743 CALL pslaschk( 'Symm', 'Diag', n, nrhs,
744 $ mem( ipb ), 1, 1, descb,
745 $ iaseed, 1, 1, desca,
746 $ ibseed, anorm, sresid2,
747 $ mem( ipw ) )
748*
749 IF( iam.EQ.0 .AND. sresid2.GT.thresh )
750 $ WRITE( nout, fmt = 9985 ) sresid2
751*
752* check for memory overwrite
753*
754 CALL pschekpad( ictxt, 'PSLASCHK', np,
755 $ myrhs, mem( ipb-iprepad ),
756 $ descb( lld_ ), iprepad,
757 $ ipostpad, padval )
758 CALL pschekpad( ictxt, 'PSLASCHK',
759 $ worksiz-ipostpad, 1,
760 $ mem( ipw-iprepad ),
761 $ worksiz-ipostpad,
762 $ iprepad, ipostpad,
763 $ padval )
764 END IF
765 END IF
766*
767* Gather maximum of all CPU and WALL clock timings
768*
769 CALL slcombine( ictxt, 'All', '>', 'W', 2, 1,
770 $ wtime )
771 CALL slcombine( ictxt, 'All', '>', 'C', 2, 1,
772 $ ctime )
773*
774* Print results
775*
776 IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
777*
778* 1/3 N^3 + 1/2 N^2 flops for LLt factorization
779*
780 nops = (dble(n)**3)/3.0d+0 +
781 $ (dble(n)**2)/2.0d+0
782*
783* nrhs * 2 N^2 flops for LLt solve.
784*
785 nops = nops + 2.0d+0*(dble(n)**2)*dble(nrhs)
786*
787* Calculate total megaflops -- factorization and
788* solve -- for WALL and CPU time, and print output
789*
790* Print WALL time if machine supports it
791*
792 IF( wtime( 1 ) + wtime( 2 ) .GT. 0.0d+0 ) THEN
793 tmflops = nops /
794 $ ( ( wtime( 1 )+wtime( 2 ) ) * 1.0d+6 )
795 ELSE
796 tmflops = 0.0d+0
797 END IF
798*
799 IF( wtime( 2 ).GE.0.0d+0 )
800 $ WRITE( nout, fmt = 9993 ) 'WALL', uplo, n,
801 $ nb, nrhs, nbrhs, nprow, npcol,
802 $ wtime( 1 ), wtime( 2 ), tmflops,
803 $ passed
804*
805* Print CPU time if machine supports it
806*
807 IF( ctime( 1 )+ctime( 2 ).GT.0.0d+0 ) THEN
808 tmflops = nops /
809 $ ( ( ctime( 1 )+ctime( 2 ) ) * 1.0d+6 )
810 ELSE
811 tmflops = 0.0d+0
812 END IF
813*
814 IF( ctime( 2 ).GE.0.0d+0 )
815 $ WRITE( nout, fmt = 9993 ) 'CPU ', uplo, n,
816 $ nb, nrhs, nbrhs, nprow, npcol,
817 $ ctime( 1 ), ctime( 2 ), tmflops,
818 $ passed
819*
820 END IF
821 10 CONTINUE
822 20 CONTINUE
823*
824 IF( check .AND. sresid.GT.thresh ) THEN
825*
826* Compute FRESID = ||A - LL'|| / (||A|| * N * eps)
827*
828 CALL pspotrrv( uplo, n, mem( ipa ), 1, 1, desca,
829 $ mem( ipw ) )
830 CALL pslafchk( 'Symm', 'Diag', n, n, mem( ipa ), 1, 1,
831 $ desca, iaseed, anorm, fresid,
832 $ mem( ipw ) )
833*
834* Check for memory overwrite
835*
836 CALL pschekpad( ictxt, 'PSPOTRRV', np, nq,
837 $ mem( ipa-iprepad ), desca( lld_ ),
838 $ iprepad, ipostpad, padval )
839 CALL pschekpad( ictxt, 'PSGETRRV',
840 $ worksiz-ipostpad, 1,
841 $ mem( ipw-iprepad ), worksiz-ipostpad,
842 $ iprepad, ipostpad, padval )
843*
844 IF( iam.EQ.0 ) THEN
845 IF( lsame( uplo, 'L' ) ) THEN
846 WRITE( nout, fmt = 9986 ) 'L*L''', fresid
847 ELSE
848 WRITE( nout, fmt = 9986 ) 'U''*U', fresid
849 END IF
850 END IF
851 END IF
852*
853 30 CONTINUE
854 40 CONTINUE
855 CALL blacs_gridexit( ictxt )
856 50 CONTINUE
857*
858* Print ending messages and close output file
859*
860 60 CONTINUE
861 IF( iam.EQ.0 ) THEN
862 ktests = kpass + kfail + kskip
863 WRITE( nout, fmt = * )
864 WRITE( nout, fmt = 9992 ) ktests
865 IF( check ) THEN
866 WRITE( nout, fmt = 9991 ) kpass
867 WRITE( nout, fmt = 9989 ) kfail
868 ELSE
869 WRITE( nout, fmt = 9990 ) kpass
870 END IF
871 WRITE( nout, fmt = 9988 ) kskip
872 WRITE( nout, fmt = * )
873 WRITE( nout, fmt = * )
874 WRITE( nout, fmt = 9987 )
875 IF( nout.NE.6 .AND. nout.NE.0 )
876 $ CLOSE ( nout )
877 END IF
878*
879 CALL blacs_exit( 0 )
880*
881 9999 FORMAT( 'ILLEGAL ', a6, ': ', a5, ' = ', i3,
882 $ '; It should be at least 1' )
883 9998 FORMAT( 'ILLEGAL GRID: nprow*npcol = ', i4, '. It can be at most',
884 $ i4 )
885 9997 FORMAT( 'Bad ', a6, ' parameters: going on to next test case.' )
886 9996 FORMAT( 'Unable to perform ', a, ': need TOTMEM of at least',
887 $ i11 )
888 9995 FORMAT( 'TIME UPLO N NB NRHS NBRHS P Q LLt Time ',
889 $ 'Slv Time MFLOPS CHECK' )
890 9994 FORMAT( '---- ---- ----- --- ---- ----- ---- ---- -------- ',
891 $ '-------- -------- ------' )
892 9993 FORMAT( a4, 4x, a1, 1x, i5, 1x, i3, 1x, i4, 1x, i5, 1x, i4, 1x,
893 $ i4, 1x, f8.2, 1x, f8.2, 1x, f8.2, 1x, a6 )
894 9992 FORMAT( 'Finished ', i6, ' tests, with the following results:' )
895 9991 FORMAT( i5, ' tests completed and passed residual checks.' )
896 9990 FORMAT( i5, ' tests completed without checking.' )
897 9989 FORMAT( i5, ' tests completed and failed residual checks.' )
898 9988 FORMAT( i5, ' tests skipped because of illegal input values.' )
899 9987 FORMAT( 'END OF TESTS.' )
900 9986 FORMAT( '||A - ', a4, '|| / (||A|| * N * eps) = ', g25.7 )
901 9985 FORMAT( '||Ax-b||/(||x||*||A||*eps*N) ', f25.7 )
902*
903 stop
904*
905* End of PSLLTDRIVER
906*
907 END
pslafchk
subroutine pslafchk(aform, diag, m, n, a, ia, ja, desca, iaseed, anorm, fresid, work)
Definition pslafchk.f:3
psmatgen
subroutine psmatgen(ictxt, aform, diag, m, n, mb, nb, a, lda, iarow, iacol, iseed, iroff, irnum, icoff, icnum, myrow, mycol, nprow, npcol)
Definition psmatgen.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
pschekpad
subroutine pschekpad(ictxt, mess, m, n, a, lda, ipre, ipost, chkval)
Definition pschekpad.f:3
psfillpad
subroutine psfillpad(ictxt, m, n, a, lda, ipre, ipost, chkval)
Definition psfillpad.f:2
pslansy
real function pslansy(norm, uplo, n, a, ia, ja, desca, work)
Definition pslansy.f:3
pslaschk
subroutine pslaschk(symm, diag, n, nrhs, x, ix, jx, descx, iaseed, ia, ja, desca, ibseed, anorm, resid, work)
Definition pslaschk.f:4
pslltdriver
program pslltdriver
Definition pslltdriver.f:1
pslltinfo
subroutine pslltinfo(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 pslltinfo.f:6
pspocon
subroutine pspocon(uplo, n, a, ia, ja, desca, anorm, rcond, work, lwork, iwork, liwork, info)
Definition pspocon.f:3
psporfs
subroutine psporfs(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 psporfs.f:4
pspotrf
subroutine pspotrf(uplo, n, a, ia, ja, desca, info)
Definition pspotrf.f:2
pspotrrv
subroutine pspotrrv(uplo, n, a, ia, ja, desca, work)
Definition pspotrrv.f:2
pspotrs
subroutine pspotrs(uplo, n, nrhs, a, ia, ja, desca, b, ib, jb, descb, info)
Definition pspotrs.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