PROGRAM dotProductMPI
!
! This program computes the dot product of two vectors X,Y 
! (each of size N) with component i having value i
! in parallel using P processes.
! Vectors are initialized in the code by the root process,
! then statically distributed in blocks to all processes.
! It is not assumed N is divisible by P.
!
  INCLUDE 'mpif.h'
  ! variable declarations

  INTEGER, PARAMETER :: N = 100
  REAL, PARAMETER :: ROOT = 0

  INTEGER :: P, NBAR 
  INTEGER :: RANK, I, EXTRA, INDEX, OFFSET = 0
  INTEGER :: IERR
  REAL :: X(N), Y(N)
  REAL :: DOT, DOT_LOC = 0.0

  ! initialize MPI
  CALL MPI_INIT(IERR)
  IF (IERR.NE.MPI_SUCCESS) THEN
     PRINT*, "ERROR: MPI not initialized."
     STOP
  ENDIF
  
  ! Get the number of processors:
  CALL MPI_COMM_SIZE(MPI_COMM_WORLD, P, IERR)
  IF (IERR.NE.MPI_SUCCESS) THEN
     PRINT*, "ERROR: MPI processors not established."
     STOP
  ENDIF

  ! Get ranks of processors:
  CALL MPI_COMM_RANK(MPI_COMM_WORLD, RANK, IERR)
  IF (IERR.NE.MPI_SUCCESS) THEN
     PRINT*, "ERROR: MPI ranks not established."
     STOP
  ENDIF

  ! Root process initializes vectors x,y and distributes them
  IF (RANK.EQ.ROOT) THEN
     DO 10 I=1,N
        X(I) = I
        Y(I) = I
10   END DO
  ENDIF
  ! this could probably be done more efficiently by packing x and y
  ! into one entity and broadcasting it
  CALL MPI_BCAST(X, N, MPI_REAL, ROOT, MPI_COMM_WORLD, IERR)
  IF (IERR.NE.MPI_SUCCESS) THEN
     PRINT*, "ERROR: MPI_BCAST not successful."
     STOP
  ENDIF
  CALL MPI_BCAST(Y, N, MPI_REAL, ROOT, MPI_COMM_WORLD, IERR)
  IF (IERR.NE.MPI_SUCCESS) THEN
     PRINT*, "ERROR: MPI_BCAST not successful."
     STOP
  ENDIF

  ! determine which block of data to work on and compute dot product
  NBAR  = N/P
  EXTRA = MOD(N,P)
  IF (RANK < EXTRA) OFFSET = 1
  DO 20 I=1,NBAR+OFFSET
     INDEX = RANK*NBAR + I + MIN(EXTRA,RANK)
     DOT_LOC = DOT_LOC + X(INDEX)*Y(INDEX)
20 END DO

  ! gather and reduce the data and print the result
  CALL MPI_REDUCE(DOT_LOC, DOT, 1, MPI_REAL, MPI_SUM, ROOT, &
       MPI_COMM_WORLD, IERR)
  IF (RANK.EQ.ROOT) THEN
     IF (IERR.NE.MPI_SUCCESS) THEN
        PRINT*, "ERROR: MPI_REDUCE not successful."
        STOP
     ENDIF
     PRINT*, 'The dot product is: ', DOT
     PRINT*, 'The answer should be: ', N*(N+1)*(2*N+1)/6
  ENDIF
  
  ! Finalize MPI:
  CALL MPI_FINALIZE(IERR)
  IF (IERR.NE.MPI_SUCCESS) THEN
     PRINT*, "ERROR: MPI not finalized."
     STOP
  ENDIF
  
END PROGRAM dotProductMPI
