We recalled the composite trapezoidal rule. A very simple implementation in C is in traprule0.c . With 100 million function evaluations we get pi correct up to its first 8 decimal places. Here is the output of a run:
prompt]$ /tmp/traprule0 Approximation for pi = 3.141592648260279e+00 with error = -5.330e-09 prompt]$While this goes fast on our computer, we mentioned Romberg integration (extrapolation on the results of the composite trapezoidal rule) as a recommended algorithm for this slowly converging process. A first parallel program is obtained via data partitioning: in the embarrassingly parallel program compute_pi , every processor computes its own integration interval. There is only communication at the very end, when the integrals at the subintervals are gathered at the manager node. Instead of an MPI_Gather, followed by a summation, the program compute_pi uses MPI_Reduce, a command thus far left as an exercise. In this program we keep the number of subintervals modest (one thousand), we will see the accuracy improve as the number of processors increases. Below are the results of this experiment:
Executing the script "run_compute_pi" with content
mpirun -np 2 compute_pi mpirun -np 3 compute_pi mpirun -np 4 compute_pi mpirun -np 5 compute_pi mpirun -np 6 compute_pi mpirun -np 7 compute_pi mpirun -np 8 compute_pi mpirun -np 9 compute_pi mpirun -np 10 compute_pi mpirun -np 11 compute_pi mpirun -np 12 compute_pi mpirun -np 13 compute_pi mpirun -np 14 compute_piproduces the following output:
Approximation for pi = 3.141579505912095e+00 with error = -1.315e-05 Total wall time : 0.000306 seconds. Approximation for pi = 3.141585496863981e+00 with error = -7.157e-06 Total wall time : 0.000293 seconds. Approximation for pi = 3.141588005148171e+00 with error = -4.648e-06 Total wall time : 0.004175 seconds. Approximation for pi = 3.141589327430631e+00 with error = -3.326e-06 Total wall time : 0.004124 seconds. Approximation for pi = 3.141590123292210e+00 with error = -2.530e-06 Total wall time : 0.078507 seconds. Approximation for pi = 3.141590645645186e+00 with error = -2.008e-06 Total wall time : 0.074936 seconds. Approximation for pi = 3.141591010111094e+00 with error = -1.643e-06 Total wall time : 0.008084 seconds. Approximation for pi = 3.141591276267971e+00 with error = -1.377e-06 Total wall time : 0.062965 seconds. Approximation for pi = 3.141591477611349e+00 with error = -1.176e-06 Total wall time : 0.105379 seconds. Approximation for pi = 3.141591634269885e+00 with error = -1.019e-06 Total wall time : 0.008072 seconds. Approximation for pi = 3.141591758992224e+00 with error = -8.946e-07 Total wall time : 0.108264 seconds. Approximation for pi = 3.141591860203561e+00 with error = -7.934e-07 Total wall time : 0.206518 seconds. Approximation for pi = 3.141591943672518e+00 with error = -7.099e-07 Total wall time : 0.110734 seconds.
We see indeed the error decrease as we use more and more processors. By the way, the computational speedup is optimal. As the computational time is negligible, the wall time equals the startup time and fluctuates. Except for the illustration of data partitioning, this is not a realistic (or interesting) parallel computation.
To increase the amount of computational work, we sketched recursive integration to evaluate multiple integrals: a k-dimensional integral would require (n+1)^k function evaluations. Even for a modest number n of subintervals, the amount of computational work grows very fast.
More interesting parallel algorithms are obtained by adaptive integration, for which we need an error estimator. A simple principle is to compare two consecutive approximations. The program traprule1 modifies our first implementation of the composite trapezoidal rule to give an error estimate basically for free, without little extra work. Its output is
prompt] /tmp/traprule1 Integral 3.141592648262026e+00 with error estimate 4.12e-13 and error -5.33e-09 prompt]For slowly convergent processes as ours, we have to be aware that the error estimator will underestimate the error. Still, as an adaptive integration method takes the error estimate as a stop criterion, the addition of tiny numbers is not so meaningful anyway; and again, we better switch to extrapolation.
The lecture ended by sketching the bisection approach for automatic integration. On p processors, the manager would, based on a tree with subdivisions of the integration interval and error estimates, decide which interval to split next in p equal parts. We will return to this application area in a later part of the course.
Bibliography