Code

      program main
c**********************************************************************
c   pi.f - compute pi by integrating f(x) = 4/(1 + x**2)     
c     
c   Each node: 
c    1) receives the number of rectangles used in the approximation.
c    2) calculates the areas of it's rectangles.
c    3) Synchronizes for a global summation.
c   Node 0 prints the result.
c
c  Variables:
c
c    pi  the calculated result
c    n   number of points of integration.  
c    x           midpoint of each rectangle's interval
c    f           function to integrate
c    sum,pi      area of rectangles
c    tmp         temporary scratch space for global summation
c    i           do loop index
c****************************************************************************
      implicit  none
      include 'mpif.h'
      double precision  mypi, pi 
      integer             ::  nintvl
     $                      , nerr       ! flag for error/abort
     $                      , myid
     $                      , numprocs
     $                      , ierr
     $                      , rc
      integer, dimension(2)  :: nbcast
      integer, parameter  ::  lunerr = 0
c------------------------------------------------------
      call MPI_INIT( ierr )
      call MPI_COMM_RANK( MPI_COMM_WORLD, myid, ierr )
      call MPI_COMM_SIZE( MPI_COMM_WORLD, numprocs, ierr )
cd    print *, "Process ", myid, " of ", numprocs, " is alive"
 10   continue
      if ( myid .eq. 0 )  
     $   call get( nbcast(1), nbcast(2))
      write(lunerr,*) myid,' > MPI_BCAST'
      call MPI_BCAST(nbcast, 2, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr )
cd    print *,myid, ' bcast ', nbcast
      nintvl = nbcast(1)
      nerr   = nbcast(2)
c==                                 check for quit signal
      if( nintvl < 1 )            then
          if( nintvl == 0  .or.
     $        myid   == 0      )  goto  30
                                  goto  10   ! slaves > BCAST
                                  endif
cd    print *,'Proc',myid,'Processing ',nintvl,' intervals '      
      if( abs(nerr) > 0 .and. abs(nerr) <= nintvl )
     $  print *,'Proc',myid,' error value ',nerr
c==                               the calculation..
      mypi = myfpi( nintvl, myid, numprocs, nerr )
c==                               collect all the partial sums
      call MPI_REDUCE(mypi,pi,1,MPI_DOUBLE_PRECISION,MPI_SUM,0,
     $     MPI_COMM_WORLD,ierr)
c                                 node 0 prints the answer.
cd    print *,'Proc ',myid,' mypi = ',mypi
      if (myid .eq. 0) call myres(pi)
      goto 10
 30   continue
c==   if nintvl != 0, then slaves are waiting in MPI_BCAST
      if( nintvl >= -1 )                          then
        write(lunerr,*) myid, '> MPI_FINALIZE'
        call MPI_FINALIZE(rc)                     
        write(lunerr,*) myid, '  MPI_FINALIZE > ',rc
                                                  else 
        write(lunerr,*) myid, '> MPI_ABORT'
        call MPI_ABORT(MPI_COMM_WORLD, nintvl, rc)
        write(lunerr,*) myid, '  MPI_ABORT > ',rc
                                                  endif
      stop
      contains
      real*8 function myfpi( nint, myid, numprocs, nerr )
*..    Calculate pi..
*
      implicit  none
      integer, intent(in)    :: nint        !  number of intervals
     &                         ,myid        !  processor number
     &                         ,numprocs
     &                         ,nerr        !  flag for error
*
      real, dimension(10)    :: anarray     !  .. for bounds violation
      double precision       :: h           ! interval
     &                         ,sum, x 
      integer                :: i
c                                 calculate the interval size
cd    print *,'myfpi',myid,nerr
      h = 1.0d0/nint
      anarray(:)  = h
      sum  = 0.0d0
      do i = myid+1, nint, numprocs
         if ( i == nerr )                     then
            write(lunerr,*)  myid,' > Segv '
            x = sqrt(anarray(-nint))
         elseif( i == -nerr )                 then
           if( myid /= 0 )                  then
             write(lunerr,*)  myid,' > MPI_ABORT '
             call MPI_ABORT( MPI_COMM_WORLD, 
     $                      1, nerr )
             write(lunerr,*)  myid,'   MPI_ABORT >',nerr
           else
             write(lunerr,*)  myid,' > MPI_FINALIZE '
             call MPI_FINALIZE(nerr )
             write(lunerr,*)  myid,'   MPI_FINALIZE >',nerr
           endif
         else
            x = h * (dble(i) - 0.5d0)
                                              endif
         sum = sum + ( 4.d0 / (1.d0 + x*x))
cd       print *,myid,i,sum
      enddo
      myfpi = h * sum
*
      return
      end  function myfpi
      subroutine get( nintv, nerrv )
      integer, intent(out)     ::  nintv    ! # of intervals
     $                           , nerrv    ! flag to provoke error
*
      write(6,90)
 90   format(
     $  '=========================================================')
      write(6,98)
 98   format(
     $  'Enter the number of intervals: '/
     $  '      or  Quit :               '/
     $  '          Master   / Slaves    '/
     $  '      0 > FINALIZE / FINALIZE  '/
     $  '     -1 > FINALIZE / BCAST     '/
     $  '     -2 > ABORT    / BCAST     '  )
      read(5,*) nintv
      if( nintv > 0 )      then
           write(6,99) nintv, -nintv
           read(5,*)  nerrv
                           endif
 99   format(
     $ 'Enter 0 < nerr < ',i12,' to provoke segv '/
     $ '   or 0 > nerr > ',i12,' to call mpi_abort    ( slave  )'/
     $ '                 ',12x,'         mpi_finalize ( master ) ')
      write(6,90)
      end subroutine get
      subroutine myres(pi)
*..   print the result
      real*8  ::     pi
         write(6, 6000) pi, abs(pi - 4*datan(1.d0))
 6000    format('  pi is approximately: ', F18.16 /
     +          '  Error is           : ', F18.16)      
      end subroutine myres
      end