program picom
!
!*******************************************************************************
!
!! PICOM uses MPI routines to multiprocess a computational task.
!
!
!  Discussion:
!
!    The antiderivative of Arctan(X) is 1 / ( 1 + X**2 ).
!
!    Therefore, the integral from 0 to 1 of 1 / ( 1 + X**2 ) is
!    Arctan(1) - Arctan(0) = PI/4 - 0.
!
!    If we approximate the integral from 0 to 1 of 4 / ( 1 + X**2 ),
!    then the exact value will be PI.
!
!    The integral is estimated as the sum of N terms, each of which
!    is (1/n) * f(x).
!
!    The I-th term is evaluated at the center of the I-th interval, 
!    so X(I) = ( I - 1/2) / N.
!
!    Processor MYID is to compute the sum of the terms
!    MYID+1, MYID+1+NUMPROCS, MYID+1+2*NUMPROCS, ...
!
!  Reference:
!
!    Gropp, Lusk, Skjellum,
!    Using MPI,
!    Portable Parallel Programming with the Message-Passing Interface,
!    MIT Press, 1997.
!
!    Snir, Otto, Huss-Lederman, Walker, Dongarra,
!    MPI - The Complete Reference,
!    Volume 1, The MPI Core,
!    second edition,
!    MIT Press, 1998.
!
!  Modified:
!
!    02 September 1998
!
  implicit none
!
  include 'mpif.h'
!
  double precision, parameter :: pi25dt = 3.141592653589793238462643D+00
!
  double precision f
  double precision h
  integer i
  integer ierr
  integer myid
  double precision mypi
  integer n
  integer numprocs
  double precision pi
  double precision sum2
  double precision x
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'PICOM'
  write ( *, '(a)' ) '  A program using MPI, to compute PI.'
!
!  Establish the MPI environment.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Call MPI_INIT to estabish the environment:'

  call mpi_init ( ierr )

  if ( ierr /= mpi_success ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'PICOM: Warning!'
    write ( *, '(a,i6)' ) '  MPI_INIT returns IERR = ', ierr
    call mpi_finalize ( ierr )
    stop
  end if
!
!  Get this process's ID.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Call MPI_COMM_RANK to get this process''s ID:'

  call mpi_comm_rank ( mpi_comm_world, myid, ierr )

  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) 'Process ', myid, ' is active.'
!
!  Find out how many processes are available.
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Call MPI_COMM_SIZE to get the number of processes:'

  call mpi_comm_size ( mpi_comm_world, numprocs, ierr )

  write ( *, '(a,i6)' ) 'Number of processes is ', numprocs
!
!  Assume that we just got the value of N from the user.
!
  n = 100
!
!  Broadcast the value of N.
!
  call mpi_bcast ( n, 1, mpi_integer, 0, mpi_comm_world, ierr )
!
!  Check for quit signal.
!
  if ( n <= 0 ) then
    call mpi_finalize ( ierr )
    stop
  end if
!
!  Process MYID now adds up its terms.
!
  h = 1.0D+00 / dble ( n )
  sum2 = 0.0D+00
  do i = myid+1, n, numprocs
    x = h * ( dble ( i ) - 0.5D+00 )
    sum2 = sum2 + f ( x )
  end do

  mypi = h * sum2
!
!  All the partial sums are collected.
!
  call mpi_REDUCE ( mypi, pi, 1, mpi_DOUBLE_PRECISION, mpi_sum, 0, &
    mpi_COMM_WORLD, ierr )
!
!  Process 0 prints the answer.
!
  if ( myid == 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a,g14.6)' ) 'Estimate for PI is ', pi
    write ( *, '(a,g14.6)' ) 'Error is           ', pi - PI25DT
  end if
!
!  Finish up.
!
  call mpi_finalize ( ierr )

  stop
end
function f ( x )
!
!***********************************************************************
!
!! F is the function we are integrating.
!
!
!  Modified:
!
!    10 February 2000
!
!  Parameters:
!
!    Input, double precision X, the argument of the function.
!
!    Output, double precision F, the value of the function.
!
  implicit none
!
  double precision f
  double precision x
!
  f = 4.0D+00 / ( 1.0D+00 + x**2 )

  return
end