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