psscaex.f

Go to the documentation of this file.

      PROGRAM psscaex
*
*  -- 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 PSGESV. The input matrix and right-and-sides are
*     read from a file. The solution is written to a file.
*
*     .. Parameters ..
      INTEGER            realsz, intgsz, memsiz, totmem
      parameter( realsz = 4, intgsz = 4, totmem = 2000000,
     $                     memsiz = totmem / realsz )
      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 )
      REAL               one
      parameter( one = 1.0e+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
      REAL               anorm, bnorm, eps, xnorm, resid
*     ..
*     .. Local Arrays ..
      INTEGER            desca( dlen_ ), descb( dlen_ ), descx( dlen_ )
      REAL               mem( memsiz )
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_exit, blacs_get, blacs_gridexit,
     $                   blacs_gridinfo, blacs_gridinit, blacs_pinfo,
     $                   descinit, igsum2d, pdscaexinfo, psgesv,
     $                   psgemm, pslacpy, pslaprnt, pslaread, pslawrite
*     ..
*     .. External Functions ..
      INTEGER            iceil, numroc
      REAL               pslamch, pslange
      EXTERNAL           iceil, numroc, pslamch, pslange
*     ..
*     .. 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 ), realsz )
      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 )*realsz
         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 pslaread( 'SSCAEXMAT.dat', mem( ipa ), desca, 0, 0,
     $               mem( ipw ) )
      CALL pslaread( 'SSCAEXRHS.dat', mem( ipb ), descb, 0, 0,
     $               mem( ipw ) )
*
*     Make a copy of A and the rhs for checking purposes
*
      CALL pslacpy( 'All', n, n, mem( ipa ), 1, 1, desca,
     $              mem( ipacpy ), 1, 1, desca )
      CALL pslacpy( 'All', n, nrhs, mem( ipb ), 1, 1, descb,
     $              mem( ipx ), 1, 1, descx )
*
**********************************************************************
*     Call ScaLAPACK PSGESV routine
**********************************************************************
*
      IF( iam.EQ.0 ) THEN
         WRITE( nout, fmt = * )
         WRITE( nout, fmt = * )
     $         '***********************************************'
         WRITE( nout, fmt = * )
     $         'Example of ScaLAPACK routine call: (PSGESV)'
         WRITE( nout, fmt = * )
     $         '***********************************************'
         WRITE( nout, fmt = * )
         WRITE( nout, fmt = * ) 'A * X = B, Matrix A:'
         WRITE( nout, fmt = * )
      END IF
      CALL pslaprnt( 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 pslaprnt( n, nrhs, mem( ipb ), 1, 1, descb, 0, 0,
     $               'B', nout, mem( ipw ) )
*
      CALL psgesv( 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 PSGESV = ', info
         WRITE( nout, fmt = * )
         WRITE( nout, fmt = * ) 'Matrix X = A^{-1} * B'
         WRITE( nout, fmt = * )
      END IF
      CALL pslaprnt( n, nrhs, mem( ipb ), 1, 1, descb, 0, 0, 'X', nout,
     $               mem( ipw ) )
      CALL pslawrite( 'SSCAEXSOL.dat', n, nrhs, mem( ipb ), 1, 1, descb,
     $                0, 0, mem( ipw ) )
*
*     Compute residual ||A * X  - B|| / ( ||X|| * ||A|| * eps * N )
*
      eps = pslamch( ictxt, 'Epsilon' )
      anorm = pslange( 'I', n, n, mem( ipa ), 1, 1, desca, mem( ipw ) )
      bnorm = pslange( 'I', n, nrhs, mem( ipb ), 1, 1, descb,
     $                 mem( ipw ) )
      CALL psgemm( '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 = pslange( '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 PSSCAEX
*
      END