MODULE jacobi_module
  IMPLICIT NONE
  INTEGER, PARAMETER :: real4 = selected_real_kind(6,37)
  INTEGER, PARAMETER :: real8 = selected_real_kind(15,307)
  REAL(real8), DIMENSION(:,:), ALLOCATABLE         :: vnew
  REAL(real8), DIMENSION(:,:), ALLOCATABLE, TARGET :: v  ! solution array
  REAL(real8) :: tol=1.d-4, del, gdel=1.0d0
  REAL(real4) :: start_time, end_time
  INCLUDE "mpif.h"  !! This brings in pre-defined MPI constants, ...
  INTEGER :: p, ierr, below, above, k, m, mp, iter=0
  PUBLIC

CONTAINS
  SUBROUTINE bc(v, m, mp, k, p)
! PDE: Laplacian u = 0;      0<=x<=1;  0<=y<=1
! B.C.: u(x,0)=sin(pi*x); u(x,1)=sin(pi*x)*exp(-pi); u(0,y)=u(1,y)=0
! SOLUTION: u(x,y)=sin(pi*x)*exp(-pi*y)
    IMPLICIT NONE
    INTEGER m, mp, k, p, j
    REAL(real8), DIMENSION(0:m+1,0:mp+1) :: v
    REAL(real8), DIMENSION(:,:), POINTER :: c
    REAL(real8), DIMENSION(0:m+1) :: y0

    y0 = sin(3.141593*(/(j,j=0,m+1)/)/(m+1))

    IF(p > 1) THEN
      u = 0.0d0
      IF (k == 0  ) v(:,   0) = y0
      IF (k == p-1) v(:,mp+1) = y0*exp(-3.141593)
    ELSE
      v = 0.0d0
      v(:,0) = y0
      v(:,m+1) = y0*exp(-3.141593)
    END IF
    RETURN
  END SUBROUTINE bc

  SUBROUTINE borders(k, below, above, p)
    IMPLICIT NONE
    INTEGER :: k, below, above, p

    IF(k == 0) THEN
      below  = -1        ! tells MPI not to perform send/recv
      above = k+1
    ELSE IF(k == p-1) THEN
      below  = k-1
      above = -1        ! tells MPI not to perform send/recv
    ELSE
      below  = k-1
      above = k+1
    ENDIF

    RETURN
  END SUBROUTINE borders

  SUBROUTINE update_bc_2( v, m, mp, k, below, above )
    INCLUDE "mpif.h"
    INTEGER :: m, mp, k, below, above, ierr
    REAL(real8), dimension(0:m+1,0:mp+1) :: v
    INTEGER status(MPI_STATUS_SIZE)
 
    CALL MPI_SENDRECV(                                       &
              v(1,mp  ), m, MPI_DOUBLE_PRECISION, above, 0,  &
              v(1,   0), m, MPI_DOUBLE_PRECISION,  below, 0,  &
              MPI_COMM_WORLD, status, ierr )
    CALL MPI_SENDRECV(                                       &
              v(1,   1), m, MPI_DOUBLE_PRECISION,  below, 1,  &
              v(1,mp+1), m, MPI_DOUBLE_PRECISION, above, 1,  &
              MPI_COMM_WORLD, status, ierr )
    RETURN
  END SUBROUTINE update_bc_2

  SUBROUTINE update_bc_1(v, m, mp, k, below, above)
    IMPLICIT NONE
    INCLUDE 'mpif.h'
    INTEGER :: m, mp, k, ierr, below, above
    REAL(real8), DIMENSION(0:m+1,0:mp+1) :: v
    INTEGER status(MPI_STATUS_SIZE)

! Select 2nd index for domain decomposition to have stride 1
! Use odd/even scheme to reduce contention in message passing
    IF(mod(k,2) == 0) THEN     ! even numbered processes
      CALL MPI_Send( v(1,mp  ), m, MPI_DOUBLE_PRECISION, above, 0,  &
                     MPI_COMM_WORLD, ierr)
      CALL MPI_Recv( v(1,0   ), m, MPI_DOUBLE_PRECISION,  below, 0,  &
                     MPI_COMM_WORLD, status, ierr)
      CALL MPI_Send( v(1,1   ), m, MPI_DOUBLE_PRECISION,  below, 1,  &
                     MPI_COMM_WORLD, ierr)
      CALL MPI_Recv( v(1,mp+1), m, MPI_DOUBLE_PRECISION, above, 1,  &
                     MPI_COMM_WORLD, status, ierr)
    ELSE                         ! odd numbered processes
      CALL MPI_Recv( v(1,0   ), m, MPI_DOUBLE_PRECISION,  below, 0,  &
                     MPI_COMM_WORLD, status, ierr)
      CALL MPI_Send( v(1,mp  ), m, MPI_DOUBLE_PRECISION, above, 0,  &
                     MPI_COMM_WORLD, ierr)
      CALL MPI_Recv( v(1,mp+1), m, MPI_DOUBLE_PRECISION, above, 1,  &
                     MPI_COMM_WORLD, status, ierr)
      CALL MPI_Send( v(1,1   ), m, MPI_DOUBLE_PRECISION,  below, 1,  &
                     MPI_COMM_WORLD, ierr)
    ENDIF
    RETURN
  END SUBROUTINE update_bc_1

  SUBROUTINE print_mesh(v,m,mp,k,iter)
  IMPLICIT NONE
  INTEGER :: m, mp, k, iter, i, out
  REAL(real8), DIMENSION(0:m+1,0:mp+1) :: v

    out = 20 + k
    do i=0,m+1
      write(out,"(2i3,i5,' => ',4f10.3)")k,i,iter,v(i,:)
    enddo
    write(out,*)'+++++++++++++++++++++++++++++++++++++++++++++++++++++++'

    RETURN
  END SUBROUTINE print_mesh
END MODULE jacobi_module