      program gauss
c
c   gauss.f 9.2 92/07/01 17:53:07
c
c   Example program that solves LaPlace Equation in two dimensions using
c   an iterative Gauss-Seidel method. It uses column-wise data decomposition
c   of a two dimensional array, and on each cycle of the iteration it
c   averages each array element with its 4 nearest neighbors. To obtain the
c   nearest neighbors from columns that reside on different nodes, the
c   "boundary columns" are exchanged on each cycle, using "shadow buffers"
c   to receive copies of the neighbor columns from the previous cycle.
c
c
      double precision a(34,34), resid, residtemp, aold,eps
      double precision residtest, amax
c
c   totdim is the dimension of the square global grid.  It must be
c   a multiple of the number of nodes (a power of 2)
c
c   ndim is the dimension of a.  It must be at least totdim+2
c   but may be more
c
c   range is the dimension of the data owned by a node
c
c   n is the dimension of the node's grid, which includes shadow
c   buffers
c
      integer          totdim, ndim, n, iam, nbrnodes, range
      integer          xmin, xmax, ymin, ymax, maxiter,startup

c
c  -- Initialize --
c
      ndim     = 34
      totdim   = 16
      iam      = mynode()
      nbrnodes = numnodes()
      range    = totdim/nbrnodes 
      n        = range + 2

      xmin     = 2
      xmax     = range + 1
      ymin     = 2
      ymax     = totdim + 1

      maxiter = 2000
      eps   = 1.0D-05
      resid = 0.0D0
      residtest = 0.0D0
      amax      = 0.0D0

      if(range .lt. 1) then
         if(iam .eq. 0) then
            print *,''
            print *,'Try again with ',totdim, ' or less nodes'
         endif
         call exit
      endif
      call inita(iam,nbrnodes,totdim,ndim,n,a)
c
c  --  Perform the Calculation --
c
      do 101 k = 1,maxiter
         call shadow(ndim,n,totdim, iam,nbrnodes,range,a)

         do 100 j = xmin, xmax 
            do 100 i = ymin, ymax 
               aold      = a(i,j)
               a(i,j)    = 0.25 * (a(i+1,j) + a(i-1,j) + 
     x                            a(i,j+1) + a(i,j-1)) 
               residtemp  = dabs(a(i,j) - aold) 
               if(residtemp .gt. resid) then
                  resid = residtemp
                  amax  = aold 
               endif
100      continue

         if(amax .eq. 0.0) then
            residtest = 1.0 + eps
         else
            residtest = dabs( (resid)/amax )
         endif

         call gdhigh(residtest,1,residtemp)
         if( residtest .lt. eps) then
            if(iam .eq. 0) print *,''
            write(*,1003)iam,k
1003        format(I3,3x, I5,3x, 'breaking out of loop')
            goto 102
         else 
            resid = 0.0D0
         endif
101   continue
102   continue

c
c   -- Print out the solution --
c
c   First write out the left boundary column and the data columns owned by
c   node 0.  The data columns include the top and bottom boundaries
c
      if(iam .eq. 0) then
         write(*,1001)
         do 300 j = 1,n-1 
            write(*,1000) j-1, (a(i,j), i=1,totdim+2)
300      continue
1000     format('C',I2, 20(1x,F3.0))
1001     format(' ')

c
c Then node 0 receives the rest of the columns up to totdim
c Node 0 puts received data into its first data column, which 
c is already printed
c

         do 400 j = range+1, totdim+1
            call crecv(j,a(1,2),8*(totdim+2))
            write(*,1000) j, (a(i,2), i=1,totdim+2)
400      continue

c Other nodes send local columns 2 through range+1.
c The message type is the global column number 

      else
         do 500 k = 1,range
            call csend(iam*range + k, a(1,k+1),8*(totdim+2),0,0)
500      continue

c If I am the rightmost node, then also send the right boundary

         if(iam .eq. nbrnodes-1) then
           call csend(iam*range +range+1,a(1,range+2),8*(totdim+2),0,0)
         endif
      endif
      end

C 
      subroutine inita(iam,nbrnodes,totdim,ndim,n,a)
      integer          iam,nbrnodes,totdim,ndim,n
      double precision a(ndim,n) 
      double precision topbound, botbound, leftbound, rightbound

      topbound    = 0.0
      botbound    = 0.0
      leftbound   = 99.0
      rightbound  = 99.0

c Zero out a

      do 100 i=1,totdim+2
         do 100 j=1,n
            a(i,j) = 0.0
100   continue
c
c  -- Put in boundaries --
c
      do 200 j=1,n
         a(1,j)           = topbound
         a(totdim+2,j)    = botbound
200   continue

      if(iam .eq. 0) then
         do 300 i=1,ndim
            a(i,1)    = leftbound
300      continue
      endif

      if(iam .eq. nbrnodes-1) then
         do 400 i=1,ndim
            a(i,n)    = rightbound
400      continue
      endif
      end

C 

      subroutine shadow(ndim,n,totdim, iam,nbrnodes,range,a)
c
c  Shadow performs the exchange of neighbor columns into the shadow buffers
c
      integer          fromleft, fromright
      parameter (fromleft = 200, fromright = 300)
      integer          ndim,n,totdim
      integer          iam, nbrnodes, range
      double precision a(ndim,range+2)

      integer          length,leftnode,rightnode,leftid,rightid

      leftnode  = iam - 1
      rightnode = iam + 1
      length    = (totdim+2)*8

      if(iam.eq.0) then
c
c   If I am the leftmost node of the array (node 0) then only exchange
c   with the right (to the left is a boundary of the array)
c
         rightid = irecv(fromright, a(1,range+2), length)
         call csend(fromleft, a(1,range+1), length, rightnode,0)
         call msgwait(rightid)

      else if (iam .eq. nbrnodes-1) then
c
c   If I am the rightmost node of the array (highest numbered node) then
c   only exchange with the node to the left.
c
         leftid = irecv(fromleft, a(1,1), length)
         call csend(fromright, a(1,2), length, leftnode,0)
         call msgwait(leftid)

      else
c
c   Otherwise I am a node in the middle, so exchange with nodes to either side.
c
         leftid  = irecv(fromleft, a(1,1), length)
         rightid = irecv(fromright,a(1,range+2), length)

         call csend(fromright, a(1,2), length, leftnode,0)
         call csend(fromleft, a(1,range+1), length, rightnode,0)

         call msgwait(leftid)
         call msgwait(rightid)

      endif
      end
