HPC - Monte Carlo estimation of the area of the union of circles

The ile mpi-circles.c contains a serial implementation of a Monte Carlo algorithm that estimates the area of the union of \(N\) circles. Let cx[i], cy[i], and r[i] be the coordinates of the center of circle \(i\) and its radius. All circles are entirely contained within the bounding square with opposites corners \((0, 0)\) and \((1000, 1000)\).

Circles may overlap in whole or in part; therefore, it is not easy to compute the area of their union. We implement a Monte Carlo algorithm to estimate the area; the idea is similar to the estimation of the value of \(\pi\) by generating random points, and is as follows:

• Generate \(K\) random points uniformly distributed inside the bounding square \((0, 0)\), \((1000, 1000)\). Let \(c\) be the number of points that fall within at least one circle.

• The area \(A\) of the union of the circles can be estimated as \(A \approx 1000 \times 1000 \times c/K\). In other words, the area \(A\) is the product of the area of the bounding square and the fraction of points \(c/K\) that falls within at least one circle.

Figure 1 illustrates the idea.

The file mpi-circles.c uses a serial algorithm where process 0 performs the whole computation. The purpose of this exercise is to distribute the computation among all MPI processes. The input file that contains the coordinates of the circles must be read by process 0 only; therefore, only process 0 knows the number \(N\) of circles and their coordinates, so it must send all required information to the other processes. The program must work correctly for any value of \(N\), \(K\) and number \(P\) of MPI processes.

To do so, each process \(p\) generates \(K/P\) points and test each point with all the circles; at the end, process \(p\) knows number \(C_p\) of points that fall inside at least one circle. The master computes \(C = \sum_{p=0}^{P-1} C_p\) using a reduction, and estimates the area using the formula given above. Therefore, each process must receive a full copy of the arrays cx[], cy[] and r[].

Note: You might be tempted to distribute the circles among the MPI processes, so that each process handles \(N/P\) circles. This would not work: why?

To compile:

    mpicc -std = c99 -Wall -Wpedantic mpi-circles.c -o mpi-circles

To execute:

    mpirun -n P ./mpi-circles N input_file_name

Example:

    mpirun -n 4 ./mpi-circles 10000 circles-1000.in

Files

• mpi-circles.c
• circles-gen.c (to generate random inputs)
• circles-1000.in
• circles-10000.in