SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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
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
descinit
subroutine descinit(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld, info)
Definition descinit.f:3
iceil
integer function iceil(inum, idenom)
Definition iceil.f:2
ilcm
integer function ilcm(m, n)
Definition ilcm.f:2
numroc
integer function numroc(n, nb, iproc, isrcproc, nprocs)
Definition numroc.f:2
lsamen
logical function lsamen(n, ca, cb)
Definition pblastst.f:1457
max
#define max(A, B)
Definition pcgemr.c:180
pdlaset
subroutine pdlaset(uplo, m, n, alpha, beta, a, ia, ja, desca)
Definition pdblastst.f:6862
pdchekpad
subroutine pdchekpad(ictxt, mess, m, n, a, lda, ipre, ipost, chkval)
Definition pdchekpad.f:3
pdfillpad
subroutine pdfillpad(ictxt, m, n, a, lda, ipre, ipost, chkval)
Definition pdfillpad.f:2
pdgetrf
subroutine pdgetrf(m, n, a, ia, ja, desca, ipiv, info)
Definition pdgetrf.f:2
pdgetri
subroutine pdgetri(n, a, ia, ja, desca, ipiv, work, lwork, iwork, liwork, info)
Definition pdgetri.f:3
pdinvchk
subroutine pdinvchk(mattyp, n, a, ia, ja, desca, iaseed, anorm, fresid, rcond, work)
Definition pdinvchk.f:3
pdinvdriver
program pdinvdriver
Definition pdinvdriver.f:1
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
pdlange
double precision function pdlange(norm, m, n, a, ia, ja, desca, work)
Definition pdlange.f:3
pdlansy
double precision function pdlansy(norm, uplo, n, a, ia, ja, desca, work)
Definition pdlansy.f:3
pdlantr
double precision function pdlantr(norm, uplo, diag, m, n, a, ia, ja, desca, work)
Definition pdlantr.f:3
pdpotrf
subroutine pdpotrf(uplo, n, a, ia, ja, desca, info)
Definition pdpotrf.f:2
pdpotri
subroutine pdpotri(uplo, n, a, ia, ja, desca, info)
Definition pdpotri.f:2
pdtrtri
subroutine pdtrtri(uplo, diag, n, a, ia, ja, desca, info)
Definition pdtrtri.f:2
slboot
subroutine slboot()
Definition sltimer.f:2
sltimer
subroutine sltimer(i)
Definition sltimer.f:47
slcombine
subroutine slcombine(ictxt, scope, op, timetype, n, ibeg, times)
Definition sltimer.f:267