• API
• FAQ
• Tools
• Trends
• Archive
SHARE
TWEET

# CIS_mandelMPI.py

calcpage Feb 16th, 2013 1,573 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. def mandelbrot(x, y, maxit):
2.     c = x + y*1j
3.     z = 0 + 0j
4.     it = 0
5.     while abs(z) < 2 and it < maxit:
6.         z = z**2 + c
7.         it += 1
8.     return it
9.
10. x1, x2 = -2.0, 1.0
11. y1, y2 = -1.0, 1.0
12. w, h = 150, 100
13. maxit = 127
14.
15. from mpi4py import MPI
16. import numpy
17.
18. comm = MPI.COMM_WORLD
19. size = comm.Get_size()
20. rank = comm.Get_rank()
21.
22. # number of rows to compute here
23. N = h // size + (h % size > rank)
24.
25. # first row to compute here
26. start = comm.scan(N)-N
27.
28. # array to store local result
29. Cl = numpy.zeros([N, w], dtype='i')
30.
31. # compute owned rows
32.
33. dx = (x2 - x1) / w
34. dy = (y2 - y1) / h
35. for i in range(N):
36.     y = y1 + (i + start) * dy
37.     for j in range(w):
38.         x = x1 + j * dx
39.         Cl[i, j] = mandelbrot(x, y, maxit)
40.
41. # gather results at root (process 0)
42.
43. counts = comm.gather(N, root=0)
44. C = None
45. if rank == 0:
46.     C = numpy.zeros([h, w], dtype='i')
47.
48. rowtype = MPI.INT.Create_contiguous(w)
49. rowtype.Commit()
50.
51. comm.Gatherv(sendbuf=[Cl, MPI.INT],recvbuf=[C, (counts, None), rowtype],root=0)
52.
53. rowtype.Free()
54.
55. if comm.rank == 0:
56.     from matplotlib import pyplot
57.     pyplot.imshow(C, aspect='equal')
58.     pyplot.spectral()
59.     pyplot.show()
RAW Paste Data
Top