2D Laplace SOR Solver
c
c This program solves Laplace's equation on a two-dimensional rectangle
c using the successive over-relaxation method with red/black ordering.
c
c The number of particles per processor in each direction must be even.
c
c This is the MPI version.
c
c Author: David W. Walker
c Date: April 7, 1995
c
program sor
implicit none
include 'mpif.h'
integer myrank, ierr, nprocs, nprocx, nprocy, pseudohost
integer nptx, npty
integer MAXX, MAXY
parameter (pseudohost = 0)
parameter (MAXX = 299, MAXY=299)
real phi(0:MAXY+2,0:MAXX+2)
logical mask(0:MAXY+2,0:MAXX+2)
integer color(0:MAXY+2,0:MAXX+2)
integer i, j, k, ix, iy, ipx, ipy, jpx, jpy
integer RED, BLACK
parameter (RED = 0, BLACK = 1)
integer maxiter
parameter (maxiter = 100)
double precision timebegin, timeend
integer insize, pos, src, m, n, status(MPI_STATUS_SIZE)
parameter (insize = 1000)
integer input(insize)
real rnput(insize)
integer nfixed, indx(1000), indy(1000)
real omega, omega4, omegam, fixval(1000), buf(MAXY)
integer up, down, left, right
integer dims(2), coords(2), from(2)
logical reorder, periods(2)
integer comm2d
integer typex, typey
c
c Initialize MPI, find rank of each process, and the number of processes
c
call mpi_init (ierr)
call mpi_comm_rank (MPI_COMM_WORLD, myrank, ierr)
call mpi_comm_size (MPI_COMM_WORLD, nprocs, ierr)
c
c One process acts as the host and reads in the input
c
if (myrank .eq. pseudohost) call get_input (input, rnput)
c
c The input is broadcast to all processes
c
call mpi_bcast (input, 5, mpi_integer, pseudohost,
# MPI_COMM_WORLD, ierr)
nprocx = input(1)
nprocy = input(2)
nptx = input(3)
npty = input(4)
nfixed = input(5)
call mpi_bcast (input(6), 2*nfixed, mpi_integer, pseudohost,
# MPI_COMM_WORLD, ierr)
call mpi_bcast (rnput, 1+nfixed, mpi_real, pseudohost,
# MPI_COMM_WORLD, ierr)
omega = rnput(1)
do i=1,nfixed
indx(i) = input(4+2*i)
indy(i) = input(5+2*i)
fixval(i) = rnput(1+i)
end do
c
c Abort if incorrect.
c
if (nprocx .lt. 0 .or. nprocy .lt. 0) then
print *, 'Number of processes in each direction must be >0'
goto 999
else if (nprocx*nprocy .ne. nprocs) then
print *, 'nprocx*nprocy must = ',nprocs
goto 999
endif
if (nptx .gt. MAXX) then
print *, 'Number of points per process must be <= ',MAXX
goto 999
end if
if (npty .gt. MAXY) then
print *, 'Number of points per process must be <= ',MAXY
goto 999
end if
if (mod(nptx,2) .ne. 0 .or. mod(npty,2) .ne. 0) then
print *, 'Number of points per process must be even'
goto 999
end if
c
c Start timer
c
timebegin = mpi_wtime ()
c
c Set up 2D application topology
c
periods(1) = .false.
periods(2) = .false.
dims(1) = nprocx
dims(2) = nprocy
reorder = .false.
call mpi_cart_create (MPI_COMM_WORLD, 2, dims, periods, reorder,
# comm2d, ierr)
call mpi_cart_coords (comm2d, myrank, 2, coords, ierr)
call mpi_cart_shift (comm2d, 0, 1, left, right, ierr)
call mpi_cart_shift (comm2d, 1, 1, down, up, ierr)
c
c Set up the grid mask and grid arrays
c
call setup_grid (phi, fixval, color, nptx, npty, nfixed,
# indx, indy, coords, nprocx, nprocy,
# mask, MAXY+2)
c
c Set up general datatypes
c
call mpi_type_contiguous (npty/2, mpi_real, typex, ierr)
call mpi_type_vector (nptx/2, 1, MAXY+2, mpi_real, typey, ierr)
call mpi_type_commit (typex, ierr)
call mpi_type_commit (typey, ierr)
c
c Start main loop
c
omega4 = 0.25*omega
omegam = 1.0 - omega
do k=1,maxiter
c
c Exchange black boundary data
c
call mpi_sendrecv (phi(1,1), 1, typex, left, 14,
# phi(1,nptx+1), 1, typex, right, 14,
# comm2d, status, ierr)
call mpi_sendrecv (phi(npty/2+1,nptx), 1, typex, right, 15,
# phi(npty/2+1,0), 1, typex, left, 15,
# comm2d, status, ierr)
call mpi_sendrecv (phi(npty,nptx/2+1), 1, typey, up, 17,
# phi(0,nptx/2+1), 1, typey, down, 17,
# comm2d, status, ierr)
call mpi_sendrecv (phi(1,1), 1, typey, down, 16,
# phi(npty+1,1), 1, typey, up, 16,
# comm2d, status, ierr)
c
c Update red points
c
do j=1,nptx
do i=1,npty
if (mask(i,j) .and. color(i,j) .eq. RED) phi(i,j) =
# omega4*(phi(i,j-1)+phi(i-1,j)+phi(i,j+1)+phi(i+1,j))
# + omegam*phi(i,j)
end do
end do
c
c Exchange red boundary data
c
call mpi_sendrecv (phi(npty/2+1,1), 1, typex, left, 10,
# phi(npty/2+1,nptx+1), 1, typex, right, 10,
# comm2d, status, ierr)
call mpi_sendrecv (phi(1,nptx), 1, typex, right, 11,
# phi(1,0), 1, typex, left, 11,
# comm2d, status, ierr)
call mpi_sendrecv (phi(npty,1), 1, typey, up, 13,
# phi(0,1), 1, typey, down, 13,
# comm2d, status, ierr)
call mpi_sendrecv (phi(1,nptx/2+1), 1, typey, down, 12,
# phi(npty+1,nptx/2+1), 1, typey, up, 12,
# comm2d, status, ierr)
c
c Update black points
c
do j=1,nptx
do i=1,npty
if (mask(i,j) .and. color(i,j) .eq. BLACK) phi(i,j) =
# omega4*(phi(i,j-1)+phi(i-1,j)+phi(i,j+1)+phi(i+1,j))
# + omegam*phi(i,j)
end do
end do
end do
c
c End main loop
c
call mpi_type_free (typex, ierr)
call mpi_type_free (typey, ierr)
c
c Stop clock and write out timings
c
timeend = mpi_wtime ()
print *,'Node', myrank,' Elapsed time: ',
# timeend-timebegin,' seconds'
c
c Print out results
c
if (myrank .eq. pseudohost) open (7,file='sor.output')
do m=0,nprocx-1
do i=1,nptx
do n=0,nprocy-1
if (coords(1) .eq. m .and. coords(2) .eq. n) then
if (myrank .eq. pseudohost) then
write (7,100) (phi(j,i),j=1,npty)
else
call mpi_send (phi(1,i), npty, mpi_real,
# pseudohost, 20,
# MPI_COMM_WORLD, ierr)
end if
else if (myrank .eq. pseudohost) then
from(1) = m
from(2) = n
call mpi_cart_rank (comm2d, from, src, ierr)
call mpi_recv (buf, npty, mpi_real, src, 20,
# MPI_COMM_WORLD, ierr)
write (7,100) (buf(j),j=1,npty)
end if
end do
end do
end do
call mpi_comm_free (comm2d, ierr)
999 call mpi_finalize (ierr)
stop
100 format(5e15.6)
end
subroutine get_input (input, rnput)
implicit none
integer input(*)
real rnput(*)
integer i
open(4,file='sor.input')
c
c Read in number of processes in each direction of process mesh
c
read (4,*) input(1)
read (4,*) input(2)
c
c Read in number of points per process in each direction
c
read (4,*) input(3)
read (4,*) input(4)
c
c Enter the value of the SOR parameter omega
c
read (4,*) rnput(1)
c
c Input information on fixed interior points
c
read (4,*) input(5)
do i=1,input(5)
read (4,*) input(4+2*i)
read (4,*) input(5+2*i)
read (4,*) rnput(1+i)
end do
return
end
subroutine setup_grid (phi, fixval, color, nptx, npty, nfixed,
# indx, indy, coords, nprocx, nprocy,
# mask, ldim)
implicit none
integer ldim
real phi(0:ldim,0:*), fixval(*)
integer color(0:ldim,0:*), nptx, npty, nfixed
integer indx(*), indy(*), coords(2), nprocx, nprocy
logical mask(0:ldim,0:*)
integer i, j, ix, iy, ipx, ipy, jpx, jpy
integer RED, BLACK
parameter (RED = 0, BLACK = 1)
do j=1,nptx
do i=1,npty
mask(i,j) = .true.
phi(i,j) = 0.0
if ((j .le. nptx/2 .and. i .gt. npty/2) .or.
# (j .gt. nptx/2 .and. i .le. npty/2)) then
color(i,j) = RED
else
color(i,j) = BLACK
end if
end do
end do
do i=1,nfixed
ix = indx(i) - 1
iy = indy(i) - 1
ipx = ix/nptx
ipy = iy/npty
if (coords(1) .eq. ipx .and. coords(2) .eq. ipy) then
jpx = mod(ix, nptx) + 1
jpy = mod(iy, npty) + 1
mask(jpy,jpx) = .false.
phi(jpy,jpx) = fixval(i)
end if
end do
if (coords(1) .eq. 0) then
do i=1,npty
mask(i,1) = .false.
phi(i,1) = 0.0
end do
end if
if (coords(1) .eq. nprocx-1) then
do i=1,npty
mask(i,nptx) = .false.
phi(i,nptx) = 0.0
end do
end if
if (coords(2) .eq. 0) then
do i=1,nptx
mask(1,i) = .false.
phi(1,i) = 0.0
end do
end if
if (coords(2) .eq. nprocy-1) then
do i=1,nptx
mask(npty,i) = .false.
phi(npty,i) = 0.0
end do
end if
return
end
Click here to return to the MPI home page.
http://www.epm.ornl.gov/~walker/mpi/examples/sor.html
David W. Walker,
Oak Ridge National Laboratory
/ (walkerdw@ornl.gov)
Last Modified April 7, 1995