      program jacobi 
*
*     Jacobi iterative method for solving Uxx + Uyy = 1 on [0,1]x[0,1] area
*     Global Arrays are used to store and access U and exact solution vectors 
*     J.N., 06.01.94, works with GA 1.2
*
      implicit none
#include "mafdecls.fh"
#include "global.fh"
      integer heap, stack
      data heap,stack /2*40000/
c
      call pbeginf                                           ! fork processes
      if (.not. ma_init(MT_DBL, heap, stack)) call parerr(1) ! initialize MA
      call ga_initialize()                                   ! initialize GAs 
#ifdef GA_TRACE
      call trace_init(10000)				     ! initialize trace
#endif
      call iterate()                                         ! do the work      
#ifdef GA_TRACE
      call trace_end(ga_nodeid())			     ! end trace
#endif
      call ga_terminate()                                    ! terminate GAs 
      call pend()                                            ! terminate TCGMSG
      end


      subroutine iterate()
      implicit none
#include "mafdecls.fh"
#include "global.fh"
      integer m                               ! grid size
      integer g_u, g_ux, g_diff               ! handles to global arrays
      parameter (m = 10)
      integer ilo, ihi, jlo, jhi, ld,  ni, nj
      integer i, k
      integer indexl, index, index_first, indexd    ! MA addressing indices
      integer h_l, h_d                        ! handles to MA objects
      integer me, nproc                       ! my processor & number of procs 
      double precision eps                    ! stopping number
      double precision S, W, C, N, E          ! coef. of G matrix and d
      logical have_data	                      ! do I "own" any data ?
      double precision norm, product, uxnorm
      data eps /1d-5/
c
c***  check parallel environment
      me = ga_nodeid()
      nproc = ga_nnodes()
c
c***  create global arrays: g_u - approx. solution,  g_ux - exact solution
c                           g_diff will store g_u - g_ux 
      if (.not. ga_create(MT_DBL, m, m, 'u', 1, 1, g_u))
     $     call ga_error(' ga_create failed ',0)
      if (.not. ga_create(MT_DBL, m, m, 'ux', 1, 1, g_ux))
     $     call ga_error(' ga_create failed ',0)
      if (.not. ga_create(MT_DBL, m, m, 'diff', 1, 1, g_diff))
     $     call ga_error(' ga_create failed ',0)
c
c***  check which piece, if any, of g_ux (and g_u) I own
      call ga_distribution(g_ux, me, ilo, ihi, jlo, jhi)
      ni = ihi - ilo + 1                 ! number of 'local' rows 
      nj = jhi - jlo + 1                 ! number of 'local' columns 
      have_data = ni .gt. 0 .and. nj .gt. 0
c
      if(have_data)then
c
c***     allocate memory for D 
         if(.not. ma_push_get(MT_DBL,ni*nj,'d',h_d, indexd))
     &            call ga_error('memory allocation failed',0)
c
c***     access local data
         call ga_access(g_ux, ilo, ihi, jlo, jhi, index, ld)
c
c***     create coefficients of the jacobi iteration matrix G and exact solution
         call jacinit(m, S, W, C, N, E, DBL_MB(indexd), ni, 
     &                DBL_MB(index), ilo, ihi, jlo, jhi, ld)
c
c***     release access to the local data (the data were updated !)
         call ga_release_update(g_ux, ilo, ihi, jlo, jhi)
c
c***     allocate memory for a copy of my piece of g_u and neighb. grid points
         if(.not. ma_push_get(MT_DBL,(ni+2)*(nj+2),'loc',h_l, indexl))
     &            call ga_error('memory allocation failed',0)
c
c***     zero the allocated memory
         do i = 1, (ni+2)*(nj+2)
            DBL_MB(indexl - 1 + i) = 0d0        ! indexl points to first element
         enddo
      endif
c
c***  initial guess for u -- zero
      call ga_zero(g_u)
c
c***  the stopping test is ||u-ux||/||ux||, ||.|| - second norm
      uxnorm = sqrt(ga_ddot(g_ux, g_ux))                   ! ||g_ux||
c
c***  Now, synchronize and then iterate 
      k = 0
10    call ga_sync()
           k = k + 1
c
           if(have_data)then
c
c...          first, get access to the local piece of g_u
              call ga_access(g_u, ilo, ihi, jlo, jhi, index, ld)
c
c...          compute next approximation of u
              call nextu(ni, nj, dbl_mb(index), ld, dbl_mb(indexl),
     &                   ni+1, S, W, C, N, E, dbl_mb(indexd) ) 
c
c...          release access to the local data (the data were updated !)
              call ga_release_update(g_u, ilo, ihi, jlo, jhi)
           endif
c
c...       compute the stopping number -- second norm
           call ga_add(1d0, g_u, -1d0, g_ux, g_diff)    ! g_diff = g_u - g_ux 
           product = ga_ddot(g_diff, g_diff)             ! <g_diff, g_diff>
           norm    = sqrt(product)/uxnorm
           if(me.eq.0 .and. Mod(k,10).eq.1) print *,k,' error= ',norm
c
      if(norm .gt. eps) then
         if(.not. have_data) goto 10
c
c...     not done yet, get a copy of my piece of g_u and neighboring grid points
c        ... determine where from we should copy data -- consider topology 
         index_first = indexl
         if(jlo.eq.1) index_first = index_first + ni+2 
         if(ilo.eq.1) index_first = index_first + 1 
c
         call ga_get(g_u, max(1,ilo-1), min(m,ihi+1),
     &                    max(1,jlo-1), min(m,jhi+1),
     &                    dbl_mb(index_first), ni+2 ) 
         goto 10
      endif
c
      if(me.eq.0) then
         print *,' converged in ',k, ' iterations, error= ',norm
      endif 
c
c***  deallocate MA memory and destroy global arrays
      if(have_data)then
        if(.not. ma_pop_stack(h_l))call ga_error('invalid handle ?',0)
        if(.not. ma_pop_stack(h_d))call ga_error('invalid handle ?',0)
      endif
      if(.not. ga_destroy(g_u)) call ga_error('invalid handle ?',0)
      if(.not. ga_destroy(g_ux)) call ga_error('invalid handle ?',0)
      if(.not. ga_destroy(g_diff)) call ga_error('invalid handle ?',0)
c
      call ga_sync()
      end


      subroutine nextu(ni, nj, u, ld, u0, ld0, S, W, C, N, E, D ) 
      implicit none
      integer ni, nj, ld, ld0
      double precision u(1:ld,1:nj), u0(0:ld0,0:nj), D(1:ni,1:nj)
      double precision S, W, C, N, E
      integer i, j      
c
      do j = 1, nj
         do i = 1, ni
            u(i,j) = S*u0(i+1,j) + W*u0(i,j-1) +  
     &               N*u0(i-1,j) + E*u0(i,j+1) + D(i,j)
         enddo
      enddo
      end


      subroutine jacinit(m, S, W, C, N, E, D, ldd, ux,
     &                   ilo, ihi, jlo, jhi, ldu )
      implicit none
#include "mafdecls.fh"
#include "global.fh"
      double precision S, W, C, N, E, temp
      integer ilo, ihi, jlo, jhi, ldu, ldd, i, j
      integer m, ni, nj, jj, ii 
      double precision h, ux(ldu,m), D(ldd,m)
      if(m.le.1)call ga_error('jacinit: wrong value of m',0) 
      S = .25
      W = .25
      N = .25
      E = .25
      h = 1d0/(m+1)
      ni = ihi - ilo + 1                 ! number of 'local' rows 
      nj = jhi - jlo + 1                 ! number of 'local' columns 
      do j = 1, nj 
         jj = jlo + j -1
         do i = 1, ni 
            temp = 0d0
            ii = ilo + i -1
            if(jj .eq. 1) temp = temp + 1 
            if(ii .eq. 1) temp = temp + 1 
            if(ii .eq. m) temp = temp + jj*h + 1
            if(jj .eq. m) temp = temp + ii*h + 1
            D(i,j) = temp/4d0
            ux(i,j)= 1d0 + h*h*ii*jj
         enddo
      enddo
      end
