Trapezoidal Integration

In this example, we want to parallelize the computation of a simple numerical integral via the trapezoidal rule. The latter is given by

\[\int_{a}^{b}f(x)\,dx \approx h \sum_{i=1}^{N}\frac{f(x_{i-1})+f(x_{i})}{2}.\]

The function to be integrated is the following.

f(x) = 4 * √(1 - x^2)

f (generic function with 1 method)

The analytic result of the definite integral (from 0 to 1) is known to be $\pi$.

Sequential

Naturally, we implement the trapezoidal rule as a straightforward, sequential for loop.

function trapezoidal(a, b, n; h = (b - a) / n)
    y = (f(a) + f(b)) / 2.0
    for i in 1:(n - 1)
        x = a + i * h
        y = y + f(x)
    end
    return y * h
end

trapezoidal (generic function with 1 method)

Let's compute the integral of f above and see if we get the expected result. For simplicity, we choose N, the number of panels used to discretize the integration interval, as a multiple of the number of available Julia threads.

using Base.Threads: nthreads

N = nthreads() * 1_000_000

10000000

Calling trapezoidal we do indeed find the (approximate) value of $\pi$.

trapezoidal(0, 1, N) ≈ π

true

Parallel

Our strategy is the following: Divide the integration interval among the available Julia threads. On each thread, use the sequential trapezoidal rule to compute the partial integral. It is straightforward to implement this strategy with tmapreduce. The map part is, essentially, the application of trapezoidal and the reduction operator is chosen to be + to sum up the local integrals.

using OhMyThreads

function trapezoidal_parallel(a, b, N)
    n = N ÷ nthreads()
    h = (b - a) / N
    return tmapreduce(+, 1:nthreads()) do i
        local α = a + (i - 1) * n * h # the local keywords aren't necessary but good practice
        local β = α + n * h
        trapezoidal(α, β, n; h)
    end
end

# or equivalently
#
# function trapezoidal_parallel(a, b, N)
#     n = N ÷ nthreads()
#     h = (b - a) / N
#     @tasks for i in 1:nthreads()
#         @set reducer=+
#         local α = a + (i - 1) * n * h
#         local β = α + n * h
#         trapezoidal(α, β, n; h)
#     end
# end

trapezoidal_parallel (generic function with 1 method)

First, we check the correctness of our parallel implementation.

trapezoidal_parallel(0, 1, N) ≈ π

true

Then, we benchmark and compare the performance of the sequential and parallel versions.

using BenchmarkTools
@btime trapezoidal(0, 1, $N);
@btime trapezoidal_parallel(0, 1, $N);

  24.348 ms (0 allocations: 0 bytes)
  2.457 ms (69 allocations: 6.05 KiB)

Because the problem is trivially parallel - all threads to the same thing and don't need to communicate - we expect an ideal speedup of (close to) the number of available threads.

nthreads()

10

This page was generated using Literate.jl.