       program hypersum
       implicit none
       include 'mpif.h'
       integer myProc, nProc
       common /mpidat/myProc, nProc
       integer pvar, dim, isum, isize, temp, ierror, cubedim
       integer itype, log2

       itype = MPI_INTEGER
       isize = 1 
       call initmpi
       isum = myProc
       cubedim = log2(nProc)
c      print*,' ok ',myProc, cubedim, isum

c book algorithm starts here
       pvar=1
       do 10 dim = 1, cubedim
           call snd(dim,isum,itype,xor(myProc,pvar), isize)
           call rcv(dim,temp,itype,xor(myProc,pvar), isize)
           isum=isum+temp
           pvar=pvar *2
10     continue
c and it ends here

       if (myProc.EQ.0) print*,'isum= ', isum
       call MPI_FINALIZE(ierror)
       end


      subroutine snd(tag,value,type,dst,size)
      implicit none
      integer tag, value, type, dst, size, ierror
      include 'mpif.h'
      integer myProc, nProc
      common /mpidat/myProc, nProc
  
      call MPI_Send(value, size, type, dst, tag,  
     +             MPI_COMM_WORLD, ierror)

      return
      end

      subroutine rcv(tag,value,type,src,size)
      implicit none
      integer tag, value, type, src, size
      include 'mpif.h'
      integer myProc, nProc
      common /mpidat/myProc, nProc
      integer istatus(2), ierror(2)
      call MPI_Recv(value, size, type, src, tag, 
     +             MPI_COMM_WORLD, istatus, ierror)

      return
      end

      subroutine initmpi
      implicit none
      include 'mpif.h'
      integer myProc, nProc
      common /mpidat/myProc, nProc
      integer ierror
      call MPI_INIT(ierror)
      call MPI_COMM_RANK(MPI_COMM_WORLD, myProc, ierror)
      call MPI_COMM_SIZE(MPI_COMM_WORLD, nProc, ierror)
      return
      end

      integer function log2(n)
      implicit none
      integer n
      double precision dn
      dn = n
      log2 = log10(dn) / log10(2.D0)
      return
      end
 
