pdscaex.f

Go to the documentation of this file.

      PROGRAM pdscaex
*
*  -- ScaLAPACK example code --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*
*     Written by Antoine Petitet, August 1995 (petitet@cs.utk.edu)
*
*     This program solves a linear system by calling the ScaLAPACK 
*     routine PDGESV. The input matrix and right-and-sides are
*     read from a file. The solution is written to a file.
*
*     .. Parameters ..
      INTEGER            dblesz, intgsz, memsiz, totmem
      parameter( dblesz = 8, intgsz = 4, totmem = 2000000,
     $                     memsiz = totmem / dblesz )
      INTEGER            block_cyclic_2d, csrc_, ctxt_, dlen_, dt_,
     $                   lld_, mb_, m_, nb_, n_, rsrc_
      parameter( block_cyclic_2d = 1, dlen_ = 9, dt_ = 1,
     $                     ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
     $                     rsrc_ = 7, csrc_ = 8, lld_ = 9 )
      DOUBLE PRECISION   one
      parameter( one = 1.0d+0 )
*     ..
*     .. Local Scalars ..
      CHARACTER*80       outfile
      INTEGER            iam, ictxt, info, ipa, ipacpy, ipb, ippiv, ipx,
     $                   ipw, lipiv, mycol, myrow, n, nb, nout, npcol,
     $                   nprocs, nprow, np, nq, nqrhs, nrhs, worksiz
      DOUBLE PRECISION   anorm, bnorm, eps, xnorm, resid
*     ..
*     .. Local Arrays ..
      INTEGER            desca( dlen_ ), descb( dlen_ ), descx( dlen_ )
      DOUBLE PRECISION   mem( memsiz )
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_exit, blacs_get, blacs_gridexit,
     $                   blacs_gridinfo, blacs_gridinit, blacs_pinfo,
     $                   descinit, igsum2d, pdscaexinfo, pdgesv,
     $                   pdgemm, pdlacpy, pdlaprnt, pdlaread, pdlawrite
*     ..
*     .. External Functions ..
      INTEGER            iceil, numroc
      DOUBLE PRECISION   pdlamch, pdlange
      EXTERNAL           iceil, numroc, pdlamch, pdlange
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          dble, max
*     ..
*     .. Executable Statements ..
*
*     Get starting information
*
      CALL blacs_pinfo( iam, nprocs )
      CALL pdscaexinfo( outfile, nout, n, nrhs, nb, nprow, npcol, mem,
     $                  iam, nprocs )
*
*     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 process grid loop if this case doesn't use my
*     process
*
      IF( myrow.GE.nprow .OR. mycol.GE.npcol )
     $   GO TO 20
*
      np    = numroc( n, nb, myrow, 0, nprow )
      nq    = numroc( n, nb, mycol, 0, npcol )
      nqrhs = numroc( nrhs, nb, mycol, 0, npcol )
*
*     Initialize the array descriptor for the matrix A and B
*
      CALL descinit( desca, n, n, nb, nb, 0, 0, ictxt, max( 1, np ),
     $               info )
      CALL descinit( descb, n, nrhs, nb, nb, 0, 0, ictxt, max( 1, np ),
     $               info )
      CALL descinit( descx, n, nrhs, nb, nb, 0, 0, ictxt, max( 1, np ),
     $               info )
*
*     Assign pointers into MEM for SCALAPACK arrays, A is
*     allocated starting at position MEM( 1 )
*
      ipa = 1
      ipacpy = ipa + desca( lld_ )*nq
      ipb = ipacpy + desca( lld_ )*nq
      ipx = ipb + descb( lld_ )*nqrhs
      ippiv = ipx + descb( lld_ )*nqrhs
      lipiv = iceil( intgsz*( np+nb ), dblesz )
      ipw = ippiv + max( np, lipiv )
*
      worksiz = nb
*
*     Check for adequate memory for problem size
*
      info = 0
      IF( ipw+worksiz.GT.memsiz ) THEN
         IF( iam.EQ.0 )
     $      WRITE( nout, fmt = 9998 ) 'test', ( ipw+worksiz )*dblesz
         info = 1
      END IF
*
*     Check all processes for an error
*
      CALL igsum2d( ictxt, 'All', ' ', 1, 1, info, 1, -1, 0 )
      IF( info.GT.0 ) THEN
         IF( iam.EQ.0 )
     $      WRITE( nout, fmt = 9999 ) 'MEMORY'
         GO TO 10
      END IF
*
*     Read from file and distribute matrices A and B
*
      CALL pdlaread( 'DSCAEXMAT.dat', mem( ipa ), desca, 0, 0,
     $               mem( ipw ) )
      CALL pdlaread( 'DSCAEXRHS.dat', mem( ipb ), descb, 0, 0,
     $               mem( ipw ) )
*
*     Make a copy of A and the rhs for checking purposes
*
      CALL pdlacpy( 'All', n, n, mem( ipa ), 1, 1, desca,
     $              mem( ipacpy ), 1, 1, desca )
      CALL pdlacpy( 'All', n, nrhs, mem( ipb ), 1, 1, descb,
     $              mem( ipx ), 1, 1, descx )
*
**********************************************************************
*     Call ScaLAPACK PDGESV routine
**********************************************************************
*
      IF( iam.EQ.0 ) THEN
         WRITE( nout, fmt = * )
         WRITE( nout, fmt = * )
     $         '***********************************************'
         WRITE( nout, fmt = * )
     $         'Example of ScaLAPACK routine call: (PDGESV)'
         WRITE( nout, fmt = * )
     $         '***********************************************'
         WRITE( nout, fmt = * )
         WRITE( nout, fmt = * ) 'A * X = B, Matrix A:'
         WRITE( nout, fmt = * )
      END IF
      CALL pdlaprnt( n, n, mem( ipa ), 1, 1, desca, 0, 0,
     $               'A', nout, mem( ipw ) )
      IF( iam.EQ.0 ) THEN
         WRITE( nout, fmt = * )
         WRITE( nout, fmt = * ) 'Matrix B:'
         WRITE( nout, fmt = * )
      END IF
      CALL pdlaprnt( n, nrhs, mem( ipb ), 1, 1, descb, 0, 0,
     $               'B', nout, mem( ipw ) )
*
      CALL pdgesv( n, nrhs, mem( ipa ), 1, 1, desca, mem( ippiv ),
     $             mem( ipb ), 1, 1, descb, info )
*
      IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
         WRITE( nout, fmt = * )
         WRITE( nout, fmt = * ) 'INFO code returned by PDGESV = ', info
         WRITE( nout, fmt = * )
         WRITE( nout, fmt = * ) 'Matrix X = A^{-1} * B'
         WRITE( nout, fmt = * )
      END IF
      CALL pdlaprnt( n, nrhs, mem( ipb ), 1, 1, descb, 0, 0, 'X', nout,
     $               mem( ipw ) )
      CALL pdlawrite( 'DSCAEXSOL.dat', n, nrhs, mem( ipb ), 1, 1, descb,
     $                0, 0, mem( ipw ) )
*
*     Compute residual ||A * X  - B|| / ( ||X|| * ||A|| * eps * N )
*
      eps = pdlamch( ictxt, 'Epsilon' )
      anorm = pdlange( 'I', n, n, mem( ipa ), 1, 1, desca, mem( ipw ) )
      bnorm = pdlange( 'I', n, nrhs, mem( ipb ), 1, 1, descb,
     $                 mem( ipw ) )
      CALL pdgemm( 'No transpose', 'No transpose', n, nrhs, n, one,
     $             mem( ipacpy ), 1, 1, desca, mem( ipb ), 1, 1, descb,
     $             -one, mem( ipx ), 1, 1, descx )
      xnorm = pdlange( 'I', n, nrhs, mem( ipx ), 1, 1, descx,
     $                 mem( ipw ) )
      resid = xnorm / ( anorm * bnorm * eps * dble( n ) )
*
      IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
         WRITE( nout, fmt = * )
         WRITE( nout, fmt = * )
     $     '||A * X  - B|| / ( ||X|| * ||A|| * eps * N ) = ', resid
         WRITE( nout, fmt = * )
         IF( resid.LT.10.0d+0 ) THEN
            WRITE( nout, fmt = * ) 'The answer is correct.'
         ELSE
            WRITE( nout, fmt = * ) 'The answer is suspicious.'
         END IF
      END IF
*
   10 CONTINUE
*
      CALL blacs_gridexit( ictxt )
*
   20 CONTINUE
*
*     Print ending messages and close output file
*
      IF( iam.EQ.0 ) THEN
         WRITE( nout, fmt = * )
         WRITE( nout, fmt = * )
         WRITE( nout, fmt = 9997 )
         WRITE( nout, fmt = * )
         IF( nout.NE.6 .AND. nout.NE.0 )
     $      CLOSE ( nout )
      END IF
*
      CALL blacs_exit( 0 )
*
 9999 FORMAT( 'Bad ', a6, ' parameters: going on to next test case.' )
 9998 FORMAT( 'Unable to perform ', a, ': need TOTMEM of at least',
     $        i11 )
 9997 FORMAT( 'END OF TESTS.' )
*
      stop
*
*     End of PDSCAEX
*
      END