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