c -----------------------------------------------------------
c Matrix\times vector routine:  y = A*x
c -----------------------------------------------------------
      subroutine matmul(y, m,n,s,e,w, x, wrk, buf_out,buf_in,
     >                  me,north,south,east,west,lx,ly,bufsize)
      implicit none
      include 'mpif.h'
      integer i,j,lx,ly,me,north,south,east,west,bufsize
      integer ierr,status
      double precision m(lx,ly), 
     >                 n(lx,ly), s(lx,ly), 
     >                 e(lx,ly), w(lx,ly), 
     >                 x(lx,ly), y(lx,ly),
     >                 wrk(0:lx+1,0:ly+1),
     >                 buf_out(bufsize),buf_in(bufsize)

       call vimbed(wrk,x,lx,ly)
c ----------------------------------------------------------------
c ------------- send to north, receive from south, barrier
       do i=1,lx
          buf_out(i)=x(i,ly)
       enddo
       if (north.ne.-1) then
          call MPI_SEND( buf_out,lx,MPI_DOUBLE_PRECISION,north,
     >                   me,MPI_COMM_WORLD,ierr)
       endif
       if (south.ne.-1) then
          call MPI_RECV( buf_in,lx,MPI_DOUBLE_PRECISION,south,
     >                   south,MPI_COMM_WORLD,status,ierr)
       endif
       call MPI_BARRIER(MPI_COMM_WORLD,ierr)
       if (south.ne.-1) then
          do i=1,lx
             wrk(i,0)=buf_in(i)
          enddo
       endif
c ----------------------------------------------------------------
c ------------- send to south, receive from north, barrier
       do i=1,lx
          buf_out(i)=x(i,1)
       enddo
       if (south.ne.-1) then
          call MPI_SEND( buf_out,lx,MPI_DOUBLE_PRECISION,south,
     >                   me,MPI_COMM_WORLD,ierr)
       endif
       if (north.ne.-1) then
          call MPI_RECV( buf_in,lx,MPI_DOUBLE_PRECISION,north,
     >                   north,MPI_COMM_WORLD,status,ierr)
       endif
       call MPI_BARRIER(MPI_COMM_WORLD,ierr)
       if (north.ne.-1) then
          do i=1,lx
             wrk(i,ly+1)=buf_in(i)
          enddo
       endif
c ----------------------------------------------------------------
c ------------- send to east, receive from west, barrier
       do j=1,ly
          buf_out(j)=x(lx,j)
       enddo
       if (east.ne.-1) then
          call MPI_SEND( buf_out,ly,MPI_DOUBLE_PRECISION,east,
     >                   me,MPI_COMM_WORLD,ierr)
       endif
       if (west.ne.-1) then
          call MPI_RECV( buf_in,ly,MPI_DOUBLE_PRECISION,west,
     >                   west,MPI_COMM_WORLD,status,ierr)
       endif
       call MPI_BARRIER(MPI_COMM_WORLD,ierr)
       if (west.ne.-1) then
          do j=1,ly
             wrk(0,j)=buf_in(j)
          enddo
       endif
c ----------------------------------------------------------------
c ------------- send to west, receive from east, barrier
       do j=1,ly
          buf_out(j)=x(1,j)
       enddo
       if (west.ne.-1) then
          call MPI_SEND( buf_out,ly,MPI_DOUBLE_PRECISION,west,
     >                   me,MPI_COMM_WORLD,ierr)
       endif
       if (east.ne.-1) then
          call MPI_RECV( buf_in,ly,MPI_DOUBLE_PRECISION,east,
     >                   east,MPI_COMM_WORLD,status,ierr)
       endif
       call MPI_BARRIER(MPI_COMM_WORLD,ierr)
       if (east.ne.-1) then
          do j=1,ly
             wrk(lx+1,j)=buf_in(j)
          enddo
       endif
c ----------------------------------------------------------------
       do i=1,lx
          do j=1,ly
             y(i,j) = m(i,j)*wrk(i,j)   + 
     >                e(i,j)*wrk(i+1,j) +
     >                w(i,j)*wrk(i-1,j) +
     >                n(i,j)*wrk(i,j+1) +
     >                s(i,j)*wrk(i,j-1)
          enddo
       enddo

      return
      end

c -----------------------------------------------------------
c Matrix generator routine:  
c 5-point stencil, central differences, homogeneous Dirichlet b.c.
c -----------------------------------------------------------
      subroutine matgen(m,n,s,e,w, 
     >           lx,ly,me,north,south,east,west)
      implicit none
      include 'mpif.h'
      integer i,j,lx,ly,north,south,east,west,me
      double precision m(lx,ly), 
     >                 n(lx,ly), s(lx,ly), 
     >                 e(lx,ly), w(lx,ly)

       do i=1,lx
         do j=1,ly
            m(i,j) =  4.0d0
            n(i,j) = -1.0d0
            s(i,j) = -1.0d0
            e(i,j) = -1.0d0
            w(i,j) = -1.0d0
         enddo
      enddo

c -------- find the global position of the subregion
      if (west .eq. -1) then
        do j=1,ly
            w(1,j) = 0.0d0
        enddo
      endif

      if (south .eq. -1) then
        do i=1,lx
            s(i,1) = 0.0d0
        enddo
      endif

      if (east .eq. -1) then
        do j=1,ly
            e(lx,j) = 0.0d0
        enddo
      endif

      if (north .eq. -1) then
        do i=1,lx
            n(i,ly) = 0.0d0
        enddo
      endif

      return
      end
c -----------------------------------------------------------
c Right-hand-side routine:
c -----------------------------------------------------------
      subroutine rhsgen(rhs,u, m,n,s,e,w,
     >                  wrk, buf_out, buf_in,
     >                  me,north,south,east,west,
     >                  lx, ly, bufsize)

      implicit none
      integer i,j,lx,ly,bufsize
      integer north,south,east,west,me
      double precision rhs(lx,ly),u(lx,ly)
      real rand
      double precision m(lx,ly),
     >                 n(lx,ly), s(lx,ly),
     >                 e(lx,ly), w(lx,ly),
     >                 wrk(0:lx+1,0:ly+1),
     >                 buf_out(bufsize),buf_in(bufsize)
      do i=1,lx
         do j=1,ly
            u(i,j) = 1.0d0 !rand(0)
         enddo
      enddo
      call matmul(rhs, m,n,s,e,w, u, wrk,buf_out,buf_in,
     >            me,north,south,east,west,lx,ly,bufsize)

      return
      end