OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
pclltdriver.f
Go to the documentation of this file.
1 PROGRAM pclltdriver
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* PCLLTDRIVER is the main test program for the COMPLEX
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* CPLXSZ INTEGER, default = 8 bytes.
56* INTGSZ and CPLXSZ indicate the length in bytes on the
57* given platform for an integer and a single precision
58* complex.
59* MEM COMPLEX array, dimension ( TOTMEM / CPLXSZ )
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 cplxsz, memsiz, ntests, realsz, totmem
75 REAL zero
76 COMPLEX padval
77 parameter( cplxsz = 8, realsz = 4, totmem = 2000000,
78 $ memsiz = totmem / cplxsz, ntests = 20,
79 $ padval = ( -9923.0e+0, -9923.0e+0 ),
80 $ zero = 0.0e+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 anorm, anorm1, fresid, rcond, sresid, sresid2,
95 $ thresh
96 DOUBLE PRECISION nops, 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 mem( memsiz )
105* ..
106* .. External Subroutines ..
107 EXTERNAL blacs_barrier, blacs_exit, blacs_gridexit,
108 $ blacs_gridinfo, blacs_gridinit, descinit,
109 $ igsum2d, blacs_pinfo, pcchekpad, pcfillpad,
110 $ pclafchk, pclaschk, pclltinfo,
111 $ pcmatgen, pcpocon, pcporfs,
112 $ pcpotrf, pcpotrrv, pcpotrs, slboot,
113 $ slcombine, sltimer
114* ..
115* .. External Functions ..
116 LOGICAL lsame
117 INTEGER iceil, ilcm, numroc
118 REAL pclanhe
119 EXTERNAL iceil, ilcm, lsame, numroc, pclanhe
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 pclltinfo( 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 PCLAFCHK, PCPOTRRV, and
290* PCLANHE
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( realsz * itemp, cplxsz ) )
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 )*cplxsz
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 IF( ierr( 1 ).GT.0 ) THEN
327 IF( iam.EQ.0 )
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 pcmatgen( 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 pcfillpad( ictxt, np, nq, mem( ipa-iprepad ),
345 $ desca( lld_ ), iprepad, ipostpad,
346 $ padval )
347 CALL pcfillpad( ictxt, worksiz-ipostpad, 1,
348 $ mem( ipw-iprepad ), worksiz-ipostpad,
349 $ iprepad, ipostpad, padval )
350 anorm = pclanhe( 'I', uplo, n, mem( ipa ), 1, 1,
351 $ desca, mem( ipw ) )
352 anorm1 = pclanhe( '1', uplo, n, mem( ipa ), 1, 1,
353 $ desca, mem( ipw ) )
354 CALL pcchekpad( ictxt, 'PCLANHE', np, nq,
355 $ mem( ipa-iprepad ), desca( lld_ ),
356 $ iprepad, ipostpad, padval )
357 CALL pcchekpad( ictxt, 'PCLANHE', worksiz-ipostpad,
358 $ 1, mem( ipw-iprepad ),
359 $ worksiz-ipostpad, iprepad, ipostpad,
360 $ padval )
361 END IF
362*
363 IF( est ) THEN
364 CALL pcmatgen( 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 pcfillpad( 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 pcpotrf( uplo, n, mem( ipa ), 1, 1, desca, info )
384*
385 CALL sltimer( 1 )
386*
387 IF( info.NE.0 ) THEN
388 IF( iam.EQ.0 )
389 $ WRITE( nout, fmt = * ) 'PCPOTRF 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 pcchekpad( ictxt, 'PCPOTRF', 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 PCPOCON
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*realsz, cplxsz ) + 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 )*cplxsz
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 pcfillpad( ictxt, lwork, 1,
439 $ mem( ipw-iprepad ), lwork,
440 $ iprepad, ipostpad, padval )
441 CALL pcfillpad( 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 pcpocon( 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 pcchekpad( ictxt, 'PCPOCON', np, nq,
455 $ mem( ipa-iprepad ), desca( lld_ ),
456 $ iprepad, ipostpad, padval )
457 CALL pcchekpad( ictxt, 'PCPOCON',
458 $ lwork, 1, mem( ipw-iprepad ),
459 $ lwork, iprepad, ipostpad,
460 $ padval )
461 CALL pcchekpad( ictxt, 'PCPOCON',
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 PCLASCHK
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 )*cplxsz
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 pcmatgen( 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 pcfillpad( ictxt, np, myrhs,
549 $ mem( ipb-iprepad ),
550 $ descb( lld_ ),
551 $ iprepad, ipostpad, padval )
552*
553 IF( est ) THEN
554 CALL pcmatgen( 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 pcfillpad( ictxt, np, myrhs,
564 $ mem( ipb0-iprepad ),
565 $ descb( lld_ ), iprepad,
566 $ ipostpad, padval )
567 CALL pcfillpad( ictxt, 1, myrhs,
568 $ mem( ipferr-iprepad ), 1,
569 $ iprepad, ipostpad,
570 $ padval )
571 CALL pcfillpad( 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 pcpotrs( 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 pcchekpad( ictxt, 'PCPOTRS', np, nq,
594 $ mem( ipa-iprepad ),
595 $ desca( lld_ ),
596 $ iprepad, ipostpad, padval )
597 CALL pcchekpad( ictxt, 'PCPOTRS', np,
598 $ myrhs, mem( ipb-iprepad ),
599 $ descb( lld_ ), iprepad,
600 $ ipostpad, padval )
601*
602 CALL pcfillpad( 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 pclaschk( '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 pcchekpad( ictxt, 'PCLASCHK', np,
620 $ myrhs, mem( ipb-iprepad ),
621 $ descb( lld_ ), iprepad,
622 $ ipostpad, padval )
623 CALL pcchekpad( ictxt, 'PCLASCHK',
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 IF( ( sresid.LE.thresh ).AND.
632 $ ( (sresid-sresid).EQ.0.0e+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 PCPORFS
648*
649 lwork = max( 1, 2*np )
650 ipw2 = ipw + lwork + ipostpad + iprepad
651 lrwork = max( 1, np )
652 lw2 = iceil( lrwork*realsz, cplxsz ) +
653 $ ipostpad
654*
655 ierr( 1 ) = 0
656 IF( ipw2+lw2.GT.memsiz ) THEN
657 IF( iam.EQ.0 )
658 $ WRITE( nout, fmt = 9996 )
659 $ 'iter ref', ( ipw2+lw2 )*cplxsz
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 IF( ierr( 1 ).GT.0 ) THEN
669 IF( iam.EQ.0 )
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 pcfillpad( ictxt, lwork, 1,
678 $ mem( ipw-iprepad ),
679 $ lwork, iprepad, ipostpad,
680 $ padval )
681 CALL pcfillpad( 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 pcporfs( 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 pcchekpad( ictxt, 'PCPORFS', np,
703 $ nq, mem( ipa0-iprepad ),
704 $ desca( lld_ ), iprepad,
705 $ ipostpad, padval )
706 CALL pcchekpad( ictxt, 'PCPORFS', np,
707 $ nq, mem( ipa-iprepad ),
708 $ desca( lld_ ), iprepad,
709 $ ipostpad, padval )
710 CALL pcchekpad( ictxt, 'PCPORFS', np,
711 $ myrhs, mem( ipb-iprepad ),
712 $ descb( lld_ ), iprepad,
713 $ ipostpad, padval )
714 CALL pcchekpad( ictxt, 'PCPORFS', np,
715 $ myrhs,
716 $ mem( ipb0-iprepad ),
717 $ descb( lld_ ), iprepad,
718 $ ipostpad, padval )
719 CALL pcchekpad( ictxt, 'PCPORFS', 1,
720 $ myrhs,
721 $ mem( ipferr-iprepad ), 1,
722 $ iprepad, ipostpad,
723 $ padval )
724 CALL pcchekpad( ictxt, 'PCPORFS', 1,
725 $ myrhs,
726 $ mem( ipberr-iprepad ), 1,
727 $ iprepad, ipostpad,
728 $ padval )
729 CALL pcchekpad( ictxt, 'PCPORFS', lwork,
730 $ 1, mem( ipw-iprepad ),
731 $ lwork, iprepad, ipostpad,
732 $ padval )
733 CALL pcchekpad( ictxt, 'PCPORFS',
734 $ lw2-ipostpad, 1,
735 $ mem( ipw2-iprepad ),
736 $ lw2-ipostpad,
737 $ iprepad, ipostpad,
738 $ padval )
739*
740 CALL pcfillpad( 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 pclaschk( 'Herm', 'Diag', n, nrhs,
748 $ mem( ipb ), 1, 1, descb,
749 $ iaseed, 1, 1, desca,
750 $ ibseed, anorm, sresid2,
751 $ mem( ipw ) )
752*
753 IF( iam.EQ.0 .AND. sresid2.GT.thresh )
754 $ WRITE( nout, fmt = 9985 ) sresid2
755*
756* check for memory overwrite
757*
758 CALL pcchekpad( ictxt, 'PCLASCHK', np,
759 $ myrhs, mem( ipb-iprepad ),
760 $ descb( lld_ ), iprepad,
761 $ ipostpad, padval )
762 CALL pcchekpad( ictxt, 'PCLASCHK',
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 PCPOTRRV( UPLO, N, MEM( IPA ), 1, 1, DESCA,
833 $ MEM( IPW ) )
834 CALL PCLAFCHK( '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 PCCHEKPAD( ICTXT, 'pcpotrrv', NP, NQ,
841 $ MEM( IPA-IPREPAD ), DESCA( LLD_ ),
842 $ IPREPAD, IPOSTPAD, PADVAL )
843 CALL PCCHEKPAD( ICTXT, 'pcgetrrv',
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 PCLLTDRIVER
910*
911 END
pclafchk
subroutine pclafchk(aform, diag, m, n, a, ia, ja, desca, iaseed, anorm, fresid, work)
Definition pclafchk.f:3
pcmatgen
subroutine pcmatgen(ictxt, aform, diag, m, n, mb, nb, a, lda, iarow, iacol, iseed, iroff, irnum, icoff, icnum, myrow, mycol, nprow, npcol)
Definition pcmatgen.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
pcpotrf
subroutine pcpotrf(uplo, n, a, ia, ja, desca, info)
Definition mpi.f:801
pcpotrs
subroutine pcpotrs(uplo, n, nrhs, a, ia, ja, desca, b, ib, jb, descb, info)
Definition mpi.f:1183
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
pcchekpad
subroutine pcchekpad(ictxt, mess, m, n, a, lda, ipre, ipost, chkval)
Definition pcchekpad.f:3
pcfillpad
subroutine pcfillpad(ictxt, m, n, a, lda, ipre, ipost, chkval)
Definition pcfillpad.f:2
pcgetrrv
subroutine pcgetrrv(m, n, a, ia, ja, desca, ipiv, work)
Definition pcgetrrv.f:2
pclanhe
real function pclanhe(norm, uplo, n, a, ia, ja, desca, work)
Definition pclanhe.f:3
pclaschk
subroutine pclaschk(symm, diag, n, nrhs, x, ix, jx, descx, iaseed, ia, ja, desca, ibseed, anorm, resid, work)
Definition pclaschk.f:4
pclltdriver
program pclltdriver
Definition pclltdriver.f:1
pclltinfo
subroutine pclltinfo(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 pclltinfo.f:6
pcpocon
subroutine pcpocon(uplo, n, a, ia, ja, desca, anorm, rcond, work, lwork, rwork, lrwork, info)
Definition pcpocon.f:3
pcporfs
subroutine pcporfs(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 pcporfs.f:4
pcpotrrv
subroutine pcpotrrv(uplo, n, a, ia, ja, desca, work)
Definition pcpotrrv.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