SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
All Classes Files Functions Variables Typedefs Macros
pztrddriver.f
Go to the documentation of this file.
1 PROGRAM pztrddriver
2*
3* -- ScaLAPACK testing driver (version 1.7) --
4* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
5* and University of California, Berkeley.
6* October 15, 1999
7*
8* Purpose
9* ========
10*
11* PZTRDDRIVER is the main test program for the COMPLEX*16
12* SCALAPACK TRD (symmetric tridiagonal reduction) routines.
13*
14* The program must be driven by a short data file. An annotated
15* example of a data file can be obtained by deleting the first 3
16* characters from the following 13 lines:
17* 'ScaLAPACK TRD computation input file'
18* 'PVM machine'
19* 'TRD.out' output file name
20* 6 device out
21* 'L' define Lower or Upper
22* 3 number of problems sizes
23* 5 31 201 values of N
24* 3 number of NB's
25* 2 3 5 values of NB
26* 7 number of process grids (ordered pairs of P & Q)
27* 1 2 1 4 2 3 8 values of P
28* 1 2 4 1 3 2 1 values of Q
29* 1.0 threshold
30*
31* Internal Parameters
32* ===================
33*
34* TOTMEM INTEGER, default = 2000000
35* TOTMEM is a machine-specific parameter indicating the
36* maximum amount of available memory in bytes.
37* The user should customize TOTMEM to his platform. Remember
38* to leave room in memory for the operating system, the BLACS
39* buffer, etc. For example, on a system with 8 MB of memory
40* per process (e.g., one processor on an Intel iPSC/860), the
41* parameters we use are TOTMEM=6200000 (leaving 1.8 MB for OS,
42* code, BLACS buffer, etc). However, for PVM, we usually set
43* TOTMEM = 2000000. Some experimenting with the maximum value
44* of TOTMEM may be required.
45*
46* INTGSZ INTEGER, default = 4 bytes.
47* ZPLXSZ INTEGER, default = 16 bytes.
48* INTGSZ and ZPLXSZ indicate the length in bytes on the
49* given platform for an integer and a double precision
50* complex.
51* MEM COMPLEX*16 array, dimension ( TOTMEM / ZPLXSZ )
52*
53* All arrays used by SCALAPACK routines are allocated from
54* this array and referenced by pointers. The integer IPA,
55* for example, is a pointer to the starting element of MEM for
56* the matrix A.
57*
58* =====================================================================
59*
60* .. Parameters ..
61 INTEGER block_cyclic_2d, csrc_, ctxt_, dlen_, dtype_,
62 $ lld_, mb_, m_, nb_, n_, rsrc_
63 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
64 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
65 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
66 INTEGER dblesz, totmem, zplxsz, memsiz, ntests
67 COMPLEX*16 padval
68 parameter( dblesz = 8, totmem = 2000000, zplxsz = 16,
69 $ memsiz = totmem / zplxsz, ntests = 20,
70 $ padval = ( -9923.0d+0, -9923.0d+0 ) )
71* ..
72* .. Local Scalars ..
73 LOGICAL check
74 CHARACTER uplo
75 CHARACTER*6 passed
76 CHARACTER*80 outfile
77 INTEGER i, iam, iaseed, ictxt, imidpad, info, ipa, ipd,
78 $ ipe, ipostpad, iprepad, ipt, ipw, itemp, j, k,
79 $ kfail, kpass, kskip, ktests, lcm, lwork, mycol,
80 $ myrow, n, nb, ndiag, ngrids, nmat, nnb, noffd,
81 $ nout, np, npcol, nprocs, nprow, nq, worksiz,
82 $ worktrd
83 REAL thresh
84 DOUBLE PRECISION anorm, fresid, nops, tmflops
85* ..
86* .. Local Arrays ..
87 INTEGER desca( dlen_ ), ierr( 1 ), nbval( ntests ),
88 $ nval( ntests ), pval( ntests ), qval( ntests )
89 DOUBLE PRECISION ctime( 1 ), wtime( 1 )
90 COMPLEX*16 mem( memsiz )
91* ..
92* .. External Subroutines ..
93 EXTERNAL blacs_barrier, blacs_exit, blacs_get,
94 $ blacs_gridexit, blacs_gridinfo, blacs_gridinit,
95 $ blacs_pinfo, descinit, igsum2d, pzchekpad,
99* ..
100* .. External Functions ..
101 LOGICAL lsame
102 INTEGER iceil, ilcm, numroc
103 DOUBLE PRECISION pzlanhe
104 EXTERNAL lsame, iceil, ilcm, numroc, pzlanhe
105* ..
106* .. Intrinsic Functions ..
107 INTRINSIC dble, max
108* ..
109* .. Data statements ..
110 DATA ktests, kpass, kfail, kskip / 4*0 /
111* ..
112* .. Executable Statements ..
113*
114 IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
115 $ rsrc_.LT.0 )stop
116* Get starting information
117*
118 CALL blacs_pinfo( iam, nprocs )
119 iaseed = 100
120 CALL pztrdinfo( outfile, nout, uplo, nmat, nval, ntests, nnb,
121 $ nbval, ntests, ngrids, pval, ntests, qval, ntests,
122 $ thresh, mem, iam, nprocs )
123 check = ( thresh.GE.0.0e+0 )
124*
125* Print headings
126*
127 IF( iam.EQ.0 ) THEN
128 WRITE( nout, fmt = * )
129 WRITE( nout, fmt = 9995 )
130 WRITE( nout, fmt = 9994 )
131 WRITE( nout, fmt = * )
132 END IF
133*
134* Loop over different process grids
135*
136 DO 30 i = 1, ngrids
137*
138 nprow = pval( i )
139 npcol = qval( i )
140*
141* Make sure grid information is correct
142*
143 ierr( 1 ) = 0
144 IF( nprow.LT.1 ) THEN
145 IF( iam.EQ.0 )
146 $ WRITE( nout, fmt = 9999 )'GRID', 'nprow', nprow
147 ierr( 1 ) = 1
148 ELSE IF( npcol.LT.1 ) THEN
149 IF( iam.EQ.0 )
150 $ WRITE( nout, fmt = 9999 )'GRID', 'npcol', npcol
151 ierr( 1 ) = 1
152 ELSE IF( nprow*npcol.GT.nprocs ) THEN
153 IF( iam.EQ.0 )
154 $ WRITE( nout, fmt = 9998 )nprow*npcol, nprocs
155 ierr( 1 ) = 1
156 END IF
157*
158 IF( ierr( 1 ).GT.0 ) THEN
159 IF( iam.EQ.0 )
160 $ WRITE( nout, fmt = 9997 )'grid'
161 kskip = kskip + 1
162 GO TO 30
163 END IF
164*
165* Define process grid
166*
167 CALL blacs_get( -1, 0, ictxt )
168 CALL blacs_gridinit( ictxt, 'Row-major', nprow, npcol )
169 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
170*
171* Go to bottom of loop if this case doesn't use my process
172*
173 IF( myrow.GE.nprow .OR. mycol.GE.npcol )
174 $ GO TO 30
175*
176 DO 20 j = 1, nmat
177*
178 n = nval( j )
179*
180* Make sure matrix information is correct
181*
182 ierr( 1 ) = 0
183 IF( n.LT.1 ) THEN
184 IF( iam.EQ.0 )
185 $ WRITE( nout, fmt = 9999 )'MATRIX', 'N', n
186 ierr( 1 ) = 1
187 END IF
188*
189* Make sure no one had error
190*
191 CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1, 0 )
192*
193 IF( ierr( 1 ).GT.0 ) THEN
194 IF( iam.EQ.0 )
195 $ WRITE( nout, fmt = 9997 )'matrix'
196 kskip = kskip + 1
197 GO TO 20
198 END IF
199*
200* Loop over different blocking sizes
201*
202 DO 10 k = 1, nnb
203*
204 nb = nbval( k )
205*
206* Make sure nb is legal
207*
208 ierr( 1 ) = 0
209 IF( nb.LT.1 ) THEN
210 ierr( 1 ) = 1
211 IF( iam.EQ.0 )
212 $ WRITE( nout, fmt = 9999 )'NB', 'NB', nb
213 END IF
214*
215* Check all processes for an error
216*
217 CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1, 0 )
218*
219 IF( ierr( 1 ).GT.0 ) THEN
220 IF( iam.EQ.0 )
221 $ WRITE( nout, fmt = 9997 )'NB'
222 kskip = kskip + 1
223 GO TO 10
224 END IF
225*
226* Padding constants
227*
228 np = numroc( n, nb, myrow, 0, nprow )
229 nq = numroc( n, nb, mycol, 0, npcol )
230 IF( check ) THEN
231 iprepad = max( nb, np )
232 imidpad = nb
233 ipostpad = max( nb, nq )
234 ELSE
235 iprepad = 0
236 imidpad = 0
237 ipostpad = 0
238 END IF
239*
240* Initialize the array descriptor for the matrix A
241*
242 CALL descinit( desca, n, n, nb, nb, 0, 0, ictxt,
243 $ max( 1, np )+imidpad, ierr( 1 ) )
244*
245* Check all processes for an error
246*
247 CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1, 0 )
248*
249 IF( ierr( 1 ).LT.0 ) THEN
250 IF( iam.EQ.0 )
251 $ WRITE( nout, fmt = 9997 )'descriptor'
252 kskip = kskip + 1
253 GO TO 10
254 END IF
255*
256* Assign pointers into MEM for SCALAPACK arrays, A is
257* allocated starting at position MEM( IPREPAD+1 )
258*
259 ndiag = nq
260 IF( lsame( uplo, 'U' ) ) THEN
261 noffd = nq
262 ELSE
263 noffd = numroc( n-1, nb, mycol, 0, npcol )
264 END IF
265 ndiag = iceil( dblesz*ndiag, zplxsz )
266 noffd = iceil( dblesz*noffd, zplxsz )
267*
268 ipa = iprepad + 1
269 ipd = ipa + desca( lld_ )*nq + ipostpad + iprepad
270 ipe = ipd + ndiag + ipostpad + iprepad
271 ipt = ipe + noffd + ipostpad + iprepad
272 ipw = ipt + nq + ipostpad + iprepad
273*
274* Calculate the amount of workspace required for the
275* reduction
276*
277 lwork = max( nb*( np+1 ), 3*nb )
278 worktrd = lwork + ipostpad
279 worksiz = worktrd
280*
281* Figure the amount of workspace required by the check
282*
283 IF( check ) THEN
284 itemp = 2*nq + np
285 IF( nprow.NE.npcol ) THEN
286 lcm = ilcm( nprow, npcol )
287 itemp = nb*iceil( iceil( np, nb ), lcm / nprow ) +
288 $ itemp
289 END IF
290 itemp = max( iceil( dblesz*itemp, zplxsz ),
291 $ 2*( nb+np )*nb )
292 worksiz = max( lwork, itemp ) + ipostpad
293 END IF
294*
295* Check for adequate memory for problem size
296*
297 ierr( 1 ) = 0
298 IF( ipw+worksiz.GT.memsiz ) THEN
299 IF( iam.EQ.0 )
300 $ WRITE( nout, fmt = 9996 )'Tridiagonal reduction',
301 $ ( ipw+worksiz )*zplxsz
302 ierr( 1 ) = 1
303 END IF
304*
305* Check all processes for an error
306*
307 CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1, 0 )
308*
309 IF( ierr( 1 ).GT.0 ) THEN
310 IF( iam.EQ.0 )
311 $ WRITE( nout, fmt = 9997 )'MEMORY'
312 kskip = kskip + 1
313 GO TO 10
314 END IF
315*
316* Generate the matrix A
317*
318 CALL pzmatgen( ictxt, 'Hemm', 'N', desca( m_ ),
319 $ desca( n_ ), desca( mb_ ), desca( nb_ ),
320 $ mem( ipa ), desca( lld_ ), desca( rsrc_ ),
321 $ desca( csrc_ ), iaseed, 0, np, 0, nq,
322 $ myrow, mycol, nprow, npcol )
323*
324* Need Infinity-norm of A for checking
325*
326 IF( check ) THEN
327 CALL pzfillpad( ictxt, np, nq, mem( ipa-iprepad ),
328 $ desca( lld_ ), iprepad, ipostpad,
329 $ padval )
330 CALL pzfillpad( ictxt, ndiag, 1, mem( ipd-iprepad ),
331 $ ndiag, iprepad, ipostpad, padval )
332 CALL pzfillpad( ictxt, noffd, 1, mem( ipe-iprepad ),
333 $ noffd, iprepad, ipostpad, padval )
334 CALL pzfillpad( ictxt, nq, 1, mem( ipt-iprepad ), nq,
335 $ iprepad, ipostpad, padval )
336 CALL pzfillpad( ictxt, worksiz-ipostpad, 1,
337 $ mem( ipw-iprepad ), worksiz-ipostpad,
338 $ iprepad, ipostpad, padval )
339 anorm = pzlanhe( 'I', uplo, n, mem( ipa ), 1, 1,
340 $ desca, mem( ipw ) )
341 CALL pzchekpad( ictxt, 'PZLANHE', np, nq,
342 $ mem( ipa-iprepad ), desca( lld_ ),
343 $ iprepad, ipostpad, padval )
344 CALL pzchekpad( ictxt, 'PZLANHE', worksiz-ipostpad, 1,
345 $ mem( ipw-iprepad ), worksiz-ipostpad,
346 $ iprepad, ipostpad, padval )
347 CALL pzfillpad( ictxt, worktrd-ipostpad, 1,
348 $ mem( ipw-iprepad ), worktrd-ipostpad,
349 $ iprepad, ipostpad, padval )
350 END IF
351*
352 CALL slboot
353 CALL blacs_barrier( ictxt, 'All' )
354 CALL sltimer( 1 )
355*
356* Reduce to symmetric tridiagonal form
357*
358 CALL pzhetrd( uplo, n, mem( ipa ), 1, 1, desca,
359 $ mem( ipd ), mem( ipe ), mem( ipt ),
360 $ mem( ipw ), lwork, info )
361*
362 CALL sltimer( 1 )
363*
364 IF( check ) THEN
365*
366* Check for memory overwrite
367*
368 CALL pzchekpad( ictxt, 'PZHETRD', np, nq,
369 $ mem( ipa-iprepad ), desca( lld_ ),
370 $ iprepad, ipostpad, padval )
371 CALL pzchekpad( ictxt, 'PZHETRD', ndiag, 1,
372 $ mem( ipd-iprepad ), ndiag, iprepad,
373 $ ipostpad, padval )
374 CALL pzchekpad( ictxt, 'PZHETRD', noffd, 1,
375 $ mem( ipe-iprepad ), noffd, iprepad,
376 $ ipostpad, padval )
377 CALL pzchekpad( ictxt, 'PZHETRD', nq, 1,
378 $ mem( ipt-iprepad ), nq, iprepad,
379 $ ipostpad, padval )
380 CALL pzchekpad( ictxt, 'PZHETRD', worktrd-ipostpad, 1,
381 $ mem( ipw-iprepad ), worktrd-ipostpad,
382 $ iprepad, ipostpad, padval )
383 CALL pzfillpad( ictxt, worksiz-ipostpad, 1,
384 $ mem( ipw-iprepad ), worksiz-ipostpad,
385 $ iprepad, ipostpad, padval )
386*
387* Compute fctres = ||A - QTQ'|| / (||A|| * N * eps)
388*
389 CALL pzhetdrv( uplo, n, mem( ipa ), 1, 1, desca,
390 $ mem( ipd ), mem( ipe ), mem( ipt ),
391 $ mem( ipw ), ierr( 1 ) )
392 CALL pzlafchk( 'Hemm', 'No', n, n, mem( ipa ), 1, 1,
393 $ desca, iaseed, anorm, fresid,
394 $ mem( ipw ) )
395*
396* Check for memory overwrite
397*
398 CALL pzchekpad( ictxt, 'PZHETDRV', np, nq,
399 $ mem( ipa-iprepad ), desca( lld_ ),
400 $ iprepad, ipostpad, padval )
401 CALL pzchekpad( ictxt, 'PZHETDRV', ndiag, 1,
402 $ mem( ipd-iprepad ), ndiag, iprepad,
403 $ ipostpad, padval )
404 CALL pzchekpad( ictxt, 'PZHETDRV', noffd, 1,
405 $ mem( ipe-iprepad ), noffd, iprepad,
406 $ ipostpad, padval )
407 CALL pzchekpad( ictxt, 'PZHETDRV', worksiz-ipostpad,
408 $ 1, mem( ipw-iprepad ),
409 $ worksiz-ipostpad, iprepad, ipostpad,
410 $ padval )
411*
412* Test residual and detect NaN result
413*
414 IF( fresid.LE.thresh .AND. fresid-fresid.EQ.
415 $ 0.0d+0 .AND. ierr( 1 ).EQ.0 ) THEN
416 kpass = kpass + 1
417 passed = 'PASSED'
418 ELSE
419 IF( myrow.EQ.0 .AND. mycol.EQ.0 )
420 $ WRITE( nout, fmt = 9986 )fresid
421 kfail = kfail + 1
422 passed = 'FAILED'
423 END IF
424*
425 IF( myrow.EQ.0 .AND. mycol.EQ.0 .AND. ierr( 1 ).NE.0 )
426 $ WRITE( nout, fmt = * )'D or E copies incorrect ...'
427 ELSE
428*
429* Don't perform the checking, only the timing operation
430*
431 kpass = kpass + 1
432 fresid = fresid - fresid
433 passed = 'BYPASS'
434 END IF
435*
436* Gather maximum of all CPU and WALL clock timings
437*
438 CALL slcombine( ictxt, 'All', '>', 'W', 1, 1, wtime )
439 CALL slcombine( ictxt, 'All', '>', 'C', 1, 1, ctime )
440*
441* Print results
442*
443 IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
444*
445* TRD requires 16/3 N^3 floating point operations
446*
447 nops = dble( n )
448*
449 nops = ( 4.0d+0 / 3.0d+0 )*nops**3
450 nops = nops / 1.0d+6
451*
452* Print WALL time
453*
454 IF( wtime( 1 ).GT.0.0d+0 ) THEN
455 tmflops = nops / wtime( 1 )
456 ELSE
457 tmflops = 0.0d+0
458 END IF
459 IF( wtime( 1 ).GE.0.0d+0 )
460 $ WRITE( nout, fmt = 9993 )'WALL', uplo, n, nb,
461 $ nprow, npcol, wtime( 1 ), tmflops, fresid, passed
462*
463* Print CPU time
464*
465 IF( ctime( 1 ).GT.0.0d+0 ) THEN
466 tmflops = nops / ctime( 1 )
467 ELSE
468 tmflops = 0.0d+0
469 END IF
470 IF( ctime( 1 ).GE.0.0d+0 )
471 $ WRITE( nout, fmt = 9993 )'CPU ', uplo, n, nb,
472 $ nprow, npcol, ctime( 1 ), tmflops, fresid, passed
473 END IF
474 10 CONTINUE
475 20 CONTINUE
476*
477 CALL blacs_gridexit( ictxt )
478 30 CONTINUE
479*
480 CALL pzttrdtester( iam, nprocs, check, nout, thresh, nval, nmat,
481 $ mem, totmem, kpass, kfail, kskip )
482*
483* Print ending messages and close output file
484*
485 IF( iam.EQ.0 ) THEN
486 ktests = kpass + kfail + kskip
487 WRITE( nout, fmt = * )
488 WRITE( nout, fmt = 9992 )ktests
489 IF( check ) THEN
490 WRITE( nout, fmt = 9991 )kpass
491 WRITE( nout, fmt = 9989 )kfail
492 ELSE
493 WRITE( nout, fmt = 9990 )kpass
494 END IF
495 WRITE( nout, fmt = 9988 )kskip
496 WRITE( nout, fmt = * )
497 WRITE( nout, fmt = * )
498 WRITE( nout, fmt = 9987 )
499 IF( nout.NE.6 .AND. nout.NE.0 )
500 $ CLOSE ( nout )
501 END IF
502*
503 CALL blacs_exit( 0 )
504*
505 9999 FORMAT( 'ILLEGAL ', a6, ': ', a5, ' = ', i3,
506 $ '; It should be at least 1' )
507 9998 FORMAT( 'ILLEGAL GRID: nprow*npcol = ', i4, '. It can be at most',
508 $ i4 )
509 9997 FORMAT( 'Bad ', a6, ' parameters: going on to next test case.' )
510 9996 FORMAT( 'Unable to perform ', a, ': need TOTMEM of at least',
511 $ i11 )
512 9995 FORMAT( 'TIME UPLO N NB P Q TRD Time ',
513 $ ' MFLOPS Residual CHECK' )
514 9994 FORMAT( '---- ---- ------ --- ----- ----- --------- ',
515 $ '----------- -------- ------' )
516 9993 FORMAT( a4, 1x, a4, 1x, i6, 1x, i3, 1x, i5, 1x, i5, 1x, f9.2, 1x,
517 $ f11.2, 1x, f8.2, 1x, a6 )
518 9992 FORMAT( 'Finished', i4, ' tests, with the following results:' )
519 9991 FORMAT( i5, ' tests completed and passed residual checks.' )
520 9990 FORMAT( i5, ' tests completed without checking.' )
521 9989 FORMAT( i5, ' tests completed and failed residual checks.' )
522 9988 FORMAT( i5, ' tests skipped because of illegal input values.' )
523 9987 FORMAT( 'END OF TESTS.' )
524 9986 FORMAT( '||A - Q*T*Q''|| / (||A|| * N * eps) = ', g25.7 )
525*
526 stop
527*
528* End of PZTRDDRIVER
529*
530 END
subroutine pzlafchk(aform, diag, m, n, a, ia, ja, desca, iaseed, anorm, fresid, work)
Definition pzlafchk.f:3
subroutine pzmatgen(ictxt, aform, diag, m, n, mb, nb, a, lda, iarow, iacol, iseed, iroff, irnum, icoff, icnum, myrow, mycol, nprow, npcol)
Definition pzmatgen.f:4
subroutine descinit(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld, info)
Definition descinit.f:3
integer function iceil(inum, idenom)
Definition iceil.f:2
integer function ilcm(m, n)
Definition ilcm.f:2
integer function numroc(n, nb, iproc, isrcproc, nprocs)
Definition numroc.f:2
#define max(A, B)
Definition pcgemr.c:180
subroutine pzchekpad(ictxt, mess, m, n, a, lda, ipre, ipost, chkval)
Definition pzchekpad.f:3
subroutine pzfillpad(ictxt, m, n, a, lda, ipre, ipost, chkval)
Definition pzfillpad.f:2
subroutine pzhetdrv(uplo, n, a, ia, ja, desca, d, e, tau, work, info)
Definition pzhetdrv.f:3
subroutine pzhetrd(uplo, n, a, ia, ja, desca, d, e, tau, work, lwork, info)
Definition pzhetrd.f:3
double precision function pzlanhe(norm, uplo, n, a, ia, ja, desca, work)
Definition pzlanhe.f:3
program pztrddriver
Definition pztrddriver.f:1
subroutine pztrdinfo(summry, nout, uplo, nmat, nval, ldnval, nnb, nbval, ldnbval, ngrids, pval, ldpval, qval, ldqval, thresh, work, iam, nprocs)
Definition pztrdinfo.f:4
subroutine pzttrdtester(iam, nprocs, check, nout, thresh, nval, nmat, mem, totmem, kpass, kfail, kskip)
Definition pzttrdtester.f:3
subroutine slboot()
Definition sltimer.f:2
subroutine sltimer(i)
Definition sltimer.f:47
subroutine slcombine(ictxt, scope, op, timetype, n, ibeg, times)
Definition sltimer.f:267
logical function lsame(ca, cb)
Definition tools.f:1724