program jacobi_parallel
implicit none
include "mpif.h"
real, dimension(0:10001) :: x,newx, z
real :: dx2
integer :: n, noiters, i, k
integer :: p, me, ln, lm, tag, ierr
integer, dimension(MPI_STATUS_SIZE) :: status
! Ввод исходных данных
 open(unit = 12, file = "laplace1d.in")
! Number of cells
 read(12,  *)  n
! Number of iterations
 read(12,  *)  noiters
 close(12)

dx2 = (1. / n)**2

call MPI_Init(ierr)

call MPI_Comm_rank(MPI_COMM_WORLD, me, ierr)
call MPI_Comm_size(MPI_COMM_WORLD, p, ierr)

tag = 0
ln = n / p

do i = 1, ln
 x(i) = 1.0
enddo

if(me == 0) then
x(0) = 0.0; lm = 0
else
lm = ln * me
endif

if(me == p - 1) then
x(ln + 1) = 0.0
endif
do k = 1, noiters
! Пересылки граничных значений
 if(me - 1 >= 0) call MPI_Send(newx(1), 1, MPI_REAL, me - 1, tag, MPI_COMM_WORLD, ierr)
 if(me + 1 < p)  call MPI_Recv(x(ln + 1), 1, MPI_REAL, me + 1, tag, MPI_COMM_WORLD, status, ierr)
 tag = tag + 1
 if(me + 1 < p)  call MPI_Send(newx(ln), 1, MPI_REAL, me + 1, tag, MPI_COMM_WORLD, ierr)
 if(me - 1 >= 0) call MPI_Recv(x(0), 1, MPI_REAL, me - 1, tag, MPI_COMM_WORLD, status, ierr)
 tag = tag + 1
! Итерации Якоби
 do i = 1, ln
  newx(i) = 0.5 * (x(i - 1) + x(i + 1))
  !newx(i) = 0.5 * (x(i - 1) + x(i + 1) - dx2 * x(i))
 enddo

 do i = 1, ln
  x(i) = newx(i)
 enddo
enddo
! Собираем решение
if(me == 0) then
  do i = 1, ln
   z(i) = x(i)
  enddo
  do k = 1, p - 1
   lm = ln * k
   call MPI_Recv(z(lm), ln, MPI_REAL, k, k, MPI_COMM_WORLD, status, ierr)  
  enddo
 else
   call MPI_Send(x(1), ln, MPI_REAL, 0, me, MPI_COMM_WORLD, ierr)
 endif

call MPI_Finalize(ierr)
! Запись результата в файл
if(me == 0) then
 open(unit = 11, file = "laplace1d_parallel.dat", status = "NEW")
  write(11, "(2x, e8.3)") (z(i),   i = 1, n)
 close(11)
endif
end



