pdtrddriver.f

Go to the documentation of this file.

      PROGRAM pdtrddriver
*
*  -- ScaLAPACK testing driver (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     October 15, 1999
*
*  Purpose
*  ========
*
*  PDTRDDRIVER is the main test program for the DOUBLE PRECISION
*  SCALAPACK TRD (symmetric tridiagonal reduction) routines.
*
*  The program must be driven by a short data file.  An annotated
*  example of a data file can be obtained by deleting the first 3
*  characters from the following 13 lines:
*  'ScaLAPACK TRD computation input file'
*  'PVM machine'
*  'TRD.out'       output file name
*  6               device out
*  'L'             define Lower or Upper
*  3               number of problems sizes
*  5 31 201        values of N
*  3               number of NB's
*  2 3 5           values of NB
*  7               number of process grids (ordered pairs of P & Q)
*  1 2 1 4 2 3 8   values of P
*  1 2 4 1 3 2 1   values of Q
*  1.0             threshold
*
*  Internal Parameters
*  ===================
*
*  TOTMEM   INTEGER, default = 2000000
*           TOTMEM is a machine-specific parameter indicating the
*           maximum amount of available memory in bytes.
*           The user should customize TOTMEM to his platform.  Remember
*           to leave room in memory for the operating system, the BLACS
*           buffer, etc.  For example, on a system with 8 MB of memory
*           per process (e.g., one processor on an Intel iPSC/860), the
*           parameters we use are TOTMEM=6200000 (leaving 1.8 MB for OS,
*           code, BLACS buffer, etc).  However, for PVM, we usually set
*           TOTMEM = 2000000.  Some experimenting with the maximum value
*           of TOTMEM may be required.
*
*  INTGSZ   INTEGER, default = 4 bytes.
*  DBLESZ   INTEGER, default = 8 bytes.
*           INTGSZ and DBLESZ indicate the length in bytes on the
*           given platform for an integer and a double precision real.
*  MEM      DOUBLE PRECISION array, dimension ( TOTMEM / DBLESZ )
*
*           All arrays used by SCALAPACK routines are allocated from
*           this array and referenced by pointers.  The integer IPA,
*           for example, is a pointer to the starting element of MEM for
*           the matrix A.
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            block_cyclic_2d, csrc_, ctxt_, dlen_, dtype_,
     $                   lld_, mb_, m_, nb_, n_, rsrc_
      parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
     $                     ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
     $                     rsrc_ = 7, csrc_ = 8, lld_ = 9 )
      INTEGER            dblesz, totmem, memsiz, ntests
      DOUBLE PRECISION   padval
      parameter( dblesz = 8, totmem = 2000000,
     $                   memsiz = totmem / dblesz, ntests = 20,
     $                   padval = -9923.0d+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            check
      CHARACTER          uplo
      CHARACTER*6        passed
      CHARACTER*80       outfile
      INTEGER            i, iam, iaseed, ictxt, imidpad, info, ipa, ipd,
     $                   ipe, ipostpad, iprepad, ipt, ipw, itemp, j, k,
     $                   kfail, kpass, kskip, ktests, lcm, lwork, mycol,
     $                   myrow, n, nb, ndiag, ngrids, nmat, nnb, noffd,
     $                   nout, np, npcol, nprocs, nprow, nq, worksiz,
     $                   worktrd
      REAL               thresh
      DOUBLE PRECISION   anorm, fresid, nops, tmflops
*     ..
*     .. Local Arrays ..
      INTEGER            desca( dlen_ ), ierr( 1 ), nbval( ntests ),
     $                   nval( ntests ), pval( ntests ), qval( ntests )
      DOUBLE PRECISION   ctime( 1 ), mem( memsiz ), wtime( 1 )
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_barrier, blacs_exit, blacs_get,
     $                   blacs_gridexit, blacs_gridinfo, blacs_gridinit,
     $                   blacs_pinfo, descinit, igsum2d, pdchekpad,
     $                   pdfillpad, pdlafchk, pdmatgen, pdsytdrv,
     $                   pdsytrd, pdtrdinfo, pdttrdtester, slboot,
     $                   slcombine, sltimer
*     ..
*     .. External Functions ..
      LOGICAL            lsame
      INTEGER            iceil, ilcm, numroc
      DOUBLE PRECISION   pdlansy
      EXTERNAL           lsame, iceil, ilcm, numroc, pdlansy
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          dble, max
*     ..
*     .. Data statements ..
      DATA               ktests, kpass, kfail, kskip / 4*0 /
*     ..
*     .. Executable Statements ..
*
      IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
     $    rsrc_.LT.0 )stop
*     Get starting information
*
      CALL blacs_pinfo( iam, nprocs )
      iaseed = 100
      CALL pdtrdinfo( outfile, nout, uplo, nmat, nval, ntests, nnb,
     $                nbval, ntests, ngrids, pval, ntests, qval, ntests,
     $                thresh, mem, iam, nprocs )
      check = ( thresh.GE.0.0e+0 )
*
*     Print headings
*
      IF( iam.EQ.0 ) THEN
         WRITE( nout, fmt = * )
         WRITE( nout, fmt = 9995 )
         WRITE( nout, fmt = 9994 )
         WRITE( nout, fmt = * )
      END IF
*
*     Loop over different process grids
*
      DO 30 i = 1, ngrids
*
         nprow = pval( i )
         npcol = qval( i )
*
*        Make sure grid information is correct
*
         ierr( 1 ) = 0
         IF( nprow.LT.1 ) THEN
            IF( iam.EQ.0 )
     $         WRITE( nout, fmt = 9999 )'GRID', 'nprow', nprow
            ierr( 1 ) = 1
         ELSE IF( npcol.LT.1 ) THEN
            IF( iam.EQ.0 )
     $         WRITE( nout, fmt = 9999 )'GRID', 'npcol', npcol
            ierr( 1 ) = 1
         ELSE IF( nprow*npcol.GT.nprocs ) THEN
            IF( iam.EQ.0 )
     $         WRITE( nout, fmt = 9998 )nprow*npcol, nprocs
            ierr( 1 ) = 1
         END IF
*
         IF( ierr( 1 ).GT.0 ) THEN
            IF( iam.EQ.0 )
     $         WRITE( nout, fmt = 9997 )'grid'
            kskip = kskip + 1
            GO TO 30
         END IF
*
*        Define process grid
*
         CALL blacs_get( -1, 0, ictxt )
         CALL blacs_gridinit( ictxt, 'Row-major', nprow, npcol )
         CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
*
*        Go to bottom of loop if this case doesn't use my process
*
         IF( myrow.GE.nprow .OR. mycol.GE.npcol )
     $      GO TO 30
*
         DO 20 j = 1, nmat
*
            n = nval( j )
*
*           Make sure matrix information is correct
*
            ierr( 1 ) = 0
            IF( n.LT.1 ) THEN
               IF( iam.EQ.0 )
     $            WRITE( nout, fmt = 9999 )'MATRIX', 'N', n
               ierr( 1 ) = 1
            END IF
*
*           Make sure no one had error
*
            CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1, 0 )
*
            IF( ierr( 1 ).GT.0 ) THEN
               IF( iam.EQ.0 )
     $            WRITE( nout, fmt = 9997 )'matrix'
               kskip = kskip + 1
               GO TO 20
            END IF
*
*           Loop over different blocking sizes
*
            DO 10 k = 1, nnb
*
               nb = nbval( k )
*
*              Make sure nb is legal
*
               ierr( 1 ) = 0
               IF( nb.LT.1 ) THEN
                  ierr( 1 ) = 1
                  IF( iam.EQ.0 )
     $               WRITE( nout, fmt = 9999 )'NB', 'NB', nb
               END IF
*
*              Check all processes for an error
*
               CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1, 0 )
*
               IF( ierr( 1 ).GT.0 ) THEN
                  IF( iam.EQ.0 )
     $               WRITE( nout, fmt = 9997 )'NB'
                  kskip = kskip + 1
                  GO TO 10
               END IF
*
*              Padding constants
*
               np = numroc( n, nb, myrow, 0, nprow )
               nq = numroc( n, nb, mycol, 0, npcol )
               IF( check ) THEN
                  iprepad = max( nb, np )
                  imidpad = nb
                  ipostpad = max( nb, nq )
               ELSE
                  iprepad = 0
                  imidpad = 0
                  ipostpad = 0
               END IF
*
*              Initialize the array descriptor for the matrix A
*
               CALL descinit( desca, n, n, nb, nb, 0, 0, ictxt,
     $                        max( 1, np )+imidpad, ierr( 1 ) )
*
*              Check all processes for an error
*
               CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1, 0 )
*
               IF( ierr( 1 ).LT.0 ) THEN
                  IF( iam.EQ.0 )
     $               WRITE( nout, fmt = 9997 )'descriptor'
                  kskip = kskip + 1
                  GO TO 10
               END IF
*
*              Assign pointers into MEM for SCALAPACK arrays, A is
*              allocated starting at position MEM( IPREPAD+1 )
*
               ndiag = nq
               IF( lsame( uplo, 'U' ) ) THEN
                  noffd = nq
               ELSE
                  noffd = numroc( n-1, nb, mycol, 0, npcol )
               END IF
*
               ipa = iprepad + 1
               ipd = ipa + desca( lld_ )*nq + ipostpad + iprepad
               ipe = ipd + ndiag + ipostpad + iprepad
               ipt = ipe + noffd + ipostpad + iprepad
               ipw = ipt + nq + ipostpad + iprepad
*
*              Calculate the amount of workspace required for the
*              reduction
*
               lwork = max( nb*( np+1 ), 3*nb )
               worktrd = lwork + ipostpad
               worksiz = worktrd
*
*              Figure the amount of workspace required by the check
*
               IF( check ) THEN
                  itemp = 2*nq + np
                  IF( nprow.NE.npcol ) THEN
                     lcm = ilcm( nprow, npcol )
                     itemp = nb*iceil( iceil( np, nb ), lcm / nprow ) +
     $                       itemp
                  END IF
                  itemp = max( itemp, 2*( nb+np )*nb )
                  worksiz = max( lwork, itemp ) + ipostpad
               END IF
*
*              Check for adequate memory for problem size
*
               ierr( 1 ) = 0
               IF( ipw+worksiz.GT.memsiz ) THEN
                  IF( iam.EQ.0 )
     $               WRITE( nout, fmt = 9996 )'Tridiagonal reduction',
     $               ( ipw+worksiz )*dblesz
                  ierr( 1 ) = 1
               END IF
*
*              Check all processes for an error
*
               CALL igsum2d( ictxt, 'All', ' ', 1, 1, ierr, 1, -1, 0 )
*
               IF( ierr( 1 ).GT.0 ) THEN
                  IF( iam.EQ.0 )
     $               WRITE( nout, fmt = 9997 )'MEMORY'
                  kskip = kskip + 1
                  GO TO 10
               END IF
*
*              Generate the matrix A
*
               CALL pdmatgen( ictxt, 'Symm', 'N', desca( m_ ),
     $                        desca( n_ ), desca( mb_ ), desca( nb_ ),
     $                        mem( ipa ), desca( lld_ ), desca( rsrc_ ),
     $                        desca( csrc_ ), iaseed, 0, np, 0, nq,
     $                        myrow, mycol, nprow, npcol )
*
*              Need Infinity-norm of A for checking
*
               IF( check ) THEN
                  CALL pdfillpad( ictxt, np, nq, mem( ipa-iprepad ),
     $                            desca( lld_ ), iprepad, ipostpad,
     $                            padval )
                  CALL pdfillpad( ictxt, ndiag, 1, mem( ipd-iprepad ),
     $                            ndiag, iprepad, ipostpad, padval )
                  CALL pdfillpad( ictxt, noffd, 1, mem( ipe-iprepad ),
     $                            noffd, iprepad, ipostpad, padval )
                  CALL pdfillpad( ictxt, nq, 1, mem( ipt-iprepad ), nq,
     $                            iprepad, ipostpad, padval )
                  CALL pdfillpad( ictxt, worksiz-ipostpad, 1,
     $                            mem( ipw-iprepad ), worksiz-ipostpad,
     $                            iprepad, ipostpad, padval )
                  anorm = pdlansy( 'I', uplo, n, mem( ipa ), 1, 1,
     $                    desca, mem( ipw ) )
                  CALL pdchekpad( ictxt, 'PDLANSY', np, nq,
     $                            mem( ipa-iprepad ), desca( lld_ ),
     $                            iprepad, ipostpad, padval )
                  CALL pdchekpad( ictxt, 'PDLANSY', worksiz-ipostpad, 1,
     $                            mem( ipw-iprepad ), worksiz-ipostpad,
     $                            iprepad, ipostpad, padval )
                  CALL pdfillpad( ictxt, worktrd-ipostpad, 1,
     $                            mem( ipw-iprepad ), worktrd-ipostpad,
     $                            iprepad, ipostpad, padval )
               END IF
*
               CALL slboot
               CALL blacs_barrier( ictxt, 'All' )
               CALL sltimer( 1 )
*
*              Reduce to symmetric tridiagonal form
*
               CALL pdsytrd( uplo, n, mem( ipa ), 1, 1, desca,
     $                       mem( ipd ), mem( ipe ), mem( ipt ),
     $                       mem( ipw ), lwork, info )
*
               CALL sltimer( 1 )
*
               IF( check ) THEN
*
*                 Check for memory overwrite
*
                  CALL pdchekpad( ictxt, 'PDSYTRD', np, nq,
     $                            mem( ipa-iprepad ), desca( lld_ ),
     $                            iprepad, ipostpad, padval )
                  CALL pdchekpad( ictxt, 'PDSYTRD', ndiag, 1,
     $                            mem( ipd-iprepad ), ndiag, iprepad,
     $                            ipostpad, padval )
                  CALL pdchekpad( ictxt, 'PDSYTRD', noffd, 1,
     $                            mem( ipe-iprepad ), noffd, iprepad,
     $                            ipostpad, padval )
                  CALL pdchekpad( ictxt, 'PDSYTRD', nq, 1,
     $                            mem( ipt-iprepad ), nq, iprepad,
     $                            ipostpad, padval )
                  CALL pdchekpad( ictxt, 'PDSYTRD', worktrd-ipostpad, 1,
     $                            mem( ipw-iprepad ), worktrd-ipostpad,
     $                            iprepad, ipostpad, padval )
                  CALL pdfillpad( ictxt, worksiz-ipostpad, 1,
     $                            mem( ipw-iprepad ), worksiz-ipostpad,
     $                            iprepad, ipostpad, padval )
*
*                 Compute fctres = ||A - QTQ'|| / (||A|| * N * eps)
*
                  CALL pdsytdrv( uplo, n, mem( ipa ), 1, 1, desca,
     $                           mem( ipd ), mem( ipe ), mem( ipt ),
     $                           mem( ipw ), ierr( 1 ) )
                  CALL pdlafchk( 'Symm', 'No', n, n, mem( ipa ), 1, 1,
     $                           desca, iaseed, anorm, fresid,
     $                           mem( ipw ) )
*
*                 Check for memory overwrite
*
                  CALL pdchekpad( ictxt, 'PDSYTDRV', np, nq,
     $                            mem( ipa-iprepad ), desca( lld_ ),
     $                            iprepad, ipostpad, padval )
                  CALL pdchekpad( ictxt, 'PDSYTDRV', ndiag, 1,
     $                            mem( ipd-iprepad ), ndiag, iprepad,
     $                            ipostpad, padval )
                  CALL pdchekpad( ictxt, 'PDSYTDRV', noffd, 1,
     $                            mem( ipe-iprepad ), noffd, iprepad,
     $                            ipostpad, padval )
                  CALL pdchekpad( ictxt, 'PDSYTDRV', worksiz-ipostpad,
     $                            1, mem( ipw-iprepad ),
     $                            worksiz-ipostpad, iprepad, ipostpad,
     $                            padval )
*
*                 Test residual and detect NaN result
*
                  IF( fresid.LE.thresh .AND. fresid-fresid.EQ.
     $                0.0d+0 .AND. ierr( 1 ).EQ.0 ) THEN
                     kpass = kpass + 1
                     passed = 'PASSED'
                  ELSE
                     IF( myrow.EQ.0 .AND. mycol.EQ.0 )
     $                  WRITE( nout, fmt = 9986 )fresid
                     kfail = kfail + 1
                     passed = 'FAILED'
                  END IF
*
                  IF( myrow.EQ.0 .AND. mycol.EQ.0 .AND. ierr( 1 ).NE.0 )
     $               WRITE( nout, fmt = * )'D or E copies incorrect ...'
               ELSE
*
*                 Don't perform the checking, only the timing operation
*
                  kpass = kpass + 1
                  fresid = fresid - fresid
                  passed = 'BYPASS'
               END IF
*
*              Gather maximum of all CPU and WALL clock timings
*
               CALL slcombine( ictxt, 'All', '>', 'W', 1, 1, wtime )
               CALL slcombine( ictxt, 'All', '>', 'C', 1, 1, ctime )
*
*              Print results
*
               IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
*
*                 TRD requires 4/3 N^3 floating point operations
*
                  nops = dble( n )
*
                  nops = ( 4.0d+0 / 3.0d+0 )*nops**3
                  nops = nops / 1.0d+6
*
*                 Print WALL time
*
                  IF( wtime( 1 ).GT.0.0d+0 ) THEN
                     tmflops = nops / wtime( 1 )
                  ELSE
                     tmflops = 0.0d+0
                  END IF
                  IF( wtime( 1 ).GE.0.0d+0 )
     $               WRITE( nout, fmt = 9993 )'WALL', uplo, n, nb,
     $               nprow, npcol, wtime( 1 ), tmflops, fresid, passed
*
*                 Print CPU time
*
                  IF( ctime( 1 ).GT.0.0d+0 ) THEN
                     tmflops = nops / ctime( 1 )
                  ELSE
                     tmflops = 0.0d+0
                  END IF
                  IF( ctime( 1 ).GE.0.0d+0 )
     $               WRITE( nout, fmt = 9993 )'CPU ', uplo, n, nb,
     $               nprow, npcol, ctime( 1 ), tmflops, fresid, passed
               END IF
   10       CONTINUE
   20    CONTINUE
*
         CALL blacs_gridexit( ictxt )
   30 CONTINUE
*
      CALL pdttrdtester( iam, nprocs, check, nout, thresh, nval, nmat,
     $                   mem, totmem, kpass, kfail, kskip )
*
*     Print ending messages and close output file
*
      IF( iam.EQ.0 ) THEN
         ktests = kpass + kfail + kskip
         WRITE( nout, fmt = * )
         WRITE( nout, fmt = 9992 )ktests
         IF( check ) THEN
            WRITE( nout, fmt = 9991 )kpass
            WRITE( nout, fmt = 9989 )kfail
         ELSE
            WRITE( nout, fmt = 9990 )kpass
         END IF
         WRITE( nout, fmt = 9988 )kskip
         WRITE( nout, fmt = * )
         WRITE( nout, fmt = * )
         WRITE( nout, fmt = 9987 )
         IF( nout.NE.6 .AND. nout.NE.0 )
     $      CLOSE ( nout )
      END IF
*
      CALL blacs_exit( 0 )
*
 9999 FORMAT( 'ILLEGAL ', a6, ': ', a5, ' = ', i3,
     $      '; It should be at least 1' )
 9998 FORMAT( 'ILLEGAL GRID: nprow*npcol = ', i4, '. It can be at most',
     $      i4 )
 9997 FORMAT( 'Bad ', a6, ' parameters: going on to next test case.' )
 9996 FORMAT( 'Unable to perform ', a, ': need TOTMEM of at least',
     $      i11 )
 9995 FORMAT( 'TIME UPLO      N  NB     P     Q  TRD Time ',
     $      '     MFLOPS Residual  CHECK' )
 9994 FORMAT( '---- ---- ------ --- ----- ----- --------- ',
     $      '----------- -------- ------' )
 9993 FORMAT( a4, 1x, a4, 1x, i6, 1x, i3, 1x, i5, 1x, i5, 1x, f9.2, 1x,
     $      f11.2, 1x, f8.2, 1x, a6 )
 9992 FORMAT( 'Finished', i4, ' tests, with the following results:' )
 9991 FORMAT( i5, ' tests completed and passed residual checks.' )
 9990 FORMAT( i5, ' tests completed without checking.' )
 9989 FORMAT( i5, ' tests completed and failed residual checks.' )
 9988 FORMAT( i5, ' tests skipped because of illegal input values.' )
 9987 FORMAT( 'END OF TESTS.' )
 9986 FORMAT( '||A - Q*T*Q''|| / (||A|| * N * eps) = ', g25.7 )
*
      stop
*
*     End of PDTRDDRIVER
*
      END