University:
Boston University
Course:
MET CS 201 | Introduction to Programming
Academic year:

2023
Views:

170

Pages:

20
Author:

vielay3og

Compute Pi Z π= 0 1 n−1 1X 4 4 dx ≈ 2 2 1+x n 1 + ( i+0.5 n ) i=0 Compute Pi – sequential 1 import math 2 3 4 5 6 7 8 9 def compute_pi(n): h = 1.0 / n s = 0.0 for i in range(n): x = h * (i + 0.5) s += 4.0 / (1.0 + x**2) return s * h 10 11 12 13 n = 10 pi = compute_pi(n) error = abs(pi - math.pi) 14 15 16 print ("pi is approximately %.16f, " "error is %.16f" % (pi, error)) Compute Pi – parallel [1] 1 2 from mpi4py import MPI import math 3 4 5 6 7 8 9 10 def compute_pi(n, start=0, step=1): h = 1.0 / n s = 0.0 for i in range(start, n, step): x = h * (i + 0.5) s += 4.0 / (1.0 + x**2) return s * h 11 12 13 14 comm = MPI.COMM_WORLD nprocs = comm.Get_size() myrank = comm.Get_rank() Compute Pi – parallel [2] 1 2 3 4 if myrank == 0: n = 10 else: n = None 5 6 n = comm.bcast(n, root=0) 7 8 mypi = compute_pi(n, myrank, nprocs) 9 10 pi = comm.reduce(mypi, op=MPI.SUM, root=0) 11 12 13 14 15 if myrank == 0: error = abs(pi - math.pi) print ("pi is approximately %.16f, " "error is %.16f" % (pi, error)) Exercise #4 Modify Compute Pi example to employ NumPy. Tip: you can convert a Python int/float object to a NumPy scalar with x = numpy.array(x). Mandelbrot Set Mandelbrot Set – sequential [1] 1 2 3 4 5 6 7 8 def mandelbrot(x, y, maxit): c = x + y*1j z = 0 + 0j it = 0 while abs(z) < 2 and it < maxit: z = z**2 + c it += 1 return it 9 10 11 12 13 x1, x2 y1, y2 w, h maxit = = = = -2.0, 1.0 -1.0, 1.0 150, 100 127 Mandelbrot Set – sequential [2] 1 2 3 4 5 6 7 8 9 import numpy C = numpy.zeros([h, w], dtype=’i’) dx = (x2 - x1) / w dy = (y2 - y1) / h for i in range(h): y = y1 + i * dy for j in range(w): x = x1 + j * dx C[i, j] = mandelbrot(x, y, maxit) 10 11 12 13 14 from matplotlib import pyplot pyplot.imshow(C, aspect=’equal’) pyplot.spectral() pyplot.show() Mandelbrot Set – partitioning 0 1 2 Block distribution 0 0 0 1 1 1 2 2 2 Cyclic distribution Mandelbrot Set – parallel, block [1] 1 2 3 4 5 6 7 8 def mandelbrot(x, y, maxit): c = x + y*1j z = 0 + 0j it = 0 while abs(z) < 2 and it < maxit: z = z**2 + c it += 1 return it 9 10 11 12 13 x1, x2 y1, y2 w, h maxit = = = = -2.0, 1.0 -1.0, 1.0 150, 100 127 Mandelbrot Set – parallel, block [2] 1 2 from mpi4py import MPI import numpy 3 4 5 6 comm = MPI.COMM_WORLD size = comm.Get_size() rank = comm.Get_rank() 7 8 9 # number of rows to compute here N = h // size + (h % size > rank) 10 11 12 # first row to compute here start = comm.scan(N)-N 13 14 15 # array to store local result Cl = numpy.zeros([N, w], dtype=’i’) Mandelbrot Set – parallel, block [3] 1 # compute owned rows 2 3 4 5 6 7 8 9 dx = (x2 - x1) / w dy = (y2 - y1) / h for i in range(N): y = y1 + (i + start) * dy for j in range(w): x = x1 + j * dx Cl[i, j] = mandelbrot(x, y, maxit) Mandelbrot Set – parallel, block [4] 1 # gather results at root (process 0) 2 3 4 5 6 counts = comm.gather(N, root=0) C = None if rank == 0: C = numpy.zeros([h, w], dtype=’i’) 7 8 9 rowtype = MPI.INT.Create_contiguous(w) rowtype.Commit() 10 11 12 13 comm.Gatherv(sendbuf=[Cl, MPI.INT], recvbuf=[C, (counts, None), rowtype], root=0) 14 15 rowtype.Free() Mandelbrot Set – parallel, block [5] 1 2 3 4 5 if comm.rank == 0: from matplotlib import pyplot pyplot.imshow(C, aspect=’equal’) pyplot.spectral() pyplot.show() Exercise #5 Measure the wall clock time Ti of local computations at each process for the Mandelbrot Set example with block and cyclic row distributions. What row distribution is better regarding load balancing? Tip: use Wtime() to measure wall time, compute the ratio Tmax /Tmin to compare load balancing. Dynamic Process Management I Useful in assembling complex distributed applications. Can couple independent parallel codes written in different languages. I Create new processes from a running program. – Comm.Spawn() and Comm.Get parent() I Connect two running applications together. – Comm.Connect() and Comm.Accept() Dynamic Process Management – Spawning Spawning new processes is a collective operation that creates an intercommunicator. I Local group is group of spawning processes (parent). I Remote group is group of new processes (child). I Comm.Spawn() lets parent processes spawn the child processes. It returns a new intercommunicator. I Comm.Get parent() lets child processes find intercommunicator to the parent group. Child processes have own COMM WORLD. I Comm.Disconnect() ends the parent–child connection. After that, both groups can continue running. Dynamic Process Management – Compute Pi (parent) 1 2 from mpi4py import MPI import sys, numpy 3 4 5 6 comm = MPI.COMM_SELF.Spawn(sys.executable, args=[’compute_pi-child.py’], maxprocs=5) 7 8 9 10 11 12 13 N = numpy.array(10, ’i’) comm.Bcast([N, MPI.INT], root=MPI.ROOT) PI = numpy.array(0.0, ’d’) comm.Reduce(None, [PI, MPI.DOUBLE], op=MPI.SUM, root=MPI.ROOT) comm.Disconnect() 14 15 16 17 error = abs(PI - numpy.pi) print ("pi is approximately %.16f, " "error is %.16f" % (PI, error)) Dynamic Process Management – Compute Pi (child) 1 2 from mpi4py import MPI import numpy 3 4 5 6 comm = MPI.Comm.Get_parent() size = comm.Get_size() rank = comm.Get_rank() 7 8 9 10 11 12 13 14 15 16 N = numpy.array(0, dtype=’i’) comm.Bcast([N, MPI.INT], root=0) h = 1.0 / N; s = 0.0 for i in range(rank, N, size): x = h * (i + 0.5) s += 4.0 / (1.0 + x**2) PI = numpy.array(s * h, dtype=’d’) comm.Reduce([PI, MPI.DOUBLE], None, op=MPI.SUM, root=0) 17 18 comm.Disconnect() Exercise #5 a) Implement the Compute Pi child code in C or C++ . Adjust the parent code in Python to spawn the new implementation. b) Compute and plot the Mandelbrot Set using spawning with parent/child codes implemented in Python. Tip: Reuse the provided parent code in Python and translate the child code in Fortran 90 to Python.