C
C The Fortran code below solves Laplace's equation on the white part of the
C following rectangular domain.
C
C The solution on the outer boundary is 0, the the solution on the inner
C boundary (containing the grey box) is 1. See the image grid.gif given
C with this program ...
C 
C Author: David Walker
C---------------------------------------------------------------

      program laplace
      include 'mpif.h'
      parameter (nxmax=50, nymax=50)
      integer rank, ibuf(4), lsrc, ldest, rsrc, rdest
      integer dsrc, ddest, usrc, udest
      integer contig, strided, coords(2), src, comm2d
      integer rstatus(mpi_status_size)
      real phi(0:nymax+1,0:nxmax+1)
      real oldphi(0:nymax+1,0:nxmax+1)
      logical mask(0:nymax+1,0:nxmax+1)
      call mpi_init (ierr)
      call mpi_comm_rank (mpi_comm_world, rank, ierr)
      if (rank .eq. 0) then
         ibuf(1) = 2
         ibuf(2) = 2
         ibuf(3) = 16
         ibuf(4) = 16
      end if
      call mpi_bcast (ibuf, 4, mpi_integer, 0, mpi_comm_world, ierr)
      npx    = ibuf(1)
      npy    = ibuf(2)
      nptsx  = ibuf(3)
      nptsy  = ibuf(4)
      call setup_comm (rank, npx, npy, ipx, ipy,
     #                 lsrc, ldest, rsrc, rdest,
     #                 dsrc, ddest, usrc, udest, comm2d)
      call setup_grid (phi, npx, npy, nptsx, nptsy, ipx, ipy, mask)
      call mpi_type_contiguous (nptsy, mpi_real, contig, ierr)
      call mpi_type_vector (nptsx, 1, nymax+2, mpi_real, strided, ierr)
      call mpi_type_commit (contig, ierr)
      call mpi_type_commit (strided, ierr)
      do iter=1,100
         do j=0,nptsx+1
            do i=0,nptsy+1
               oldphi(i,j) = phi(i,j)
            end do
         end do
         call mpi_sendrecv (oldphi(1,1), 1, contig, ldest, 30,
     #                      oldphi(1,nptsx+1), 1, contig, lsrc, 30,
     #                      comm2d, rstatus, ierr)
         call mpi_sendrecv (oldphi(1,nptsx), 1, contig, rdest, 31,
     #                      oldphi(1,0), 1, contig, rsrc, 31,
     #                      comm2d, rstatus, ierr)
         call mpi_sendrecv (oldphi(1,1), 1, strided, udest, 40,
     #                      oldphi(nptsy+1,1), 1, strided, usrc, 40,
     #                      comm2d, rstatus, ierr)
         call mpi_sendrecv (oldphi(nptsy,1), 1, strided, ddest, 41,
     #                      oldphi(0,1), 1, strided, dsrc, 41,
     #                      comm2d, rstatus, ierr)
         do j=1,nptsx
            do i=1,nptsy
               if (mask(i,j)) phi(i,j) = 0.25*(oldphi(i-1,j) +
     #            oldphi(i+1,j) + oldphi(i,j-1) + oldphi(i,j+1))
            end do
         end do
      end do
      call mpi_type_free (contig, ierr)
      call mpi_type_free (strided, ierr)
      if (rank .eq. 0) open (7, file='jacobi.output')
      do m=0,npx-1
         do j=1,nptsx
            do n=0,npy-1
               if (ipx .eq. m .and. ipy .eq. n) then
                  if (rank .eq. 0) then
                     write(6,600)(phi(i,j),i=1,nptsy)
 600                 format(5e15.6)
                  else
                     call mpi_send (phi(1,j), nptsy, mpi_real,
     #                              0, 50, mpi_comm_world, ierr)
                  end if
               elseif (rank .eq. 0) then
                  coords(1) = n
                  coords(2) = m
                  call mpi_cart_rank (comm2d, coords, src, ierr)
                  call mpi_recv (oldphi(1,1), nptsy, mpi_real, src,
     #                           50, mpi_comm_world, rstatus, ierr)
                  write (6,500)(oldphi(i,1),i=1,nptsy)
 500              format(5e15.6)
               end if
             end do
          end do
      end do
      call mpi_finalize (ierr)
      stop
      end

      subroutine setup_comm (rank, npx, npy, ipx, ipy,
     #                       lsrc, ldest, rsrc, rdest,
     #                       dsrc, ddest, usrc, udest, comm2d)
      include 'mpif.h'
      integer rank, comm2d
      integer lsrc, ldest, rsrc, rdest, dsrc, ddest, usrc, udest
C      integer contig, strided
      integer coords(2), dims(2)
      logical reorder, periods(2)
      periods(1) = .false.
      periods(2) = .false.
      reorder    = .false.
      dims(1)    = npy
      dims(2)    = npx
      call mpi_cart_create (mpi_comm_world, 2, dims, periods, reorder,
     #                      comm2d, ierr)
      call mpi_cart_coords (comm2d, rank, 2, coords, ierr)
      ipy = coords(1)
      ipx = coords(2)
      call mpi_cart_shift (comm2d, 0,  1, dsrc, ddest, ierr)
      call mpi_cart_shift (comm2d, 0, -1, usrc, udest, ierr)
      call mpi_cart_shift (comm2d, 1,  1, rsrc, rdest, ierr)
      call mpi_cart_shift (comm2d, 1, -1, lsrc, ldest, ierr)
      return
      end

      subroutine setup_grid (phi, npx, npy, nptsx, nptsy,
     #                       ipx, ipy, mask)
      implicit none
      integer nxmax, nymax
      parameter (nxmax=50, nymax=50)
      real phi(0:nymax+1,0:nxmax+1)
      logical mask(0:nymax+1,0:nxmax+1)
      integer npx, npy, nptsx, nptsy, ipx, ipy
      integer i, j, ntotx, ntoty, ntotx4, ntoty4
      integer nglobx, ngloby

      do j=0,nptsx+1
         do i=0,nptsy+1
            phi(i,j) = 0.0
            mask(i,j)   = .true.
         end do
      end do

      if (ipx .eq. 0) then
         do i=1,nptsy
            mask(i,1)     = .false.
         end do
      end if

      if (ipx .eq. npx-1) then
         do i=1,nptsy
            mask(i,nptsx) = .false.
         end do
      end if

      if (ipy .eq. 0) then
         do j=1,nptsx
            mask(1,j) = .false.
         end do
      end if

      if (ipy .eq. npy-1) then
         do j=1,nptsx
            mask(nptsy,j) = .false.
         end do
      end if

      ntotx  = npx*nptsx
      ntoty  = npy*nptsy
      ntotx4 = ntotx/4
      ntoty4 = ntoty/4
      do j=1,nptsx
         nglobx = ipx*nptsx + j
         if (nglobx .gt. ntotx4 .and. nglobx .le. 3*ntotx4) then
            do i=1,nptsy
               ngloby = ipy*nptsy + i
               if (ngloby .gt. ntoty4 .and. ngloby .le. 3*ntoty4) then
                  phi(i,j) = 1.0
                  mask(i,j) = .false.
               end if
            end do
         end if
      end do

      return
      end
