BLACS Examples

Example (sub)routines:
----------------------------------------------------------------------------

HELLO WORLD

The following routine takes the available processes, forms them into a
process grid, and then has each process check in with the process at {0,0}
in the process grid.

      PROGRAM HELLO
*     -- BLACS example code --
*     Written by Clint Whaley 7/26/94
*     Performs a simple check-in type hello world
*     ..
*     .. External Functions ..
      INTEGER KPNUM
      EXTERNAL KPNUM
*     ..
*     .. Variable Declaration ..
      INTEGER IAM, NPROCS, NPROW, NPCOL, MYPROW, MYPCOL
      INTEGER ICALLER, I, J, HISROW, HISCOL
*
*     Determine my process number and the number of processes in
*     machine
*
      CALL INITPINFO(IAM, NPROCS)
*
*     If in PVM, create virtual machine if it doesn't exist
*
      IF (NPROCS .LT. 1) THEN
         IF (IAM .EQ. 0) THEN
            WRITE(*, 1000)
            READ(*, 2000) NPROCS
         END IF
         CALL AUXSETUP(IAM, NPROCS)
      END IF
*
*     Set up process grid that is as close to square as possible
*
      NPROW = INT( SQRT( REAL(NPROCS) ) )
      NPCOL = NPROCS / NPROW
      CALL BLACSINIT(NPROW, NPCOL)
      CALL GRIDINFO(NPROW, NPCOL, MYPROW, MYPCOL)
*
*     If I'm not in grid, go to end of program
*
      IF ( (MYPROW.GE.NPROW) .OR. (MYPCOL.GE.NPCOL) ) GOTO 30
*
*     Get my process ID from my grid coordinates
*
      ICALLER = KPNUM(MYPROW, MYPCOL)
*
*     If I am process {0,0}, receive check-in messages from
*     all nodes
*
      IF ( (MYPROW.EQ.0) .AND. (MYPCOL.EQ.0) ) THEN

         WRITE(*,*) ' '
         DO 20 I = 0, NPROW-1
            DO 10 J = 0, NPCOL-1

               IF ( (I.NE.0) .OR. (J.NE.0) ) THEN
                  CALL IGERV2D(1, 1, ICALLER, 1, I, J)
               ENDIF
*
*              Make sure ICALLER is where we think in process grid
*
               CALL PCOORD(ICALLER, HISROW, HISCOL)
               IF ( (HISROW.NE.I) .OR. (HISCOL.NE.J) ) THEN
                  WRITE(*,*) 'Grid error!  Halting . . .'
                  STOP
               END IF
               WRITE(*, 3000) I, J, ICALLER

10          CONTINUE
20       CONTINUE
         WRITE(*,*) ' '
         WRITE(*,*) 'All processes checked in.  Run finished.'
*
*     All processes but {0,0} send process ID as a check-in
*
      ELSE
         CALL IGESD2D(1, 1, ICALLER, 1, 0, 0)
      END IF

30    CONTINUE

      CALL BLACSEXIT(0)

1000  FORMAT('How many processes should be spawned?')
2000  FORMAT(I)
3000  FORMAT('Process {',i2,',',i2,'} (node number =',I,
     $       ') has checked in.')

      STOP
      END

----------------------------------------------------------------------------

PARALLEL DOT PRODUCT

This routine does a bone-headed parallel double precision dot product of two
vectors. Arguments are input on process {0,0}, and output everywhere else.

      DOUBLE PRECISION FUNCTION PDDOT( N, X, Y )
*
*     -- BLACS example code --
*     Written by Clint Whaley 7/26/94.
*     ..
*     .. Scalar Arguments ..
      INTEGER N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION X(*), Y(*)
*     ..
*
*  Purpose
*  =======
*  PDDOT is a restricted parallel version of the BLAS routine
*  DDOT.  It assumes that the increment on both vectors is one,
*  and that process {0,0} starts out owning the vectors and
*  has N.  It returns the dot product of the two N-length vectors
*  X and Y, i.e. PDDOT = X' Y.
*
*  Arguments
*  =========
*  N            (input/ouput) INTEGER
*               The length of the vectors X and Y.  Input
*               for {0,0}, output for everyone else.
*
*  X            (input/output) DOUBLE PRECISION array of dimension (N)
*               The vector X of PDDOT = X' Y.  Input for {0,0},
*               output for everyone else.
*
*  Y            (input/output) DOUBLE PRECISION array of dimension (N)
*               The vector Y of PDDOT = X' Y.  Input for {0,0},
*               output for everyone else.
*
*  =====================================================================
*     ..
*     .. External Functions ..
      DOUBLE PRECISION DDOT
      EXTERNAL DDOT
*     ..
*     .. External Subroutines ..
      EXTERNAL GRIDINFO, DGEBS2D, DGEBR2D, DGSUM2D
*     ..
*     .. Local Scalars ..
      INTEGER IAM, NPROCS, NPROW, NPCOL, MYPROW, MYPCOL, I, LN
      DOUBLE PRECISION LDDOT
*     ..
*     .. Executable Statements ..
*
*     Find out what grid has been set up, and pretend it's 1-D
*
      CALL GRIDINFO( NPROW, NPCOL, MYPROW, MYPCOL )
      IAM = MYPROW*NPCOL + MYPCOL
      NPROCS = NPROW * NPCOL
*
*     Do bone-headed thing, and just send entire X and Y to
*     everyone
*
      IF ( (MYPROW.EQ.0) .AND. (MYPCOL.EQ.0) ) THEN
         CALL IGEBS2D('All', 'i-ring', 1, 1, N, 1 )
         CALL DGEBS2D('All', 'i-ring', N, 1, X, N )
         CALL DGEBS2D('All', 'i-ring', N, 1, Y, N )
      ELSE
         CALL IGEBR2D('All', 'i-ring', 1, 1, N, 1, 0, 0 )
         CALL DGEBR2D('All', 'i-ring', N, 1, X, N, 0, 0 )
         CALL DGEBR2D('All', 'i-ring', N, 1, Y, N, 0, 0 )
      ENDIF
*
*     Find out the number of local rows to multiply (LN), and
*     where in vectors to start (I)
*
      LN = N / NPROCS
      I = 1 + IAM * LN
*
*     Last process does any extra rows
*
      IF (IAM .EQ. NPROCS-1) LN = LN + MOD(N, NPROCS)
*
*     Figure dot product of my piece of X and Y
*
      LDDOT = DDOT( LN, X(I), 1, Y(I), 1 )
*
*     Add local dot products to get global dot product;
*     give all procs the answer
*
      CALL DGSUM2D( 'All', '1-tree', 1, 1, LDDOT, 1, -1, 0 )

      PDDOT = LDDOT

      RETURN
      END

----------------------------------------------------------------------------

PARALLEL MATRIX INFINITY NORM

The routine does a parallel infinity norm on a distributed double precision
matrix. Unlike the PDDOT example, this routine assumes the matrix has
already been distributed.

      DOUBLE PRECISION FUNCTION PDINFNRM(LM, LN, A, LDA, WORK)
*
*     -- BLACS example routine --
*     Written by Clint Whaley.
*
*     .. Scalar Arguments ..
      INTEGER LM, LN, LDA
*
*     .. Array Arguments ..
      DOUBLE PRECISION A(LDA, *), WORK(*)
*
*  PURPOSE:
*  ========
*
*  Compute the infinity norm of a distributed matrix, where
*  the matrix is spread across a 2D process grid.  The result is
*  left on all processes.
*
*  Arguments
*  =========
*
*  LM        (input) INTEGER
*            Number of rows of the global matrix owned by this process.
*
*  LN        (input) INTEGER
*            Number of columns of the global matrix owned by this process.
*
*  A         (input) DOUBLE PRECISION, dimension (LDA,N)
*            The matrix who's norm you wish to compute.
*
*  LDA       (input) INTEGER
*            Leading Dimension of A.
*
*  WORK      (temporary) DOUBLE PRECISION array, dimension (LM)
*            Temporary work space used for summing rows.
*
*     .. External Subroutines ..
      EXTERNAL DGSUM2D, DGMAX2D, GRIDINFO
*
*     .. External Functions ..
      INTEGER IDAMAX
      DOUBLE PRECISION DASUM
*
*     .. Local Scalars ..
      INTEGER NPROW, NPCOL, MYROW, MYCOL,  I, J
      DOUBLE PRECISION MAX
*
*     .. Executable Statements ..
*
*     Get process grid information
*
      CALL GRIDINFO(NPROW, NPCOL, MYROW, MYCOL)
*
*     Add all local rows together
*
      DO 20 I = 1, LM
         WORK(I) = DASUM(LN, A(I,1), LDA)
20    CONTINUE
*
*     Find sum of global matrix rows and store on column 0 of
*     process grid
*
      CALL DGSUM2D('Row', '1-tree', LM, 1, WORK, LM, MYROW, 0)
*
*     Find maximum sum of rows for supnorm
*
      IF (MYCOL .EQ. 0) THEN
         MAX = WORK(IDAMAX(LM,WORK,1))
         IF (LM .LT. 1) MAX = 0.0D0
         CALL DGMAX2D('Col', 'h', 1, 1, MAX, 1, I, J, 1, -1, 0)
      END IF
*
*     Process column 0 has answer; send answer to all nodes
*
      IF (MYCOL.EQ.0) THEN
         CALL DGEBS2D('Row', 'hyp', 1, 1, MAX, 1)
      ELSE
         CALL DGEBR2D('Row', 'hyp', 1, 1, MAX, 1, 0, 0)
      END IF
*
      PDINFNRM = MAX
*
      RETURN
*
*     End of PDINFNRM
*
      END

----------------------------------------------------------------------------

PLAYING WITH MESSAGE IDS

The following piece of code reserves message IDs for the programmer's use.
Please note that to do this we first call IDRANGES to get the initial range,
and then shift the range upwards by the correct amount. We do this because
the initial range may vary between platforms. Thus, if we hardcoded the
reserved range to be from 0 to 99, for instance, this might be valid on one
platform, but not on others. The way it is done here is safer.

THIS ROUTINE NOT TESTED

      SUBROUTINE GETIDS(NIDS, MYBEGRNG, MYENDRNG)
*
*     -- BLACS example code --
*     Written by Clint Whaley 7/27/94.
*     ..
*     .. Scalar Arguments ..
      INTEGER NIDS, MYBEGRNG, MYENDRNG
*     ..
*
*  Purpose
*  =======
*  GETIDS reserves NIDS message IDs for the user.  The reserved
*  IDs are on a continuous range, who's beginning is returned in
*  MYBEGRNG, and who's ending is returned in MYENDRNG.  Note:
*  this routine should be called _before_ BLACSINIT or GRIDMAP.
*
*  Arguments
*  =========
*  NIDS         (input) INTEGER
*               The number if IDs to reserve.
*
*  MYBEGRNG     (output) INTEGER
*               The beginning of the reserved range.
*
*  MYENDRNG     (output) INTEGER
*               The ending of the reserved range.
*
*  =====================================================================
*     ..
*     .. External Subroutines ..
      EXTERNAL IDRANGES, SHIFT_RANGE
*     ..
*     .. Local Scalars ..
      INTEGER BEGRNG1, ENDRNG1, BEGRNG1, ENDRNG2
*     ..
*     .. Executable Statements ..
*
*     Get system-dependant initial ID ranges
*
      CALL IDRANGES( BEGRNG1, ENDRNG1, BEGRNG2, ENDRNG2 )
*
*     Reserve IDs for my own use
*
      MYBEGRNG = BEGRNG1
      MYENDRNG = MYBEGRNG + NIDS
      BEGRNG1 = MYENDRNG + 1
      CALL SHIFT_RANGE( BEGRNG1, ENDRNG1, BEGRNG2, ENDRNG2 )

      RETURN
      END

----------------------------------------------------------------------------

 PROCMAP 

This routine maps processes to a grid using gridmap.

THIS ROUTINE NOT TESTED

      SUBROUTINE PROCMAP( MAPPING, NPROW, NPCOL, IMAP )
*
*     -- BLACS example code --
*     Written by Clint Whaley 7/26/94.
*     ..
*     .. Scalar Arguments ..
      INTEGER MAPPING, NPROW, NPCOL
*     ..
*     .. Array Arguments ..
      INTEGER IMAP(NPROW, *)
*     ..
*
*  Purpose
*  =======
*  PROCMAP maps the first NPROW*NPCOL processes to the grid in
*  a variety of ways depending on the parameter MAPPING.
*
*  Arguments
*  =========
*  MAPPING      (input) INTEGER
*               Way to map processes to grid.  Choices are:
*               1 : row-major natural ordering
*               2 : column-major natural ordering
*
*  NPROW        (input) INTEGER
*               The number of process rows the created grid
*               should have.
*
*  NPCOL        (input) INTEGER
*               The number of process columns the created grid
*               should have.
*
*  IMAP         (workspace) INTEGER array of dimension (NPROW, NPCOL)
*               Workspace where the array which maps the
*               processes to the grid will be stored for the
*               call to GRIDMAP.
*
*  =====================================================================
*
*     ..
*     .. External Functions ..
      INTEGER  KPNUM
      EXTERNAL KPNUM
*     ..
*     .. External Subroutines ..
      EXTERNAL INITPINFO, BLACSINIT, GRIDMAP
*     ..
*     .. Local Scalars ..
      INTEGER NPROCS, I, J, K
*     ..
*     .. Executable Statements ..
*
*     See how many processes there are in the system
*
      CALL INITPINFO( I, NPROCS )
      IF (NPROCS .LT. NPROW*NPCOL) THEN
         WRITE(*,*) 'Not enough processes for grid'
         STOP
      ENDIF
*
*     Temporarily map all processes into 1 x NPROCS grid
*
      CALL BLACSINIT( 1, NPROCS )
*
*     If we want a row-major natural ordering
*
      IF (MAPPING .EQ. 1) THEN
         K = 1
         DO I = 1, NPROW
            DO J = 1, NPCOL
               IMAP(I, J) = KPNUM(0, K)
               K = K + 1
            END DO
         END DO
*
*     If we want a column-major natural ordering
*
      ELSE IF (MAPPING .EQ. 2) THEN
         DO J = 1, NPCOL
            DO I = 1, NPROW
               IMAP(I, J) = KPNUM(0, K)
               K = K + 1
            END DO
         END DO
      ELSE
         WRITE(*,*) 'Uknown mapping.'
         STOP
      END IF
*
*     Apply the new mapping
*
      CALL GRIDMAP( IMAP, NPROW, NPROW, NPCOL )

      RETURN
      END

This material was reproduced for educational use only from the original document located at: http://www.netlib.org/blacs/BLACS/Examples.html.