# How can I reduce the communication bottleneck of a parallel explicit finite difference scheme?

Suppose i was trying to solve a parabolic PDE (heat equation) on a rectangular domain using an explicit finite difference scheme. I am storing my solution vector in a matrix form (because it closely resembles the shape of the domain). I know that I can break the matrix into rectangular blocks and send information across the boundaries of each block (as needed) using ghost points. However, when I use a fine grid, the processors have to send a lot of information between them, and this seems to slow down the overall performance. Is there anything that I can do to increase the scalability of the algorithm? Is there a theoretical limit to how much improvement I can expect through the use of ghost points?

Jonathan's answer gives a useful performance model for a simple explicit method. In practice, memory bandwidth and network latency are more common performance limiters than network bandwidth.

To improve the efficiency, it is worth considering:

• high order methods reach a given accuracy with fewer grid cells while also doing more computational work per cell, choose variants with small vertex or edge separators (like spectral elements or DG) to keep communication volume low
• implicit methods solve the problem without the painful stability criteria for explicit methods, but tend to require more communication per step and are more involved to program
• Chebyshev Runge-Kutta methods with many stages greatly extend the stability region relative to low-stage methods for pure diffusion, but these are fragile if the eigenvalues drift from the negative real axis
• $s$-step methods (computing multiple time steps in one sweep) further hide latency and reduce memory bandwidth limitations, but need to communicate multiple ghost layers, so the constants are poor for higher order and higher dimension problems, especially with smallish subdomains
• What methods, then, would you recommend for an implicit method? – Paul Jan 16 '12 at 14:22
• For long time steps or steady-state, it is desperately hard to beat multigrid. On regular grids, it is simple to implement and one F-cycle will provably reach discretization error (about 50 flops per vertex in 3D). You can use a lighter weight cycle for moderate length time steps. For the constant coefficient heat equation (or piecewise constant with simple coefficients), the Green's function solution can be evaluated with FMM, giving you a (highly parallel) $\mathcal O(N / P) + \mathcal O(\log_{2^d} P)$ to evaluate the solution at any time $t$ (with no time discretization error). – Jed Brown Jan 16 '12 at 14:37
• Thanks for the tip... i somewhat suspected that multigrid would be a likely candidate.. – Paul Jan 16 '12 at 14:48
• Implicit solves will avoid stringent stability costs, but often by transporting the problem of stability from the user (i.e. in choosing the timestep) to the linear solver, which ends up doing the same amount of work as an explicit method with appropriately chosen timesteps. That is hidden behind the linear solver black-box, but it's all effectively polynomial iterations. Unless of course you can factorize the system, then implicit wins hands down. – Reid.Atcheson Mar 23 '13 at 5:52
• Also for the second-to-last point, the closely related DUMKA schemes can be used (say you're negative definite but non-normal). dumkaland.org These pull the stability region a little away from the neagative real axis. RK-Chebyshev methods have dubious intersections with real axis at single points, leading to the mentioned instability. DUMKA fixes this, to an extent. – Reid.Atcheson Mar 23 '13 at 5:58

Tying this question into a previous one ( When analyzing a parallel algorithm, how do you take communication costs into account? ) on modeling parallel efficiency, let's create a simple model where we take a look at how the computation behaves as we change the number of processors.

Let's say we have a grid of size $N^d$, and $P$ processors -- we'll assume that $P$ is a perfect square/cube/whatever and that $N$ is evenly divided by $P^{1/d}$ (if it isn't, it just introduces some modest load imbalance).

Let's say the actual explicit computation involves $n_\mathrm{ops}$ operations (memory, or floating point) operations per zone, each of which has (time) cost $c_\mathrm{ops}$. Further, let's model the time spent to perform a communication with a latency and a bandwidth; it takes a message of even zero size some time $l$ to go fron node to node, and there is in addition a finite bandwidth $b$ to the network, so that a message of size $s$ takes time $l + s/b$. (In practice, things are usually a bit more complicated than this; local neighbours might be on-node, for instance, and splitting bandwidth between cores; and distant neighbours might not even be on the same switch. So you might need to use something other than the nominal network parameters. Still, this is a sensible start.)

The total computation time taken by the decomposed problem per timestep is $$t_\mathrm{comp} = \frac{N^d}{P} n_\mathrm{ops} c_\mathrm{ops}$$ and the total communication time is $$t_\mathrm{comm} = 2^d \left ( l + \frac{n_g N^{d-1}}{b P^{(d-1)/d}} \right )$$ where at each timestep, one has to send to each of one's $2^d$ neighbours slice of $n_g \times (N/P^{1/d})^{(d-1)}$ guard cells. (eg, in 2d, $n_g \times N/\sqrt{P}$).

(Further, explicit parabolic methods generally for stability require timesteps that scale as $\Delta x^{-2}$, so that we would need to take a number of timesteps that scales as $N^2$ to reach a given physical elapsed time. That's relevant for how the computational cost scales wtih $N$, but not parallel scalability, so we won't need this here.)

The parallel speedup here would be $$\begin{eqnarray*} S & = & \frac{t(N,P=1)}{t(N,P)} \\ & = & \frac{N^d n_\mathrm{ops} c_\mathrm{ops}}{\frac{N^d}{P} n_\mathrm{ops} c_\mathrm{ops} + 2^d \left ( l + \frac{n_g N^{d-1}}{b P^{(d-1)/d}} \right )} \\ & = & \frac{P}{1 + \frac{P}{N^d n_\mathrm{ops} c_\mathrm{ops}} 2^d \left ( l + \frac{n_g N^{d-1}}{b P^{(d-1)/d}} \right )} \\ & = & \frac{P}{1 + \frac{t_\mathrm{comm}}{t_\mathrm{comp}}} \\ \end{eqnarray*}$$ -- that is, the speedup would be exactly $P$ if the communications time were zero, which isn't all that surprising. Looking a bit more closely at the communcations to computation ratio: $$\begin{eqnarray*} \frac{t_\mathrm{comm}}{t_\mathrm{comp}} & = & \frac{P}{N^d n_\mathrm{ops} c_\mathrm{ops}} 2^d \left ( l + \frac{n_g N^{d-1}}{b P^{(d-1)/d}} \right ) \\ \end{eqnarray*}$$

Let's look at the latency-dominated scenario; that is, you're sending small messages, so that the time for communications is dominated just by the number of messages you're sending. (You would always get to this point eventually, as the message size, as we'll see, goes down with number of processors.) In this case, $$\begin{eqnarray*} \frac{t_\mathrm{comm}}{t_\mathrm{comp}} & = & \frac{P}{N^d n_\mathrm{ops} c_\mathrm{ops}} 2^d l \\ & \propto & P / N^d \end{eqnarray*}$$ Here the relative parallel overhead goes linearly as you add more processors -- this is not because communications is getting more expensive, but because the computational cost is going down linearly in P as we're adding more processors and each only has to do $N^d/P$ work.

In the bandwidth dominated case -- for very large messages or very slow networks -- you'd instead have: $$\begin{eqnarray*} \frac{t_\mathrm{comm}}{t_\mathrm{comp}} & = & \frac{P}{N^d n_\mathrm{ops} c_\mathrm{ops}} 2^d \frac{n_g N^{d-1}}{b P^{(d-1)/d}} \\ & \propto & P^{1/d} / N \end{eqnarray*}$$ In this regime -- which you are likely in at the beginning (especially in 3d) with comparitively large messages, especially on ethernet -- the overhead is much more modest, scaling as the cube root of P (in 3d), because as you increase the number of processes, you are also reducing the size of the messages you have to send because the grid is more finely divided and you have less overlap with your neighbouring processor.

It's fun to put some numbers in; let's say we do a weak scaling study in 3D, with $N^3/P$ fixed at $256^3$, and with 1 guardcell of double precision numbers, with gigabit ethernet so say latency is 100 microseconds and 1/bandwidth is about 0.06 microseconds/double. Let's say you have to do 100 operations per zone, with each operation taking 1 nanoseconds.

Then you have $$\begin{eqnarray*} \frac{t_\mathrm{comm}}{t_\mathrm{comp}} & = & \frac{P}{N^d n_\mathrm{ops} c_\mathrm{ops}} 2^d \left( l + \frac{n_g N^{d-1}}{b P^{(d-1)/d}} \right) \\ & = & \frac{8}{256^3 (100 \times 10^{-9} \mathrm{sec})} \left ( 10^{-4} \mathrm{sec} + 256^2 \cdot 0.06 \cdot 10^{-6} \mathrm{sec} \right ) \\ & = & 0.02 \\ \end{eqnarray*}$$

Note that this is constant in the weak scaling regime, since the work per processor -- even the communication work! -- is constant. Since this is constant, you would have your speedup be

$$\begin{eqnarray*} S & = & \frac{P}{1 + 0.02} \\ & = & 0.98 P \end{eqnarray*}$$ or an efficiency of 98%, which is quite good. Note too that in this regime you are definitely bandwidth, rather than latency, dominated.

You can play with other scaling regimes, or play with other network transports. Let's play with strong scaling, say a constant $N^3 = 1024^3$: here you'd have $$\begin{eqnarray*} \frac{t_\mathrm{comm}}{t_\mathrm{comp}} & = & \frac{P}{N^d n_\mathrm{ops} c_\mathrm{ops}} 2^d \left( l + \frac{n_g N^{d-1}}{b P^{(d-1)/d}} \right) \\ & = & \frac{8 P}{1024^3 (10 \times 10^{-9} \mathrm{sec})} \left ( 10^{-4} \mathrm{sec} + \frac{1024^2}{P^{2/3}} \cdot 0.06 \cdot 10^{-6} \mathrm{sec} \right ) \\ & = & 0.0745 P \, \mathrm{sec}^{-1} \, \left ( 10^{-4} \,\mathrm{sec} + 6.29 \times 10^{-2} \,\mathrm{sec}\, P^{-2/3} \right ) \\ \end{eqnarray*}$$

This is plotted above; obviously, treating things continuously when we've assumed that P is a perfect cube and divides N is problematic, but it still gives you a picture of what's going on. For the first while, bandwidth dominates your communication overhead, but eventualy you switch over to latency. Here, with the relative overhead increasing with $P$ (again, due to the fact that computation is decreasing), you definitely do not get a constant efficiency, but one which falls with $P$.

Now, can you do better than this -- improve things with different approaches to coding? Yes, there's a couple classic things you can do. Probably most important is to overlap communications with computation. Rather than communicating the guard cells, then doing all the processor's computation, so that the total time taken is $t_\mathrm{comm} + t_\mathrm{comp}$, you can start the process of communicating the guardcells, then do the inner zones computation while communications is taking place (with infiniband, or with communications threads), and then once the receive is done, do the outer zones calculation. This allows you to hide some of the communications time by doing real computation while it is occuring. Then the time is, instead, $\max{t_\mathrm{comm}, t_\mathrm{comp}}$ -- or more precisely, $\max{t_\mathrm{comm}, t_\mathrm{comp,inner}} + t_\mathrm{comp,outer}$. This can be a significant improvement!

Note too here that in such a code there is often some kind of reduction operation at the end of each step; that isn't being modelled here. Neither is I/O of any sort.

• Ok, so let's nip this in the bud before it gets too entrenched. I've seen a few times on scicomp people taking it upon themselves to just edit an answer to actually change what is said. That is not cool; on Stack Overflow, for instance, this would be considered extremely rude. The "How To Edit" box that appears when you edit makes this fairly explicit; don't change meaning. If you think something should be changed, make a comment, don't just decide to rewrite what was said -- even if it is something clearly wrong like a math mistake (as I made, above, last night). Ok? – user389 Jan 16 '12 at 15:13
• In particular, the high opcount was intended to include amortized memory costs. If you don't like that assumption, say so, don't just erase it. – user389 Jan 16 '12 at 15:18
• My apologies. Note that an op count of 10 per cell is almost spot-on for low-order finite difference methods on Blue Gene/P (it is between our measured values for 7- and 27-point stencils). – Jed Brown Jan 16 '12 at 16:30