ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
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  IF( myrow.EQ.0 .AND. mycol.EQ.0 ) 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  IF( wtime( 1 ) + wtime( 2 ) .GT. 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  IF( wtime( 2 ).GE.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  IF( ctime( 1 )+ctime( 2 ).GT.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  IF( ctime( 2 ).GE.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  IF( check .AND. sresid.GT.thresh ) 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  IF( iam.EQ.0 ) 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  IF( iam.EQ.0 ) 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  IF( nout.NE.6 .AND. nout.NE.0 )
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
max
#define max(A, B)
Definition: pcgemr.c:180
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
ilcm
integer function ilcm(M, N)
Definition: ilcm.f:2
pclanhe
real function pclanhe(NORM, UPLO, N, A, IA, JA, DESCA, WORK)
Definition: pclanhe.f:3
sltimer
subroutine sltimer(I)
Definition: sltimer.f:47
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
pcpotrf
subroutine pcpotrf(UPLO, N, A, IA, JA, DESCA, INFO)
Definition: pcpotrf.f:2
lsame
logical function lsame(CA, CB)
Definition: tools.f:1724
pclltdriver
program pclltdriver
Definition: pclltdriver.f:1
pcpocon
subroutine pcpocon(UPLO, N, A, IA, JA, DESCA, ANORM, RCOND, WORK, LWORK, RWORK, LRWORK, INFO)
Definition: pcpocon.f:3
pcchekpad
subroutine pcchekpad(ICTXT, MESS, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: pcchekpad.f:3
pcpotrrv
subroutine pcpotrrv(UPLO, N, A, IA, JA, DESCA, WORK)
Definition: pcpotrrv.f:2
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
descinit
subroutine descinit(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD, INFO)
Definition: descinit.f:3
slboot
subroutine slboot()
Definition: sltimer.f:2
pcfillpad
subroutine pcfillpad(ICTXT, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: pcfillpad.f:2
numroc
integer function numroc(N, NB, IPROC, ISRCPROC, NPROCS)
Definition: numroc.f:2
pclaschk
subroutine pclaschk(SYMM, DIAG, N, NRHS, X, IX, JX, DESCX, IASEED, IA, JA, DESCA, IBSEED, ANORM, RESID, WORK)
Definition: pclaschk.f:4
slcombine
subroutine slcombine(ICTXT, SCOPE, OP, TIMETYPE, N, IBEG, TIMES)
Definition: sltimer.f:267
min
#define min(A, B)
Definition: pcgemr.c:181
pcpotrs
subroutine pcpotrs(UPLO, N, NRHS, A, IA, JA, DESCA, B, IB, JB, DESCB, INFO)
Definition: pcpotrs.f:3
iceil
integer function iceil(INUM, IDENOM)
Definition: iceil.f:2