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