module mysub

use mpi

implicit none

real(kind(0.d0)) :: pi=acos(-1d0)
integer :: numprocs, proc_num

 contains 

subroutine log_likelihood(y, theta, lli, ll)
 double precision, allocatable, intent(in) :: y(:) 
 double precision, intent(in) :: theta(2)
 double precision, allocatable, intent(out) :: lli(:)
 double precision, intent(out) :: ll
 double precision, allocatable :: proc_contrib(:)
 integer :: i, size_y, diff,ierr, points_per_proc, istart, iend
 
 size_y=size(y)
 
 ALLOCATE (lli(size_y))

 !Broadcasting
 call MPI_BCAST(theta, 2, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
 call MPI_BCAST(y, size_y, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)

 ! Determine how many points to handle with each proc
 points_per_proc = (size_y + numprocs - 1)/numprocs
 print *, proc_num, points_per_proc, numprocs
 
 
 ! Determine start and end index for this proc's points
 istart = proc_num * points_per_proc + 1
 iend = min((proc_num + 1)*points_per_proc, size_y)
 diff = iend-istart+1

 ALLOCATE(proc_contrib(diff))

 do i = istart, iend
   
    proc_contrib(i-istart+1) = -log(sqrt(theta(2))) - 0.5*log(2.0d0*pi) - (1.0d0/(2.0d0*theta(2)))*((y(i)-theta(1))**2)
    
 end do
 
 call MPI_ALLGATHER(proc_contrib, diff, MPI_DOUBLE_PRECISION, lli, diff, MPI_DOUBLE_PRECISION, MPI_COMM_WORLD, ierr)
 
 ll = sum(lli)
 
 call MPI_BARRIER(MPI_COMM_WORLD,ierr)
 
end subroutine log_likelihood


end module mysub



program log_like

use mpi
use mysub

implicit none 

double precision, allocatable :: y(:) 
double precision :: theta(2)
double precision, allocatable :: lli(:)
double precision :: ll
integer :: i, ierr

call MPI_INIT(ierr)
 
call MPI_COMM_SIZE(MPI_COMM_WORLD,numprocs,ierr)
call MPI_COMM_RANK(MPI_COMM_WORLD,proc_num,ierr)

print *, numprocs
print *, proc_num

allocate(y(10))

if (proc_num == 0) then
  theta(1)=0d0
  theta(2)=10d0
  y=(/(i*1d0,i=1,10)/)
  print *, 'ok init'
end if

call log_likelihood(y,theta,lli,ll)

print *, proc_num, lli, ll
 
 call MPI_FINALIZE(ierr)
 
end program log_like