      program driver
#include "realdef.inc"
      include "mpif.h"
c
c     /***************************************************************/
c     /* driver routine for testing the communication modules needed */
c     /* for the solution of a 2d poisson problem in a 2nd order     */
c     /* finite-difference representation and using jacobi iteration */
c     /*                                                             */
c     /* the basic algorithm is as simple as follows:                */
c     /*                                                             */
c     /* (1) initialize the arrays & set physical b.c.'s             */
c     /* (2) loop while not_converged:                               */
c     /*     (a) iterate one step                                    */
c     /*     (b) if(converged): break                                */
c     /*     (c) exchange ghost-cell data                            */
c     /* (3) output & termination                                    */
c     /*                                                             */
c     /* compile:                                                    */
c     /* make MAINP                                                  */
c     /* run:                                                        */
c     /* mpirun -np x ./MAINP   [where 'x' is no. of procs]          */
c     /* clean:                                                      */
c     /* make clean                                                  */
c     /*-------------------------------------------------------------*/
c     /* CHANGELOG                                                   */
c     /* v0.0 (nov.03) .call exchange_ghostline1/2 several times     */
c     /*                with uninitialized data; perform timing.     */
c     /*                                                             */
c     /* v1.0 04.12.03 .initialize data; enable check of correct     */
c     /*                communication by output to file of arrays    */
c     /*                by each processor.                           */
c     /*                added "realdef.inc" file, defining the float */
c     /*                type to be used in the code. this also means */
c     /*                that we now have an additional argument to   */
c     /*                be passed to the communication routines,     */
c     /*                called 'MY_MPI_REAL'.                        */
c     /*                removed the second communication call. for   */
c     /*                checking different comm-routines, edit the   */
c     /*                driver source code.                          */
c     /* v1.1 10.12.03  corrected small bugs in 'assign' and 'check_ */
c     /*                array'.                                      */
c     /* v2.0 16.12.03  modified such that the "real" problem is     */
c     /*                solved, i.e. jacobi iteration of poisson.    */
c     /*                this affects: 'assign','check_array',        */
c     /*                'jacobi_step','lconverged','write2d_p',      */
c     /*                'do_wrtie2d'. results are written to file in */
c     /*                ascii-gnuplot format; timing/convergence is  */
c     /*                written to unit 55.                          */
c     /***************************************************************/
c
c     /* (global) problem size and number of processors */
      parameter(m=10,n=10,nproc=1)
c
c     /* max. number of iterations allowed & tolerance*/
      parameter(nitmx=100000,tol=1.e-5)
c
c     /* derived dimensions: the second index 'n' is decomposed */
      parameter(np=n/nproc+1)
c
c     /* arrays */
      dimension aloc(0:m+1,np+2),adebug(0:m+1,np+2)
      dimension jbeg(0:nproc-1),jend(0:nproc-1)
c
c     /* mpi related */
      integer istat(MPI_STATUS_SIZE)
      real*8 t0,t1,t2           !note: result from MPI_WTIME() is real*8
c
c     /* other stuff */
      character*80 filename
c-----|---------------------------------------------------------------|
c     /* initialize MPI & define the float data-type 'MY_MPI_REAL' */
      call MPI_INIT(ierr)
      t0=MPI_WTIME()
      call MPI_COMM_RANK(MPI_COMM_WORLD,myid,ierr)
      call MPI_COMM_SIZE(MPI_COMM_WORLD,numprocs,ierr)
#if double_precision
      MY_MPI_REAL=MPI_DOUBLE_PRECISION
#else
      MY_MPI_REAL=MPI_REAL
#endif
      if(numprocs.ne.nproc)then !thats not what we compiled for
         write(*,*)'difference in no of procs: ',numprocs,nproc
         goto 111
      endif
c
c     /* write a welcome message to screen */
      write(*,*)'hello. I am ',myid,' out of ',numprocs
c
c     /* compute the pointers into local arrays */
      call MPE_DECOMP1D_all(n,nproc,jbeg,jend,'_')
c
#ifdef DEBUG_MAIN
c     /* write decomposition to screen and check data lengths */
      write(*,*)myid,' jbeg,jend: ',(jbeg(i),jend(i),i=0,nproc-1)
      do ip=0,nproc-1
         if((jend(ip)-jbeg(ip)+1).gt.np)
     $        write(*,*)'problem with excessive data length ',
     $        ip,jbeg(ip),jend(ip)
      enddo
#endif /* DEBUG_MAIN */
c
#ifdef DEBUG_COMM
c     /* initialize the arrays with some data; set ghosts to zero */
      call assign(aloc,1,m,jbeg(myid),jend(myid),m,n,myid,nproc,'_',1)
      call write_array(aloc,1,m,jbeg(myid),jend(myid),myid,20+myid)
c
c     /* synchronize for timing */
      call MPI_BARRIER(MPI_COMM_WORLD,ierr)
      t1=MPI_WTIME()
      call exchange_ghostline1(aloc,m,jbeg(myid),jend(myid),
     $     myid,nproc,MPI_COMM_WORLD,MY_MPI_REAL)
c     
      write(*,*)myid,' after 1st comm'
      call MPI_BARRIER(MPI_COMM_WORLD,ierr)
      t2=MPI_WTIME()
      if(myid.eq.0)then
         write(*,*)'communication took: ',t2-t1,' sec'
      endif
c     /* check the content of the arrays after comm */
      call assign(adebug,1,m,jbeg(myid),jend(myid),m,n,myid,nproc,
     $     'g',1)
      call check_array(aloc,adebug,1,m,jbeg(myid),jend(myid),
     $     myid,nproc,MPI_COMM_WORLD,MY_MPI_REAL,'a')
      call write_array(aloc  ,1,m,jbeg(myid),jend(myid),myid,40+myid)
      call write_array(adebug,1,m,jbeg(myid),jend(myid),myid,50+myid)
#else /* DEBUG_COMM */
c     /* initialize the physical boundaries; set interior to zero*/
      call assign(aloc,1,m,jbeg(myid),jend(myid),m,n,myid,nproc,'b',2)
c
      filename='init.dat'
      call write2d_p(aloc,m,n,jbeg(myid),jend(myid),jbeg,jend,
     $     MPI_COMM_WORLD,MY_MPI_REAL,myid,nproc,filename)
c     
      t1=MPI_WTIME()
c     
      do 100 iter=1,nitmx       !this is the main iteration loop
c     
         call jacobi_step(aloc,adebug,1,m,jbeg(myid),jend(myid),m,n,
     $        myid,nproc,diffl,valmxl) 
c     
         if(lconverged(diffl,valmxl,tol,iter,myid,nproc,
     $        MPI_COMM_WORLD,MY_MPI_REAL))goto 110
c     
         call exchange_ghostline1(aloc,m,jbeg(myid),jend(myid),
     $        myid,nproc,MPI_COMM_WORLD,MY_MPI_REAL)
c
 100  continue
 110  continue
c
      t2=MPI_WTIME()
      if(myid.eq.0)then
         write(*,*)'time for iterative loop: ',t2-t1,' sec ',
     $        '(unit 55)'
         write(55,*)m,n,nproc,t2-t1,iter
      endif
c
c     /* measure the error */
      call assign(adebug,1,m,jbeg(myid),jend(myid),m,n,myid,nproc,
     $     'a',2)
      call check_array(aloc,adebug,1,m,jbeg(myid),jend(myid),
     $     myid,nproc,MPI_COMM_WORLD,MY_MPI_REAL,'a')
c
c     /* output of the result */
      filename='numeric.dat'
      call write2d_p(aloc,m,n,jbeg(myid),jend(myid),jbeg,jend,
     $     MPI_COMM_WORLD,MY_MPI_REAL,myid,nproc,filename)
c
      filename='analytic.dat'
      call write2d_p(adebug,m,n,jbeg(myid),jend(myid),jbeg,jend,
     $     MPI_COMM_WORLD,MY_MPI_REAL,myid,nproc,filename)
c     
#endif /* DEBUG_COMM */
c     
c     /* say 'bye' before finalizing MPI */
 111  write(*,*)myid,' about to finalize MPI'
      call MPI_FINALIZE(ierr)
      write(*,*)'bye. ',myid      
c     
      stop
      end
c-----|---------------------------------------------------------------|
      subroutine MPE_DECOMP1D_all(n,numprocs,s,e,chi)
c
c     /* input                                                       */
c     /* n             data length to be decomposed (1:n)            */
c     /* numprocs      no. of processors                             */
c     /* chi           'p' - shift the resulting pointers by 1       */
c     /*               [else] nothing                                */
c     /*                                                             */
c     /* output                                                      */
c     /* s(0:numprocs-1) pointer to start of memory local to procs   */
c     /* e(0:numprocs-1) pointer to end of memory local to procs     */
c
      integer n, numprocs, myid, s(0:numprocs-1), e(0:numprocs-1)
      integer nlocal
      integer deficit
      character*1 chi
c     /* chi  'p' augment count by one (to account for ghost)        */
c     
      do myid=0,numprocs-1
         nlocal  = n / numprocs
         s(myid) = myid * nlocal + 1
         deficit = mod(n,numprocs)
         s(myid) = s(myid) + min(myid,deficit)
         if (myid .lt. deficit) then
            nlocal = nlocal + 1
         endif
         e(myid) = s(myid) + nlocal - 1
         if (e(myid) .gt. n .or. myid .eq. numprocs-1) e(myid) = n
         if(chi.eq.'p')then
c     /* finally: augment all counts by one */
            s(myid) = s(myid) + 1
            e(myid) = e(myid) + 1
         endif
      enddo
      return
      end
c----+----------------------------------------------------------------|
      subroutine assign(a,if,il,jf,jl,m,n,my_rank,nproc,cghost,itype)
#include "realdef.inc"
c
c     /***************************************************************/
c     /* initializes the matrix a to some values.                    */
c     /*                                                             */
c     /* input:                                                      */
c     /* a(if-1:il+1,jf-1:jl+1)  matrix in local memory              */
c     /* my_rank                 processor id                        */
c     /* cghost                  character*1 flag which determines   */
c     /*                         whether the ghost cells (in the 2nd */
c     /*                         direction!) should be set to zero   */
c     /*                         or also filled with data.           */
c     /*                         ='g' or 'a': fill the ghosts w. data*/
c     /*                         ='b': fill the phys. boundary only  */
c     /*                         [else] set ghosts to zero, else data*/
c     /* itype                   type of the requested values.       */
c     /*                         =1: a(i,j)=i*j                      */
c     /*                         =2: a(i,j)=sin(pi*xi)*exp(-pi*yj)   */
c     /*                                                             */
c     /* output:                                                     */
c     /* a(if-1:il+1,jf-1:jl+1)  newly assigned matrix               */
c     /***************************************************************/
c
      dimension a(if-1:il+1,jf-1:jl+1)
      character*1 cghost
c----+----------------------------------------------------------------|
      pi=4.*atan(1.)
      if(cghost.eq.'g'.or.cghost.eq.'a')then !fill ghost-cells in 'j' with data
         jff=jf-1
         jll=jl+1
      else                      !set j-ghosts to zero
         jff=jf
         jll=jl
         do i=if-1,il+1
            a(i,jf-1)=0.0
            a(i,jl+1)=0.0
         enddo
      endif
      if(itype.eq.1)then
c     /* initialize the matrix as: a(i,j)=i*j */
         do j=jff,jll
            do i=if-1,il+1
               a(i,j)=float(i*j)
            enddo
         enddo
      elseif(itype.eq.2)then
c     /* initialize the matrix as zero, except physical boundaries:   */
c     /* aij=sin(pi*xi)*exp(-pi*yj) with xi=i/(m+1) \forall 0<=i<=m+1 */
c     /*                            with yj=j/(n+1) \forall 0<=j<=n+1 */
         do j=jf-1,jl+1
            do i=if-1,il+1
               a(i,j)=0.0       !all to zero
            enddo
         enddo
         if(cghost.ne.'b')then  !fill interior (and ghosts in case of 'g'/'a')
            do j=jff,jll
               yy=float(j)/float(n+1)
               do i=if-1,il+1
                  xx=float(i)/float(m+1)
                  a(i,j)=solution(xx,yy,itype)
               enddo
            enddo
         endif
         do j=jff,jll
            yy=float(j)/float(n+1)
            i=if-1
            xx=float(i)/float(m+1)
            a(i,j)=solution(xx,yy,itype)
            i=il+1
            xx=float(i)/float(m+1)
            a(i,j)=solution(xx,yy,itype)
         enddo
         if(my_rank.eq.0)then   !i am the first processor
            j=jf-1
            yy=float(j)/float(n+1)
            do i=if-1,il+1
               xx=float(i)/float(m+1)
               a(i,j)=solution(xx,yy,itype)
            enddo
         endif
         if(my_rank.eq.nproc-1)then !i am the last proc.
            j=jl+1
            yy=float(j)/float(n+1)
            do i=if-1,il+1
               xx=float(i)/float(m+1)
               a(i,j)=solution(xx,yy,itype)
            enddo
         endif
      else
         write(*,*)'assign: this type not implemented: ',itype
         stop
      endif
c
      return
      end subroutine assign
c----+----------------------------------------------------------------|
      subroutine check_array(a,aref,if,il,jf,jl,my_rank,nproc,mpi_comm,
     $     irealtype,cghost)
#include "realdef.inc"
      include "mpif.h"
c
c     /***************************************************************/
c     /* compares the matrix with some reference values.             */
c     /*                                                             */
c     /* input:                                                      */
c     /* a(if-1:il+1,jf-1:jl+1)    matrix in local memory            */
c     /* aref(if-1:il+1,jf-1:jl+1) matrix containing reference vals  */
c     /* my_rank                   processor id                      */
c     /* cghost                  character*1 flag which determines   */
c     /*                         whether the ghost cells (in the 2nd */
c     /*                         direction!) should be checked.      */
c     /*                         ='g': check _only_ the ghost-cells  */
c     /*                         ='a': check _all_ the array         */
c     /*                         [else] check only the interior      */
c     /*                                                             */
c     /* output:                                                     */
c     /* -                                                           */
c     /***************************************************************/
c
      dimension a(if-1:il+1,jf-1:jl+1),aref(if-1:il+1,jf-1:jl+1)
      character*1 cghost
c----+----------------------------------------------------------------|
      if(cghost.eq.'g')then     !check _only_ the ghost cells
                                !end-procs: don't check boundaries
         if(my_rank.gt.0)then
            jff=jf-1
         else
            jff=jf
         endif
         if(my_rank.lt.nproc-1)then
            jll=jl+1
         else
            jll=jl
         endif
         jskip=jll-jff
      elseif(cghost.eq.'a')then !check all of the array
                                !end-procs: don't check boundaries
         if(my_rank.gt.0)then
            jff=jf-1
         else
            jff=jf
         endif
         if(my_rank.lt.nproc-1)then
            jll=jl+1
         else
            jll=jl
         endif
         jskip=1
      else                      !check only the "interior"
         jff=jf
         jll=jl
         jskip=1
      endif
c
      errmx=0.
      valmx=0.
      do j=jff,jll,jskip
         do i=if-1,il+1
            err=abs(a(i,j)-aref(i,j))
            errmx=max(errmx,err)
            valmx=max(valmx,aref(i,j))
         enddo
      enddo
c
      CALL MPI_REDUCE(errmx,errmxg,1,irealtype,
     &     MPI_MAX,0,mpi_comm,ierr)
      CALL MPI_REDUCE(valmx,valmxg,1,irealtype,
     &     MPI_MAX,0,mpi_comm,ierr)
      
      if(my_rank.eq.0)write(*,*)'maximum relative error: ',
     $     errmxg/valmxg,errmxg,valmxg
c
      return
      end subroutine check_array
c----+----------------------------------------------------------------|
      subroutine write_array(a,if,il,jf,jl,my_rank,iunit)
#include "realdef.inc"
c
c     /***************************************************************/
c     /* writes the matrix a to a unit.                              */
c     /*                                                             */
c     /* input:                                                      */
c     /* a(if-1:il+1,jf-1:jl+1)    matrix in local memory            */
c     /* my_rank                   processor id                      */
c     /* iunit                     unit for writing (='6': screen)   */
c     /*                                                             */
c     /* output:                                                     */
c     /* -                                                           */
c     /***************************************************************/
c
      dimension a(if-1:il+1,jf-1:jl+1)
c----+----------------------------------------------------------------|
      write(iunit,110)my_rank
      do i=if-1,il+1
         do j=jf-1,jl+1
            write(iunit,111)my_rank,i,j,a(i,j)
         enddo
      enddo
c
 110  format('write_array: myid=',i5)
 111  format(1x,3(i5,1x),1(e15.8,1x))
c
      return
      end subroutine write_array
c----+----------------------------------------------------------------|
      REALN function solution(x,y,itype)
#include "realdef.inc"
      pi=4.*atan(1.)
c
      if(itype.eq.2)then
         solution=sin(pi*x)*exp(-pi*y)
      else
         write(*,*)'solution: this type not implmented! ',itype
         stop
      endif
c     
      return
      end function solution
c----+----------------------------------------------------------------|
      subroutine jacobi_step(a,tmp,if,il,jf,jl,m,n,my_rank,nproc,
     $     diff,valmx)
#include "realdef.inc"
c
c     /***************************************************************/
c     /* performs one step of jacobi iteration for 2d matrix in 1d   */
c     /* cartesian domain decomposition along 2nd index.             */
c     /*                                                             */
c     /* input:                                                      */
c     /* a(if-1:il+1,jf-1:jl+1)  matrix in local memory              */
c     /* tmp(if-1:il+1,jf-1:jl+1) work array                         */
c     /* my_rank                 processor id                        */
c     /*                                                             */
c     /* output:                                                     */
c     /* a(if-1:il+1,jf-1:jl+1)  updated matrix                      */
c     /*                         (updated only for if:il,jf:jl)      */
c     /* diff                    the r.m.s. of the increment         */
c     /*                         relative to the maximum value of the*/
c     /*                         function (for convergence measure)  */
c     /***************************************************************/
c
      dimension a(if-1:il+1,jf-1:jl+1),tmp(if-1:il+1,jf-1:jl+1)
c----+----------------------------------------------------------------|
      pi=4.*atan(1.)
c
      diff=0.0
      valmx=0.0
      do j=jf,jl
         do i=if,il
            tmp(i,j)=(a(i,j+1)+a(i+1,j)+a(i-1,j)+a(i,j-1))*.25
            diff=max(diff,abs(tmp(i,j)-a(i,j)))
            valmx=max(valmx,abs(tmp(i,j)))
         enddo
      enddo
      do j=jf,jl
         do i=if,il
            a(i,j)=tmp(i,j)
         enddo
      enddo
c
      return
      end subroutine jacobi_step
c----+----------------------------------------------------------------|
      logical function lconverged(diffl,valmxl,tol,iter,
     $     my_rank,nproc,mpi_comm,irealtype)
#include "realdef.inc"
      include "mpif.h"
c     
      lconverged=.false.
c     
      CALL MPI_Allreduce(diffl,diffg,1,irealtype,
     &     MPI_MAX,mpi_comm,ierr)
      CALL MPI_Allreduce(valmxl,valmxg,1,irealtype,
     &     MPI_MAX,mpi_comm,ierr)
c     
      diffl=diffl/valmxg
      diffg=diffg/valmxg
      if(diffg.le.tol)then
         lconverged=.true.
         write(*,*)my_rank,' we have converged!',diffl,diffg,iter
      endif
      diffl=diffg
#ifdef DEBUG_CONVG
      if(my_rank.eq.0)write(*,*)'increment: ',diffl
#endif
c     
      return
      end function lconverged
c---- +----------------------------------------------------------------|
      subroutine write2d_p(a,m,n,jf,jl,jbeg,jend,
     $     icomm,irealtype,my_rank,nproc,filename)
#include "realdef.inc"
      include "mpif.h"
      integer istat(MPI_STATUS_SIZE)
      dimension a(0:m+1,jf-1:jl+1)
      dimension jbeg(0:nproc-1),jend(0:nproc-1)
      character*80 filename
      parameter(iout=30)
c-----|---------------------------------------------------------------|
      nline=m+2
#ifdef DEBUG_WRITE
      write(*,*)my_rank,' hello inside write2d_p'
#endif
      if(my_rank.eq.0)then
         write(*,*)'writing ascii to ',filename
         open(iout,file=filename)
         jll=jl
         if(nproc.eq.1)jll=jll+1
         call do_write2d(a,m,jf,jl,jf-1,jll,iout)
         do iproc=1,nproc-1
#ifdef DEBUG_WRITE
            write(*,*)my_rank,' about to recv from ',iproc
#endif
            call MPI_RECV(a,(jend(iproc)-jbeg(iproc)+3)*nline,
     $           irealtype,iproc,1,icomm,istat,ierr)
#ifdef DEBUG_WRITE
            write(*,*)my_rank,' recvd from ',iproc
#endif
            jll=jend(iproc)
            if(iproc.eq.nproc-1)jll=jend(iproc)+1
            call do_write2d(a,m,jbeg(iproc),jend(iproc),jbeg(iproc),
     $           jll,iout)
         enddo
         close(iout)
      else
#ifdef DEBUG_WRITE
         write(*,*)my_rank,' about to send to 0'
#endif
         call MPI_SEND(a,(jend(my_rank)-jbeg(my_rank)+3)*nline,
     $        irealtype,0,1,icomm,istat,ierr) !includes ghosts!!!
#ifdef DEBUG_WRITE
         write(*,*)my_rank,' sent to 0'
#endif
      endif
c
      call MPI_BARRIER(icomm,ierr)
c
      return
      end subroutine write2d_p
c----+----------------------------------------------------------------|
      subroutine do_write2d(a,m,jf,jl,jff,jll,iunit)
#include "realdef.inc"
      dimension a(0:m+1,jf-1:jl+1)
c
      do j=jff,jll
         do i=0,m+1
            write(iunit,111)i,j,a(i,j)
         enddo
         write(iunit,112)
      enddo
c     
 111  format(1x,2(i5,1x),1(e15.8,1x))
 112  format(' ')
c
      return
      end subroutine do_write2d
c----+----------------------------------------------------------------|
