Next: HPF Interface to ScaLAPACK Up: Example Programs Previous: Example Programs

# Example Program #2

In Chapter 2, we presented an example program using ScaLAPACK. Here we present a second example--a more flexible and memory efficient program to solve a system of linear equations using the ScaLAPACK driver routine PDGESV. In this example we will read the input matrices from a file, distribute these matrices to the processes in the grid. After calling the ScaLAPACK routine, we will write the output solution matrix to a file. The input data files for the program are SCAEX.dat, SCAEXMAT.dat, and SCAEXRHS.dat.

This program is also available in the scalapack directory on netlib
(http://www.netlib.org/scalapack/examples/scaex.shar).

SCAEX.dat is:

    'ScaLAPACK Example Program 2'
    'May 1997'
    'SCAEX.out'             output file name (if any)
    6                       device out
    6                       value of N
    1                       value of NRHS
    2                       values of NB
    2                       values of NPROW
    2                       values of NPCOL

SCAEXMAT.dat is:

    6 6
     6.0000D+0
     3.0000D+0
     0.0000D+0
     0.0000D+0
     3.0000D+0
     0.0000D+0
     0.0000D+0
    -3.0000D+0
    -1.0000D+0
     1.0000D+0
     1.0000D+0
     0.0000D+0
    -1.0000D+0
     0.0000D+0
    11.0000D+0
     0.0000D+0
     0.0000D+0
    10.0000D+0
     0.0000D+0
     0.0000D+0
     0.0000D+0
    -11.0000D+0
     0.0000D+0
     0.0000D+0
     0.0000D+0
     0.0000D+0
     0.0000D+0
     2.0000D+0
    -4.0000D+0
     0.0000D+0
     0.0000D+0
     0.0000D+0
     0.0000D+0
     8.0000D+0
     0.0000D+0
    -10.0000D+0

SCAEXRHS.dat is:

       6  1
         72.000000000000000000D+00
          0.000000000000000000D+00
        160.000000000000000000D+00
          0.000000000000000000D+00
          0.000000000000000000D+00
          0.000000000000000000D+00

          PROGRAM PDSCAEX
    *
    *  -- ScaLAPACK example code --
    *     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
    *     and University of California, Berkeley.
    *
    *     Written by Antoine Petitet, (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_, 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 )
          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 = MAX( NB, NP )
    *
    *     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( 'SCAEXMAT.dat', MEM( IPA ), DESCA, 0, 0,
         $               MEM( IPW ) )
          CALL PDLAREAD( 'SCAEXRHS.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( 'SCAEXSOL.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

          SUBROUTINE PDSCAEXINFO( SUMMRY, NOUT, N, NRHS, NB, NPROW, NPCOL,
         $                        WORK, IAM, NPROCS )
    *
    *  -- ScaLAPACK example code --
    *     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
    *     and University of California, Berkeley.
    *
    *     Written by Antoine Petitet, (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.
    *
    *     .. Scalar Arguments ..
          CHARACTER*( * )    SUMMRY
          INTEGER            IAM, N, NRHS, NB, NOUT, NPCOL, NPROCS, NPROW
    *     ..
    *     .. Array Arguments ..
          INTEGER            WORK( * )
    *     ..
    *
    * ======================================================================
    *
    *     .. Parameters ..
          INTEGER            NIN
          PARAMETER          ( NIN = 11 )
    *     ..
    *     .. Local Scalars ..
          CHARACTER*79       USRINFO
          INTEGER            ICTXT
    *     ..
    *     .. External Subroutines ..
          EXTERNAL           BLACS_ABORT, BLACS_GET, BLACS_GRIDEXIT,
         $                   BLACS_GRIDINIT, BLACS_SETUP, IGEBR2D, IGEBS2D
    *     ..
    *     .. Intrinsic Functions ..
          INTRINSIC          MAX, MIN
    *     ..
    *     .. Executable Statements ..
    *
    *     Process 0 reads the input data, broadcasts to other processes and
    *     writes needed information to NOUT
    *
          IF( IAM.EQ.0 ) THEN
    *
    *        Open file and skip data file header
    *
             OPEN( NIN, FILE='SCAEX.dat', STATUS='OLD' )
             READ( NIN, FMT = * ) SUMMRY
             SUMMRY = ' '
    *
    *        Read in user-supplied info about machine type, compiler, etc.
    *
             READ( NIN, FMT = 9999 ) USRINFO
    *
    *        Read name and unit number for summary output file
    *
             READ( NIN, FMT = * ) SUMMRY
             READ( NIN, FMT = * ) NOUT
             IF( NOUT.NE.0 .AND. NOUT.NE.6 )
         $      OPEN( NOUT, FILE = SUMMRY, STATUS = 'UNKNOWN' )
    *
    *        Read and check the parameter values for the tests.
    *
    *        Get matrix dimensions
    *
             READ( NIN, FMT = * ) N
             READ( NIN, FMT = * ) NRHS
    *
    *        Get value of NB
    *
             READ( NIN, FMT = * ) NB
    *
    *        Get grid shape
    *
             READ( NIN, FMT = * ) NPROW
             READ( NIN, FMT = * ) NPCOL
    *
    *        Close input file
    *
             CLOSE( NIN )
    *
    *        If underlying system needs additional set up, do it now
    *
             IF( NPROCS.LT.1 ) THEN
                NPROCS = NPROW * NPCOL
                CALL BLACS_SETUP( IAM, NPROCS )
             END IF
    *
    *        Temporarily define blacs grid to include all processes so
    *        information can be broadcast to all processes
    *
             CALL BLACS_GET( -1, 0, ICTXT )
             CALL BLACS_GRIDINIT( ICTXT, 'Row-major', 1, NPROCS )
    *
    *        Pack information arrays and broadcast
    *
             WORK( 1 ) = N
             WORK( 2 ) = NRHS
             WORK( 3 ) = NB
             WORK( 4 ) = NPROW
             WORK( 5 ) = NPCOL
             CALL IGEBS2D( ICTXT, 'All', ' ', 5, 1, WORK, 5 )
    *
    *        regurgitate input
    *
             WRITE( NOUT, FMT = 9999 )
         $               'SCALAPACK example driver.'
             WRITE( NOUT, FMT = 9999 ) USRINFO
             WRITE( NOUT, FMT = * )
             WRITE( NOUT, FMT = 9999 )
         $               'The matrices A and B are read from '//
         $               'a file.'
             WRITE( NOUT, FMT = * )
             WRITE( NOUT, FMT = 9999 )
         $               'An explanation of the input/output '//
         $               'parameters follows:'
    *
             WRITE( NOUT, FMT = 9999 )
         $               'N       : The order of the matrix A.'
             WRITE( NOUT, FMT = 9999 )
         $               'NRHS    : The number of right and sides.'
             WRITE( NOUT, FMT = 9999 )
         $               'NB      : The size of the square blocks the'//
         $               ' matrices A and B are split into.'
             WRITE( NOUT, FMT = 9999 )
         $               'P       : The number of process rows.'
             WRITE( NOUT, FMT = 9999 )
         $               'Q       : The number of process columns.'
             WRITE( NOUT, FMT = * )
             WRITE( NOUT, FMT = 9999 )
         $               'The following parameter values will be used:'
             WRITE( NOUT, FMT = 9998 ) 'N    ', N
             WRITE( NOUT, FMT = 9998 ) 'NRHS ', NRHS
             WRITE( NOUT, FMT = 9998 ) 'NB   ', NB
             WRITE( NOUT, FMT = 9998 ) 'P    ', NPROW
             WRITE( NOUT, FMT = 9998 ) 'Q    ', NPCOL
             WRITE( NOUT, FMT = * )
    *
          ELSE
    *
    *        If underlying system needs additional set up, do it now
    *
             IF( NPROCS.LT.1 )
         $      CALL BLACS_SETUP( IAM, NPROCS )
    *
    *        Temporarily define blacs grid to include all processes so
    *        information can be broadcast to all processes
    *
             CALL BLACS_GET( -1, 0, ICTXT )
             CALL BLACS_GRIDINIT( ICTXT, 'Row-major', 1, NPROCS )
    *
             CALL IGEBR2D( ICTXT, 'All', ' ', 5, 1, WORK, 5, 0, 0 )
             N     = WORK( 1 )
             NRHS  = WORK( 2 )
             NB    = WORK( 3 )
             NPROW = WORK( 4 )
             NPCOL = WORK( 5 )
    *
          END IF
    *
          CALL BLACS_GRIDEXIT( ICTXT )
    *
          RETURN
    *
       20 WRITE( NOUT, FMT = 9997 )
          CLOSE( NIN )
          IF( NOUT.NE.6 .AND. NOUT.NE.0 )
         $   CLOSE( NOUT )
          CALL BLACS_ABORT( ICTXT, 1 )
    *
          STOP
    *
     9999 FORMAT( A )
     9998 FORMAT( 2X, A5, '   :        ', I6 )
     9997 FORMAT( ' Illegal input in file ',40A,'.  Aborting run.' )
    *
    *     End of PDSCAEXINFO
    *
          END

          SUBROUTINE PDLAREAD( FILNAM, A, DESCA, IRREAD, ICREAD, WORK )
    *
    *  -- ScaLAPACK example --
    *     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
    *     and University of California, Berkeley.
    *
    *     written by Antoine Petitet, (petitet@cs.utk.edu)
    *
    *     .. Scalar Arguments ..
          INTEGER            ICREAD, IRREAD
    *     ..
    *     .. Array Arguments ..
          CHARACTER*(*)      FILNAM
          INTEGER            DESCA( * )
          DOUBLE PRECISION   A( * ), WORK( * )
    *     ..
    *
    *  Purpose
    *  =======
    *
    *  PDLAREAD reads from a file named FILNAM a matrix and distribute
    *  it to the process grid.
    *
    *  Only the process of coordinates {IRREAD, ICREAD} read the file.
    *
    *  WORK must be of size >= MB_ = DESCA( MB_ ).
    *
    *  =====================================================================
    *
    *     .. Parameters ..
          INTEGER            NIN
          PARAMETER          ( NIN = 11 )
          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 )
    *     ..
    *     .. Local Scalars ..
          INTEGER            H, I, IB, ICTXT, ICURCOL, ICURROW, II, J, JB,
         $                   JJ, K, LDA, M, MYCOL, MYROW, N, NPCOL, NPROW
    *     ..
    *     .. Local Arrays ..
          INTEGER            IWORK( 2 )
    *     ..
    *     .. External Subroutines ..
          EXTERNAL           BLACS_GRIDINFO, INFOG2L, DGERV2D, DGESD2D,
         $                   IGEBS2D, IGEBR2D
    *     ..
    *     .. External Functions ..
          INTEGER            ICEIL
          EXTERNAL           ICEIL
    *     ..
    *     .. Intrinsic Functions ..
          INTRINSIC          MIN
    *     ..
    *     .. Executable Statements ..
    *
    *     Get grid parameters
    *
          ICTXT = DESCA( CTXT_ )
          CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
    *
          IF( MYROW.EQ.IRREAD .AND. MYCOL.EQ.ICREAD ) THEN
             OPEN( NIN, FILE=FILNAM, STATUS='OLD' )
             READ( NIN, FMT = * ) ( IWORK( I ), I = 1, 2 )
             CALL IGEBS2D( ICTXT, 'All', ' ', 2, 1, IWORK, 2 )
          ELSE
             CALL IGEBR2D( ICTXT, 'All', ' ', 2, 1, IWORK, 2, IRREAD,
         $                 ICREAD )
          END IF
          M = IWORK( 1 )
          N = IWORK( 2 )
    *
          IF( M.LE.0 .OR. N.LE.0 )
         $   RETURN
    *
          IF( M.GT.DESCA( M_ ).OR. N.GT.DESCA( N_ ) ) THEN
             IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
                WRITE( *, FMT = * ) 'PDLAREAD: Matrix too big to fit in'
                WRITE( *, FMT = * ) 'Abort ...'
             END IF
             CALL BLACS_ABORT( ICTXT, 0 )
          END IF
    *
          II = 1
          JJ = 1
          ICURROW = DESCA( RSRC_ )
          ICURCOL = DESCA( CSRC_ )
          LDA = DESCA( LLD_ )
    *
    *     Loop over column blocks
    *
          DO 50 J = 1, N, DESCA( NB_ )
             JB = MIN(  DESCA( NB_ ), N-J+1 )
             DO 40 H = 0, JB-1
    *
    *           Loop over block of rows
    *
                DO 30 I = 1, M, DESCA( MB_ )
                   IB = MIN( DESCA( MB_ ), M-I+1 )
                   IF( ICURROW.EQ.IRREAD .AND. ICURCOL.EQ.ICREAD ) THEN
                      IF( MYROW.EQ.IRREAD .AND. MYCOL.EQ.ICREAD ) THEN
                         DO 10 K = 0, IB-1
                            READ( NIN, FMT = * ) A( II+K+(JJ+H-1)*LDA )
       10                CONTINUE
                      END IF
                   ELSE
                      IF( MYROW.EQ.ICURROW .AND. MYCOL.EQ.ICURCOL ) THEN
                         CALL DGERV2D( ICTXT, IB, 1, A( II+(JJ+H-1)*LDA ),
         $                             LDA, IRREAD, ICREAD )
                       ELSE IF( MYROW.EQ.IRREAD .AND. MYCOL.EQ.ICREAD ) THEN
                         DO 20 K = 1, IB
                            READ( NIN, FMT = * ) WORK( K )
       20                CONTINUE
                         CALL DGESD2D( ICTXT, IB, 1, WORK, DESCA( MB_ ),
         $                             ICURROW, ICURCOL )
                      END IF
                   END IF
                   IF( MYROW.EQ.ICURROW )
         $            II = II + IB
                   ICURROW = MOD( ICURROW+1, NPROW )
       30       CONTINUE
    *
                II = 1
                ICURROW = DESCA( RSRC_ )
       40    CONTINUE
    *
             IF( MYCOL.EQ.ICURCOL )
         $      JJ = JJ + JB
             ICURCOL = MOD( ICURCOL+1, NPCOL )
    *
       50 CONTINUE
    *
          IF( MYROW.EQ.IRREAD .AND. MYCOL.EQ.ICREAD ) THEN
             CLOSE( NIN )
          END IF
    *
          RETURN
    *
    *     End of PDLAREAD
    *
          END

          SUBROUTINE PDLAWRITE( FILNAM, M, N, A, IA, JA, DESCA, IRWRIT,
         $                      ICWRIT, WORK )
    *
    *  -- ScaLAPACK example --
    *     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
    *     and University of California, Berkeley.
    *
    *     written by Antoine Petitet, (petitet@cs.utk.edu)
    *
    *     .. Scalar Arguments ..
          INTEGER            IA, ICWRIT, IRWRIT, JA, M, N
    *     ..
    *     .. Array Arguments ..
          CHARACTER*(*)      FILNAM
          INTEGER            DESCA( * )
          DOUBLE PRECISION   A( * ), WORK( * )
    *     ..
    *
    *  Purpose
    *  =======
    *
    *  PDLAWRITE writes to a file named FILNAMa distributed matrix sub( A )
    *  denoting A(IA:IA+M-1,JA:JA+N-1). The local pieces are sent to and
    *  written by the process of coordinates (IRWWRITE, ICWRIT).
    *
    *  WORK must be of size >= MB_ = DESCA( MB_ ).
    *
    *  =====================================================================
    *
    *     .. Parameters ..
          INTEGER            NOUT
          PARAMETER          ( NOUT = 13 )
          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 )
    *     ..
    *     .. Local Scalars ..
          INTEGER            H, I, IACOL, IAROW, IB, ICTXT, ICURCOL,
         $                   ICURROW, II, IIA, IN, J, JB, JJ, JJA, JN, K,
         $                   LDA, MYCOL, MYROW, NPCOL, NPROW
    *     ..
    *     .. External Subroutines ..
          EXTERNAL           BLACS_BARRIER, BLACS_GRIDINFO, INFOG2L,
         $                   DGERV2D, DGESD2D
    *     ..
    *     .. External Functions ..
          INTEGER            ICEIL
          EXTERNAL           ICEIL
    *     ..
    *     .. Intrinsic Functions ..
          INTRINSIC          MIN
    *     ..
    *     .. Executable Statements ..
    *
    *     Get grid parameters
    *
          ICTXT = DESCA( CTXT_ )
          CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
    *
          IF( MYROW.EQ.IRWRIT .AND. MYCOL.EQ.ICWRIT ) THEN
             OPEN( NOUT, FILE=FILNAM, STATUS='UNKNOWN' )
             WRITE( NOUT, FMT = * ) M, N
          END IF
    *
          CALL INFOG2L( IA, JA, DESCA, NPROW, NPCOL, MYROW, MYCOL,
         $              IIA, JJA, IAROW, IACOL )
          ICURROW = IAROW
          ICURCOL = IACOL
          II = IIA
          JJ = JJA
          LDA = DESCA( LLD_ )
    *
    *     Handle the first block of column separately
    *
          JN = MIN( ICEIL( JA, DESCA( NB_ ) ) * DESCA( NB_ ), JA+N-1 )
          JB = JN-JA+1
          DO 60 H = 0, JB-1
             IN = MIN( ICEIL( IA, DESCA( MB_ ) ) * DESCA( MB_ ), IA+M-1 )
             IB = IN-IA+1
             IF( ICURROW.EQ.IRWRIT .AND. ICURCOL.EQ.ICWRIT ) THEN
                IF( MYROW.EQ.IRWRIT .AND. MYCOL.EQ.ICWRIT ) THEN
                   DO 10 K = 0, IB-1
                      WRITE( NOUT, FMT = 9999 ) A( II+K+(JJ+H-1)*LDA )
       10          CONTINUE
                END IF
             ELSE
                IF( MYROW.EQ.ICURROW .AND. MYCOL.EQ.ICURCOL ) THEN
                   CALL DGESD2D( ICTXT, IB, 1, A( II+(JJ+H-1)*LDA ), LDA,
         $                       IRWRIT, ICWRIT )
                ELSE IF( MYROW.EQ.IRWRIT .AND. MYCOL.EQ.ICWRIT ) THEN
                   CALL DGERV2D( ICTXT, IB, 1, WORK, DESCA( MB_ ),
         $                       ICURROW, ICURCOL )
                   DO 20 K = 1, IB
                      WRITE( NOUT, FMT = 9999 ) WORK( K )
       20          CONTINUE
                END IF
             END IF
             IF( MYROW.EQ.ICURROW )
         $      II = II + IB
             ICURROW = MOD( ICURROW+1, NPROW )
             CALL BLACS_BARRIER( ICTXT, 'All' )
    *
    *        Loop over remaining block of rows
    *
             DO 50 I = IN+1, IA+M-1, DESCA( MB_ )
                IB = MIN( DESCA( MB_ ), IA+M-I )
                IF( ICURROW.EQ.IRWRIT .AND. ICURCOL.EQ.ICWRIT ) THEN
                   IF( MYROW.EQ.IRWRIT .AND. MYCOL.EQ.ICWRIT ) THEN
                      DO 30 K = 0, IB-1
                         WRITE( NOUT, FMT = 9999 ) A( II+K+(JJ+H-1)*LDA )
       30             CONTINUE
                   END IF
                ELSE
                   IF( MYROW.EQ.ICURROW .AND. MYCOL.EQ.ICURCOL ) THEN
                      CALL DGESD2D( ICTXT, IB, 1, A( II+(JJ+H-1)*LDA ),
         $                          LDA, IRWRIT, ICWRIT )
                   ELSE IF( MYROW.EQ.IRWRIT .AND. MYCOL.EQ.ICWRIT ) THEN
                      CALL DGERV2D( ICTXT, IB, 1, WORK, DESCA( MB_ ),
         $                          ICURROW, ICURCOL )
                      DO 40 K = 1, IB
                         WRITE( NOUT, FMT = 9999 ) WORK( K )
       40             CONTINUE
                   END IF
                END IF
                IF( MYROW.EQ.ICURROW )
         $         II = II + IB
                ICURROW = MOD( ICURROW+1, NPROW )
                CALL BLACS_BARRIER( ICTXT, 'All' )
       50    CONTINUE
    *
            II = IIA
            ICURROW = IAROW
       60 CONTINUE
    *
          IF( MYCOL.EQ.ICURCOL )
         $   JJ = JJ + JB
          ICURCOL = MOD( ICURCOL+1, NPCOL )
          CALL BLACS_BARRIER( ICTXT, 'All' )
    *
    *     Loop over remaining column blocks
    *
          DO 130 J = JN+1, JA+N-1, DESCA( NB_ )
             JB = MIN(  DESCA( NB_ ), JA+N-J )
             DO 120 H = 0, JB-1
                IN = MIN( ICEIL( IA, DESCA( MB_ ) ) * DESCA( MB_ ), IA+M-1 )
                IB = IN-IA+1
                IF( ICURROW.EQ.IRWRIT .AND. ICURCOL.EQ.ICWRIT ) THEN
                   IF( MYROW.EQ.IRWRIT .AND. MYCOL.EQ.ICWRIT ) THEN
                      DO 70 K = 0, IB-1
                         WRITE( NOUT, FMT = 9999 ) A( II+K+(JJ+H-1)*LDA )
       70             CONTINUE
                   END IF
                ELSE
                   IF( MYROW.EQ.ICURROW .AND. MYCOL.EQ.ICURCOL ) THEN
                      CALL DGESD2D( ICTXT, IB, 1, A( II+(JJ+H-1)*LDA ),
         $                          LDA, IRWRIT, ICWRIT )
                   ELSE IF( MYROW.EQ.IRWRIT .AND. MYCOL.EQ.ICWRIT ) THEN
                      CALL DGERV2D( ICTXT, IB, 1, WORK, DESCA( MB_ ),
         $                          ICURROW, ICURCOL )
                      DO 80 K = 1, IB
                         WRITE( NOUT, FMT = 9999 ) WORK( K )
       80             CONTINUE
                   END IF
                END IF
                IF( MYROW.EQ.ICURROW )
         $         II = II + IB
                ICURROW = MOD( ICURROW+1, NPROW )
                CALL BLACS_BARRIER( ICTXT, 'All' )
    *
    *           Loop over remaining block of rows
    *
                DO 110 I = IN+1, IA+M-1, DESCA( MB_ )
                   IB = MIN( DESCA( MB_ ), IA+M-I )
                   IF( ICURROW.EQ.IRWRIT .AND. ICURCOL.EQ.ICWRIT ) THEN
                      IF( MYROW.EQ.IRWRIT .AND. MYCOL.EQ.ICWRIT ) THEN
                         DO 90 K = 0, IB-1
                            WRITE( NOUT, FMT = 9999 ) A( II+K+(JJ+H-1)*LDA )
       90                CONTINUE
                      END IF
                   ELSE
                      IF( MYROW.EQ.ICURROW .AND. MYCOL.EQ.ICURCOL ) THEN
                         CALL DGESD2D( ICTXT, IB, 1, A( II+(JJ+H-1)*LDA ),
         $                             LDA, IRWRIT, ICWRIT )
                       ELSE IF( MYROW.EQ.IRWRIT .AND. MYCOL.EQ.ICWRIT ) THEN
                         CALL DGERV2D( ICTXT, IB, 1, WORK, DESCA( MB_ ),
         $                             ICURROW, ICURCOL )
                         DO 100 K = 1, IB
                            WRITE( NOUT, FMT = 9999 ) WORK( K )
      100                CONTINUE
                      END IF
                   END IF
                   IF( MYROW.EQ.ICURROW )
         $            II = II + IB
                   ICURROW = MOD( ICURROW+1, NPROW )
                   CALL BLACS_BARRIER( ICTXT, 'All' )
      110       CONTINUE
    *
                II = IIA
                ICURROW = IAROW
      120    CONTINUE
    *
             IF( MYCOL.EQ.ICURCOL )
         $      JJ = JJ + JB
             ICURCOL = MOD( ICURCOL+1, NPCOL )
             CALL BLACS_BARRIER( ICTXT, 'All' )
    *
      130 CONTINUE
    *
          IF( MYROW.EQ.IRWRIT .AND. MYCOL.EQ.ICWRIT ) THEN
             CLOSE( NOUT )
          END IF
    *
     9999 FORMAT( D30.18 )
    *
          RETURN
    *
    *     End of PDLAWRITE
    *
          END