ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
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
max
#define max(A, B)
Definition: pcgemr.c:180
pdlansy
double precision function pdlansy(NORM, UPLO, N, A, IA, JA, DESCA, WORK)
Definition: pdlansy.f:3
pdchekpad
subroutine pdchekpad(ICTXT, MESS, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: pdchekpad.f:3
ilcm
integer function ilcm(M, N)
Definition: ilcm.f:2
pdlltdriver
program pdlltdriver
Definition: pdlltdriver.f:1
sltimer
subroutine sltimer(I)
Definition: sltimer.f:47
pdlaschk
subroutine pdlaschk(SYMM, DIAG, N, NRHS, X, IX, JX, DESCX, IASEED, IA, JA, DESCA, IBSEED, ANORM, RESID, WORK)
Definition: pdlaschk.f:4
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
lsame
logical function lsame(CA, CB)
Definition: tools.f:1724
descinit
subroutine descinit(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD, INFO)
Definition: descinit.f:3
slboot
subroutine slboot()
Definition: sltimer.f:2
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
pdpocon
subroutine pdpocon(UPLO, N, A, IA, JA, DESCA, ANORM, RCOND, WORK, LWORK, IWORK, LIWORK, INFO)
Definition: pdpocon.f:3
pdlafchk
subroutine pdlafchk(AFORM, DIAG, M, N, A, IA, JA, DESCA, IASEED, ANORM, FRESID, WORK)
Definition: pdlafchk.f:3
pdpotrs
subroutine pdpotrs(UPLO, N, NRHS, A, IA, JA, DESCA, B, IB, JB, DESCB, INFO)
Definition: pdpotrs.f:3
numroc
integer function numroc(N, NB, IPROC, ISRCPROC, NPROCS)
Definition: numroc.f:2
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
pdfillpad
subroutine pdfillpad(ICTXT, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: pdfillpad.f:2
pdpotrrv
subroutine pdpotrrv(UPLO, N, A, IA, JA, DESCA, WORK)
Definition: pdpotrrv.f:2
pdpotrf
subroutine pdpotrf(UPLO, N, A, IA, JA, DESCA, INFO)
Definition: pdpotrf.f:2
slcombine
subroutine slcombine(ICTXT, SCOPE, OP, TIMETYPE, N, IBEG, TIMES)
Definition: sltimer.f:267
min
#define min(A, B)
Definition: pcgemr.c:181
iceil
integer function iceil(INUM, IDENOM)
Definition: iceil.f:2