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