program main
c This program calculates the variance of the first 10**9 integers
c using MPI with multiple dual processor Macintoshes.
      implicit none
c get definition of MPI constants
      include 'mpif.h'
      integer TickCount
      external TickCount
      integer NINT
      parameter(NINT=1000000000)
      integer n, n1, n2, start1, start2
      integer start_program, end_program
      double precision sum11, sum12, sum21, sum22, mean, mesq
      external adder
      integer ierror, nproc, idproc, root, idtask, nargs
      double precision partial(2), result(2)
      root = 0
      nargs = 4
c initialize the MPI environment
      call MPI_INIT(ierror)
      call MPI_COMM_SIZE(MPI_COMM_WORLD,nproc,ierror)
      call MPI_COMM_RANK(MPI_COMM_WORLD,idproc,ierror)
      start_program = TickCount()
c set arguments for each node, last node might have more work
      n = NINT/nproc
      start1 = n*idproc
      n1 = n + start1
      if (idproc.eq.(nproc-1)) n1 = NINT

c initialize the MP environment
      call MP_INIT(nproc)
c set arguments for each task
      n = n/2
      n2 = n1
      n1 = n1 - n
      start2 = start1 + n

      call MP_TASKSTART(idtask,adder,nargs,n2,start2,sum21,sum22)
      call adder(n1,start1,sum11,sum12)
      call MP_TASKWAIT(idtask)

c terminate the MP environment
      call MP_END()

c add up all the partial sums from each node
      partial(1) = sum11 + sum21
      partial(2) = sum12 + sum22
      call MPI_REDUCE(partial,result,2,MPI_DOUBLE_PRECISION,MPI_SUM,root
     &,MPI_COMM_WORLD,ierror)

c terminate MPI environment
      call MPI_BARRIER(MPI_COMM_WORLD,ierror)
      call MPI_FINALIZE(ierror)

c node 0 has the results
      if (idproc.eq.root) then
         n = NINT
         mean = result(1) / n
         mesq = result(2) / n
         write (*,'(i15,1x,f20.3,1x,f21.3)') n, mean, mesq-(mean * mean)
      endif

      end_program = TickCount()

      write (*,*)
      write (*,'("time in sec: ",f7.2)') (end_program-start_program)/60.
c node 0 waits for carriage return
      if (idproc.eq.root) then
         write (*,*) 'hit carriage return to quit'
         pause
      endif
      stop
      end


      subroutine adder(n,start,s1,s2)
      integer n, start
      double precision s1, s2
      integer i
      double precision ls1, ls2
      ls1 = 0
      ls2 = 0
      do 10 i = start, n-1
      ls1 = ls1 + dble(i)
      ls2 = ls2 + dble(i)*i
   10 continue
      s1 = ls1
      s2 = ls2
      return
      end