1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
|
# Jacobi computation, using persistent requests
from mpi4py import MPI
try:
import numpy
except ImportError:
raise SystemExit
n = 5 * MPI.COMM_WORLD.Get_size()
# compute number of processes and myrank
p = MPI.COMM_WORLD.Get_size()
myrank = MPI.COMM_WORLD.Get_rank()
# compute size of local block
m = n//p
if myrank < (n - p * m):
m = m + 1
#compute neighbors
if myrank == 0:
left = MPI.PROC_NULL
else:
left = myrank - 1
if myrank == p - 1:
right = MPI.PROC_NULL
else:
right = myrank + 1
# allocate local arrays
A = numpy.empty((n+2, m+2), dtype=float, order='f')
B = numpy.empty((n, m), dtype=float, order='f')
A.fill(1)
A[0, :] = A[-1, :] = 0
A[:, 0] = A[:, -1] = 0
# create persistent requests
tag = 0
sreq1 = MPI.COMM_WORLD.Send_init((B[:, 0], MPI.DOUBLE), left, tag)
sreq2 = MPI.COMM_WORLD.Send_init((B[:, -1], MPI.DOUBLE), right, tag)
rreq1 = MPI.COMM_WORLD.Recv_init((A[:, 0], MPI.DOUBLE), left, tag)
rreq2 = MPI.COMM_WORLD.Recv_init((A[:, -1], MPI.DOUBLE), right, tag)
reqlist = [sreq1, sreq2, rreq1, rreq2]
for req in reqlist:
assert req != MPI.REQUEST_NULL
# main loop
converged = False
while not converged:
# compute boundary columns
N, S = A[ :-2, 1], A[2:, 1]
E, W = A[1:-1, 0], A[1:-1, 2]
C = B[:, 0]
numpy.add(N, S, C)
numpy.add(C, E, C)
numpy.add(C, W, C)
C *= 0.25
N, S = A[ :-2, -2], A[2:, -2]
E, W = A[1:-1, -3], A[1:-1, -1]
C = B[:, -1]
numpy.add(N, S, C)
numpy.add(C, E, C)
numpy.add(C, W, C)
C *= 0.25
# start communication
#MPI.Prequest.Startall(reqlist)
for r in reqlist:
r.Start()
# compute interior
N, S = A[ :-2, 2:-2], A[2, 2:-2]
E, W = A[1:-1, 2:-2], A[1:-1, 2:-2]
C = B[:, 1:-1]
numpy.add(N, S, C)
numpy.add(E, C, C)
numpy.add(W, C, C)
C *= 0.25
A[1:-1, 1:-1] = B
# complete communication
MPI.Prequest.Waitall(reqlist)
# convergence
myconv = numpy.allclose(B, 0)
loc_conv = numpy.asarray(myconv, dtype='i')
glb_conv = numpy.asarray(0, dtype='i')
MPI.COMM_WORLD.Allreduce([loc_conv, MPI.INT],
[glb_conv, MPI.INT],
op=MPI.LAND)
converged = bool(glb_conv)
# free persistent requests
for req in reqlist:
req.Free()
|
