C ---------------------------------------------------------------------------
C HEAT2D Example - Parallelized Fortran Version
C FILE: pvm_heat2D.f
C OTHER FILES: make.pvm_heat2D.f  draw_heat.c pvm_heat2D.h
C DESCRIPTIONS:  This example is based on a simplified two-dimensional heat
C   equation domain decomposition.  The initial temperature is computed to be
C   high in the middle of the domain and zero at the boundaries.  The
C   boundaries are held at zero throughout the simulation.  During the
C   time-stepping, an array containing two domains is used; these domains
C   alternate between old data and new data.
C 
C   In this parallelized version, the grid is decomposed by the parent
C   process and then distributed by cols to children processes.  At each
C   time step, children processes must exchange border data with neighbors,
C   because a grid point's current temperature depends upon it's previous
C   time step value plus the values of the neighboring grid points.  Upon
C   completion of all time steps, the child processes return their results
C   to the parent process.
C 
C   Two data files are produced: an initial data set and a final data set.
C   An X graphic of these two states displays after all calculations have
C   completed.
C 
C AUTHOR: Blaise Barney - adapted from D. Turner's serial C version
C REVISED: 7/24/94 Blaise Barney - created
C ---------------------------------------------------------------------------
C Explanation of constants and variables
C   NXPROB                         =  x dimension of problem grid 
C   NYPROB                         =  y dimension of problem grid
C   STEPS                          =  number of time steps
C   AOUT                           =  name of PVM executable 
C   DONTCARE                       =  accept message from anytask
C   MAXCHILD                       =  maximum number of children tasks
C   MINCHILD                       =  minimum number of children tasks
C   BEGIN, NGHBOR1, NGHBOR2, DONE  =  message types 
C   NONE                           =  indicates no neighbor
C   CX, CY                         =  used in heat equation
C   u                              =  array for grids
C   taskid,parent,tids()           =  PVM taskids 
C   nproc,narch                    =  number of PVM machines & architect
C   nchild                         =  number of children processes
C   avecol,cols,offset,extra       =  for sending cols of data
C   neighbor1,neighbor2            =  neighbor tasks
C   rc,start,end                   =  misc
C   dtid,name,arch,speed           =  misc PVM
C   i,ix,iy,iz,it                  =  loop variables
C ---------------------------------------------------------------------------

      program heat2D
C     PVM version 3 include file
      (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
      include 'pvm_heat2D.h'

C     Routine for creating the X graph of the wave
      external draw_heat

      integer STEPS,DONTCARE,MAXCHILD,MINCHILD,BEGIN,
     &        NGHBOR1,NGHBOR2,DONE,NONE
      parameter(STEPS=50)
      parameter(DONTCARE=-1)
      parameter(MAXCHILD=8)
      parameter(MINCHILD=3)
      parameter(BEGIN=1)
      parameter(NGHBOR1=2)
      parameter(NGHBOR2=3)
      parameter(DONE=4)
      parameter(NONE=0)
      integer taskid,parent,tids(MAXCHILD),nproc,narch,nchild,avecol,
     &        cols,offset,extra,neighbor1,neighbor2,rc,start,end,
     &        i,ix,iy,iz,it,dtid,speed,utemp
      character*6 AOUT/'heat2D'/
      character*20 name,arch

C     Enroll process and find parent task id
      (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
      (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)

      if (parent .eq. PvmNoParent) then
C     ****************************** parent code *****************************
C     Determine the number of machines in this configuration. Then spawn one
C     child process for each processor excluding the parent.  If number of
C     children exceed MAXCHILD, then just spawn the maximum. If less than
C     MINCHILD then just spawn MINCHILD child processes. 

      (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
      nchild=nproc-1
      if (nchild .gt. MAXCHILD) then
        nchild=MAXCHILD
      end if
      if (nchild .lt. MINCHILD) then
        nchild=MINCHILD
      end if
      print *,'Spawning',nchild,'children processes'
      do i=1,nchild
        (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
        print *,'...child task id=',tids(i)
      end do

C     Initialize grid 
      print *,'Grid size: X=',NXPROB,' Y=',NYPROB,' Time steps=',STEPS
      print *,'Initializing grid and writing initial.dat file...'
      call inidat
      call prtdat(1)

C     Distribute work to children.  Must first figure out how many cols to
C     send and what to do with extra cols.  
      avecol=NYPROB/nchild
      extra=mod(NYPROB,nchild)
      offset=1
      do i=1, nchild
        if (i .le. extra) then
          cols=avecol+1
        else
          cols=avecol 
        end if
C       Tell each child which other children are its neighbors, since
C       they must exchange data with each other later. 
        if (i .eq. 1) then
          neighbor1=NONE
        else  
          neighbor1=tids(i-1)
        end if
        if (i .eq. nchild) then
          neighbor2=NONE
        else
          neighbor2=tids(i+1)
        end if 
C       Now send startup information to each child
        (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
        (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
        (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
        (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
        (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
        (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
        (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
        print *,'Sent to=',tids(i),'offset=',offset,'cols=',cols,
     .    'neighbor1=',neighbor1,'neighbor2=',neighbor2
        offset = offset + cols
      end do
          
C     Now wait for results from all children tasks 
      do i=1, nchild
        (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
        (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
        (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
        (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
      end do
 
C     Exit PVM
      (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
 
C     Print and show results 
      print *,'Creating final.dat file and generating graph...'
      call prtdat(2)
      call draw_heat(%VAL(NXPROB),%VAL(NYPROB))
 
C     End of parent code
      end if
 
 
      if (parent .ne. PvmNoParent)  then
C     ****************************** child code ******************************
C     Initialize everything - including the borders - to zero 
      do ix=1,NXPROB
        do iy=1,NYPROB
          do iz=1,2
            u(ix,iy,iz) = 0.0
          end do
        end do
      end do

C     Now receive my offset, cols, neighbors and grid partition from parent 
      (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
      (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
      (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
      (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
      (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
      (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
 
C     Determine border elements.  Need to consider first and last columns.
C     Obviously, col 1 can't exchange with col 1-1.  Likewise, the last
C     col can't exchange with last+1. 
      if (offset .eq. 1) then
        start=2
      else
        start=offset
      end if
      if (offset + cols .gt. NYPROB) then
        end=NYPROB-1
      else
        end=offset+cols-1
      end if

C     Begin doing STEPS iterations.  Must communicate border cols with
C     neighbors.  If I have the first or last grid col, then I only need to 
C     communicate with one neighbor. 
      iz=1
      do it=1, STEPS
        if (neighbor1 .ne. NONE) then
          (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
          (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
          (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
 
          (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
          (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
        end if
 
        if (neighbor2 .ne. NONE) then
          (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
          (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
          (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
 
          (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
          (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
        end if
      
C       Now call update to update the value of grid points
        call update(start,end,u(1,1,iz),u(1,1,3-iz))
        
        iz=3-iz
      end do
    
C     Send my portion of final results back to parent 
      (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
      (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
      (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
      (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)
      (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)

C     Exit PVM  
      (XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX)

C     End of child code
      end if

      end


C****************************************************************************
      subroutine update (start, end, u1, u2)
C****************************************************************************
      include 'pvm_heat2D.h'
      integer start, end, ix, iy
      real*4 u1, u2
      dimension u1(NXPROB,NYPROB),u2(NXPROB,NYPROB)

      do iy=start, end
        do ix=2, NXPROB-1
            u2(ix,iy) = u1(ix,iy)
     .       + CX * ( u1(ix+1,iy) + u1(ix-1,iy) - 2.0 * u1(ix,iy))
     .       + CY * ( u1(ix,iy+1) + u1(ix,iy-1) - 2.0 * u1(ix,iy))
        end do
      end do
	  end

C*****************************************************************************
      subroutine inidat
C*****************************************************************************
      include 'pvm_heat2D.h'
      integer ix,iy

      do ix=0,NXPROB-1
        do iy=0,NYPROB-1
          u(ix+1,iy+1,1) = float(ix*(NXPROB-ix-1) * iy*(NYPROB-iy-1))
        end do
      end do
      end


C**************************************************************************
      subroutine prtdat(i)
C**************************************************************************
      include 'pvm_heat2D.h'
      integer i,ix, iy
      character*11 fname

      if (i .eq. 1) then
        fname ='initial.dat'
      else if (i .eq. 2) then
        fname = 'final.dat'
      end if

      open(21, file=fname, form='formatted')
      do ix=1,NXPROB
        do iy=1,NYPROB
          write(21,'(f6.1,1x,$)')u(ix,iy,1)
        end do
        write(21,'(1x)')
      end do
      close(21)
      end