This page imported from: /afs/bu.edu/cwis/webuser/web/s/c/scv/documentation/tutorials/MPI/alliance/apply/solvers/jacobi_module.html

MODULE jacobi_module

MODULE jacobi_module


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) THEN
v = 0.0d0
IF (k == 0 ) u(:, 0) = y0
IF (k == p-1) u(:,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