1 c**********************************************************************
2 c twod.f - a solution to the Poisson problem by using Jacobi
3 c interation on a 2-d decomposition
5 c .... the rest of this is from pi3.f to show the style ...
8 c 1) receives the number of rectangles used in the approximation.
9 c 2) calculates the areas of it's rectangles.
10 c 3) Synchronizes for a global summation.
11 c Node 0 prints the result.
15 c pi the calculated result
16 c n number of points of integration.
17 c x midpoint of each rectangle's interval
18 c f function to integrate
19 c sum,pi area of rectangles
20 c tmp temporary scratch space for global summation
23 c This code is included (without the prints) because one version of
24 c MPICH SEGV'ed (probably because of errors in handling send/recv of
25 c MPI_PROC_NULL source/destination).
27 c****************************************************************************
31 parameter (maxn = 128)
32 double precision a(maxn,maxn), b(maxn,maxn), f(maxn,maxn)
34 integer myid, numprocs, it, rc, comm2d, ierr, stride
35 integer nbrleft, nbrright, nbrtop, nbrbottom
36 integer sx, ex, sy, ey
39 double precision diff2d, diffnorm, dwork
40 double precision t1, t2
42 data periods/2*.false./
45 call MPI_COMM_RANK( MPI_COMM_WORLD, myid, ierr )
46 call MPI_COMM_SIZE( MPI_COMM_WORLD, numprocs, ierr )
47 c print *, "Process ", myid, " of ", numprocs, " is alive"
50 c Get the size of the problem
56 c print *, 'About to do bcast on ', myid
57 call MPI_BCAST(nx,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
60 c Get a new communicator for a decomposition of the domain. Let MPI
61 c find a "good" decomposition
65 call MPI_DIMS_CREATE( numprocs, 2, dims, ierr )
66 call MPI_CART_CREATE( MPI_COMM_WORLD, 2, dims, periods, .true.,
69 c Get my position in this communicator
71 call MPI_COMM_RANK( comm2d, myid, ierr )
72 c print *, "Process ", myid, " of ", numprocs, " is alive"
74 c My neighbors are now +/- 1 with my rank. Handle the case of the
75 c boundaries by using MPI_PROCNULL.
76 call fnd2dnbrs( comm2d, nbrleft, nbrright, nbrtop, nbrbottom )
77 c print *, "Process ", myid, ":",
78 c * nbrleft, nbrright, nbrtop, nbrbottom
80 c Compute the decomposition
82 call fnd2ddecomp( comm2d, nx, sx, ex, sy, ey )
83 c print *, "Process ", myid, ":", sx, ex, sy, ey
85 c Create a new, "strided" datatype for the exchange in the "non-contiguous"
88 call mpi_Type_vector( ey-sy+1, 1, ex-sx+3,
89 $ MPI_DOUBLE_PRECISION, stride, ierr )
90 call mpi_Type_commit( stride, ierr )
93 c Initialize the right-hand-side (f) and the initial solution guess (a)
95 call twodinit( a, b, f, nx, sx, ex, sy, ey )
97 c Actually do the computation. Note the use of a collective operation to
98 c check for convergence, and a do-loop to bound the number of iterations.
100 call MPI_BARRIER( MPI_COMM_WORLD, ierr )
103 call exchng2( b, sx, ex, sy, ey, comm2d, stride,
104 $ nbrleft, nbrright, nbrtop, nbrbottom )
105 call sweep2d( b, f, nx, sx, ex, sy, ey, a )
106 call exchng2( a, sx, ex, sy, ey, comm2d, stride,
107 $ nbrleft, nbrright, nbrtop, nbrbottom )
108 call sweep2d( a, f, nx, sx, ex, sy, ey, b )
109 dwork = diff2d( a, b, nx, sx, ex, sy, ey )
110 call MPI_Allreduce( dwork, diffnorm, 1, MPI_DOUBLE_PRECISION,
111 $ MPI_SUM, comm2d, ierr )
112 if (diffnorm .lt. 1.0e-5) goto 20
113 if (myid .eq. 0) print *, 2*it, ' Difference is ', diffnorm
115 if (myid .eq. 0) print *, 'Failed to converge'
118 c if (myid .eq. 0) then
119 c print *, 'Converged after ', 2*it, ' Iterations in ', t2 - t1,
124 call MPI_Type_free( stride, ierr )
125 call MPI_Comm_free( comm2d, ierr )
126 call MPI_FINALIZE(rc)
129 c Perform a Jacobi sweep for a 2-d decomposition
131 subroutine sweep2d( a, f, n, sx, ex, sy, ey, b )
132 integer n, sx, ex, sy, ey
133 double precision a(sx-1:ex+1, sy-1:ey+1), f(sx-1:ex+1, sy-1:ey+1),
134 + b(sx-1:ex+1, sy-1:ey+1)
139 h = 1.0d0 / dble(n+1)
142 b(i,j) = 0.25 * (a(i-1,j)+a(i,j+1)+a(i,j-1)+a(i+1,j)) -
148 subroutine exchng2( a, sx, ex, sy, ey,
150 $ nbrleft, nbrright, nbrtop, nbrbottom )
152 integer sx, ex, sy, ey, stride
153 double precision a(sx-1:ex+1, sy-1:ey+1)
154 integer nbrleft, nbrright, nbrtop, nbrbottom, comm2d
155 integer status(MPI_STATUS_SIZE), ierr, nx
158 c These are just like the 1-d versions, except for less data
159 call MPI_SENDRECV( a(sx,ey), nx, MPI_DOUBLE_PRECISION,
161 $ a(sx,sy-1), nx, MPI_DOUBLE_PRECISION,
162 $ nbrbottom, 0, comm2d, status, ierr )
163 call MPI_SENDRECV( a(sx,sy), nx, MPI_DOUBLE_PRECISION,
165 $ a(sx,ey+1), nx, MPI_DOUBLE_PRECISION,
166 $ nbrtop, 1, comm2d, status, ierr )
168 c This uses the "strided" datatype
169 call MPI_SENDRECV( a(ex,sy), 1, stride, nbrright, 0,
170 $ a(sx-1,sy), 1, stride, nbrleft, 0,
171 $ comm2d, status, ierr )
172 call MPI_SENDRECV( a(sx,sy), 1, stride, nbrleft, 1,
173 $ a(ex+1,sy), 1, stride, nbrright, 1,
174 $ comm2d, status, ierr )
179 c The rest of the 2-d program
181 double precision function diff2d( a, b, nx, sx, ex, sy, ey )
182 integer nx, sx, ex, sy, ey
183 double precision a(sx-1:ex+1, sy-1:ey+1), b(sx-1:ex+1, sy-1:ey+1)
191 sum = sum + (a(i,j) - b(i,j)) ** 2
197 subroutine twodinit( a, b, f, nx, sx, ex, sy, ey )
198 integer nx, sx, ex, sy, ey
199 double precision a(sx-1:ex+1, sy-1:ey+1), b(sx-1:ex+1, sy-1:ey+1),
200 & f(sx-1:ex+1, sy-1:ey+1)
211 c Handle boundary conditions
236 c This file contains a routine for producing a decomposition of a 1-d array
237 c when given a number of processors. It may be used in "direct" product
238 c decomposition. The values returned assume a "global" domain in [1:n]
240 subroutine MPE_DECOMP1D( n, numprocs, myid, s, e )
241 integer n, numprocs, myid, s, e
245 nlocal = n / numprocs
246 s = myid * nlocal + 1
247 deficit = mod(n,numprocs)
248 s = s + min(myid,deficit)
249 if (myid .lt. deficit) then
253 if (e .gt. n .or. myid .eq. numprocs-1) e = n
257 c This routine show how to determine the neighbors in a 2-d decomposition of
258 c the domain. This assumes that MPI_Cart_create has already been called
260 subroutine fnd2dnbrs( comm2d,
261 $ nbrleft, nbrright, nbrtop, nbrbottom )
262 integer comm2d, nbrleft, nbrright, nbrtop, nbrbottom
266 call MPI_Cart_shift( comm2d, 0, 1, nbrleft, nbrright, ierr )
267 call MPI_Cart_shift( comm2d, 1, 1, nbrbottom, nbrtop, ierr )
272 c Note: THIS IS A TEST PROGRAM. THE ACTUAL VALUES MOVED ARE NOT
273 c CORRECT FOR A POISSON SOLVER.
275 subroutine fnd2ddecomp( comm2d, n, sx, ex, sy, ey )
277 integer n, sx, ex, sy, ey
278 integer dims(2), coords(2), ierr
281 call MPI_Cart_get( comm2d, 2, dims, periods, coords, ierr )
283 call MPE_DECOMP1D( n, dims(1), coords(1), sx, ex )
284 call MPE_DECOMP1D( n, dims(2), coords(2), sy, ey )