*********************************************************************
*   Jacobi Method by Vertical Slices for 2Dim Laplace's Equation    *
*********************************************************************
* 
* Laplace Eqn 
* T is initially 0.0 everywhere except at boundaries where T=100.
* Number of Jacobi iterations is fixed by NITER.
* Ghost boundary exchange depends on Column-Wise Fortran Environment
*
*        T=400.0
*       _________               ____ ____ ____ ____    Y,I
*       |*******|              |    |    |    |    |   |
*       |Initial|              |    |    |    |    |   |
* T=200.|T=100.0| T=300.0      | 0  | 1  | 2  | 3  |   |
*       |*Start*|              |    |    |    |    |   |
*       |*******|              |    |    |    |    |   |_______X,J
*       |_______|              |____|____|____|____|
*        T=100.0                         
* Use Central Differencing Method and 4 Slice Domain Decomposition
* Each process only has subgrid.
* Each Processor works on a sub grid and then sends its
* Boundaries to neighbours
*FH  Submit  MCS572 class nqs script "qsub pgm.job" after "cp laplace.f pgm.f" 
*FH  Converted from SPMD to Master/slave after SPMD_EX example
*
* Susheel Chitre PSC .. Mar 1994; Revised FBH/UIC Nov 1995,Nov 1996
***********************************************************************
*
      Program Laplace
      Parameter(Nrows=1002,Ncols=1002)
      Parameter(Nrint=Nrows-2,Ncint=Ncols-2)
      Parameter(Nwrows=Nrows,Nwcols=Ncint/4)
      Parameter(NRC=Nrint*Nwcols)  ! Number Points per Slice
      Parameter(Maxproc=8)
      Parameter(Niter=10)
      Real T(Nrows,0:Nwcols+1),Told(Nrows,0:Nwcols+1)
      Real D(Nrint,Nwcols)   ! T interation change: T-Told
      Integer Msgtype,Info,MyTID,MyPE,MsgPE,Lnbor,Rnbor,Nnbors
c  Note changed required include from old original.
      include '/usr/include/mpp/fpvm3.h'
      character*12 nodename
*
c  Enroll into PVM
*
c --- Obtain initial PVM data.
      call pvmfmytid(MyTID)
      call pvmfgetpe(MyTID,MyPE)
      call pvmfgsize(PVMALL,Nproc)
*
      write(6,*) 'MyTID,MyPE=',MyTID,MyPE
*
c  Master Spawns Nproc-1 slaves
*
            nodename = "Master"
*
*  Number of Neighbors Defined
      If(Nproc .gt. 1) then
         Nnbors=2
         If(MyPE .eq. 0 .or. MyPE .eq. nproc-1)then
            Nnbors=1
         Endif
      Endif
*
c  Initial Starting Values Set:
      Do jj = 0,Nwcols+1
         Do ii = 1,Nrows
            Told(ii,jj) = 100.0
         Enddo
      Enddo
*
c  Left Boundary Condition Set:
      If(MyPE .eq. 0) then
         Do ii=1,Nrows
            Told(ii,0)=200.0
         Enddo
      Endif
*
c  Right Boundary Condtion Set:
      If(MyPE .eq. Nproc-1) then
         Do ii=1,Nrows
            Told(ii,Nwcols+1)=300.0
         Enddo
      Endif
*
*  Top and Bottom Boundary Conditions Set:
      Do jj=0,Nwcols+1
         Told(1,jj) = 100.0
         Told(Nrows,jj) = 400.0
      Enddo
*
c  Do Computation on Sub-grid for Niter Iiterations
      write(6,*) 'Starting Iteration Big Loop 100'
*
      Do 100 Iter=1,Niter
*
         Do j=1,Nwcols
            Do i=2,Nrows-1
Comment:  Take average of 4 neighbors:
               T(i,j) = 0.25 * (Told(i+1,j)+Told(i-1,j)
     &                  + Told(i,j+1)+Told(i,j-1))
            Enddo
         Enddo
*
c  Save Current T into Told, but get difference first:
         Do j=1,Nwcols
            Do i=2,Nrows-1
               D(i-1,j) = abs(T(i,j)-Told(i,j))
               Told(i,j) = T(i,j)
            Enddo
         Enddo
*
c  Broadcast values of T to other Processors:
         If(Nproc  .gt. 1) then
c  Send to Neighbors Lnbor on Left Only:
            If(MyPE .ne. 0 ) then
               msgtype=Iter
               call pvmfinitsend(pvmraW,info)
               call pvmfgettid(PVMALL,MyPE-1,Lnbor)
               call pvmfpack(Integer4,MyPE,1,1,info)
               call pvmfpack(Real8,T(1,1),Nrows,1,info)
               call pvmfsend(Lnbor,msgtype,info)
            Endif
c  Send to Neighbors Rnbor on Right Only:
            If(MyPE .ne. nproc-1) then
               msgtype=Iter
               call pvmfinitsend(pvmraw,info)
               call pvmfgettid(PVMALL,MyPE+1,Rnbor)
               call pvmfpack(Integer4,MyPE,1,1,info)
               call pvmfpack(Real8,T(1,Nwcols),Nrows,1,info)
               call pvmfsend(Rnbor,msgtype,info)
            Endif
*
c  Receive T from Neigbours
            Do jj = 1,Nnbors
               msgtype=Iter
               call pvmfrecv(-1,msgtype,info)
               call pvmfunpack(Integer4,MsgPE,1,1,info)
               If(MsgPE .lt. MyPE)then
C  Unpack T Column as Right Neighbor: MsgPE:T(1,Nwcols) -> MyPE:Told(1,0)
                  call pvmfunpack(Real8,Told(1,0),Nrows,1,info)
               Else
C  Unpack T Column as Left Neighbor: MsgPE:T(1,1) -> MyPE:Told(Nwcols+1,0)
                  call pvmfunpack(Real8,Told(1,Nwcols+1),Nrows,1,info)
               Endif
            Enddo
         Endif
*
c  Print some Values (Caution: diagnostic only, else not wise):
*
         If(MyPE .eq. 0 .and. mod(Iter,10) .eq. 0) then
            write(6,66) Iter, MyPE, T(10,10),D(9,10)
         Endif
         If(MyPE .eq. nproc-1 .and. mod(Iter,10) .eq. 0) then
            write(6,66) Iter, MyPE, T(10,nwcols-9),D(9,nwcols-9)
            call flush (6)     ! write buffered output to file for unit 6
         Endif
66    Format(1x,'Iter,MyPE,Temp,DT= ',i3,i3,F10.3,E10.2)
*
c  Go Back for Next Jacobi Iteration
Caution:  Should use Convergence Stopping Criterion instead.
*
100   CONTINUE
*
c  Exit PVM
      call pvmfexit(info)
*
      STOP
      END