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