SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
pchrddriver.f
Go to the documentation of this file.
1 PROGRAM pchrddriver
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* PCHRDDRIVER is the main test program for the COMPLEX
12* ScaLAPACK HRD (Hessenberg 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 14 lines:
17* 'ScaLAPACK HRD input file'
18* 'PVM machine'
19* 'HRD.out' output file name (if any)
20* 6 device out
21* 2 number of problems sizes
22* 100 101 values of N
23* 2 1 values of ILO
24* 99 101 values of IHI
25* 3 number of NB's
26* 2 3 5 values of NB
27* 7 number of process grids (ordered pairs of P & Q)
28* 1 2 1 4 2 3 8 values of P
29* 1 2 4 1 3 2 1 values of Q
30* 3.0 threshold
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
68 COMPLEX padval
69 parameter( cplxsz = 8, totmem = 2000000,
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, ihi, ihip, ihlp, ihlq,
78 $ ilcol, ilo, ilrow, info, inlq, imidpad, ipa,
79 $ ipt, ipw, ipostpad, iprepad, itemp, j, k,
80 $ kfail, kpass, kskip, ktests, lcm, lcmq, loff,
81 $ lwork, mycol, myrow, n, nb, ngrids, nmat, nnb,
82 $ nprocs, nout, np, npcol, nprow, nq, workhrd,
83 $ worksiz
84 REAL anorm, fresid, thresh
85 DOUBLE PRECISION nops, tmflops
86* ..
87* .. Local Arrays ..
88 INTEGER desca( dlen_ ), ierr( 1 ), nbval( ntests ),
89 $ nval( ntests ), nvhi( ntests ), nvlo( ntests ),
90 $ pval( ntests ), qval( ntests )
91 DOUBLE PRECISION ctime( 1 ), wtime( 1 )
92 COMPLEX mem( memsiz )
93* ..
94* .. External Subroutines ..
95 EXTERNAL blacs_barrier, blacs_exit, blacs_get,
96 $ blacs_gridexit, blacs_gridinit, blacs_gridinfo,
97 $ descinit, igsum2d, blacs_pinfo, pcfillpad,
98 $ pclafchk, pcgehdrv, pcgehrd,
99 $ pchrdinfo, pcmatgen, slboot,
100 $ slcombine, sltimer
101* ..
102* .. External Functions ..
103 INTEGER ilcm, indxg2p, numroc
104 REAL pclange
105 EXTERNAL ilcm, indxg2p, numroc, pclange
106* ..
107* .. Intrinsic Functions ..
108 INTRINSIC dble, max
109* ..
110* .. Data statements ..
111 DATA ktests, kpass, kfail, kskip / 4*0 /
112* ..
113* .. Executable Statements ..
114*
115* Get starting information
116*
117 CALL blacs_pinfo( iam, nprocs )
118 iaseed = 100
119 CALL pchrdinfo( outfile, nout, nmat, nval, nvlo, nvhi, ntests,
120 $ nnb, nbval, ntests, ngrids, pval, ntests, qval,
121 $ ntests, 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 $ GOTO 30
174*
175 DO 20 j = 1, nmat
176*
177 n = nval( j )
178 ilo = nvlo( j )
179 ihi = nvhi( j )
180*
181* Make sure matrix information is correct
182*
183 ierr( 1 ) = 0
184 IF( n.LT.1 ) THEN
185 IF( iam.EQ.0 )
186 $ WRITE( nout, fmt = 9999 ) 'MATRIX', 'N', n
187 ierr( 1 ) = 1
188 END IF
189*
190* Check all processes for an error
191*
192 CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1, 0 )
193*
194 IF( ierr( 1 ).GT.0 ) THEN
195 IF( iam.EQ.0 )
196 $ WRITE( nout, fmt = 9997 ) 'matrix'
197 kskip = kskip + 1
198 GO TO 20
199 END IF
200*
201 DO 10 k = 1, nnb
202 nb = nbval( k )
203*
204* Make sure nb is legal
205*
206 ierr( 1 ) = 0
207 IF( nb.LT.1 ) THEN
208 ierr( 1 ) = 1
209 IF( iam.EQ.0 )
210 $ WRITE( nout, fmt = 9999 ) 'NB', 'NB', nb
211 END IF
212*
213* Check all processes for an error
214*
215 CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1, 0 )
216*
217 IF( ierr( 1 ).GT.0 ) THEN
218 IF( iam.EQ.0 )
219 $ WRITE( nout, fmt = 9997 ) 'NB'
220 kskip = kskip + 1
221 GO TO 10
222 END IF
223*
224 np = numroc( n, nb, myrow, 0, nprow )
225 nq = numroc( n, nb, mycol, 0, npcol )
226 IF( check ) THEN
227 iprepad = max( nb, np )
228 imidpad = nb
229 ipostpad = max( nb, nq )
230 ELSE
231 iprepad = 0
232 imidpad = 0
233 ipostpad = 0
234 END IF
235*
236* Initialize the array descriptor for the matrix A
237*
238 CALL descinit( desca, n, n, nb, nb, 0, 0, ictxt,
239 $ max( 1, np ) + imidpad, info )
240*
241* Check all processes for an error
242*
243 CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1, 0 )
244*
245 IF( ierr( 1 ).LT.0 ) THEN
246 IF( iam.EQ.0 )
247 $ WRITE( nout, fmt = 9997 ) 'descriptor'
248 kskip = kskip + 1
249 GO TO 10
250 END IF
251*
252* Assign pointers into MEM for SCALAPACK arrays, A is
253* allocated starting at position MEM( IPREPAD+1 )
254*
255 ipa = iprepad + 1
256 ipt = ipa + desca( lld_ )*nq + ipostpad + iprepad
257 ipw = ipt + nq + ipostpad + iprepad
258*
259* Calculate the amount of workspace required for the
260* reduction
261*
262 ihip = numroc( ihi, nb, myrow, desca( rsrc_ ), nprow )
263 loff = mod( ilo-1, nb )
264 ilrow = indxg2p( ilo, nb, myrow, desca( rsrc_ ), nprow )
265 ilcol = indxg2p( ilo, nb, mycol, desca( csrc_ ), npcol )
266 ihlp = numroc( ihi-ilo+loff+1, nb, myrow, ilrow, nprow )
267 inlq = numroc( n-ilo+loff+1, nb, mycol, ilcol, npcol )
268 lwork = nb*( nb + max( ihip+1, ihlp+inlq ) )
269 workhrd = lwork + ipostpad
270 worksiz = workhrd
271*
272* Figure the amount of workspace required by the check
273*
274 IF( check ) THEN
275 lcm = ilcm( nprow, npcol )
276 lcmq = lcm / npcol
277 ihlq = numroc( ihi-ilo+loff+1, nb, mycol, ilcol,
278 $ npcol )
279 itemp = nb*max( ihlp+inlq, ihlq+max( ihip,
280 $ ihlp+numroc( numroc( ihi-ilo+loff+1, nb, 0, 0,
281 $ npcol ), nb, 0, 0, lcmq ) ) )
282 worksiz = max( nb*nb + nb*ihlp + itemp, nb * np ) +
283 $ ipostpad
284 END IF
285*
286* Check for adequate memory for problem size
287*
288 ierr( 1 ) = 0
289 IF( ipw+worksiz.GT.memsiz ) THEN
290 IF( iam.EQ.0 )
291 $ WRITE( nout, fmt = 9996 ) 'Hessenberg reduction',
292 $ ( ipw+worksiz )*cplxsz
293 ierr( 1 ) = 1
294 END IF
295*
296* Check all processes for an error
297*
298 CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1, 0 )
299*
300 IF( ierr( 1 ).GT.0 ) THEN
301 IF( iam.EQ.0 )
302 $ WRITE( nout, fmt = 9997 ) 'MEMORY'
303 kskip = kskip + 1
304 GO TO 10
305 END IF
306*
307* Generate A
308*
309 CALL pcmatgen( ictxt, 'No', 'No', desca( m_ ),
310 $ desca( n_ ), desca( mb_ ), desca( nb_ ),
311 $ mem( ipa ), desca( lld_ ), desca( rsrc_ ),
312 $ desca( csrc_ ),
313 $ iaseed, 0, np, 0, nq, myrow, mycol,
314 $ nprow, npcol )
315*
316* Need Infinity-norm of A for checking
317*
318 IF( check ) THEN
319 CALL pcfillpad( ictxt, np, nq, mem( ipa-iprepad ),
320 $ desca( lld_ ), iprepad, ipostpad,
321 $ padval )
322 CALL pcfillpad( ictxt, nq, 1, mem( ipt-iprepad ),
323 $ nq, iprepad, ipostpad, padval )
324 CALL pcfillpad( ictxt, worksiz-ipostpad, 1,
325 $ mem( ipw-iprepad ), worksiz-ipostpad,
326 $ iprepad, ipostpad, padval )
327 anorm = pclange( 'I', n, n, mem( ipa ), 1, 1, desca,
328 $ mem( ipw ) )
329 CALL pcchekpad( ictxt, 'PCLANGE', np, nq,
330 $ mem( ipa-iprepad ), desca( lld_ ),
331 $ iprepad, ipostpad, padval )
332 CALL pcchekpad( ictxt, 'PCLANGE',
333 $ worksiz-ipostpad, 1,
334 $ mem( ipw-iprepad ), worksiz-ipostpad,
335 $ iprepad, ipostpad, padval )
336 CALL pcfillpad( ictxt, workhrd-ipostpad, 1,
337 $ mem( ipw-iprepad ), workhrd-ipostpad,
338 $ iprepad, ipostpad, padval )
339 END IF
340*
341 CALL slboot()
342 CALL blacs_barrier( ictxt, 'All' )
343 CALL sltimer( 1 )
344*
345* Reduce Hessenberg form
346*
347 CALL pcgehrd( n, ilo, ihi, mem( ipa ), 1, 1, desca,
348 $ mem( ipt ), mem( ipw ), lwork, info )
349 CALL sltimer( 1 )
350*
351 IF( check ) THEN
352*
353* Check for memory overwrite
354*
355 CALL pcchekpad( ictxt, 'PCGEHRD', np, nq,
356 $ mem( ipa-iprepad ), desca( lld_ ),
357 $ iprepad, ipostpad, padval )
358 CALL pcchekpad( ictxt, 'PCGEHRD', nq, 1,
359 $ mem( ipt-iprepad ), nq, iprepad,
360 $ ipostpad, padval )
361 CALL pcchekpad( ictxt, 'PCGEHRD', workhrd-ipostpad,
362 $ 1, mem( ipw-iprepad ),
363 $ workhrd-ipostpad, iprepad,
364 $ ipostpad, padval )
365 CALL pcfillpad( ictxt, worksiz-ipostpad, 1,
366 $ mem( ipw-iprepad ), worksiz-ipostpad,
367 $ iprepad, ipostpad, padval )
368*
369* Compute fctres = ||A - Q H Q'|| / (||A||*N*eps)
370*
371 CALL pcgehdrv( n, ilo, ihi, mem( ipa ), 1, 1, desca,
372 $ mem( ipt ), mem( ipw ) )
373 CALL pclafchk( 'No', 'No', n, n, mem( ipa ), 1, 1,
374 $ desca, iaseed, anorm, fresid,
375 $ mem( ipw ) )
376*
377* Check for memory overwrite
378*
379 CALL pcchekpad( ictxt, 'PCGEHDRV', np, nq,
380 $ mem( ipa-iprepad ), desca( lld_ ),
381 $ iprepad, ipostpad, padval )
382 CALL pcchekpad( ictxt, 'PCGEHDRV', nq, 1,
383 $ mem( ipt-iprepad ), nq, iprepad,
384 $ ipostpad, padval )
385 CALL pcchekpad( ictxt, 'PCGEHDRV',
386 $ worksiz-ipostpad, 1,
387 $ mem( ipw-iprepad ), worksiz-ipostpad,
388 $ iprepad, ipostpad, padval )
389*
390* Test residual and detect NaN result
391*
392 IF( fresid.LE.thresh .AND. fresid-fresid.EQ.0.0e+0 )
393 $ THEN
394 kpass = kpass + 1
395 passed = 'PASSED'
396 ELSE
397 IF( myrow.EQ.0 .AND. mycol.EQ.0 )
398 $ WRITE( nout, fmt = 9986 ) fresid
399 kfail = kfail + 1
400 passed = 'FAILED'
401 END IF
402 ELSE
403*
404* Don't perform the checking, only the timing operation
405*
406 kpass = kpass + 1
407 fresid = fresid - fresid
408 passed = 'BYPASS'
409 END IF
410*
411* Gather max. of all CPU and WALL clock timings
412*
413 CALL slcombine( ictxt, 'All', '>', 'W', 1, 1, wtime )
414 CALL slcombine( ictxt, 'All', '>', 'C', 1, 1, ctime )
415*
416* Print results
417*
418 IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
419*
420* HRD requires 40/3 * N^3 floating point ops. (flops)
421* more precisely,
422* HRD requires 16/3*(IHI-ILO)^3+8*IHI*(IHI-ILO)^2 flops
423*
424 nops = dble( ihi-ilo )
425 nops = nops * nops *
426 $ ( 8.0d0*dble( ihi ) + (16.0d0/3.0d0)*nops )
427 nops = nops / 1.0d+6
428*
429* Print WALL time
430*
431 IF( wtime( 1 ).GT.0.0d+0 ) THEN
432 tmflops = nops / wtime( 1 )
433 ELSE
434 tmflops = 0.0d+0
435 END IF
436 IF( wtime( 1 ).GE.0.0d+0 )
437 $ WRITE( nout, fmt = 9993 ) 'WALL', n, ilo, ihi, nb,
438 $ nprow, npcol, wtime( 1 ), tmflops, fresid,
439 $ passed
440*
441* Print CPU time
442*
443 IF( ctime( 1 ).GT.0.0d+0 ) THEN
444 tmflops = nops / ctime( 1 )
445 ELSE
446 tmflops = 0.0d+0
447 END IF
448 IF( ctime( 1 ).GE.0.0d+0 )
449 $ WRITE( nout, fmt = 9993 ) 'CPU ', n, ilo, ihi, nb,
450 $ nprow, npcol, ctime( 1 ), tmflops, fresid,
451 $ passed
452 END IF
453 10 CONTINUE
454 20 CONTINUE
455*
456 CALL blacs_gridexit( ictxt )
457 30 CONTINUE
458*
459* Print ending messages and close output file
460*
461 IF( iam.EQ.0 ) THEN
462 ktests = kpass + kfail + kskip
463 WRITE( nout, fmt = * )
464 WRITE( nout, fmt = 9992 ) ktests
465 IF( check ) THEN
466 WRITE( nout, fmt = 9991 ) kpass
467 WRITE( nout, fmt = 9989 ) kfail
468 ELSE
469 WRITE( nout, fmt = 9990 ) kpass
470 END IF
471 WRITE( nout, fmt = 9988 ) kskip
472 WRITE( nout, fmt = * )
473 WRITE( nout, fmt = * )
474 WRITE( nout, fmt = 9987 )
475 IF( nout.NE.6 .AND. nout.NE.0 )
476 $ CLOSE( nout )
477 END IF
478*
479 CALL blacs_exit( 0 )
480*
481 9999 FORMAT( 'ILLEGAL ', a6, ': ', a5, ' = ', i3,
482 $ '; It should be at least 1' )
483 9998 FORMAT( 'ILLEGAL GRID: nprow*npcol = ', i4, '. It can be at most',
484 $ i4 )
485 9997 FORMAT( 'Bad ', a6, ' parameters: going on to next test case.' )
486 9996 FORMAT( 'Unable to perform ', a, ': need TOTMEM of at least',
487 $ i11 )
488 9995 FORMAT( 'TIME N ILO IHI NB P Q HRD Time ',
489 $ ' MFLOPS Residual CHECK' )
490 9994 FORMAT( '---- ------ ------ ------ --- ----- ----- --------- ',
491 $ '----------- -------- ------' )
492 9993 FORMAT( a4, 1x, i6, 1x, i6, 1x, i6, 1x, i3, 1x, i5, 1x, i5, 1x,
493 $ f9.2, 1x, f11.2, 1x, f8.2, 1x, a6 )
494 9992 FORMAT( 'Finished', i4, ' tests, with the following results:' )
495 9991 FORMAT( i5, ' tests completed and passed residual checks.' )
496 9990 FORMAT( i5, ' tests completed without checking.' )
497 9989 FORMAT( i5, ' tests completed and failed residual checks.' )
498 9988 FORMAT( i5, ' tests skipped because of illegal input values.' )
499 9987 FORMAT( 'END OF TESTS.' )
500 9986 FORMAT( '||A - Q*H*Q''|| / (||A|| * N * eps) = ', g25.7 )
501*
502 stop
503*
504* End of PCHRDDRIVER
505*
506 END
pclafchk
subroutine pclafchk(aform, diag, m, n, a, ia, ja, desca, iaseed, anorm, fresid, work)
Definition pclafchk.f:3
pcmatgen
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
descinit
subroutine descinit(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld, info)
Definition descinit.f:3
ilcm
integer function ilcm(m, n)
Definition ilcm.f:2
indxg2p
integer function indxg2p(indxglob, nb, iproc, isrcproc, nprocs)
Definition indxg2p.f:2
numroc
integer function numroc(n, nb, iproc, isrcproc, nprocs)
Definition numroc.f:2
pcchekpad
subroutine pcchekpad(ictxt, mess, m, n, a, lda, ipre, ipost, chkval)
Definition pcchekpad.f:3
pcfillpad
subroutine pcfillpad(ictxt, m, n, a, lda, ipre, ipost, chkval)
Definition pcfillpad.f:2
pcgehdrv
subroutine pcgehdrv(n, ilo, ihi, a, ia, ja, desca, tau, work)
Definition pcgehdrv.f:2
pcgehrd
subroutine pcgehrd(n, ilo, ihi, a, ia, ja, desca, tau, work, lwork, info)
Definition pcgehrd.f:3
max
#define max(A, B)
Definition pcgemr.c:180
pchrddriver
program pchrddriver
Definition pchrddriver.f:1
pchrdinfo
subroutine pchrdinfo(summry, nout, nmat, nval, nvlo, nvhi, ldnval, nnb, nbval, ldnbval, ngrids, pval, ldpval, qval, ldqval, thresh, work, iam, nprocs)
Definition pchrdinfo.f:5
pclange
real function pclange(norm, m, n, a, ia, ja, desca, work)
Definition pclange.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