PROGRAM PI
*
*     -- BLACS example code --
*     University of Tennessee, Knoxville
*
*     Written by Antoine Petitet, August 1995 (petitet@cs.utk.edu)
*     Edited for use in instructional labs for JICS by Chris Hastings
*     7/15/96 (hastings@cs.utk.edu)
*
*     This program calculates the value of pi, using numerical
*     integration with parallel processing, and clocks the solution
*     time. The user selects the number of processors and the number
*     of points of integration.  By selecting and timing different
*     grid sizes you obtain a measure of the speedup with perfectly
*     parallel problems.
*
*     .. Local Scalars ..
      INTEGER            ICTXT, IAM, NPROCS, NPROW, NPCOL, MYPROW,
     $                   MYPCOL, I, NPOINTS, EXTRAPOINTS, MYPOINTS,
     $                   BASICPOINTS
      REAL               A, B, X, SLICESIZE, PARTIALINT, MYA, MYB,
     $                   CTIME, WTIME
*     ..
*     .. External Subroutines ..
      EXTERNAL           BLACS_GET, BLACS_PINFO, BLACS_SETUP,
     $                   BLACS_GRIDINIT, BLACS_GRIDINFO, SLTIMER,
     $                   SLBOOT, SLCOMBINE, SGSUM2D
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MOD
*     ..
*     .. External Functions ..
      REAL               F
      EXTERNAL           F
*
*     Determine my process number and the number of processes in
*     machine
*
      CALL BLACS_PINFO( IAM, NPROCS )
*
      A = 0.0E+0
      B = 1.0E+0
*
      IF( IAM.EQ.0 ) THEN
*
         WRITE( *, FMT = * )
         WRITE( *, FMT = * ) '******************** ',
     $            'INTEGRATION EXAMPLE ***********************'
         WRITE( *, FMT = * )
         WRITE( *, FMT = * ) 'This Example calculates pi by ',
     $            'integrating the function:'
         WRITE( *, FMT = * ) 'f(x) = 4 / (1 + x**2) between x=0 ',
     $            'and x=1.'
         WRITE( *, FMT = * ) 'The method used is the n-point ',
     $            'rectangle quadrature rule.'
         WRITE( *, FMT = * )
*
         WRITE( *, FMT = * ) 'How many points do you want (0 or ',
     $            'neg. value quits)? '
         READ( *, FMT = * ) NPOINTS
*
      END IF
*
*     If underlying system needs additional set up, do it now
*
      IF( NPROCS.LT.1 ) THEN
         IF( IAM.EQ.0 ) THEN
            WRITE( *, FMT = * ) ' How many processes should I use ? '
            READ( *, FMT = * ) NPROCS
         END IF
         CALL BLACS_SETUP( IAM, NPROCS )
      END IF
*
*     Get default system context, and define grid
*
      NPROW = 1
      NPCOL = NPROCS
      CALL BLACS_GET( -1, 0, ICTXT )
      CALL BLACS_GRIDINIT( ICTXT, 'Row-major', NPROW, NPCOL )
      CALL BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYPROW, MYPCOL )
*
*     If I'm not in grid, go to end of program
*
      IF( MYPROW.GE.NPROW .OR. MYPCOL.GE.NPCOL )
     $  GO TO 30
*
*     Broadcast problem parameters
*
      IF( MYPROW.EQ.0 .AND. MYPCOL.EQ.0 ) THEN
*
         CALL IGEBS2D( ICTXT, 'All', ' ', 1, 1, NPOINTS, 1 )
         CALL SGEBS2D( ICTXT, 'All', ' ', 1, 1, A, 1 )
         CALL SGEBS2D( ICTXT, 'All', ' ', 1, 1, B, 1 )
*
      ELSE
*
         CALL IGEBR2D( ICTXT, 'All', ' ', 1, 1, NPOINTS, 1, 0, 0 )
         CALL SGEBR2D( ICTXT, 'All', ' ', 1, 1, A, 1, 0, 0 )
         CALL SGEBR2D( ICTXT, 'All', ' ', 1, 1, B, 1, 0, 0 )
*
      END IF
*
      IF( NPOINTS.LE.0 )
     $  GO TO 30
*
      PARTIALINT = 0.0E+0
*
*     start the clock running
*
      CALL SLBOOT()
      CALL SLTIMER( 1 )
*
*     Calculate slice size:
*
      SLICESIZE = ( B - A ) / NPOINTS
*
*     Calculate # of points per process & local sub-interval:
*
      BASICPOINTS = NPOINTS / ( NPROW * NPCOL )
      EXTRAPOINTS = MOD( NPOINTS, NPROW * NPCOL )
*
      IF( IAM.LT.EXTRAPOINTS ) THEN 
         MYPOINTS = BASICPOINTS + 1
         MYA = A + SLICESIZE * IAM * MYPOINTS
      ELSE 
         MYPOINTS = BASICPOINTS
         MYA = A + SLICESIZE * (IAM * MYPOINTS + EXTRAPOINTS)
      END IF
*
      MYB = MYA + SLICESIZE * MYPOINTS
*
*     calculate the partial integral of the function over my
*     sub-interval:
*
      DO 20 I = 1, MYPOINTS
         X = MYA + ( I - 0.5D+0 )*SLICESIZE
         PARTIALINT = PARTIALINT + F( X ) * SLICESIZE
   20 CONTINUE
*
*     if we're process (0,0), gather all the contributions from all
*     the other processes, add them up, and figure elapsed time
*
      CALL SGSUM2D( ICTXT, 'All', ' ', 1, 1, PARTIALINT, 1, -1, -1 )
*
      CALL SLTIMER( 1 )
      CALL SLCOMBINE( ICTXT, 'All', '>', 'W', 1, 1, WTIME )
      CALL SLCOMBINE( ICTXT, 'All', '>', 'C', 1, 1, CTIME )
*
      IF( MYPROW.EQ.0 .AND. MYPCOL.EQ.0 ) THEN
         WRITE( *, FMT = * )
         WRITE( *, FMT = 9999 ) PARTIALINT
         WRITE( *, FMT = 9998 ) CTIME
         WRITE( *, FMT = 9997 ) WTIME
         WRITE( *, FMT = * )
      END IF
*
   30 CONTINUE
*   
      CALL BLACS_GRIDEXIT( ICTXT )
      CALL BLACS_EXIT( 0 )
*
 9999 FORMAT( ' PI is approximately : ', F18.16 )
 9998 FORMAT( ' CPU  time to solution in seconds: ', 1PE8.2 )
 9997 FORMAT( ' WALL time to solution in seconds: ', 1PE8.2 )
*
      STOP
*
      END
*
      REAL             FUNCTION F( X )
*
*     Function to integrate in the integration example.
*
      REAL             X
*
      F = 4.0E+0 / ( 1.0E+0 + X*X )
*
      RETURN 
*
      END