! quadrature_ex.f90
include "mpiSim.f90"
program main
include "useSim.f90"

!********************************************************************72
! Code from: http://people.scs.fsu.edu/~burkardt/f_src/mpi/mpi.html
!  MAIN is the main program for QUADRATURE.
!
!  Discussion:
!
!    QUADRATURE uses MPI routines to multiprocess a computational task.
!
!    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 I is to compute the sum of the terms
!    I+1, I+1+NUM_PROCS, I+1+2*NUM_PROCS, ...
!
!  Modified:
!
!    13 September 2002
!
!  Reference:
!
!    William Gropp, Ewing Lusk, Anthony Skjellum,
!    Using MPI: Portable Parallel Programming with the
!    Message-Passing Interface,
!    Second Edition,
!    MIT Press, 1999,
!    ISBN: 0262571323.
!
!    Snir, Otto, Huss-Lederman, Walker, Dongarra,
!    MPI - The Complete Reference,
!    Volume 1, The MPI Core,
!    second edition,
!    MIT Press, 1998.
!

!
!  Fortran77 include file:
!
! include 'mpif.h'
!
!  Fortran90 module:
!
!  use mpi
!
! implicit none
!
  integer dest
  Double Precision end_time
  Double Precision f
  Double Precision h
  integer i
  integer ierr
  integer, parameter :: master = 0
  integer my_id
  Double Precision mypi
  integer n
  integer num_procs
  Double Precision pi
  Double Precision, parameter :: pi25dt = 3.141592653589793238462643D+00
  integer source
  Double Precision start_time
  Double Precision sum2
  Double Precision x
!
!  Establish the MPI environment.
!
  call MPI_Init ( ierr )
!
!  Get this process's ID.
!
  call MPI_Comm_rank ( MPI_COMM_WORLD, my_id, ierr )
!
!  Find out how many processes are available.
!
  call MPI_Comm_size ( MPI_COMM_WORLD, num_procs, ierr )

  if ( my_id == master ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'QUADRATURE - Master process:'
    write ( *, '(a)' ) '  FORTRAN90 version'
    write ( *, '(a)' ) '  An MPI example program to compute PI.'
    write ( *, '(a,i8)' ) '  The number of processes is ', num_procs

    start_time = MPI_Wtime ( )

  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8,a)' ) 'Process ', my_id, ' is active.'
!
!  Assume that the master process just got the value of N from the user.
!  Here, we'll use an assignment statement.
!
  if ( my_id == master ) then
    n = 100
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'QUADRATURE - Master process:'
    write ( *, '(a,i8,a)' ) 'Number of intervals being used is ', n
  end if
!
!  The master process must broadcast the value of N to all processes.
!
  source = master

  call MPI_Bcast ( n, 1, MPI_INTEGER, source, MPI_COMM_WORLD, ierr )
!
!  Process MY_ID now adds up its terms.
!
  h = 1.0D+00 / real ( n, kind(0.D0))
  sum2 = 0.0D+00
  do i = my_id+1, n, num_procs
    x = h * ( real ( i, kind(0.0D0)) - 0.5D+00 )
    sum2 = sum2 + f ( x )
  end do

  mypi = h * sum2
!
!  All the partial sums are collected and summed by the master process.
!
  dest = master

  call MPI_Reduce ( mypi, pi, 1, MPI_DOUBLE_PRECISION, MPI_SUM, dest, &
    MPI_COMM_WORLD, ierr )
!
!  The master process prints the answer.
!
  if ( my_id == master ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'QUADRATURE - Master process:'
    write ( *, '(a,g14.6)' ) '  Estimate for PI is ', pi
    write ( *, '(a,g14.6)' ) '  Error is           ', pi - PI25DT

    end_time = MPI_Wtime ( )

    write ( *, '(a)' ) ' '
    write ( *, '(a,f14.6)' ) '  Elapsed wall clock seconds = ', &
      end_time - start_time

  end if
!
!  Finish up.
!
  call MPI_Finalize ( ierr )

  if ( my_id == master ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'QUADRATURE - Master process:'
    write ( *, '(a)' ) '  Normal end of execution.'
  end if

  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