c     domain-dcmp/poisson/p-2d
c     2-D domain decomposition 
c
      include "mpif.h"
      parameter (NDIM=2)
      parameter (MAX_X=200,MAX_Y=200)
      integer my_rank,nproc,ierr,stride
      real*8 hx,hy,xl,yl,u1(MAX_X,MAX_Y),u2(MAX_X,MAX_Y)
      real*8 array(400),diff1,diffnorm,uu(400)
      integer comm2d,my_coord(2)
      integer l_nbr,r_nbr
      integer t_nbr,b_nbr,my_crank,my_coord(NDIM)
c     
      call init(nproc,my_rank)
      call nbr2d(comm2d,my_crank,my_coord,l_nbr,r_nbr,t_nbr,b_nbr)
c
c     define number of computational points for each process in x-dir
      nx=3
      ny=3
c
c     total number of points: grid pts for each process + 2 bndry pts
c
      ntotx=nx*nproc+2
      ntoty=ny*nproc+2
      xl=1.0
      yl=1.0
c     spacing hx and hy
      hx=xl/dfloat(ntotx-1)
      hy=yl/dfloat(ntoty-1)
      call init_con(my_coord,nx,MAX_X,ny,u1)
c     MPI_Type_vector:
c     arg1: number of block
c     arg2: number of elements in each block
c     arg3: number of elements between start of each block
c     arg4: type
c     arg5: name of new type
       call MPI_Type_vector(ny+2,1,MAX_X,
     *     MPI_DOUBLE_PRECISION,stride,ierr)
      call MPI_type_commit(stride,ierr)
      diffnorm=10.
      iter=0
      do while(dabs(diffnorm).gt.1.e-4)
         call comput(my_crank,nx,MAX_X,ny,hx,hy,u1,u2)
         call exchange(l_nbr,r_nbr,t_nbr,b_nbr,nx,MAX_X,ny,stride,u2)
         call comput(my_crank,nx,MAX_X,ny,hx,hy,u2,u1)
         call exchange(l_nbr,r_nbr,t_nbr,b_nbr,nx,MAX_X,ny,stride,u1)
         call diff(nx,MAX_X,ny,u1,u2,diff1)
         call MPI_Allreduce(diff1,diffnorm,1,MPI_DOUBLE_PRECISION,
     *        MPI_SUM,MPI_COMM_WORLD,ierr)
         iter=iter+2
         if(my_rank.eq.1) then
            write(*,*) 'diffnorm=',diffnorm,'   iter=',iter
         endif
      enddo
c
c     collocting u1 of each proc to global array
c
      iroot=0
      if(my_crank.eq.0) then
         do i=2,nx+1
            write(*,*) (u1(i,j),j=2,ny+1)
         enddo
      endif
      call MPI_Finalize(ierr)
      stop
      end

      subroutine init(nproc,my_rank)

      include "mpif.h"
c
c     initialize the MPI
c
      integer my_rank,nproc
      character*30 nodename
      integer nchar
      call MPI_init(ierr)
      call MPI_Comm_rank(MPI_COMM_WORLD,my_rank,ierr)
      call MPI_Comm_size(MPI_COMM_WORLD,nproc,ierr)
      call MPI_Get_processor_name(nodename,nchar,ierr)
      write(*,*) 'my_rank=',my_rank,' processor=',nodename
      return
      end
c
      subroutine nbr2d(comm2d,my_crank,my_coord,
     *     l_nbr,r_nbr,t_nbr,b_nbr)
      include "mpif.h"
      parameter (NDIM=2)
      integer idim(2),comm2d,my_crank,l_nbr,r_nbr,t_nbr,b_nbr
      integer sideways,updown,right,up,my_coord(NDIM)
      logical periods(2),reorder
      idim(1)=2
      idim(2)=2
      periods(1)=.false.
      periods(2)=.false.
      reorder=.true.
      sideways=1
      updown=0
      right=1
      up=1
      call MPI_Cart_Create(MPI_COMM_WORLD,NDIM,idim,
     *     periods,reorder,comm2d, ierr)
      call MPI_Comm_Rank(comm2d,my_crank,ierr)
      call MPI_Cart_Coords(comm2d,my_crank,NDIM,my_coord,ierr)
      call MPI_Cart_Shift(comm2d,sideways,right,l_nbr,r_nbr,ierr)
      call MPI_Cart_Shift(comm2d,updown,up,t_nbr,b_nbr,ierr)
      write(*,*) 'my_crank=',my_crank,
     *     '  coor=(',my_coord(1),',',my_coord(2),' )'
      write(*,*) 'my_crank=',my_crank,' l_nbr=',l_nbr,
     *     ' r_nbr=',r_nbr,
     * 't_nbr=',t_nbr,'  b_nbr=',b_nbr
      return
      end
c
c--------------------------------------------------------
c     Initial Guess * Boundary Condition
c--------------------------------------------------------
      subroutine init_con(my_coord,nx,MAX_X,ny,u1)
      real*8 u1(MAX_X,1)
      integer nx,ny,my_coord(2)
      do i=1,nx+2
         do j=1,ny+2
            u1(i,j)=0.0
         enddo
      enddo
c
c     bndry u(0,y)=u(x,1)=1, u(x,0)=u(1,y)=0
c
      if(my_coord(2).eq.0) then
         do j=2,ny+1
            u1(1,j)=1.0
         enddo
      endif
      if(my_coord(1).eq.1) then
         do i=2,nx+1
            u1(i,1)=1.0
         enddo
      endif
      return
      end
c
c
c--------------------------------------------------------------
c     Computing u for the next time step
c--------------------------------------------------------------
      subroutine comput(my_crank,nx,MAX_X,ny,hx,hy,u1,u2)
      real*8 u1(MAX_X,1),u2(MAX_X,1),hx,hy,h
      real*8 f(200,400)
      integer my_crank
      h=hx
      do i=1,nx+2
         do j=1,ny+2
            u2(i,j)=u1(i,j)
         enddo
      enddo
      do i=2,nx+1
         do j=2,ny+1
            f(i,j)=0.0
            u2(i,j)=0.25*(u1(i-1,j)+u1(i,j+1)+u1(i,j-1)+u1(i+1,j)-
     *           h*h*f(i,j))
         enddo
      enddo
      return
      end

c-----------------------------------------------------------
c     Exchangeing the first and the last u to the neighbors
c-----------------------------------------------------------
      subroutine exchange(l_nbr,r_nbr,t_nbr,b_nbr,nx,MAX_X,
     *     ny,stride,u)
      include "mpif.h"
      real*8 u(MAX_X,1)
      integer status(MPI_STATUS_SIZE),stride
      integer l_nbr,r_nbr,t_nbr,b_nbr,itag

      itag=10
c     send u(2,j) to left neighbor and receive u(n+2,j) from right neighbor
c
      call MPI_Sendrecv(u(2,1),1,stride,l_nbr,itag,
     *     u(nx+2,1),1,stride,r_nbr,itag,
     *     MPI_COMM_WORLD,status,ierr)
c     
c     sending the last u to my_rank+1 for 0< my_rank < nproc-1
c     
      itag=20
c     send u(n+1,j) to right neighbor and receive u(1,j) from left neighbor
c
      call MPI_Sendrecv(u(nx+1,1),1,stride,r_nbr,itag,
     *     u(1,1),1,stride,l_nbr,itag,
     *     MPI_COMM_WORLD,status,ierr)
c
c     send u(i,2) to bottom neighbor and receive u(i,n+2) from top neighbor
c
      itag=30
      call MPI_Sendrecv(u(1,2),nx+2,MPI_DOUBLE_PRECISION,b_nbr,
     *     itag,u(1,ny+2),nx+2,MPI_DOUBLE_PRECISION,t_nbr,itag,
     *     MPI_COMM_WORLD,status,ierr)
c
c     sending the last u to my_rank+1 for 0< my_rank < nproc-1
c
      itag=40
c     send u(i,n+1) to top neighbor and receive u(i,1) from bottom neighbor
c
      call MPI_Sendrecv(u(1,ny+1),nx+2,MPI_DOUBLE_PRECISION,t_nbr,
     *     itag,u(1,1),nx+2,MPI_DOUBLE_PRECISION,b_nbr,itag,
     *     MPI_COMM_WORLD,status,ierr)
      return
      end
c
      subroutine diff(nx,MAX_X,ny,u1,u2,diff1)
      real*8 u1(MAX_X,1),u2(MAX_X,1),diff1
      diff1=0.0d0
      do i=2,nx+1
         do j=2,ny+1
            diff1=diff1+(u1(i,j)-u2(i,j))**2
         enddo
      enddo
      return
      end