ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pslltdriver.f
Go to the documentation of this file.
1  PROGRAM pslltdriver
2 *
3 * -- ScaLAPACK testing driver (version 1.7) --
4 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
5 * and University of California, Berkeley.
6 * May 1, 1997
7 *
8 * Purpose
9 * =======
10 *
11 * PSLLTDRIVER is the main test program for the REAL
12 * ScaLAPACK Cholesky routines. This test driver performs an
13 * A = L*L**T or A = U**T*U factorization and solve, and optionally
14 * performs condition estimation and iterative refinement.
15 *
16 * The program must be driven by a short data file. An annotated
17 * example of a data file can be obtained by deleting the first 3
18 * characters from the following 18 lines:
19 * 'ScaLAPACK LLt factorization input file'
20 * 'Intel iPSC/860 hypercube, gamma model.'
21 * 'LLT.out' output file name (if any)
22 * 6 device out
23 * 'U' define Lower or Upper
24 * 1 number of problems sizes
25 * 31 100 200 values of N
26 * 1 number of NB's
27 * 2 10 24 values of NB
28 * 1 number of NRHS's
29 * 1 values of NRHS
30 * 1 Number of NBRHS's
31 * 1 values of NBRHS
32 * 1 number of process grids (ordered pairs of P & Q)
33 * 2 values of P
34 * 2 values of Q
35 * 1.0 threshold
36 * T (T or F) Test Cond. Est. and Iter. Ref. Routines
37 *
38 *
39 * Internal Parameters
40 * ===================
41 *
42 * TOTMEM INTEGER, default = 2000000
43 * TOTMEM is a machine-specific parameter indicating the
44 * maximum amount of available memory in bytes.
45 * The user should customize TOTMEM to his platform. Remember
46 * to leave room in memory for the operating system, the BLACS
47 * buffer, etc. For example, on a system with 8 MB of memory
48 * per process (e.g., one processor on an Intel iPSC/860), the
49 * parameters we use are TOTMEM=6200000 (leaving 1.8 MB for OS,
50 * code, BLACS buffer, etc). However, for PVM, we usually set
51 * TOTMEM = 2000000. Some experimenting with the maximum value
52 * of TOTMEM may be required.
53 *
54 * INTGSZ INTEGER, default = 4 bytes.
55 * REALSZ INTEGER, default = 4 bytes.
56 * INTGSZ and REALSZ indicate the length in bytes on the
57 * given platform for an integer and a single precision real.
58 * MEM REAL array, dimension ( TOTMEM / REALSZ )
59 *
60 * All arrays used by SCALAPACK routines are allocated from
61 * this array and referenced by pointers. The integer IPA,
62 * for example, is a pointer to the starting element of MEM for
63 * the matrix A.
64 *
65 * =====================================================================
66 *
67 * .. Parameters ..
68  INTEGER block_cyclic_2d, csrc_, ctxt_, dlen_, dtype_,
69  $ lld_, mb_, m_, nb_, n_, rsrc_
70  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
71  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
72  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
73  INTEGER intgsz, memsiz, ntests, realsz, totmem
74  REAL padval, zero
75  parameter( intgsz = 4, realsz = 4, totmem = 2000000,
76  $ memsiz = totmem / realsz, ntests = 20,
77  $ padval = -9923.0e+0, zero = 0.0e+0 )
78 * ..
79 * .. Local Scalars ..
80  LOGICAL check, est
81  CHARACTER uplo
82  CHARACTER*6 passed
83  CHARACTER*80 outfile
84  INTEGER hh, i, iam, iaseed, ibseed, ictxt, imidpad,
85  $ info, ipa, ipa0, ipb, ipb0, ipberr, ipferr,
86  $ iprepad, ipostpad, ipw, ipw2, itemp, j, k,
87  $ kfail, kk, kpass, kskip, ktests, lcm, lcmq,
88  $ liwork, lwork, lw2, mycol, myrhs, myrow, n, nb,
89  $ nbrhs, ngrids, nmat, nnb, nnbr, nnr, nout, np,
90  $ npcol, nprocs, nprow, nq, nrhs, worksiz
91  REAL anorm, anorm1, fresid, rcond, sresid, sresid2,
92  $ thresh
93  DOUBLE PRECISION nops, tmflops
94 * ..
95 * .. Local Arrays ..
96  INTEGER desca( dlen_ ), descb( dlen_ ), ierr( 1 ),
97  $ nbrval( ntests ), nbval( ntests ),
98  $ nrval( ntests ), nval( ntests ),
99  $ pval( ntests ), qval( ntests )
100  REAL mem( memsiz )
101  DOUBLE PRECISION ctime( 2 ), wtime( 2 )
102 * ..
103 * .. External Subroutines ..
104  EXTERNAL blacs_barrier, blacs_exit, blacs_gridexit,
105  $ blacs_gridinfo, blacs_gridinit, descinit,
106  $ igsum2d, blacs_pinfo, pschekpad, psfillpad,
107  $ pslafchk, pslaschk, pslltinfo,
108  $ psmatgen, pspocon, psporfs,
109  $ pspotrf, pspotrrv, pspotrs, slboot,
110  $ slcombine, sltimer
111 * ..
112 * .. External Functions ..
113  LOGICAL lsame
114  INTEGER iceil, ilcm, numroc
115  REAL pslansy
116  EXTERNAL iceil, ilcm, lsame, numroc, pslansy
117 * ..
118 * .. Intrinsic Functions ..
119  INTRINSIC dble, max, min
120 * ..
121 * .. Data Statements ..
122  DATA kfail, kpass, kskip, ktests / 4*0 /
123 * ..
124 * .. Executable Statements ..
125 *
126 * Get starting information
127 *
128  CALL blacs_pinfo( iam, nprocs )
129  iaseed = 100
130  ibseed = 200
131  CALL pslltinfo( outfile, nout, uplo, nmat, nval, ntests, nnb,
132  $ nbval, ntests, nnr, nrval, ntests, nnbr, nbrval,
133  $ ntests, ngrids, pval, ntests, qval, ntests,
134  $ thresh, est, mem, iam, nprocs )
135  check = ( thresh.GE.0.0e+0 )
136 *
137 * Print headings
138 *
139  IF( iam.EQ.0 ) THEN
140  WRITE( nout, fmt = * )
141  WRITE( nout, fmt = 9995 )
142  WRITE( nout, fmt = 9994 )
143  WRITE( nout, fmt = * )
144  END IF
145 *
146 * Loop over different process grids
147 *
148  DO 50 i = 1, ngrids
149 *
150  nprow = pval( i )
151  npcol = qval( i )
152 *
153 * Make sure grid information is correct
154 *
155  ierr( 1 ) = 0
156  IF( nprow.LT.1 ) THEN
157  IF( iam.EQ.0 )
158  $ WRITE( nout, fmt = 9999 ) 'GRID', 'nprow', nprow
159  ierr( 1 ) = 1
160  ELSE IF( npcol.LT.1 ) THEN
161  IF( iam.EQ.0 )
162  $ WRITE( nout, fmt = 9999 ) 'GRID', 'npcol', npcol
163  ierr( 1 ) = 1
164  ELSE IF( nprow*npcol.GT.nprocs ) THEN
165  IF( iam.EQ.0 )
166  $ WRITE( nout, fmt = 9998 ) nprow*npcol, nprocs
167  ierr( 1 ) = 1
168  END IF
169 *
170  IF( ierr( 1 ).GT.0 ) THEN
171  IF( iam.EQ.0 )
172  $ WRITE( nout, fmt = 9997 ) 'grid'
173  kskip = kskip + 1
174  GO TO 50
175  END IF
176 *
177 * Define process grid
178 *
179  CALL blacs_get( -1, 0, ictxt )
180  CALL blacs_gridinit( ictxt, 'Row-major', nprow, npcol )
181  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
182 *
183 * Go to bottom of process grid loop if this case doesn't use my
184 * process
185 *
186  IF( myrow.GE.nprow .OR. mycol.GE.npcol )
187  $ GO TO 50
188 *
189  DO 40 j = 1, nmat
190 *
191  n = nval( j )
192 *
193 * Make sure matrix information is correct
194 *
195  ierr( 1 ) = 0
196  IF( n.LT.1 ) THEN
197  IF( iam.EQ.0 )
198  $ WRITE( nout, fmt = 9999 ) 'MATRIX', 'N', n
199  ierr( 1 ) = 1
200  ELSE IF( n.LT.1 ) THEN
201  IF( iam.EQ.0 )
202  $ WRITE( nout, fmt = 9999 ) 'MATRIX', 'N', n
203  ierr( 1 ) = 1
204  END IF
205 *
206 * Check all processes for an error
207 *
208  CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1, 0 )
209 *
210  IF( ierr( 1 ).GT.0 ) THEN
211  IF( iam.EQ.0 )
212  $ WRITE( nout, fmt = 9997 ) 'matrix'
213  kskip = kskip + 1
214  GO TO 40
215  END IF
216 *
217  DO 30 k = 1, nnb
218 *
219  nb = nbval( k )
220 *
221 * Make sure nb is legal
222 *
223  ierr( 1 ) = 0
224  IF( nb.LT.1 ) THEN
225  ierr( 1 ) = 1
226  IF( iam.EQ.0 )
227  $ WRITE( nout, fmt = 9999 ) 'NB', 'NB', nb
228  END IF
229 *
230 * Check all processes for an error
231 *
232  CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1, 0 )
233 *
234  IF( ierr( 1 ).GT.0 ) THEN
235  IF( iam.EQ.0 )
236  $ WRITE( nout, fmt = 9997 ) 'NB'
237  kskip = kskip + 1
238  GO TO 30
239  END IF
240 *
241 * Padding constants
242 *
243  np = numroc( n, nb, myrow, 0, nprow )
244  nq = numroc( n, nb, mycol, 0, npcol )
245  IF( check ) THEN
246  iprepad = max( nb, np )
247  imidpad = nb
248  ipostpad = max( nb, nq )
249  ELSE
250  iprepad = 0
251  imidpad = 0
252  ipostpad = 0
253  END IF
254 *
255 * Initialize the array descriptor for the matrix A
256 *
257  CALL descinit( desca, n, n, nb, nb, 0, 0, ictxt,
258  $ max( 1, np )+imidpad, ierr( 1 ) )
259 *
260 * Check all processes for an error
261 *
262  CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1, 0 )
263 *
264  IF( ierr( 1 ).LT.0 ) THEN
265  IF( iam.EQ.0 )
266  $ WRITE( nout, fmt = 9997 ) 'descriptor'
267  kskip = kskip + 1
268  GO TO 30
269  END IF
270 *
271 * Assign pointers into MEM for SCALAPACK arrays, A is
272 * allocated starting at position MEM( IPREPAD+1 )
273 *
274  ipa = iprepad+1
275  IF( est ) THEN
276  ipa0 = ipa + desca( lld_ )*nq + ipostpad + iprepad
277  ipw = ipa0 + desca( lld_ )*nq + ipostpad + iprepad
278  ELSE
279  ipw = ipa + desca( lld_ )*nq + ipostpad + iprepad
280  END IF
281 *
282 *
283  IF( check ) THEN
284 *
285 * Calculate the amount of workspace required by
286 * the checking routines PSLAFCHK, PSPOTRRV, and
287 * PSLANSY
288 *
289  worksiz = np * desca( nb_ )
290 *
291  worksiz = max( worksiz, desca( mb_ ) * desca( nb_ ) )
292 *
293  lcm = ilcm( nprow, npcol )
294  itemp = max( 2, 2 * nq ) + np
295  IF( nprow.NE.npcol ) THEN
296  itemp = itemp +
297  $ nb * iceil( iceil( np, nb ), lcm / nprow )
298  END IF
299  worksiz = max( worksiz, itemp )
300  worksiz = worksiz + ipostpad
301 *
302  ELSE
303 *
304  worksiz = ipostpad
305 *
306  END IF
307 *
308 * Check for adequate memory for problem size
309 *
310  ierr( 1 ) = 0
311  IF( ipw+worksiz.GT.memsiz ) THEN
312  IF( iam.EQ.0 )
313  $ WRITE( nout, fmt = 9996 ) 'factorization',
314  $ ( ipw+worksiz )*realsz
315  ierr( 1 ) = 1
316  END IF
317 *
318 * Check all processes for an error
319 *
320  CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1, 0 )
321 *
322  IF( ierr( 1 ).GT.0 ) THEN
323  IF( iam.EQ.0 )
324  $ WRITE( nout, fmt = 9997 ) 'MEMORY'
325  kskip = kskip + 1
326  GO TO 30
327  END IF
328 *
329 * Generate a symmetric positive definite matrix A
330 *
331  CALL psmatgen( ictxt, 'Symm', 'Diag', desca( m_ ),
332  $ desca( n_ ), desca( mb_ ), desca( nb_ ),
333  $ mem( ipa ), desca( lld_ ), desca( rsrc_ ),
334  $ desca( csrc_ ), iaseed, 0, np, 0, nq,
335  $ myrow, mycol, nprow, npcol )
336 *
337 * Calculate inf-norm of A for residual error-checking
338 *
339  IF( check ) THEN
340  CALL psfillpad( ictxt, np, nq, mem( ipa-iprepad ),
341  $ desca( lld_ ), iprepad, ipostpad,
342  $ padval )
343  CALL psfillpad( ictxt, worksiz-ipostpad, 1,
344  $ mem( ipw-iprepad ), worksiz-ipostpad,
345  $ iprepad, ipostpad, padval )
346  anorm = pslansy( 'I', uplo, n, mem( ipa ), 1, 1,
347  $ desca, mem( ipw ) )
348  anorm1 = pslansy( '1', uplo, n, mem( ipa ), 1, 1,
349  $ desca, mem( ipw ) )
350  CALL pschekpad( ictxt, 'PSLANSY', np, nq,
351  $ mem( ipa-iprepad ), desca( lld_ ),
352  $ iprepad, ipostpad, padval )
353  CALL pschekpad( ictxt, 'PSLANSY',
354  $ worksiz-ipostpad, 1,
355  $ mem( ipw-iprepad ), worksiz-ipostpad,
356  $ iprepad, ipostpad, padval )
357  END IF
358 *
359  IF( est ) THEN
360  CALL psmatgen( ictxt, 'Symm', 'Diag', desca( m_ ),
361  $ desca( n_ ), desca( mb_ ),
362  $ desca( nb_ ), mem( ipa0 ),
363  $ desca( lld_ ), desca( rsrc_ ),
364  $ desca( csrc_ ), iaseed, 0, np, 0, nq,
365  $ myrow, mycol, nprow, npcol )
366  IF( check )
367  $ CALL psfillpad( ictxt, np, nq,
368  $ mem( ipa0-iprepad ), desca( lld_ ),
369  $ iprepad, ipostpad, padval )
370  END IF
371 *
372  CALL slboot()
373  CALL blacs_barrier( ictxt, 'All' )
374 *
375 * Perform LLt factorization
376 *
377  CALL sltimer( 1 )
378 *
379  CALL pspotrf( uplo, n, mem( ipa ), 1, 1, desca, info )
380 *
381  CALL sltimer( 1 )
382 *
383  IF( info.NE.0 ) THEN
384  IF( iam.EQ.0 )
385  $ WRITE( nout, fmt = * ) 'PSPOTRF INFO=', info
386  kfail = kfail + 1
387  rcond = zero
388  GO TO 60
389  END IF
390 *
391  IF( check ) THEN
392 *
393 * Check for memory overwrite in LLt factorization
394 *
395  CALL pschekpad( ictxt, 'PSPOTRF', np, nq,
396  $ mem( ipa-iprepad ), desca( lld_ ),
397  $ iprepad, ipostpad, padval )
398  END IF
399 *
400  IF( est ) THEN
401 *
402 * Calculate workspace required for PSPOCON
403 *
404  lwork = max( 1, 2*np ) + max( 1, 2*nq ) +
405  $ max( 2, desca( nb_ )*
406  $ max( 1, iceil( nprow-1, npcol ) ),
407  $ nq + desca( nb_ )*
408  $ max( 1, iceil( npcol-1, nprow ) ) )
409  ipw2 = ipw + lwork + ipostpad + iprepad
410  liwork = max( 1, np )
411  lw2 = iceil( liwork*intgsz, realsz ) + ipostpad
412 *
413  ierr( 1 ) = 0
414  IF( ipw2+lw2.GT.memsiz ) THEN
415  IF( iam.EQ.0 )
416  $ WRITE( nout, fmt = 9996 )'cond est',
417  $ ( ipw2+lw2 )*realsz
418  ierr( 1 ) = 1
419  END IF
420 *
421 * Check all processes for an error
422 *
423  CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1,
424  $ -1, 0 )
425 *
426  IF( ierr( 1 ).GT.0 ) THEN
427  IF( iam.EQ.0 )
428  $ WRITE( nout, fmt = 9997 ) 'MEMORY'
429  kskip = kskip + 1
430  GO TO 60
431  END IF
432 *
433  IF( check ) THEN
434  CALL psfillpad( ictxt, lwork, 1,
435  $ mem( ipw-iprepad ), lwork,
436  $ iprepad, ipostpad, padval )
437  CALL psfillpad( ictxt, lw2-ipostpad, 1,
438  $ mem( ipw2-iprepad ),
439  $ lw2-ipostpad, iprepad,
440  $ ipostpad, padval )
441  END IF
442 *
443 * Compute condition number of the matrix
444 *
445  CALL pspocon( uplo, n, mem( ipa ), 1, 1, desca,
446  $ anorm1, rcond, mem( ipw ), lwork,
447  $ mem( ipw2 ), liwork, info )
448 *
449  IF( check ) THEN
450  CALL pschekpad( ictxt, 'PSPOCON', np, nq,
451  $ mem( ipa-iprepad ), desca( lld_ ),
452  $ iprepad, ipostpad, padval )
453  CALL pschekpad( ictxt, 'PSPOCON',
454  $ lwork, 1, mem( ipw-iprepad ),
455  $ lwork, iprepad, ipostpad,
456  $ padval )
457  CALL pschekpad( ictxt, 'PSPOCON',
458  $ lw2-ipostpad, 1,
459  $ mem( ipw2-iprepad ), lw2-ipostpad,
460  $ iprepad, ipostpad, padval )
461  END IF
462  END IF
463 *
464 * Loop over the different values for NRHS
465 *
466  DO 20 hh = 1, nnr
467 *
468  nrhs = nrval( hh )
469 *
470  DO 10 kk = 1, nnbr
471 *
472  nbrhs = nbrval( kk )
473 *
474 * Initialize Array Descriptor for rhs
475 *
476  CALL descinit( descb, n, nrhs, nb, nbrhs, 0, 0,
477  $ ictxt, max( 1, np )+imidpad,
478  $ ierr( 1 ) )
479 *
480 * move IPW to allow room for RHS
481 *
482  myrhs = numroc( descb( n_ ), descb( nb_ ), mycol,
483  $ descb( csrc_ ), npcol )
484  ipb = ipw
485 *
486  IF( est ) THEN
487  ipb0 = ipb + descb( lld_ )*myrhs + ipostpad +
488  $ iprepad
489  ipferr = ipb0 + descb( lld_ )*myrhs + ipostpad
490  $ + iprepad
491  ipberr = myrhs + ipferr + ipostpad + iprepad
492  ipw = myrhs + ipberr + ipostpad + iprepad
493  ELSE
494  ipw = ipb + descb( lld_ )*myrhs + ipostpad +
495  $ iprepad
496  END IF
497 *
498  IF( check ) THEN
499 *
500 * Calculate the amount of workspace required by
501 * the checking routines PSLASCHK
502 *
503  lcmq = lcm / npcol
504  worksiz = max( worksiz-ipostpad,
505  $ nq * nbrhs + np * nbrhs +
506  $ max( max( nq*nb, 2*nbrhs ),
507  $ nbrhs * numroc( numroc(n,nb,0,0,npcol),nb,
508  $ 0,0,lcmq ) ) )
509  worksiz = ipostpad + worksiz
510  ELSE
511  worksiz = ipostpad
512  END IF
513 *
514  ierr( 1 ) = 0
515  IF( ipw+worksiz.GT.memsiz ) THEN
516  IF( iam.EQ.0 )
517  $ WRITE( nout, fmt = 9996 )'solve',
518  $ ( ipw+worksiz )*realsz
519  ierr( 1 ) = 1
520  END IF
521 *
522 * Check all processes for an error
523 *
524  CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1,
525  $ -1, 0 )
526 *
527  IF( ierr( 1 ).GT.0 ) THEN
528  IF( iam.EQ.0 )
529  $ WRITE( nout, fmt = 9997 ) 'MEMORY'
530  kskip = kskip + 1
531  GO TO 10
532  END IF
533 *
534 * Generate RHS
535 *
536  CALL psmatgen( ictxt, 'No', 'No', descb( m_ ),
537  $ descb( n_ ), descb( mb_ ),
538  $ descb( nb_ ), mem( ipb ),
539  $ descb( lld_ ), descb( rsrc_ ),
540  $ descb( csrc_ ), ibseed, 0, np, 0,
541  $ myrhs, myrow, mycol, nprow, npcol )
542 *
543  IF( check )
544  $ CALL psfillpad( ictxt, np, myrhs,
545  $ mem( ipb-iprepad ),
546  $ descb( lld_ ),
547  $ iprepad, ipostpad, padval )
548 *
549  IF( est ) THEN
550  CALL psmatgen( ictxt, 'No', 'No', descb( m_ ),
551  $ descb( n_ ), descb( mb_ ),
552  $ descb( nb_ ), mem( ipb0 ),
553  $ descb( lld_ ), descb( rsrc_ ),
554  $ descb( csrc_ ), ibseed, 0, np, 0,
555  $ myrhs, myrow, mycol, nprow,
556  $ npcol )
557 *
558  IF( check ) THEN
559  CALL psfillpad( ictxt, np, myrhs,
560  $ mem( ipb0-iprepad ),
561  $ descb( lld_ ), iprepad,
562  $ ipostpad, padval )
563  CALL psfillpad( ictxt, 1, myrhs,
564  $ mem( ipferr-iprepad ), 1,
565  $ iprepad, ipostpad,
566  $ padval )
567  CALL psfillpad( ictxt, 1, myrhs,
568  $ mem( ipberr-iprepad ), 1,
569  $ iprepad, ipostpad,
570  $ padval )
571  END IF
572  END IF
573 *
574  CALL blacs_barrier( ictxt, 'All' )
575  CALL sltimer( 2 )
576 *
577 * Solve linear system via Cholesky factorization
578 *
579  CALL pspotrs( uplo, n, nrhs, mem( ipa ), 1, 1,
580  $ desca, mem( ipb ), 1, 1, descb,
581  $ info )
582 *
583  CALL sltimer( 2 )
584 *
585  IF( check ) THEN
586 *
587 * check for memory overwrite
588 *
589  CALL pschekpad( ictxt, 'PSPOTRS', np, nq,
590  $ mem( ipa-iprepad ),
591  $ desca( lld_ ),
592  $ iprepad, ipostpad, padval )
593  CALL pschekpad( ictxt, 'PSPOTRS', np,
594  $ myrhs, mem( ipb-iprepad ),
595  $ descb( lld_ ), iprepad,
596  $ ipostpad, padval )
597 *
598  CALL psfillpad( ictxt, worksiz-ipostpad, 1,
599  $ mem( ipw-iprepad ),
600  $ worksiz-ipostpad, iprepad,
601  $ ipostpad, padval )
602 *
603 * check the solution to rhs
604 *
605  CALL pslaschk( 'Symm', 'Diag', n, nrhs,
606  $ mem( ipb ), 1, 1, descb,
607  $ iaseed, 1, 1, desca, ibseed,
608  $ anorm, sresid, mem( ipw ) )
609 *
610  IF( iam.EQ.0 .AND. sresid.GT.thresh )
611  $ WRITE( nout, fmt = 9985 ) sresid
612 *
613 * check for memory overwrite
614 *
615  CALL pschekpad( ictxt, 'PSLASCHK', np,
616  $ myrhs, mem( ipb-iprepad ),
617  $ descb( lld_ ), iprepad,
618  $ ipostpad, padval )
619  CALL pschekpad( ictxt, 'PSLASCHK',
620  $ worksiz-ipostpad, 1,
621  $ mem( ipw-iprepad ),
622  $ worksiz-ipostpad, iprepad,
623  $ ipostpad, padval )
624 *
625 * The second test is a NaN trap
626 *
627  IF( ( sresid.LE.thresh ).AND.
628  $ ( (sresid-sresid).EQ.0.0e+0 ) ) THEN
629  kpass = kpass + 1
630  passed = 'PASSED'
631  ELSE
632  kfail = kfail + 1
633  passed = 'FAILED'
634  END IF
635  ELSE
636  kpass = kpass + 1
637  sresid = sresid - sresid
638  passed = 'BYPASS'
639  END IF
640 *
641  IF( est ) THEN
642 *
643 * Calculate workspace required for PSPORFS
644 *
645  lwork = max( 1, 3*np )
646  ipw2 = ipw + lwork + ipostpad + iprepad
647  liwork = max( 1, np )
648  lw2 = iceil( liwork*intgsz, realsz ) +
649  $ ipostpad
650 *
651  ierr( 1 ) = 0
652  IF( ipw2+lw2.GT.memsiz ) THEN
653  IF( iam.EQ.0 )
654  $ WRITE( nout, fmt = 9996 )
655  $ 'iter ref', ( ipw2+lw2 )*realsz
656  ierr( 1 ) = 1
657  END IF
658 *
659 * Check all processes for an error
660 *
661  CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr,
662  $ 1, -1, 0 )
663 *
664  IF( ierr( 1 ).GT.0 ) THEN
665  IF( iam.EQ.0 )
666  $ WRITE( nout, fmt = 9997 )
667  $ 'MEMORY'
668  kskip = kskip + 1
669  GO TO 10
670  END IF
671 *
672  IF( check ) THEN
673  CALL psfillpad( ictxt, lwork, 1,
674  $ mem( ipw-iprepad ),
675  $ lwork, iprepad, ipostpad,
676  $ padval )
677  CALL psfillpad( ictxt, lw2-ipostpad,
678  $ 1, mem( ipw2-iprepad ),
679  $ lw2-ipostpad,
680  $ iprepad, ipostpad,
681  $ padval )
682  END IF
683 *
684 * Use iterative refinement to improve the
685 * computed solution
686 *
687  CALL psporfs( uplo, n, nrhs, mem( ipa0 ),
688  $ 1, 1, desca, mem( ipa ), 1, 1,
689  $ desca, mem( ipb0 ), 1, 1,
690  $ descb, mem( ipb ), 1, 1, descb,
691  $ mem( ipferr ), mem( ipberr ),
692  $ mem( ipw ), lwork, mem( ipw2 ),
693  $ liwork, info )
694 *
695 * check for memory overwrite
696 *
697  IF( check ) THEN
698  CALL pschekpad( ictxt, 'PSPORFS', np,
699  $ nq, mem( ipa0-iprepad ),
700  $ desca( lld_ ), iprepad,
701  $ ipostpad, padval )
702  CALL pschekpad( ictxt, 'PSPORFS', np,
703  $ nq, mem( ipa-iprepad ),
704  $ desca( lld_ ), iprepad,
705  $ ipostpad, padval )
706  CALL pschekpad( ictxt, 'PSPORFS', np,
707  $ myrhs, mem( ipb-iprepad ),
708  $ descb( lld_ ), iprepad,
709  $ ipostpad, padval )
710  CALL pschekpad( ictxt, 'PSPORFS', np,
711  $ myrhs,
712  $ mem( ipb0-iprepad ),
713  $ descb( lld_ ), iprepad,
714  $ ipostpad, padval )
715  CALL pschekpad( ictxt, 'PSPORFS', 1,
716  $ myrhs,
717  $ mem( ipferr-iprepad ), 1,
718  $ iprepad, ipostpad,
719  $ padval )
720  CALL pschekpad( ictxt, 'PSPORFS', 1,
721  $ myrhs,
722  $ mem( ipberr-iprepad ), 1,
723  $ iprepad, ipostpad,
724  $ padval )
725  CALL pschekpad( ictxt, 'PSPORFS', lwork,
726  $ 1, mem( ipw-iprepad ),
727  $ lwork, iprepad, ipostpad,
728  $ padval )
729  CALL pschekpad( ictxt, 'PSPORFS',
730  $ lw2-ipostpad, 1,
731  $ mem( ipw2-iprepad ),
732  $ lw2-ipostpad,
733  $ iprepad, ipostpad,
734  $ padval )
735 *
736  CALL psfillpad( ictxt, worksiz-ipostpad,
737  $ 1, mem( ipw-iprepad ),
738  $ worksiz-ipostpad, iprepad,
739  $ ipostpad, padval )
740 *
741 * check the solution to rhs
742 *
743  CALL pslaschk( 'Symm', 'Diag', n, nrhs,
744  $ mem( ipb ), 1, 1, descb,
745  $ iaseed, 1, 1, desca,
746  $ ibseed, anorm, sresid2,
747  $ mem( ipw ) )
748 *
749  IF( iam.EQ.0 .AND. sresid2.GT.thresh )
750  $ WRITE( nout, fmt = 9985 ) sresid2
751 *
752 * check for memory overwrite
753 *
754  CALL pschekpad( ictxt, 'PSLASCHK', np,
755  $ myrhs, mem( ipb-iprepad ),
756  $ descb( lld_ ), iprepad,
757  $ ipostpad, padval )
758  CALL pschekpad( ictxt, 'PSLASCHK',
759  $ worksiz-ipostpad, 1,
760  $ mem( ipw-iprepad ),
761  $ worksiz-ipostpad,
762  $ iprepad, ipostpad,
763  $ padval )
764  END IF
765  END IF
766 *
767 * Gather maximum of all CPU and WALL clock timings
768 *
769  CALL slcombine( ictxt, 'All', '>', 'W', 2, 1,
770  $ wtime )
771  CALL slcombine( ictxt, 'All', '>', 'C', 2, 1,
772  $ ctime )
773 *
774 * Print results
775 *
776  IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
777 *
778 * 1/3 N^3 + 1/2 N^2 flops for LLt factorization
779 *
780  nops = (dble(n)**3)/3.0d+0 +
781  $ (dble(n)**2)/2.0d+0
782 *
783 * nrhs * 2 N^2 flops for LLt solve.
784 *
785  nops = nops + 2.0d+0*(dble(n)**2)*dble(nrhs)
786 *
787 * Calculate total megaflops -- factorization and
788 * solve -- for WALL and CPU time, and print output
789 *
790 * Print WALL time if machine supports it
791 *
792  IF( wtime( 1 ) + wtime( 2 ) .GT. 0.0d+0 ) THEN
793  tmflops = nops /
794  $ ( ( wtime( 1 )+wtime( 2 ) ) * 1.0d+6 )
795  ELSE
796  tmflops = 0.0d+0
797  END IF
798 *
799  IF( wtime( 2 ).GE.0.0d+0 )
800  $ WRITE( nout, fmt = 9993 ) 'WALL', uplo, n,
801  $ nb, nrhs, nbrhs, nprow, npcol,
802  $ wtime( 1 ), wtime( 2 ), tmflops,
803  $ passed
804 *
805 * Print CPU time if machine supports it
806 *
807  IF( ctime( 1 )+ctime( 2 ).GT.0.0d+0 ) THEN
808  tmflops = nops /
809  $ ( ( ctime( 1 )+ctime( 2 ) ) * 1.0d+6 )
810  ELSE
811  tmflops = 0.0d+0
812  END IF
813 *
814  IF( ctime( 2 ).GE.0.0d+0 )
815  $ WRITE( nout, fmt = 9993 ) 'CPU ', uplo, n,
816  $ nb, nrhs, nbrhs, nprow, npcol,
817  $ ctime( 1 ), ctime( 2 ), tmflops,
818  $ passed
819 *
820  END IF
821  10 CONTINUE
822  20 CONTINUE
823 *
824  IF( check .AND. sresid.GT.thresh ) THEN
825 *
826 * Compute FRESID = ||A - LL'|| / (||A|| * N * eps)
827 *
828  CALL pspotrrv( uplo, n, mem( ipa ), 1, 1, desca,
829  $ mem( ipw ) )
830  CALL pslafchk( 'Symm', 'Diag', n, n, mem( ipa ), 1, 1,
831  $ desca, iaseed, anorm, fresid,
832  $ mem( ipw ) )
833 *
834 * Check for memory overwrite
835 *
836  CALL pschekpad( ictxt, 'PSPOTRRV', np, nq,
837  $ mem( ipa-iprepad ), desca( lld_ ),
838  $ iprepad, ipostpad, padval )
839  CALL pschekpad( ictxt, 'PSGETRRV',
840  $ worksiz-ipostpad, 1,
841  $ mem( ipw-iprepad ), worksiz-ipostpad,
842  $ iprepad, ipostpad, padval )
843 *
844  IF( iam.EQ.0 ) THEN
845  IF( lsame( uplo, 'L' ) ) THEN
846  WRITE( nout, fmt = 9986 ) 'L*L''', fresid
847  ELSE
848  WRITE( nout, fmt = 9986 ) 'U''*U', fresid
849  END IF
850  END IF
851  END IF
852 *
853  30 CONTINUE
854  40 CONTINUE
855  CALL blacs_gridexit( ictxt )
856  50 CONTINUE
857 *
858 * Print ending messages and close output file
859 *
860  60 CONTINUE
861  IF( iam.EQ.0 ) THEN
862  ktests = kpass + kfail + kskip
863  WRITE( nout, fmt = * )
864  WRITE( nout, fmt = 9992 ) ktests
865  IF( check ) THEN
866  WRITE( nout, fmt = 9991 ) kpass
867  WRITE( nout, fmt = 9989 ) kfail
868  ELSE
869  WRITE( nout, fmt = 9990 ) kpass
870  END IF
871  WRITE( nout, fmt = 9988 ) kskip
872  WRITE( nout, fmt = * )
873  WRITE( nout, fmt = * )
874  WRITE( nout, fmt = 9987 )
875  IF( nout.NE.6 .AND. nout.NE.0 )
876  $ CLOSE ( nout )
877  END IF
878 *
879  CALL blacs_exit( 0 )
880 *
881  9999 FORMAT( 'ILLEGAL ', a6, ': ', a5, ' = ', i3,
882  $ '; It should be at least 1' )
883  9998 FORMAT( 'ILLEGAL GRID: nprow*npcol = ', i4, '. It can be at most',
884  $ i4 )
885  9997 FORMAT( 'Bad ', a6, ' parameters: going on to next test case.' )
886  9996 FORMAT( 'Unable to perform ', a, ': need TOTMEM of at least',
887  $ i11 )
888  9995 FORMAT( 'TIME UPLO N NB NRHS NBRHS P Q LLt Time ',
889  $ 'Slv Time MFLOPS CHECK' )
890  9994 FORMAT( '---- ---- ----- --- ---- ----- ---- ---- -------- ',
891  $ '-------- -------- ------' )
892  9993 FORMAT( a4, 4x, a1, 1x, i5, 1x, i3, 1x, i4, 1x, i5, 1x, i4, 1x,
893  $ i4, 1x, f8.2, 1x, f8.2, 1x, f8.2, 1x, a6 )
894  9992 FORMAT( 'Finished ', i6, ' tests, with the following results:' )
895  9991 FORMAT( i5, ' tests completed and passed residual checks.' )
896  9990 FORMAT( i5, ' tests completed without checking.' )
897  9989 FORMAT( i5, ' tests completed and failed residual checks.' )
898  9988 FORMAT( i5, ' tests skipped because of illegal input values.' )
899  9987 FORMAT( 'END OF TESTS.' )
900  9986 FORMAT( '||A - ', a4, '|| / (||A|| * N * eps) = ', g25.7 )
901  9985 FORMAT( '||Ax-b||/(||x||*||A||*eps*N) ', f25.7 )
902 *
903  stop
904 *
905 * End of PSLLTDRIVER
906 *
907  END
pslafchk
subroutine pslafchk(AFORM, DIAG, M, N, A, IA, JA, DESCA, IASEED, ANORM, FRESID, WORK)
Definition: pslafchk.f:3
max
#define max(A, B)
Definition: pcgemr.c:180
ilcm
integer function ilcm(M, N)
Definition: ilcm.f:2
sltimer
subroutine sltimer(I)
Definition: sltimer.f:47
lsame
logical function lsame(CA, CB)
Definition: tools.f:1724
pschekpad
subroutine pschekpad(ICTXT, MESS, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: pschekpad.f:3
pspocon
subroutine pspocon(UPLO, N, A, IA, JA, DESCA, ANORM, RCOND, WORK, LWORK, IWORK, LIWORK, INFO)
Definition: pspocon.f:3
pslansy
real function pslansy(NORM, UPLO, N, A, IA, JA, DESCA, WORK)
Definition: pslansy.f:3
pslltinfo
subroutine pslltinfo(SUMMRY, NOUT, UPLO, NMAT, NVAL, LDNVAL, NNB, NBVAL, LDNBVAL, NNR, NRVAL, LDNRVAL, NNBR, NBRVAL, LDNBRVAL, NGRIDS, PVAL, LDPVAL, QVAL, LDQVAL, THRESH, EST, WORK, IAM, NPROCS)
Definition: pslltinfo.f:6
pspotrs
subroutine pspotrs(UPLO, N, NRHS, A, IA, JA, DESCA, B, IB, JB, DESCB, INFO)
Definition: pspotrs.f:3
descinit
subroutine descinit(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD, INFO)
Definition: descinit.f:3
slboot
subroutine slboot()
Definition: sltimer.f:2
psmatgen
subroutine psmatgen(ICTXT, AFORM, DIAG, M, N, MB, NB, A, LDA, IAROW, IACOL, ISEED, IROFF, IRNUM, ICOFF, ICNUM, MYROW, MYCOL, NPROW, NPCOL)
Definition: psmatgen.f:4
numroc
integer function numroc(N, NB, IPROC, ISRCPROC, NPROCS)
Definition: numroc.f:2
pslltdriver
program pslltdriver
Definition: pslltdriver.f:1
pspotrf
subroutine pspotrf(UPLO, N, A, IA, JA, DESCA, INFO)
Definition: pspotrf.f:2
psfillpad
subroutine psfillpad(ICTXT, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: psfillpad.f:2
pspotrrv
subroutine pspotrrv(UPLO, N, A, IA, JA, DESCA, WORK)
Definition: pspotrrv.f:2
slcombine
subroutine slcombine(ICTXT, SCOPE, OP, TIMETYPE, N, IBEG, TIMES)
Definition: sltimer.f:267
pslaschk
subroutine pslaschk(SYMM, DIAG, N, NRHS, X, IX, JX, DESCX, IASEED, IA, JA, DESCA, IBSEED, ANORM, RESID, WORK)
Definition: pslaschk.f:4
min
#define min(A, B)
Definition: pcgemr.c:181
iceil
integer function iceil(INUM, IDENOM)
Definition: iceil.f:2
psporfs
subroutine psporfs(UPLO, N, NRHS, A, IA, JA, DESCA, AF, IAF, JAF, DESCAF, B, IB, JB, DESCB, X, IX, JX, DESCX, FERR, BERR, WORK, LWORK, IWORK, LIWORK, INFO)
Definition: psporfs.f:4