ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pdinvdriver.f
Go to the documentation of this file.
1  PROGRAM pdinvdriver
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 * PDINVDRIVER is the main test program for the DOUBLE PRECISION
12 * SCALAPACK matrix inversion routines. This test driver computes the
13 * inverse of different kind of matrix and tests the results.
14 *
15 * The program must be driven by a short data file. An annotated example
16 * of a data file can be obtained by deleting the first 3 characters
17 * from the following 14 lines:
18 * 'ScaLAPACK Matrix Inversion Testing input file'
19 * 'PVM machine.'
20 * 'INV.out' output file name (if any)
21 * 6 device out
22 * 5 number of matrix types (next line)
23 * 'GEN' 'UTR' 'LTR' 'UPD' LPD' GEN, UTR, LTR, UPD, LPD
24 * 4 number of problems sizes
25 * 1000 2000 3000 4000 values of N
26 * 3 number of NB's
27 * 4 30 35 values of NB
28 * 2 number of process grids (ordered P & Q)
29 * 4 2 values of P
30 * 4 4 values of Q
31 * 1.0 threshold
32 *
33 * Internal Parameters
34 * ===================
35 *
36 * TOTMEM INTEGER, default = 2000000
37 * TOTMEM is a machine-specific parameter indicating the
38 * maximum amount of available memory in bytes.
39 * The user should customize TOTMEM to his platform. Remember
40 * to leave room in memory for the operating system, the BLACS
41 * buffer, etc. For example, on a system with 8 MB of memory
42 * per process (e.g., one processor on an Intel iPSC/860), the
43 * parameters we use are TOTMEM=6200000 (leaving 1.8 MB for OS,
44 * code, BLACS buffer, etc). However, for PVM, we usually set
45 * TOTMEM = 2000000. Some experimenting with the maximum value
46 * of TOTMEM may be required.
47 *
48 * INTGSZ INTEGER, default = 4 bytes.
49 * DBLESZ INTEGER, default = 8 bytes.
50 * INTGSZ and DBLESZ indicate the length in bytes on the
51 * given platform for an integer and a double precision real.
52 * MEM DOUBLE PRECISION array, dimension ( TOTMEM / DBLESZ )
53 *
54 * All arrays used by SCALAPACK routines are allocated from
55 * this array and referenced by pointers. The integer IPA,
56 * for example, is a pointer to the starting element of MEM for
57 * the matrix A.
58 *
59 * =====================================================================
60 *
61 * .. Parameters ..
62  INTEGER block_cyclic_2d, csrc_, ctxt_, dlen_, dtype_,
63  $ lld_, mb_, m_, nb_, n_, rsrc_
64  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
65  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
66  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
67  INTEGER dblesz, intgsz, memsiz, ntests, totmem
68  DOUBLE PRECISION padval, zero
69  parameter( dblesz = 8, intgsz = 4, totmem = 2000000,
70  $ memsiz = totmem / dblesz, ntests = 20,
71  $ padval = -9923.0d+0, zero = 0.0d+0 )
72 * ..
73 * .. Local Scalars ..
74  CHARACTER uplo
75  CHARACTER*3 mtyp
76  CHARACTER*6 passed
77  CHARACTER*80 outfile
78  LOGICAL check
79  INTEGER i, iam, iaseed, ictxt, imidpad, info, ipa,
80  $ ippiv, iprepad, ipostpad, ipiw, ipw, itemp, j,
81  $ k, ktests, kpass, kfail, kskip, l, lcm, lipiv,
82  $ liwork, lwork, mycol, myrow, n, nb, ngrids,
83  $ nmat, nmtyp, nnb, nout, np, npcol, nprocs,
84  $ nprow, nq, workiinv, workinv, worksiz
85  REAL thresh
86  DOUBLE PRECISION anorm, fresid, nops, rcond, tmflops
87 * ..
88 * .. Local Arrays ..
89  CHARACTER*3 mattyp( ntests )
90  INTEGER desca( dlen_ ), ierr( 1 ), nbval( ntests ),
91  $ nval( ntests ), pval( ntests ),
92  $ qval( ntests )
93  DOUBLE PRECISION mem( memsiz ), ctime( 2 ), wtime( 2 )
94 * ..
95 * .. External Subroutines ..
96  EXTERNAL blacs_barrier, blacs_exit, blacs_get,
97  $ blacs_gridexit, blacs_gridinfo, blacs_gridinit,
98  $ blacs_pinfo, descinit, igsum2d, pdchekpad,
99  $ pdfillpad, pdgetrf, pdgetri,
100  $ pdinvchk, pdinvinfo, pdlaset,
101  $ pdmatgen, pdpotrf, pdpotri,
102  $ pdtrtri, slboot, slcombine, sltimer
103 * ..
104 * .. External Functions ..
105  LOGICAL lsamen
106  INTEGER iceil, ilcm, numroc
107  DOUBLE PRECISION pdlange, pdlansy, pdlantr
108  EXTERNAL iceil, ilcm, lsamen, numroc, pdlange,
109  $ pdlansy, pdlantr
110 * ..
111 * .. Intrinsic Functions ..
112  INTRINSIC dble, max, mod
113 * ..
114 * .. Data Statements ..
115  DATA ktests, kpass, kfail, kskip /4*0/
116 * ..
117 * .. Executable Statements ..
118 *
119 * Get starting information
120 *
121  CALL blacs_pinfo( iam, nprocs )
122  iaseed = 100
123  CALL pdinvinfo( outfile, nout, nmtyp, mattyp, ntests, nmat, nval,
124  $ ntests, nnb, nbval, ntests, ngrids, pval, ntests,
125  $ qval, ntests, thresh, mem, iam, nprocs )
126  check = ( thresh.GE.0.0e+0 )
127 *
128 * Loop over the different matrix types
129 *
130  DO 40 i = 1, nmtyp
131 *
132  mtyp = mattyp( i )
133 *
134 * Print headings
135 *
136  IF( iam.EQ.0 ) THEN
137  WRITE( nout, fmt = * )
138  IF( lsamen( 3, mtyp, 'GEN' ) ) THEN
139  WRITE( nout, fmt = 9986 )
140  $ 'A is a general matrix.'
141  ELSE IF( lsamen( 3, mtyp, 'UTR' ) ) THEN
142  WRITE( nout, fmt = 9986 )
143  $ 'A is an upper triangular matrix.'
144  ELSE IF( lsamen( 3, mtyp, 'LTR' ) ) THEN
145  WRITE( nout, fmt = 9986 )
146  $ 'A is a lower triangular matrix.'
147  ELSE IF( lsamen( 3, mtyp, 'UPD' ) ) THEN
148  WRITE( nout, fmt = 9986 )
149  $ 'A is a symmetric positive definite matrix.'
150  WRITE( nout, fmt = 9986 )
151  $ 'Only the upper triangular part will be '//
152  $ 'referenced.'
153  ELSE IF( lsamen( 3, mtyp, 'LPD' ) ) THEN
154  WRITE( nout, fmt = 9986 )
155  $ 'A is a symmetric positive definite matrix.'
156  WRITE( nout, fmt = 9986 )
157  $ 'Only the lower triangular part will be '//
158  $ 'referenced.'
159  END IF
160  WRITE( nout, fmt = * )
161  WRITE( nout, fmt = 9995 )
162  WRITE( nout, fmt = 9994 )
163  WRITE( nout, fmt = * )
164  END IF
165 *
166 * Loop over different process grids
167 *
168  DO 30 j = 1, ngrids
169 *
170  nprow = pval( j )
171  npcol = qval( j )
172 *
173 * Make sure grid information is correct
174 *
175  ierr( 1 ) = 0
176  IF( nprow.LT.1 ) THEN
177  IF( iam.EQ.0 )
178  $ WRITE( nout, fmt = 9999 ) 'GRID', 'nprow', nprow
179  ierr( 1 ) = 1
180  ELSE IF( npcol.LT.1 ) THEN
181  IF( iam.EQ.0 )
182  $ WRITE( nout, fmt = 9999 ) 'GRID', 'npcol', npcol
183  ierr( 1 ) = 1
184  ELSE IF( nprow*npcol.GT.nprocs ) THEN
185  IF( iam.EQ.0 )
186  $ WRITE( nout, fmt = 9998 ) nprow*npcol, nprocs
187  ierr( 1 ) = 1
188  END IF
189 *
190  IF( ierr( 1 ).GT.0 ) THEN
191  IF( iam.EQ.0 )
192  $ WRITE( nout, fmt = 9997 ) 'grid'
193  kskip = kskip + 1
194  GO TO 30
195  END IF
196 *
197 * Define process grid
198 *
199  CALL blacs_get( -1, 0, ictxt )
200  CALL blacs_gridinit( ictxt, 'Row-major', nprow, npcol )
201  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
202 *
203 * Go to bottom of loop if this case doesn't use my process
204 *
205  IF( myrow.GE.nprow .OR. mycol.GE.npcol )
206  $ GO TO 30
207 *
208  DO 20 k = 1, nmat
209 *
210  n = nval( k )
211 *
212 * Make sure matrix information is correct
213 *
214  ierr( 1 ) = 0
215  IF( n.LT.1 ) THEN
216  IF( iam.EQ.0 )
217  $ WRITE( nout, fmt = 9999 ) 'MATRIX', 'N', n
218  ierr( 1 ) = 1
219  END IF
220 *
221 * Make sure no one had error
222 *
223  CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1, 0 )
224 *
225  IF( ierr( 1 ).GT.0 ) THEN
226  IF( iam.EQ.0 )
227  $ WRITE( nout, fmt = 9997 ) 'matrix'
228  kskip = kskip + 1
229  GO TO 20
230  END IF
231 *
232 * Loop over different blocking sizes
233 *
234  DO 10 l = 1, nnb
235 *
236  nb = nbval( l )
237 *
238 * Make sure nb is legal
239 *
240  ierr( 1 ) = 0
241  IF( nb.LT.1 ) THEN
242  ierr( 1 ) = 1
243  IF( iam.EQ.0 )
244  $ WRITE( nout, fmt = 9999 ) 'NB', 'NB', nb
245  END IF
246 *
247 * Check all processes for an error
248 *
249  CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1,
250  $ 0 )
251 *
252  IF( ierr( 1 ).GT.0 ) THEN
253  IF( iam.EQ.0 )
254  $ WRITE( nout, fmt = 9997 ) 'NB'
255  kskip = kskip + 1
256  GO TO 10
257  END IF
258 *
259 * Padding constants
260 *
261  np = numroc( n, nb, myrow, 0, nprow )
262  nq = numroc( n, nb, mycol, 0, npcol )
263  IF( check ) THEN
264  iprepad = max( nb, np )
265  imidpad = nb
266  ipostpad = max( nb, nq )
267  ELSE
268  iprepad = 0
269  imidpad = 0
270  ipostpad = 0
271  END IF
272 *
273 * Initialize the array descriptor for the matrix A
274 *
275  CALL descinit( desca, n, n, nb, nb, 0, 0, ictxt,
276  $ max( 1, np ) + imidpad, ierr( 1 ) )
277 *
278 * Check all processes for an error
279 *
280  CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1,
281  $ 0 )
282 *
283  IF( ierr( 1 ).LT.0 ) THEN
284  IF( iam.EQ.0 )
285  $ WRITE( nout, fmt = 9997 ) 'descriptor'
286  kskip = kskip + 1
287  GO TO 10
288  END IF
289 *
290 * Assign pointers into MEM for ScaLAPACK arrays, A is
291 * allocated starting at position MEM( IPREPAD+1 )
292 *
293  ipa = iprepad+1
294 *
295  lcm = ilcm( nprow, npcol )
296  IF( lsamen( 3, mtyp, 'GEN' ) ) THEN
297 *
298 * Pivots are needed by LU factorization
299 *
300  ippiv = ipa + desca( lld_ ) * nq + ipostpad +
301  $ iprepad
302  lipiv = iceil( intgsz * ( np + nb ), dblesz )
303  ipw = ippiv + lipiv + ipostpad + iprepad
304 *
305  lwork = max( 1, np * desca( nb_ ) )
306  workinv = lwork + ipostpad
307 *
308 * Figure the amount of workspace required by the
309 * general matrix inversion
310 *
311  IF( nprow.EQ.npcol ) THEN
312  liwork = nq + desca( nb_ )
313  ELSE
314 *
315 * change the integer workspace needed for PDGETRI
316 * LIWORK = MAX( DESCA( NB_ ), DESCA( MB_ ) *
317 * $ ICEIL( ICEIL( DESCA( LLD_ ),
318 * $ DESCA( MB_ ) ), LCM / NPROW ) )
319 * $ + NQ
320  liwork = numroc( desca( m_ ) +
321  $ desca( mb_ ) * nprow
322  $ + mod( 1 - 1, desca( mb_ ) ), desca( nb_ ),
323  $ mycol, desca( csrc_ ), npcol ) +
324  $ max( desca( mb_ ) * iceil( iceil(
325  $ numroc( desca( m_ ) + desca( mb_ ) * nprow,
326  $ desca( mb_ ), myrow, desca( rsrc_ ), nprow ),
327  $ desca( mb_ ) ), lcm / nprow ), desca( nb_ ) )
328 *
329  END IF
330  workiinv = iceil( liwork*intgsz, dblesz ) +
331  $ ipostpad
332  ipiw = ipw + workinv + iprepad
333  worksiz = workinv + iprepad + workiinv
334 *
335  ELSE
336 *
337 * No pivots or workspace needed for triangular or
338 * symmetric positive definite matrices.
339 *
340  ipw = ipa + desca( lld_ ) * nq + ipostpad + iprepad
341  worksiz = 1 + ipostpad
342 *
343  END IF
344 *
345  IF( check ) THEN
346 *
347 * Figure amount of work space for the norm
348 * computations
349 *
350  IF( lsamen( 3, mtyp, 'GEN' ).OR.
351  $ lsamen( 2, mtyp( 2:3 ), 'TR' ) ) THEN
352  itemp = nq
353  ELSE
354  itemp = 2 * nq + np
355  IF( nprow.NE.npcol ) THEN
356  itemp = itemp +
357  $ nb * iceil( iceil( np, nb ),
358  $ lcm / nprow )
359  END IF
360  END IF
361  worksiz = max( worksiz-ipostpad, itemp )
362 *
363 * Figure the amount of workspace required by the
364 * checking routine
365 *
366  worksiz = max( worksiz, 2 * nb * max( 1, np ) ) +
367  $ ipostpad
368 *
369  END IF
370 *
371 * Check for adequate memory for problem size
372 *
373  ierr( 1 ) = 0
374  IF( ipw+worksiz.GT.memsiz ) THEN
375  IF( iam.EQ.0 )
376  $ WRITE( nout, fmt = 9996 ) 'inversion',
377  $ ( ipw + worksiz ) * dblesz
378  ierr( 1 ) = 1
379  END IF
380 *
381 * Check all processes for an error
382 *
383  CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1,
384  $ 0 )
385 *
386  IF( ierr( 1 ).GT.0 ) THEN
387  IF( iam.EQ.0 )
388  $ WRITE( nout, fmt = 9997 ) 'MEMORY'
389  kskip = kskip + 1
390  GO TO 10
391  END IF
392 *
393  IF( lsamen( 3, mtyp, 'GEN' ).OR.
394  $ lsamen( 2, mtyp( 2:3 ), 'TR' ) ) THEN
395 *
396 * Generate a general diagonally dominant matrix A
397 *
398  CALL pdmatgen( ictxt, 'N', 'D', desca( m_ ),
399  $ desca( n_ ), desca( mb_ ),
400  $ desca( nb_ ), mem( ipa ),
401  $ desca( lld_ ), desca( rsrc_ ),
402  $ desca( csrc_ ), iaseed, 0, np, 0,
403  $ nq, myrow, mycol, nprow, npcol )
404 *
405  ELSE IF( lsamen( 2, mtyp( 2:3 ), 'PD' ) ) THEN
406 *
407 * Generate a symmetric positive definite matrix
408 *
409  CALL pdmatgen( ictxt, 'S', 'D', desca( m_ ),
410  $ desca( n_ ), desca( mb_ ),
411  $ desca( nb_ ), mem( ipa ),
412  $ desca( lld_ ), desca( rsrc_ ),
413  $ desca( csrc_ ), iaseed, 0, np, 0,
414  $ nq, myrow, mycol, nprow, npcol )
415 *
416  END IF
417 *
418 * Zeros not-referenced part of A, if any.
419 *
420  IF( lsamen( 1, mtyp, 'U' ) ) THEN
421 *
422  uplo = 'U'
423  CALL pdlaset( 'Lower', n-1, n-1, zero, zero,
424  $ mem( ipa ), 2, 1, desca )
425 *
426  ELSE IF( lsamen( 1, mtyp, 'L' ) ) THEN
427 *
428  uplo = 'L'
429  CALL pdlaset( 'Upper', n-1, n-1, zero, zero,
430  $ mem( ipa ), 1, 2, desca )
431 *
432  ELSE
433 *
434  uplo = 'G'
435 *
436  END IF
437 *
438 * Need 1-norm of A for checking
439 *
440  IF( check ) THEN
441 *
442  CALL pdfillpad( ictxt, np, nq, mem( ipa-iprepad ),
443  $ desca( lld_ ), iprepad, ipostpad,
444  $ padval )
445  CALL pdfillpad( ictxt, worksiz-ipostpad, 1,
446  $ mem( ipw-iprepad ),
447  $ worksiz-ipostpad, iprepad,
448  $ ipostpad, padval )
449 *
450  IF( lsamen( 3, mtyp, 'GEN' ) ) THEN
451 *
452  CALL pdfillpad( ictxt, lipiv, 1,
453  $ mem( ippiv-iprepad ), lipiv,
454  $ iprepad, ipostpad, padval )
455  anorm = pdlange( '1', n, n, mem( ipa ), 1, 1,
456  $ desca, mem( ipw ) )
457  CALL pdchekpad( ictxt, 'PDLANGE', np, nq,
458  $ mem( ipa-iprepad ),
459  $ desca( lld_ ),
460  $ iprepad, ipostpad, padval )
461  CALL pdchekpad( ictxt, 'PDLANGE',
462  $ worksiz-ipostpad, 1,
463  $ mem( ipw-iprepad ),
464  $ worksiz-ipostpad,
465  $ iprepad, ipostpad, padval )
466  CALL pdfillpad( ictxt, workinv-ipostpad, 1,
467  $ mem( ipw-iprepad ),
468  $ workinv-ipostpad,
469  $ iprepad, ipostpad, padval )
470  CALL pdfillpad( ictxt, workiinv-ipostpad, 1,
471  $ mem( ipiw-iprepad ),
472  $ workiinv-ipostpad, iprepad,
473  $ ipostpad, padval )
474  ELSE IF( lsamen( 2, mtyp( 2:3 ), 'TR' ) ) THEN
475 *
476  anorm = pdlantr( '1', uplo, 'Non unit', n, n,
477  $ mem( ipa ), 1, 1, desca,
478  $ mem( ipw ) )
479  CALL pdchekpad( ictxt, 'PDLANTR', np, nq,
480  $ mem( ipa-iprepad ),
481  $ desca( lld_ ),
482  $ iprepad, ipostpad, padval )
483  CALL pdchekpad( ictxt, 'PDLANTR',
484  $ worksiz-ipostpad, 1,
485  $ mem( ipw-iprepad ),
486  $ worksiz-ipostpad,
487  $ iprepad, ipostpad, padval )
488 *
489  ELSE IF( lsamen( 2, mtyp( 2:3 ), 'PD' ) ) THEN
490 *
491  anorm = pdlansy( '1', uplo, n, mem( ipa ), 1, 1,
492  $ desca, mem( ipw ) )
493  CALL pdchekpad( ictxt, 'PDLANSY', np, nq,
494  $ mem( ipa-iprepad ),
495  $ desca( lld_ ),
496  $ iprepad, ipostpad, padval )
497  CALL pdchekpad( ictxt, 'PDLANSY',
498  $ worksiz-ipostpad, 1,
499  $ mem( ipw-iprepad ),
500  $ worksiz-ipostpad,
501  $ iprepad, ipostpad, padval )
502 *
503  ELSE IF( lsamen( 2, mtyp( 2:3 ), 'SY' ) ) THEN
504 *
505  CALL pdfillpad( ictxt, lipiv, 1,
506  $ mem( ippiv-iprepad ), lipiv,
507  $ iprepad, ipostpad, padval )
508  anorm = pdlansy( '1', uplo, n, mem( ipa ), 1, 1,
509  $ desca, mem( ipw ) )
510  CALL pdchekpad( ictxt, 'PDLANSY', np, nq,
511  $ mem( ipa-iprepad ),
512  $ desca( lld_ ),
513  $ iprepad, ipostpad, padval )
514  CALL pdchekpad( ictxt, 'PDLANSY',
515  $ worksiz-ipostpad, 1,
516  $ mem( ipw-iprepad ),
517  $ worksiz-ipostpad,
518  $ iprepad,ipostpad, padval )
519 *
520  END IF
521 *
522  END IF
523 *
524  CALL slboot()
525  CALL blacs_barrier( ictxt, 'All' )
526 *
527  IF( lsamen( 3, mtyp, 'GEN' ) ) THEN
528 *
529 * Perform LU factorization
530 *
531  CALL sltimer( 1 )
532  CALL pdgetrf( n, n, mem( ipa ), 1, 1, desca,
533  $ mem( ippiv ), info )
534  CALL sltimer( 1 )
535 *
536  IF( check ) THEN
537 *
538 * Check for memory overwrite
539 *
540  CALL pdchekpad( ictxt, 'PDGETRF', np, nq,
541  $ mem( ipa-iprepad ),
542  $ desca( lld_ ),
543  $ iprepad, ipostpad, padval )
544  CALL pdchekpad( ictxt, 'PDGETRF', lipiv, 1,
545  $ mem( ippiv-iprepad ), lipiv,
546  $ iprepad, ipostpad, padval )
547  END IF
548 *
549 * Perform the general matrix inversion
550 *
551  CALL sltimer( 2 )
552  CALL pdgetri( n, mem( ipa ), 1, 1, desca,
553  $ mem( ippiv ), mem( ipw ), lwork,
554  $ mem( ipiw ), liwork, info )
555  CALL sltimer( 2 )
556 *
557  IF( check ) THEN
558 *
559 * Check for memory overwrite
560 *
561  CALL pdchekpad( ictxt, 'PDGETRI', np, nq,
562  $ mem( ipa-iprepad ),
563  $ desca( lld_ ),
564  $ iprepad, ipostpad, padval )
565  CALL pdchekpad( ictxt, 'PDGETRI', lipiv, 1,
566  $ mem( ippiv-iprepad ), lipiv,
567  $ iprepad, ipostpad, padval )
568  CALL pdchekpad( ictxt, 'PDGETRI',
569  $ workiinv-ipostpad, 1,
570  $ mem( ipiw-iprepad ),
571  $ workiinv-ipostpad,
572  $ iprepad, ipostpad, padval )
573  CALL pdchekpad( ictxt, 'PDGETRI',
574  $ workinv-ipostpad, 1,
575  $ mem( ipw-iprepad ),
576  $ workinv-ipostpad,
577  $ iprepad, ipostpad, padval )
578  END IF
579 *
580  ELSE IF( lsamen( 2, mtyp( 2:3 ), 'TR' ) ) THEN
581 *
582 * Perform the general matrix inversion
583 *
584  CALL sltimer( 2 )
585  CALL pdtrtri( uplo, 'Non unit', n, mem( ipa ), 1,
586  $ 1, desca, info )
587  CALL sltimer( 2 )
588 *
589  IF( check ) THEN
590 *
591 * Check for memory overwrite
592 *
593  CALL pdchekpad( ictxt, 'PDTRTRI', np, nq,
594  $ mem( ipa-iprepad ),
595  $ desca( lld_ ),
596  $ iprepad, ipostpad, padval )
597  END IF
598 *
599  ELSE IF( lsamen( 2, mtyp( 2:3 ), 'PD' ) ) THEN
600 *
601 * Perform Cholesky factorization
602 *
603  CALL sltimer( 1 )
604  CALL pdpotrf( uplo, n, mem( ipa ), 1, 1, desca,
605  $ info )
606  CALL sltimer( 1 )
607 *
608  IF( check ) THEN
609 *
610 * Check for memory overwrite
611 *
612  CALL pdchekpad( ictxt, 'PDPOTRF', np, nq,
613  $ mem( ipa-iprepad ),
614  $ desca( lld_ ),
615  $ iprepad, ipostpad, padval )
616  END IF
617 *
618 * Perform the symmetric positive definite matrix
619 * inversion
620 *
621  CALL sltimer( 2 )
622  CALL pdpotri( uplo, n, mem( ipa ), 1, 1, desca,
623  $ info )
624  CALL sltimer( 2 )
625 *
626  IF( check ) THEN
627 *
628 * Check for memory overwrite
629 *
630  CALL pdchekpad( ictxt, 'PDPOTRI', np, nq,
631  $ mem( ipa-iprepad ),
632  $ desca( lld_ ),
633  $ iprepad, ipostpad, padval )
634  END IF
635 *
636  END IF
637 *
638  IF( check ) THEN
639 *
640  CALL pdfillpad( ictxt, worksiz-ipostpad, 1,
641  $ mem( ipw-iprepad ),
642  $ worksiz-ipostpad, iprepad,
643  $ ipostpad, padval )
644 *
645 * Compute fresid = || inv(A)*A-I ||
646 *
647  CALL pdinvchk( mtyp, n, mem( ipa ), 1, 1, desca,
648  $ iaseed, anorm, fresid, rcond,
649  $ mem( ipw ) )
650 *
651 * Check for memory overwrite
652 *
653  CALL pdchekpad( ictxt, 'PDINVCHK', np, nq,
654  $ mem( ipa-iprepad ),
655  $ desca( lld_ ),
656  $ iprepad, ipostpad, padval )
657  CALL pdchekpad( ictxt, 'PDINVCHK',
658  $ worksiz-ipostpad, 1,
659  $ mem( ipw-iprepad ),
660  $ worksiz-ipostpad, iprepad,
661  $ ipostpad, padval )
662 *
663 * Test residual and detect NaN result
664 *
665  IF( fresid.LE.thresh .AND. info.EQ.0 .AND.
666  $ ( (fresid-fresid) .EQ. 0.0d+0 ) ) THEN
667  kpass = kpass + 1
668  passed = 'PASSED'
669  ELSE
670  kfail = kfail + 1
671  IF( info.GT.0 ) THEN
672  passed = 'SINGUL'
673  ELSE
674  passed = 'FAILED'
675  END IF
676  END IF
677 *
678  ELSE
679 *
680 * Don't perform the checking, only the timing
681 * operation
682 *
683  kpass = kpass + 1
684  fresid = fresid - fresid
685  passed = 'BYPASS'
686 *
687  END IF
688 *
689 * Gather maximum of all CPU and WALL clock timings
690 *
691  CALL slcombine( ictxt, 'All', '>', 'W', 2, 1, wtime )
692  CALL slcombine( ictxt, 'All', '>', 'C', 2, 1, ctime )
693 *
694 * Print results
695 *
696  IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
697 *
698  IF( lsamen( 3, mtyp, 'GEN' ) ) THEN
699 *
700 * 2/3 N^3 - 1/2 N^2 flops for LU factorization
701 *
702  nops = ( 2.0d+0 / 3.0d+0 )*( dble( n )**3 ) -
703  $ ( 1.0d+0 / 2.0d+0 )*( dble( n )**2 )
704 *
705 * 4/3 N^3 - N^2 flops for inversion
706 *
707  nops = nops +
708  $ ( 4.0d+0 / 3.0d+0 ) * ( dble( n )**3 ) -
709  $ dble( n )**2
710 *
711  ELSE IF( lsamen( 2, mtyp( 2:3 ), 'TR' ) ) THEN
712 *
713 * 1/3 N^3 + 2/3 N flops for triangular inversion
714 *
715  ctime(1) = 0.0d+0
716  wtime(1) = 0.0d+0
717  nops = ( 1.0d+0 / 3.0d+0 ) * ( dble( n )**3 ) +
718  $ ( 2.0d+0 / 3.0d+0 ) * ( dble( n ) )
719 *
720  ELSE IF( lsamen( 2, mtyp( 2:3 ), 'PD' ) ) THEN
721 *
722 * 1/3 N^3 + 1/2 N^2 flops for Cholesky
723 * factorization
724 *
725  nops = ( 1.0d+0 / 3.0d+0 ) * ( dble( n )**3 ) +
726  $ ( 1.0d+0 / 2.0d+0 ) * ( dble( n )**2 )
727 *
728 * 2/3 N^3 + 1/2 N^2 flops for Cholesky inversion
729 *
730  nops = nops +
731  $ ( 2.0d+0 / 3.0d+0 ) * ( dble( n )**3 ) +
732  $ ( 1.0d+0 / 2.0d+0 ) * ( dble( n )**2 )
733 *
734  END IF
735 *
736 * Figure total megaflops -- factorization and
737 * inversion, for WALL and CPU time, and print
738 * output.
739 *
740 * Print WALL time if machine supports it
741 *
742  IF( wtime( 1 ) + wtime( 2 ) .GT. 0.0d+0 ) THEN
743  tmflops = nops /
744  $ ( ( wtime( 1 )+wtime( 2 ) ) * 1.0d+6 )
745  ELSE
746  tmflops = 0.0d+0
747  END IF
748 *
749  IF( wtime( 2 ) .GE. 0.0d+0 )
750  $ WRITE( nout, fmt = 9993 ) 'WALL', n, nb, nprow,
751  $ npcol, wtime( 1 ), wtime( 2 ), tmflops,
752  $ rcond, fresid, passed
753 *
754 * Print CPU time if machine supports it
755 *
756  IF( ctime( 1 ) + ctime( 2 ) .GT. 0.0d+0 ) THEN
757  tmflops = nops /
758  $ ( ( ctime( 1 )+ctime( 2 ) ) * 1.0d+6 )
759  ELSE
760  tmflops = 0.0d+0
761  END IF
762 *
763  IF( ctime( 2 ) .GE. 0.0d+0 )
764  $ WRITE( nout, fmt = 9993 ) 'CPU ', n, nb, nprow,
765  $ npcol, ctime( 1 ), ctime( 2 ), tmflops,
766  $ rcond, fresid, passed
767  END IF
768 *
769  10 CONTINUE
770 *
771  20 CONTINUE
772 *
773  CALL blacs_gridexit( ictxt )
774 *
775  30 CONTINUE
776 *
777  40 CONTINUE
778 *
779 * Print out ending messages and close output file
780 *
781  IF( iam.EQ.0 ) THEN
782  ktests = kpass + kfail + kskip
783  WRITE( nout, fmt = * )
784  WRITE( nout, fmt = 9992 ) ktests
785  IF( check ) THEN
786  WRITE( nout, fmt = 9991 ) kpass
787  WRITE( nout, fmt = 9989 ) kfail
788  ELSE
789  WRITE( nout, fmt = 9990 ) kpass
790  END IF
791  WRITE( nout, fmt = 9988 ) kskip
792  WRITE( nout, fmt = * )
793  WRITE( nout, fmt = * )
794  WRITE( nout, fmt = 9987 )
795  IF( nout.NE.6 .AND. nout.NE.0 )
796  $ CLOSE ( nout )
797  END IF
798 *
799  CALL blacs_exit( 0 )
800 *
801  9999 FORMAT( 'ILLEGAL ', a6, ': ', a5, ' = ', i3,
802  $ '; It should be at least 1' )
803  9998 FORMAT( 'ILLEGAL GRID: nprow*npcol = ', i4, '. It can be at most',
804  $ i4 )
805  9997 FORMAT( 'Bad ', a6, ' parameters: going on to next test case.' )
806  9996 FORMAT( 'Unable to perform ', a, ': need TOTMEM of at least',
807  $ i11 )
808  9995 FORMAT( 'TIME N NB P Q Fct Time Inv Time ',
809  $ ' MFLOPS Cond Resid CHECK' )
810  9994 FORMAT( '---- ----- --- ----- ----- -------- -------- ',
811  $ '----------- ------- ------- ------' )
812  9993 FORMAT( a4, 1x, i5, 1x, i3, 1x, i5, 1x, i5, 1x, f8.2, 1x, f8.2,
813  $ 1x, f11.2, 1x, f7.1, 1x, f7.2, 1x, a6 )
814  9992 FORMAT( 'Finished ', i6, ' tests, with the following results:' )
815  9991 FORMAT( i5, ' tests completed and passed residual checks.' )
816  9990 FORMAT( i5, ' tests completed without checking.' )
817  9989 FORMAT( i5, ' tests completed and failed residual checks.' )
818  9988 FORMAT( i5, ' tests skipped because of illegal input values.' )
819  9987 FORMAT( 'END OF TESTS.' )
820  9986 FORMAT( a )
821 *
822  stop
823 *
824 * End of PDINVDRIVER
825 *
826  END
max
#define max(A, B)
Definition: pcgemr.c:180
pdinvchk
subroutine pdinvchk(MATTYP, N, A, IA, JA, DESCA, IASEED, ANORM, FRESID, RCOND, WORK)
Definition: pdinvchk.f:3
pdlansy
double precision function pdlansy(NORM, UPLO, N, A, IA, JA, DESCA, WORK)
Definition: pdlansy.f:3
pdgetri
subroutine pdgetri(N, A, IA, JA, DESCA, IPIV, WORK, LWORK, IWORK, LIWORK, INFO)
Definition: pdgetri.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
sltimer
subroutine sltimer(I)
Definition: sltimer.f:47
pdinvdriver
program pdinvdriver
Definition: pdinvdriver.f:1
pdgetrf
subroutine pdgetrf(M, N, A, IA, JA, DESCA, IPIV, INFO)
Definition: pdgetrf.f:2
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
pdlange
double precision function pdlange(NORM, M, N, A, IA, JA, DESCA, WORK)
Definition: pdlange.f:3
pdlaset
subroutine pdlaset(UPLO, M, N, ALPHA, BETA, A, IA, JA, DESCA)
Definition: pdblastst.f:6862
numroc
integer function numroc(N, NB, IPROC, ISRCPROC, NPROCS)
Definition: numroc.f:2
lsamen
logical function lsamen(N, CA, CB)
Definition: pblastst.f:1457
pdfillpad
subroutine pdfillpad(ICTXT, M, N, A, LDA, IPRE, IPOST, CHKVAL)
Definition: pdfillpad.f:2
pdpotrf
subroutine pdpotrf(UPLO, N, A, IA, JA, DESCA, INFO)
Definition: pdpotrf.f:2
pdlantr
double precision function pdlantr(NORM, UPLO, DIAG, M, N, A, IA, JA, DESCA, WORK)
Definition: pdlantr.f:3
pdpotri
subroutine pdpotri(UPLO, N, A, IA, JA, DESCA, INFO)
Definition: pdpotri.f:2
pdinvinfo
subroutine pdinvinfo(SUMMRY, NOUT, NMTYP, MATTYP, LDMTYP, NMAT, NVAL, LDNVAL, NNB, NBVAL, LDNBVAL, NGRIDS, PVAL, LDPVAL, QVAL, LDQVAL, THRESH, WORK, IAM, NPROCS)
Definition: pdinvinfo.f:5
slcombine
subroutine slcombine(ICTXT, SCOPE, OP, TIMETYPE, N, IBEG, TIMES)
Definition: sltimer.f:267
pdtrtri
subroutine pdtrtri(UPLO, DIAG, N, A, IA, JA, DESCA, INFO)
Definition: pdtrtri.f:2
iceil
integer function iceil(INUM, IDENOM)
Definition: iceil.f:2