program main  

c*********************************************************************72
c
cc MAIN is the main program for POISSON_MPI.
c
c  Discussion:
c
c    POISSON_BAND is the main program for a solver of the Poisson problem.
c
c    The Poisson equation
c
c      DEL^2 U(X,Y) = F(X,Y)
c
c    is solved on the unit square [0,1] x [0,1] using a grid of [0,NX+1] by
c    [0,NY+1] evenly spaced points.  
c
c    The linear system generated by the discretization is solved
c    using Jacobi iteration.
c
c  Licensing:
c
c    This code is distributed under the MIT license. 
c
c  Modified:
c
c    15 October 2005
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    William Gropp, Ewing Lusk, Anthony Skjellum,
c    Using MPI: Portable Parallel Programming with the
c    Message-Passing Interface,
c    Second Edition,
c    MIT Press, 1999,
c    ISBN: 0262571323,
c    LC: QA76.642.G76.
c
c    Marc Snir, Steve Otto, Steven Huss-Lederman, David Walker, 
c    Jack Dongarra,
c    MPI: The Complete Reference,
c    Volume I: The MPI Core,
c    Second Edition,
c    MIT Press, 1998,
c    ISBN: 0-262-69216-3,
c     LC: QA76.642.M65.
c
      include "mpif.h" 

      integer maxn 
      parameter ( maxn = 128 ) 

      double precision a(maxn,maxn)
      double precision b(maxn,maxn)
      integer comm1d
      double precision diff
      double precision diffnorm
      double precision dwork 
      integer e
      double precision f(maxn,maxn)
      integer ierr 
      integer it
      integer my_id
      integer nbrbottom
      integer nbrtop
      integer nx
      integer ny 
      integer num_procs
      integer s
      double precision t1
      double precision t2 
c
c  Initialize MPI.
c
      call MPI_INIT ( ierr )
c
c  Get the individual process ID.
c
      call MPI_COMM_RANK ( MPI_COMM_WORLD, my_id, ierr ) 
c
c  Get the number of processes.
c
      call MPI_COMM_SIZE ( MPI_COMM_WORLD, num_procs, ierr )
c
c  Print a message.
c
      if ( my_id .eq. 0 ) then
        call timestamp ( )
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'POISSON_MPI - Master process:'
        write ( *, '(a)' ) '  FORTRAN77 version'
        write ( *, '(a)' ) '  An MPI program to solve the '
        write ( *, '(a)' ) '  Poisson equation.'
        write ( *, '(a)' ) '  The MPI message passing library is used.'
        write ( *, '(a)' ) ' '
        write ( *, '(a,i6)' ) '  The number of processes is ', num_procs
      else
        write ( *, '(a,i3,a)' ) '  Process ', my_id, ' is active'
      end if
c 
c  Get the size of the problem 
c 
      if ( my_id .eq. 0 ) then 
        nx = 50
      end if
 
      call MPI_BCAST ( nx, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr )
      ny = nx
c
c  Print a message.
c
      if ( my_id .eq. 0 ) then
        write ( *, '(a,i6)' ) 
     &    '  The number of interior X grid lines is ', nx
        write ( *, '(a,i6)' ) 
     &    '  The number of interior Y grid lines is ', ny
        write ( *, '(a)' ) ' '
      end if
c 
c  Get a new communicator for a decomposition of the domain 
c 
      call MPI_CART_CREATE ( MPI_COMM_WORLD, 1, num_procs, .false.,  
     &  .true., comm1d, ierr ) 
c 
c  Get my position in this communicator, and my neighbors 
c  
      call MPI_COMM_RANK( comm1d, my_id, ierr ) 

      call MPI_Cart_shift ( comm1d, 0, 1, nbrbottom, nbrtop, ierr ) 
c 
c  Compute the actual decomposition.
c      
      call decompose ( ny, num_procs, my_id, s, e ) 
c 
c  Initialize the right-hand-side F and the initial solution guess A. 
c 
      call onedinit ( a, b, f, nx, s, e ) 
c 
c  Actually do the computation.  Note the use of a collective operation to 
c  check for convergence, and a do-loop to bound the number of iterations. 
c 
      call MPI_BARRIER ( MPI_COMM_WORLD, ierr ) 

      t1 = MPI_WTIME ( ) 

      do it = 1, 400 

        call exchng1 ( a, nx, s, e, comm1d, nbrbottom, nbrtop ) 
        call sweep1d ( a, f, nx, s, e, b ) 

        call exchng1 ( b, nx, s, e, comm1d, nbrbottom, nbrtop ) 
        call sweep1d ( b, f, nx, s, e, a ) 

        dwork = diff ( a, b, nx, s, e ) 

        call MPI_Allreduce ( dwork, diffnorm, 1, MPI_DOUBLE_PRECISION,  
     &    MPI_SUM, comm1d, ierr ) 

        if ( my_id .eq. 0 ) then
          write ( *, '(a,i6,a,g14.6)' ) '  Step ', it, 
     &      '  Diff is ', diffnorm
        end if

        if ( diffnorm .lt. 1.0e-4 ) then
          if ( my_id .eq. 0 ) then 
            print *, '  Converged after ', it, ' iterations.'
          end if 
          go to 20 
        end if

      end do

      if ( my_id .eq. 0 ) then
        print *, '  The Jacobi iteration failed to converge.' 
      end if

20    continue 

      t2 = MPI_WTIME ( ) 

      if ( my_id .eq. 0 ) then
        print *, ' ' 
        print *, '  Terminated after ', t2 - t1, ' seconds.' 
      end if 
c
c  Terminate MPI.
c
      call MPI_FINALIZE ( ierr )
c
c  Terminate.
c
      if ( my_id .eq. 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'POISSON_BAND - Master process:'
        write ( *, '(a)' ) '  Normal end of execution.'
        write ( *, '(a)' ) ' '
        call timestamp ( )
      end if

      stop
      end 
      subroutine decompose ( n, num_procs, my_id, s, e )

c*********************************************************************72
c 
cc DECOMPOSE decomposes a 1-d array for a given number of processors.
c
c  Discussion:
c
c    This routine may be used in "direct" product decomposition.  The 
c    values returned assume a "global" domain in [1:n].
c
c  Licensing:
c
c    This code is distributed under the MIT license. 
c
c  Modified:
c
c    15 October 2005
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer deficit
      integer e
      integer my_id
      integer n
      integer num_procs
      integer nlocal
      integer s

      nlocal  = n / num_procs 
      s = my_id * nlocal + 1 
      deficit = mod ( n, num_procs ) 
      s = s + min ( my_id, deficit ) 

      if ( my_id .lt. deficit ) then 
        nlocal = nlocal + 1 
      end if 

      e = s + nlocal - 1 

      if ( e .gt. n .or. my_id .eq. num_procs - 1 ) then
        e = n 
      end if

      return 
      end 
      function diff ( a, b, nx, s, e ) 

c*********************************************************************72
c
cc DIFF computes a difference between arrays A and B.
c
c  Licensing:
c
c    This code is distributed under the MIT license. 
c
c  Modified:
c
c    15 October 2005
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer e
      integer nx
      integer s

      double precision a(0:nx+1,s-1:e+1)
      double precision b(0:nx+1,s-1:e+1)
      double precision diff
      integer i
      integer j

      diff = 0.0D+00 
      do j = s, e 
        do i = 1, nx 
          diff = diff + ( a(i,j) - b(i,j) ) ** 2 
        end do
      end do

      return 
      end 
      subroutine exchng1 ( a, nx, s, e, comm1d, nbrbottom, nbrtop ) 

c*********************************************************************72
c
cc EXCHNG1 exchanges data with "neighboring" processors.
c
c  Licensing:
c
c    This code is distributed under the MIT license. 
c
c  Modified:
c
c    15 October 2005
c
c  Author:
c
c    John Burkardt
c
      include 'mpif.h' 

      integer e
      integer nx
      integer s

      double precision a(0:nx+1,s-1:e+1) 
      integer comm1d
      integer ierr
      integer nbrbottom
      integer nbrtop 
      integer status(MPI_STATUS_SIZE)

      call MPI_SENDRECV (  
     &  a(1,e), nx, MPI_DOUBLE_PRECISION, nbrtop, 0,  
     &  a(1,s-1), nx, MPI_DOUBLE_PRECISION, nbrbottom, 0,  
     &  comm1d, status, ierr ) 

      call MPI_SENDRECV (  
     &  a(1,s), nx, MPI_DOUBLE_PRECISION, nbrbottom, 1,  
     &  a(1,e+1), nx, MPI_DOUBLE_PRECISION, nbrtop, 1,  
     &  comm1d, status, ierr ) 

      return 
      end 
      subroutine onedinit ( a, b, f, nx, s, e ) 

c*********************************************************************72
c
cc ONEDINIT initializes the arrays.
c
c  Licensing:
c
c    This code is distributed under the MIT license. 
c
c  Modified:
c
c    15 October 2005
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer e
      integer nx
      integer s

      double precision a(0:nx+1,s-1:e+1)
      double precision b(0:nx+1,s-1:e+1)
      double precision f(0:nx+1,s-1:e+1) 
      integer i
      integer j

      do j = s-1, e+1 
        do i = 0, nx+1 
          a(i,j) = 0.0D+00
          b(i,j) = 0.0D+00
          f(i,j) = 0.0D+00
        end do
      end do
c       
c  Handle boundary conditions 
c 
      do j = s, e 
        a(0,j) = 1.0D+00 
        b(0,j) = 1.0D+00 
        a(nx+1,j) = 0.0D+00 
        b(nx+1,j) = 0.0D+00 
      end do

      if ( s .eq. 1 ) then 
        do i = 1, nx 
          a(i,0) = 1.0D+00 
          b(i,0) = 1.0D+00 
        end do
      end if 

      return 
      end 
      subroutine sweep1d ( a, f, nx, s, e, b )

c*********************************************************************72
c 
cc SWEEP1D performs a Jacobi sweep for a 1-d decomposition. 
c
c  Licensing:
c
c    This code is distributed under the MIT license. 
c
c  Modified:
c
c    15 October 2005
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer e
      integer nx
      integer s

      double precision a(0:nx+1,s-1:e+1)
      double precision b(0:nx+1,s-1:e+1) 
      double precision f(0:nx+1,s-1:e+1)
      double precision h
      integer i
      integer j

      h = 1.0D+00 / dble ( nx + 1 ) 

      do j = s, e 
        do i = 1, nx 

          b(i,j) = 0.25D+00 * ( 
     &      a(i-1,j) + a(i,j+1) + a(i,j-1) + a(i+1,j) ) -  
     &      h * h * f(i,j) 

        end do
      end do 

      return 
      end 
      subroutine timestamp ( )

c*********************************************************************72
c
cc TIMESTAMP prints out the current YMDHMS date as a timestamp.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    12 January 2007
c
c  Author:
c
c    John Burkardt
c
c  Parameters:
c
c    None
c
      implicit none

      character * ( 8 ) ampm
      integer d
      character * ( 8 ) date
      integer h
      integer m
      integer mm
      character * ( 9 ) month(12)
      integer n
      integer s
      character * ( 10 ) time
      integer y

      save month

      data month /
     &  'January  ', 'February ', 'March    ', 'April    ', 
     &  'May      ', 'June     ', 'July     ', 'August   ', 
     &  'September', 'October  ', 'November ', 'December ' /

      call date_and_time ( date, time )

      read ( date, '(i4,i2,i2)' ) y, m, d
      read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm

      if ( h .lt. 12 ) then
        ampm = 'AM'
      else if ( h .eq. 12 ) then
        if ( n .eq. 0 .and. s .eq. 0 ) then
          ampm = 'Noon'
        else
          ampm = 'PM'
        end if
      else
        h = h - 12
        if ( h .lt. 12 ) then
          ampm = 'PM'
        else if ( h .eq. 12 ) then
          if ( n .eq. 0 .and. s .eq. 0 ) then
            ampm = 'Midnight'
          else
            ampm = 'AM'
          end if
        end if
      end if

      write ( *, 
     &  '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) 
     &  d, month(m), y, h, ':', n, ':', s, '.', mm, ampm

      return
      end