1#if HAVE_CONFIG_H
2#   include "config.fh"
3#endif
4c $Id: jacobi.F,v 1.9 2003-06-24 19:01:59 d3h325 Exp $
5      program jacobi
6*
7*     Jacobi iterative method for solving Uxx + Uyy = 1 on [0,1]x[0,1] area
8*     Global Arrays are used to store and access U and exact solution vectors
9*     J.N., 06.01.94, works with GA 1.2
10*
11      implicit none
12#include "mafdecls.fh"
13#include "global.fh"
14      integer heap, stack
15c
16c***  Intitialize a message passing library
17c
18#include "mp3.fh"
19c
20      heap = 80000
21      stack = heap
22      call ga_initialize()                                   ! initialize GAs
23      if (.not. ma_init(MT_DBL, heap, stack))                ! initialize MA
24     $    call ga_error("ma init failed",heap+stack)
25#ifdef GA_TRACE
26      call trace_init(10000)                     ! initialize trace
27#endif
28      call iterate()                                         ! do the work
29#ifdef GA_TRACE
30      call trace_end(ga_nodeid())                 ! end trace
31#endif
32c
33      if(ga_nodeid().eq.0) print *,'All tests successful '
34c
35      call ga_terminate()                                    ! terminate GAs
36c
37      call MP_FINALIZE()
38      end
39
40
41      subroutine iterate()
42      implicit none
43#include "mafdecls.fh"
44#include "global.fh"
45      integer m                               ! grid size
46      integer g_u, g_ux, g_diff               ! handles to global arrays
47      parameter (m = 10)
48      integer ilo, ihi, jlo, jhi, ld,  ni, nj
49      integer i, k
50      MA_ACCESS_INDEX_TYPE indexl, index, index_first, indexd    ! MA addressing indices
51      integer h_l, h_d                        ! handles to MA objects
52      integer me, nproc                       ! my processor & number of procs
53      double precision eps                    ! stopping number
54      double precision S, W, C, N, E          ! coef. of G matrix and d
55      logical have_data                          ! do I "own" any data ?
56      double precision norm, product, uxnorm
57      data eps /1d-5/
58c
59c***  check parallel environment
60      me = ga_nodeid()
61      nproc = ga_nnodes()
62c
63c***  create global arrays: g_u - approx. solution,  g_ux - exact solution
64c                           g_diff will store g_u - g_ux
65      if (.not. ga_create(MT_DBL, m, m, 'u', 1, 1, g_u))
66     $     call ga_error(' ga_create failed ',0)
67      if (.not. ga_create(MT_DBL, m, m, 'ux', 1, 1, g_ux))
68     $     call ga_error(' ga_create failed ',0)
69      if (.not. ga_create(MT_DBL, m, m, 'diff', 1, 1, g_diff))
70     $     call ga_error(' ga_create failed ',0)
71c
72c***  check which piece, if any, of g_ux (and g_u) I own
73      call ga_distribution(g_ux, me, ilo, ihi, jlo, jhi)
74      ni = ihi - ilo + 1                 ! number of 'local' rows
75      nj = jhi - jlo + 1                 ! number of 'local' columns
76      have_data = ni .gt. 0 .and. nj .gt. 0
77c
78      if(have_data)then
79c
80c***     allocate memory for D
81         if(.not. ma_push_get(MT_DBL,ni*nj,'d',h_d, indexd))
82     &            call ga_error('memory allocation failed',0)
83c
84c***     access local data
85         call ga_access(g_ux, ilo, ihi, jlo, jhi, index, ld)
86c
87c***     create coefficients of the jacobi iteration matrix G and exact solution
88         call jacinit(m, S, W, C, N, E, DBL_MB(indexd), ni,
89     &                DBL_MB(index), ilo, ihi, jlo, jhi, ld)
90c
91c***     release access to the local data (the data were updated !)
92         call ga_release_update(g_ux, ilo, ihi, jlo, jhi)
93c
94c***     allocate memory for a copy of my piece of g_u and neighb. grid points
95         if(.not. ma_push_get(MT_DBL,(ni+2)*(nj+2),'loc',h_l, indexl))
96     &            call ga_error('memory allocation failed',0)
97c
98c***     zero the allocated memory
99         do i = 1, (ni+2)*(nj+2)
100            DBL_MB(indexl - 1 + i) = 0d0        ! indexl points to first element
101         enddo
102      endif
103c
104c***  initial guess for u -- zero
105      call ga_zero(g_u)
106c
107c***  the stopping test is ||u-ux||/||ux||, ||.|| - second norm
108      uxnorm = sqrt(ga_ddot(g_ux, g_ux))                   ! ||g_ux||
109c
110c***  Now, synchronize and then iterate
111      k = 0
11210    call ga_sync()
113           k = k + 1
114c
115           if(have_data)then
116c
117c...          first, get access to the local piece of g_u
118              call ga_access(g_u, ilo, ihi, jlo, jhi, index, ld)
119c
120c...          compute next approximation of u
121              call nextu(ni, nj, dbl_mb(index), ld, dbl_mb(indexl),
122     &                   ni+1, S, W, C, N, E, dbl_mb(indexd) )
123c
124c...          release access to the local data (the data were updated !)
125              call ga_release_update(g_u, ilo, ihi, jlo, jhi)
126           endif
127c
128c...       compute the stopping number -- second norm
129           call ga_add(1d0, g_u, -1d0, g_ux, g_diff)    ! g_diff = g_u - g_ux
130           product = ga_ddot(g_diff, g_diff)             ! <g_diff, g_diff>
131           norm    = sqrt(product)/uxnorm
132           if(me.eq.0 .and. Mod(k,10).eq.1) print *,k,' error= ',norm
133c
134      if(norm .gt. eps) then
135         if(.not. have_data) goto 10
136c
137c...     not done yet, get a copy of my piece of g_u and neighboring grid points
138c        ... determine where from we should copy data -- consider topology
139         index_first = indexl
140         if(jlo.eq.1) index_first = index_first + ni+2
141         if(ilo.eq.1) index_first = index_first + 1
142c
143         call ga_get(g_u, max(1,ilo-1), min(m,ihi+1),
144     &                    max(1,jlo-1), min(m,jhi+1),
145     &                    dbl_mb(index_first), ni+2 )
146         goto 10
147      endif
148c
149      if(me.eq.0) then
150         print *,' converged in ',k, ' iterations, error= ',norm
151      endif
152c
153c***  deallocate MA memory and destroy global arrays
154      if(have_data)then
155        if(.not. ma_pop_stack(h_l))call ga_error('invalid handle ?',0)
156        if(.not. ma_pop_stack(h_d))call ga_error('invalid handle ?',0)
157      endif
158      if(.not. ga_destroy(g_u)) call ga_error('invalid handle ?',0)
159      if(.not. ga_destroy(g_ux)) call ga_error('invalid handle ?',0)
160      if(.not. ga_destroy(g_diff)) call ga_error('invalid handle ?',0)
161c
162      call ga_sync()
163      end
164
165
166      subroutine nextu(ni, nj, u, ld, u0, ld0, S, W, C, N, E, D )
167      implicit none
168      integer ni, nj, ld, ld0
169      double precision u(1:ld,1:nj), u0(0:ld0,0:nj), D(1:ni,1:nj)
170      double precision S, W, C, N, E
171      integer i, j
172c
173      do j = 1, nj
174         do i = 1, ni
175            u(i,j) = S*u0(i+1,j) + W*u0(i,j-1) +
176     &               N*u0(i-1,j) + E*u0(i,j+1) + D(i,j)
177         enddo
178      enddo
179      end
180
181
182      subroutine jacinit(m, S, W, C, N, E, D, ldd, ux,
183     &                   ilo, ihi, jlo, jhi, ldu )
184      implicit none
185#include "mafdecls.fh"
186#include "global.fh"
187      double precision S, W, C, N, E, temp
188      integer ilo, ihi, jlo, jhi, ldu, ldd, i, j
189      integer m, ni, nj, jj, ii
190      double precision h, ux(ldu,m), D(ldd,m)
191      if(m.le.1)call ga_error('jacinit: wrong value of m',0)
192      S = .25
193      W = .25
194      N = .25
195      E = .25
196      h = 1d0/(m+1)
197      ni = ihi - ilo + 1                 ! number of 'local' rows
198      nj = jhi - jlo + 1                 ! number of 'local' columns
199      do j = 1, nj
200         jj = jlo + j -1
201         do i = 1, ni
202            temp = 0d0
203            ii = ilo + i -1
204            if(jj .eq. 1) temp = temp + 1
205            if(ii .eq. 1) temp = temp + 1
206            if(ii .eq. m) temp = temp + jj*h + 1
207            if(jj .eq. m) temp = temp + ii*h + 1
208            D(i,j) = temp/4d0
209            ux(i,j)= 1d0 + h*h*ii*jj
210         enddo
211      enddo
212      end
213