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