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( a, 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( b, 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 c 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 if (myid .eq. 0) then
127 print *, ' No Errors'
129 call MPI_FINALIZE(rc)
132 c Perform a Jacobi sweep for a 2-d decomposition
134 subroutine sweep2d( a, f, n, sx, ex, sy, ey, b )
135 integer n, sx, ex, sy, ey
136 double precision a(sx-1:ex+1, sy-1:ey+1), f(sx-1:ex+1, sy-1:ey+1),
137 + b(sx-1:ex+1, sy-1:ey+1)
142 h = 1.0d0 / dble(n+1)
145 b(i,j) = 0.25 * (a(i-1,j)+a(i,j+1)+a(i,j-1)+a(i+1,j)) -
151 subroutine exchng2( a, b, sx, ex, sy, ey,
153 $ nbrleft, nbrright, nbrtop, nbrbottom )
155 integer sx, ex, sy, ey, stride
156 double precision a(sx-1:ex+1, sy-1:ey+1),
157 $ b(sx-1:ex+1, sy-1:ey+1)
158 integer nbrleft, nbrright, nbrtop, nbrbottom, comm2d
159 integer status(MPI_STATUS_SIZE), ierr, nx
162 c These are just like the 1-d versions, except for less data
163 call MPI_SENDRECV( b(ex,sy), nx, MPI_DOUBLE_PRECISION,
165 $ a(sx-1,sy), nx, MPI_DOUBLE_PRECISION,
166 $ nbrbottom, 0, comm2d, status, ierr )
167 call MPI_SENDRECV( b(sx,sy), nx, MPI_DOUBLE_PRECISION,
169 $ a(ex+1,sy), nx, MPI_DOUBLE_PRECISION,
170 $ nbrtop, 1, comm2d, status, ierr )
172 c This uses the "strided" datatype
173 call MPI_SENDRECV( b(sx,ey), 1, stride, nbrright, 0,
174 $ a(sx,sy-1), 1, stride, nbrleft, 0,
175 $ comm2d, status, ierr )
176 call MPI_SENDRECV( b(sx,sy), 1, stride, nbrleft, 1,
177 $ a(sx,ey+1), 1, stride, nbrright, 1,
178 $ comm2d, status, ierr )
183 c The rest of the 2-d program
185 double precision function diff2d( a, b, nx, sx, ex, sy, ey )
186 integer nx, sx, ex, sy, ey
187 double precision a(sx-1:ex+1, sy-1:ey+1), b(sx-1:ex+1, sy-1:ey+1)
195 sum = sum + (a(i,j) - b(i,j)) ** 2
201 subroutine twodinit( a, b, f, nx, sx, ex, sy, ey )
202 integer nx, sx, ex, sy, ey
203 double precision a(sx-1:ex+1, sy-1:ey+1), b(sx-1:ex+1, sy-1:ey+1),
204 & f(sx-1:ex+1, sy-1:ey+1)
215 c Handle boundary conditions
240 c This file contains a routine for producing a decomposition of a 1-d array
241 c when given a number of processors. It may be used in "direct" product
242 c decomposition. The values returned assume a "global" domain in [1:n]
244 subroutine MPE_DECOMP1D( n, numprocs, myid, s, e )
245 integer n, numprocs, myid, s, e
249 nlocal = n / numprocs
250 s = myid * nlocal + 1
251 deficit = mod(n,numprocs)
252 s = s + min(myid,deficit)
253 if (myid .lt. deficit) then
257 if (e .gt. n .or. myid .eq. numprocs-1) e = n
261 c This routine show how to determine the neighbors in a 2-d decomposition of
262 c the domain. This assumes that MPI_Cart_create has already been called
264 subroutine fnd2dnbrs( comm2d,
265 $ nbrleft, nbrright, nbrtop, nbrbottom )
266 integer comm2d, nbrleft, nbrright, nbrtop, nbrbottom
270 call MPI_Cart_shift( comm2d, 0, 1, nbrleft, nbrright, ierr )
271 call MPI_Cart_shift( comm2d, 1, 1, nbrbottom, nbrtop, ierr )
276 c Note: THIS IS A TEST PROGRAM. THE ACTUAL VALUES MOVED ARE NOT
277 c CORRECT FOR A POISSON SOLVER.
279 subroutine fnd2ddecomp( comm2d, n, sx, ex, sy, ey )
281 integer n, sx, ex, sy, ey
282 integer dims(2), coords(2), ierr
285 call MPI_Cart_get( comm2d, 2, dims, periods, coords, ierr )
287 call MPE_DECOMP1D( n, dims(1), coords(1), sx, ex )
288 call MPE_DECOMP1D( n, dims(2), coords(2), sy, ey )