pcscaex.f

Go to the documentation of this file.

      PROGRAM pcscaex
*
*  -- 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 PCGESV. The input matrix and right-and-sides are
*     read from a file. The solution is written to a file.
*
*     .. Parameters ..
      INTEGER            cplxsz, intgsz, memsiz, totmem
      parameter( cplxsz = 8, intgsz = 4, totmem = 2000000,
     $                     memsiz = totmem / cplxsz )
      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 )
      COMPLEX            one
      parameter( one = (1.0d+0,0.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
      REAL               anorm, bnorm, eps, xnorm, resid
*     ..
*     .. Local Arrays ..
      INTEGER            desca( dlen_ ), descb( dlen_ ), descx( dlen_ )
      COMPLEX            mem( memsiz )
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_exit, blacs_get, blacs_gridexit,
     $                   blacs_gridinfo, blacs_gridinit, blacs_pinfo,
     $                   descinit, igsum2d, pcgesv,
     $                   pcgemm, pclacpy, pclaprnt, pclaread, pclawrite
*     ..
*     .. External Functions ..
      INTEGER            iceil, numroc
      REAL               pslamch, pclange
      EXTERNAL           iceil, numroc, pslamch, pclange
*     ..
*     .. 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 ), cplxsz )
      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 )*cplxsz
         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 pclaread( 'CSCAEXMAT.dat', mem( ipa ), desca, 0, 0,
     $               mem( ipw ) )
      CALL pclaread( 'CSCAEXRHS.dat', mem( ipb ), descb, 0, 0,
     $               mem( ipw ) )
*
*     Make a copy of A and the rhs for checking purposes
*
      CALL pclacpy( 'All', n, n, mem( ipa ), 1, 1, desca,
     $              mem( ipacpy ), 1, 1, desca )
      CALL pclacpy( 'All', n, nrhs, mem( ipb ), 1, 1, descb,
     $              mem( ipx ), 1, 1, descx )
*
**********************************************************************
*     Call ScaLAPACK PCGESV routine
**********************************************************************
*
      IF( iam.EQ.0 ) THEN
         WRITE( nout, fmt = * )
         WRITE( nout, fmt = * )
     $         '***********************************************'
         WRITE( nout, fmt = * )
     $         'example of scalapack routine call: (pcgesv)'
         WRITE( NOUT, FMT = * )
     $         '***********************************************'
         WRITE( NOUT, FMT = * )
         WRITE( NOUT, FMT = * ) 'a * x = b, matrix a:'
         WRITE( NOUT, FMT = * )
      END IF
      CALL PCLAPRNT( N, N, MEM( IPA ), 1, 1, DESCA, 0, 0,
     $               'a', NOUT, MEM( IPW ) )
.EQ.      IF( IAM0 ) THEN
         WRITE( NOUT, FMT = * )
         WRITE( NOUT, FMT = * ) 'matrix b:'
         WRITE( NOUT, FMT = * )
      END IF
      CALL PCLAPRNT( N, NRHS, MEM( IPB ), 1, 1, DESCB, 0, 0,
     $               'b', NOUT, MEM( IPW ) )
*
      CALL PCGESV( N, NRHS, MEM( IPA ), 1, 1, DESCA, MEM( IPPIV ),
     $             MEM( IPB ), 1, 1, DESCB, INFO )
*
.EQ..AND..EQ.      IF( MYROW0  MYCOL0 ) THEN
         WRITE( NOUT, FMT = * )
         WRITE( NOUT, FMT = * ) 'info code returned by pcgesv = ', INFO
         WRITE( NOUT, FMT = * )
         WRITE( NOUT, FMT = * ) 'matrix x = a^{-1} * b'
         WRITE( NOUT, FMT = * )
      END IF
      CALL PCLAPRNT( N, NRHS, MEM( IPB ), 1, 1, DESCB, 0, 0, 'x', NOUT,
     $               MEM( IPW ) )
      CALL PCLAWRITE( 'cscaexsol.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 = PCLANGE( 'i', N, N, MEM( IPA ), 1, 1, DESCA, MEM( IPW ) )
      BNORM = PCLANGE( 'i', N, NRHS, MEM( IPB ), 1, 1, DESCB,
     $                 MEM( IPW ) )
      CALL PCGEMM( '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 = PCLANGE( 'i', N, NRHS, MEM( IPX ), 1, 1, DESCX,
     $                 MEM( IPW ) )
      RESID = XNORM / ( ANORM * BNORM * EPS * DBLE( N ) )
*
.EQ..AND..EQ.      IF( MYROW0  MYCOL0 ) THEN
         WRITE( NOUT, FMT = * )
         WRITE( NOUT, FMT = * )
     $     '||a * x  - b|| / ( ||x|| * ||a|| * eps * n ) = ', RESID
         WRITE( NOUT, FMT = * )
.LT.         IF( RESID10.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
*
.EQ.      IF( IAM0 ) THEN
         WRITE( NOUT, FMT = * )
         WRITE( NOUT, FMT = * )
         WRITE( NOUT, FMT = 9997 )
         WRITE( NOUT, FMT = * )
.NE..AND..NE.         IF( NOUT6  NOUT0 )
     $      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 PCSCAEX
*
      END