### 6.2.2.2 Matrix-vector product, partitioning by columns

Partition

$$A \rightarrow (A_0, A_1, \dots, A_{p-1}) \quad x \rightarrow \begin{pmatrix} x_0 \\ x_1 \\ \vdots \\ x_{p-1} \end{pmatrix}, \quad \text{and} \quad y \rightarrow \begin{pmatrix} y_0 \\ y_1 \\ \vdots \\ y_{p-1} \end{pmatrix},$$

where $A_j \in \mathbb{R}^{n_j \times n_j}$ and $x_j, y_j \in \mathbb{R}^{n_j}$ with $\sum_{j=0}^{p-1} n_j = n$ and $n_j \approx n/p$.

We will start by assuming that $A_j, x_j$, and $y_j$ are originally assigned to $\mathcal{P}_j$ (but now $A_i$ is a block of columns). In this algorithm by columns, processor $i$ can compute the length $n$ vector $A_ix_i$ without prior communication. These partial results then have to be added together

$$y \leftarrow \sum_i A_i x_i$$

in a *reduce-scatter* operation: each processor $i$ scatters a part $(A_i x_i)_j$ of its result to processor $j$. The receiving processors then perform a reduction, adding all these fragments:

$$y_j = \sum_i (A_i x_i)_j.$$

The algorithm with costs is then given by:

<table><thead><tr><th>Step</th><th>Cost (lower bound)</th></tr></thead><tbody><tr><td>Locally compute y<sup>(j)</sup> = A<sub>j</sub>x<sub>j</sub></td><td>≈ 2<sup>n<sub>j</sub>2</sup>γ</td></tr><tr><td>Reduce-scatter the y<sup>(j)</sup>s so that y<sub>i</sub> = &sum;<sup>p-1</sup><sub>j=0</sub> y<sub>i</sub><sup>(j)</sup> is on P<sub>i</sub></td><td>[log<sub>2</sub>(p)]α + <sup>p-1</sup><sub>p</sub>nβ + <sup>p-1</sup><sub>p</sub>nγ<br>≈ log<sub>2</sub>(p)α + n(β + γ)</td></tr></tbody></table>

**Cost analysis** The total cost of the algorithm is given by, approximately,

$$T_p^{\text{1D-col}}(n) = 2\frac{n^2}{p}\gamma + \underbrace{\log_2(p)\alpha + n(\beta + \gamma)}_{\text{Overhead}}.$$

Notice that this is identical to the cost $T_p^{\text{1D-row}}(n)$, except with $\beta$ replaced by $(\beta + \gamma)$. It is not hard to see that the conclusions about scalability are the same.

### 6.2.2.3 Two-dimensional partitioning

Next, partition

$$A \rightarrow \begin{pmatrix} A_{00} & A_{01} & \dots & A_{0,p-1} \\ A_{10} & A_{11} & \dots & A_{1,p-1} \\ \vdots & \vdots & \ddots & \vdots \\ A_{p-1,0} & A_{p-1,0} & \dots & A_{p-1,p-1} \end{pmatrix} \quad x \rightarrow \begin{pmatrix} x_0 \\ x_1 \\ \vdots \\ x_{p-1} \end{pmatrix}, \quad \text{and} \quad y \rightarrow \begin{pmatrix} y_0 \\ y_1 \\ \vdots \\ y_{p-1} \end{pmatrix},$$
where $A_{ij} \in \mathbb{R}^{n_i \times n_j}$ and $x_i, y_i \in \mathbb{R}^{n_i}$ with $\sum_{i=0}^{p-1} n_i = N$ and $n_i \approx N/\sqrt{P}$.

We will view the nodes as an $r \times c$ mesh, with $P = rc$, and index them as $p_{ij}$, with $i = 0, \dots, r-1$ and $j = 0, \dots, c-1$. Figure 6.4, for a $12 \times 12$ matrix on a $3 \times 4$ processor grid, illustrates the assignment of data to nodes, where the $i, j$ 'cell' shows the matrix and vector elements owned by $p_{ij}$.

<table><caption>Figure 6.4: Distribution of matrix and vector elements for a problem of size 12 on a $4 \times 3$ processor grid</caption><tbody><tr><td>x<sub>0</sub></td><td></td><td></td><td></td><td>x<sub>3</sub></td><td></td><td></td><td>x<sub>6</sub></td><td></td><td></td><td>x<sub>9</sub></td></tr><tr><td>a<sub>00</sub></td><td>a<sub>01</sub></td><td>a<sub>02</sub></td><td>y<sub>0</sub></td><td>a<sub>03</sub></td><td>a<sub>04</sub></td><td>a<sub>05</sub></td><td>a<sub>06</sub></td><td>a<sub>07</sub></td><td>a<sub>08</sub></td><td>a<sub>09</sub></td><td>a<sub>0,10</sub></td><td>a<sub>0,11</sub></td></tr><tr><td>a<sub>10</sub></td><td>a<sub>11</sub></td><td>a<sub>12</sub></td><td></td><td>a<sub>13</sub></td><td>a<sub>14</sub></td><td>a<sub>15</sub></td><td>a<sub>16</sub></td><td>a<sub>17</sub></td><td>a<sub>18</sub></td><td>a<sub>19</sub></td><td>a<sub>1,10</sub></td><td>a<sub>1,11</sub></td></tr><tr><td>a<sub>20</sub></td><td>a<sub>21</sub></td><td>a<sub>22</sub></td><td></td><td>a<sub>23</sub></td><td>a<sub>24</sub></td><td>a<sub>25</sub></td><td>a<sub>26</sub></td><td>a<sub>27</sub></td><td>a<sub>28</sub></td><td>y<sub>2</sub></td><td>a<sub>29</sub></td><td>a<sub>2,10</sub></td><td>a<sub>2,11</sub></td></tr><tr><td>a<sub>30</sub></td><td>a<sub>31</sub></td><td>a<sub>32</sub></td><td></td><td>a<sub>33</sub></td><td>a<sub>34</sub></td><td>a<sub>35</sub></td><td>a<sub>37</sub></td><td>a<sub>37</sub></td><td>a<sub>38</sub></td><td></td><td>a<sub>39</sub></td><td>a<sub>3,10</sub></td><td>a<sub>3,11</sub></td><td>y<sub>3</sub></td></tr><tr><td>x<sub>1</sub></td><td></td><td></td><td></td><td>x<sub>4</sub></td><td></td><td></td><td>x<sub>7</sub></td><td></td><td></td><td>x<sub>10</sub></td><td></td><td></td><td></td></tr><tr><td>a<sub>40</sub></td><td>a<sub>41</sub></td><td>a<sub>42</sub></td><td>y<sub>4</sub></td><td>a<sub>43</sub></td><td>a<sub>44</sub></td><td>a<sub>45</sub></td><td>a<sub>46</sub></td><td>a<sub>47</sub></td><td>a<sub>48</sub></td><td>a<sub>49</sub></td><td>a<sub>4,10</sub></td><td>a<sub>4,11</sub></td></tr><tr><td>a<sub>50</sub></td><td>a<sub>51</sub></td><td>a<sub>52</sub></td><td></td><td>a<sub>53</sub></td><td>a<sub>54</sub></td><td>a<sub>55</sub></td><td>a<sub>56</sub></td><td>a<sub>57</sub></td><td>a<sub>58</sub></td><td>a<sub>59</sub></td><td>a<sub>5,10</sub></td><td>a<sub>5,11</sub></td></tr><tr><td>a<sub>60</sub></td><td>a<sub>61</sub></td><td>a<sub>62</sub></td><td></td><td>a<sub>63</sub></td><td>a<sub>64</sub></td><td>a<sub>65</sub></td><td>a<sub>66</sub></td><td>a<sub>67</sub></td><td>a<sub>68</sub></td><td>y<sub>6</sub></td><td>a<sub>69</sub></td><td>a<sub>6,10</sub></td><td>a<sub>6,11</sub></td></tr><tr><td>a<sub>70</sub></td><td>a<sub>71</sub></td><td>a<sub>72</sub></td><td></td><td>a<sub>73</sub></td><td>a<sub>74</sub></td><td>a<sub>75</sub></td><td>a<sub>77</sub></td><td>a<sub>77</sub></td><td>a<sub>78</sub></td><td></td><td>a<sub>79</sub></td><td>a<sub>7,10</sub></td><td>a<sub>7,11</sub></td><td>y<sub>7</sub></td></tr><tr><td>x<sub>2</sub></td><td></td><td></td><td></td><td>x<sub>5</sub></td><td></td><td></td><td>x<sub>8</sub></td><td></td><td></td><td>x<sub>11</sub></td><td></td><td></td><td></td></tr><tr><td>a<sub>80</sub></td><td>a<sub>81</sub></td><td>a<sub>82</sub></td><td>y<sub>8</sub></td><td>a<sub>83</sub></td><td>a<sub>84</sub></td><td>a<sub>85</sub></td><td>a<sub>86</sub></td><td>a<sub>87</sub></td><td>a<sub>88</sub></td><td>a<sub>89</sub></td><td>a<sub>8,10</sub></td><td>a<sub>8,11</sub></td></tr><tr><td>a<sub>90</sub></td><td>a<sub>91</sub></td><td>a<sub>92</sub></td><td></td><td>a<sub>93</sub></td><td>a<sub>94</sub></td><td>a<sub>95</sub></td><td>a<sub>96</sub></td><td>a<sub>97</sub></td><td>a<sub>98</sub></td><td>a<sub>99</sub></td><td>a<sub>9,10</sub></td><td>a<sub>9,11</sub></td></tr><tr><td>a<sub>10,0</sub></td><td>a<sub>10,1</sub></td><td>a<sub>10,2</sub></td><td></td><td>a<sub>10,3</sub></td><td>a<sub>10,4</sub></td><td>a<sub>10,5</sub></td><td>a<sub>10,6</sub></td><td>a<sub>10,7</sub></td><td>a<sub>10,8</sub></td><td>y<sub>10</sub></td><td>a<sub>10,9</sub></td><td>a<sub>10,10</sub></td><td>a<sub>10,11</sub></td></tr><tr><td>a<sub>11,0</sub></td><td>a<sub>11,1</sub></td><td>a<sub>11,2</sub></td><td></td><td>a<sub>11,3</sub></td><td>a<sub>11,4</sub></td><td>a<sub>11,5</sub></td><td>a<sub>11,7</sub></td><td>a<sub>11,7</sub></td><td>a<sub>11,8</sub></td><td></td><td>a<sub>11,9</sub></td><td>a<sub>11,10</sub></td><td>a<sub>11,11</sub></td><td>y<sub>11</sub></td></tr></tbody></table>

In other words, $p_{ij}$ owns the matrix block $A_{ij}$ and parts of $x$ and $y$. This makes possible the following algorithm<fnref n="2" />:

* Since $x_j$ is distributed over the $j$th column, the algorithm starts by collecting $x_j$ on each processor $p_{ij}$ by an *allgather* inside the processor columns.

* Each processor $p_{ij}$ then computes $y_{ij} = A_{ij}x_j$. This involves no further communication.

* The result $y_i$ is then collected by gathering together the pieces $y_{ij}$ in each processor row to form $y_i$, and this is then distributed over the processor row. These two operations are in fact combined to form a *reduce-scatter*.

* If $r = c$, we can transpose the $y$ data over the processors, so that it can function as the input for a subsequent matrix-vector product. If, on the other hand, we are computing $A^tAx$, then $y$ is now correctly distributed for the $A^t$ product.

<footnote n="2">2. This figure shows a partitioning of the matrix into contiguous blocks, and the vector distribution seem to be what is necessary to work with this matrix distribution. You could also look at this story the other way: start with a distribution of input and output vector, and then decide what that implies for the matrix distribution. For instance, if you distributed $x$ and $y$ the same way, you would arrive at a different matrix distribution, but otherwise the product algorithm would be much the same; see [49].</footnote>
**Algorithm** The algorithm with cost analysis is

<table><thead><tr><th>Step</th><th>Cost (lower bound)</th></tr></thead><tbody><tr><td>Allgather x<sub>i</sub>'s within columns</td><td>[[log<sub>2</sub>(r)]α + <sup>r-1</sup>&frasl;<sub>p</sub>nβ<br>≈ log<sub>2</sub>(r)α + <sup>n</sup>&frasl;<sub>c</sub>β</td></tr><tr><td>Perform local matrix-vector multiply</td><td>≈ 2<sup>n<sup>2</sup></sup>&frasl;<sub>p</sub>γ</td></tr><tr><td>Reduce-scatter y<sub>i</sub>'s within rows</td><td>[[log<sub>2</sub>(c)]α + <sup>c-1</sup>&frasl;<sub>p</sub>nβ + <sup>c-1</sup>&frasl;<sub>p</sub>nγ<br>≈ log<sub>2</sub>(c)α + <sup>n</sup>&frasl;<sub>c</sub>β + <sup>n</sup>&frasl;<sub>c</sub>γ</td></tr></tbody></table>

**Cost analysis** The total cost of the algorithm is given by, approximately,

$$T_p^{r \times c}(n) = T_p^{r \times c}(n) = 2 \frac{n^2}{p} \gamma + \underbrace{\log_2(p) \alpha + \left(\frac{n}{c} + \frac{n}{r}\right) \beta + \frac{n}{r} \gamma}_{\text{Overhead}}$$

We will now make the simplification that $r = c = \sqrt{p}$ so that

$$T_p^{\sqrt{p} \times \sqrt{p}}(n) = T_p^{\sqrt{p} \times \sqrt{p}}(n) = 2 \frac{n^2}{p} \gamma + \underbrace{\log_2(p) \alpha + \frac{n}{\sqrt{p}} (2\beta + \gamma)}_{\text{Overhead}}$$

Since the sequential cost is $T_1(n) = 2n^2\gamma$, the speedup is given by

$$S_p^{\sqrt{p} \times \sqrt{p}}(n) = \frac{T_1(n)}{T_p^{\sqrt{p} \times \sqrt{p}}(n)} = \frac{2n^2\gamma}{2\frac{n^2}{p}\gamma + \frac{n}{\sqrt{p}}(2\beta + \gamma)} = \frac{p}{1 + \frac{p\log_2(p)}{2n^2}\frac{\alpha}{\gamma} + \frac{\sqrt{p}}{2n}\frac{(2\beta+\gamma)}{\gamma}}$$

and the parallel efficiency by

$$E_p^{\sqrt{p} \times \sqrt{p}}(n) = \frac{1}{1 + \frac{p \log_2(p)}{2n^2} \frac{\alpha}{\gamma} + \frac{\sqrt{p} (2\beta + \gamma)}{2n \gamma}}$$

We again ask the question what the parallel efficiency for the largest problem that can be stored on $p$ nodes is.

$$\begin{aligned} E_p^{\sqrt{p} \times \sqrt{p}}(n_{\max}(p)) &= \frac{1}{1 + \frac{p \log_2(p)}{2n^2} \frac{\alpha}{\gamma} + \frac{\sqrt{p} (2\beta+\gamma)}{2n \gamma}} \\ &= \frac{1}{1 + \frac{\log_2(p)}{2M} \frac{\alpha}{\gamma} + \frac{1}{2\sqrt{M}} \frac{(2\beta+\gamma)}{\gamma}} \end{aligned}$$

so that still

$$\lim_{p \to \infty} E_p^{\sqrt{p} \times \sqrt{p}}(n_{\max}(p)) = \lim_{p \to \infty} \frac{1}{1 + \frac{\log_2(p)}{2M} \frac{\alpha}{\gamma} + \frac{1}{2\sqrt{M}} \frac{(2\beta+\gamma)}{\gamma}} = 0.$$
However, $\log_2 p$ grows very slowly with $p$ and is therefore considered to act much like a constant. In this case $E_p^{\sqrt{p} \times \sqrt{p}}(n_{max}(p))$ decreases very slowly and the algorithm is considered to be scalable for practical purposes.

Note that when $r = p$ the 2D algorithm becomes the "partitioned by rows" algorithm and when $c = p$ it becomes the "partitioned by columns" algorithm. It is not hard to show that the 2D algorithm is scalable in the sense of the above analysis when $r = c$, as long as $r/c$ is kept constant.

## 6.3 LU factorization in parallel

The matrix-vector and matrix-product are easy to parallelize in one sense. The elements of the output can all be computed independently and in any order, so we have many degrees of freedom in parallelizing the algorithm. This does not hold for computing an LU factorization, or solving a linear system with the factored matrix.

### 6.3.1 Solving a triangular system

The solution of a triangular system $y = L^{-1}x$ (with $L$ is lower triangular) is a matrix-vector operation, so it has its $O(N^2)$ complexity in common with the matrix-vector product. However, unlike the product operation, this solution process contains a recurrence relation between the output elements:

$$y_i = \ell_{ii}^{-1} (x_i - \sum_{j<i} \ell_{ij} x_j).$$

This means that parallelization is not trivial. In the case of a sparse matrix special strategies may be possible; see section 6.10. Here we will make a few remarks about general, dense case.

Let us assume for simplicity that communication takes no time, and that all arithmetic operations take the same unit time. First we consider the matrix distribution by rows, meaning that processor $p$ stores the elements $\ell_{p*}$. With this we can implement the triangular solution as:

* Processor 1 solves $y_1 = \ell_{11}^{-1} x_1$ and sends its value to the next processor.

* In general, processor $p$ gets the values $y_1, \dots, y_{p-1}$ from processor $p-1$, and computes $y_p$;

* Each processor $p$ then sends $y_1, \dots, y_p$ to $p+1$.

**Exercise 6.3.** Show that this algorithm takes time $2N^2$, just like the sequential algorithm.

This algorithm has each processor passing all computed $y_i$ values to its successor, in a pipeline fashion. However, this means that processor $p$ receives $y_1$ only at the last moment, whereas that value was computed already in the first step. We can formulate the solution algorithm in such a way that computed elements are made available as soon as possible:

* Processor 1 solve $y_1$, and sends it to all later processors.

* In general, processor $p$ waits for individual messages with values $y_q$ for $q < p$.

* Processor $p$ then computes $y_p$ and sends it to processors $q$ with $q > p$.

Under the assumption that communication time is negligible, this algorithm can be much faster. For instance, all processors $p > 1$ receive $y_1$ simultaneously, and can compute $\ell_{p1} y_1$ simultaneously.
Exercise 6.4. Show that this algorithm variant takes a time $O(N)$.

Exercise 6.5. Now consider the matrix distribution by columns: processor $p$ stores $\ell_{*p}$. Outline the triangular solution algorithm with this distribution, and show that the parallel solve time is $O(N)$.

### 6.3.2 Factorization

A full analysis of the scalability of dense LU factorization is quite involved, so we will state without further proof that as in the matrix-vector case a two-dimensional distribution is needed. However, we can identify a further complication. Since factorizations of any type<fnref n="3" /> progress through a matrix, processors will be inactive for part of the time.

Exercise 6.6. Consider the regular right-looking Gaussian elimination

```
for k=1..n
    p = 1/a(k,k)
    for i=k+1,n
        for j=k+1,n
            a(i,j) = a(i,j) - a(i,k) * p * a(k,j)
```

Analyze the running time, speedup, and efficiency as a function of $N$, if we assume a one-dimensional distribution, and enough processors to store one column per processor. Show that speedup is limited.

Also perform this analysis for a two-dimensional decomposition where each processor stores one element.

For this reason, an *overdecomposition* is used, where the matrix is divided in more blocks than there are processors, and each processor stores several, non-contiguous, sub-matrices. We illustrate this in figure 6.5

<figure><img src="image_88.png" /><figcaption>Figure 6.5: One-dimensional cyclic distribution: assignment of four matrix columns to two processors, and the resulting mapping of storage to matrixcolumns</figcaption></figure>

<footnote n="3">3. Gaussian elimination can be performed in right-looking, left-looking and Crout variants; see [157].</footnote>
where we divide four block columns of a matrix to two processors: each processor stores in a contiguous
block of memory two non-contiguous matrix columns.

Next, we illustrate in figure 6.6 that a matrix-vector product with such a matrix can be performed without

<figure><img src="image_89.png" /><figcaption>Figure 6.6: Matrix-vector multiplication with a cyclically distributed matrix</figcaption></figure>

knowing that the processors store non-contiguous parts of the matrix. All that is needed is that the input
vector is also cyclicly distributed.

Exercise 6.7. Now consider a 4 × 4 matrix and a 2 × 2 processor grid. Distribute the matrix cyclically both by rows and columns. Show how the matrix-vector product can again be performed using the contiguous matrix storage, as long as the input is distributed correctly. How is the output distributed? Show that more communication is needed than the reduction of the one-dimensional example.

Specifically, with $P < N$ processors, and assuming for simplicity $N = cP$, we let processor 0 store rows 0, $c$, $2c$, $3c$, ...; processor 1 stores rows 1, $c+1$, $2c+1$, ... , et cetera. This scheme can be generalized two a two-dimensional distribution, if $N = c_1P_1 = c_2P_2$ and $P = P_1P_2$. This is called a 2D *cyclic distribution*. This scheme can be further extended by considering block rows and columns (with a small block size), and assigning to processor 0 the *block* rows 0, $c$, $2c$, ... .

Exercise 6.8. Consider a square $n \times n$ matrix, and a square $p \times p$ processor grid, where $p$ divides $n$ without remainder. Consider the overdecomposition outlined above, and make a sketch of matrix element assignment for the specific case $n = 6, p = 2$. That is, draw an $n \times n$ table where location $(i, j)$ contains the processor number that stores the corresponding matrix element. Also make a table for each of the processors describing the local to global mapping, that is, giving the global $(i, j)$ coordinates of the elements in the local matrix. (You will find this task facilitated by using zero-based numbering.)

Now write functions $P, Q, I, J$ of $i, j$ that describe the global to local mapping, that is, matrix element $a_{ij}$ is stored in location $(I(i, j), J(i, j))$ on processor $(P(i, j), Q(i, j))$.

## 6.4 Matrix-matrix product

In section 1.6.11 we considered the high performance implementation of the single-processor *matrix-matrix product*. We will now briefly consider the distributed memory version of this operation.
One algorithm for this operation is known as Cannon's algorithm. It assumes a square processor grid where processor (i, j) gradually accumulates the (i, j) block $C_{i,j} = \sum_k A_{i,k}B_{k,j}$; see figure 6.7.

If you start top-left, you see that processor (0, 0) has both $A_{00}$ and $B_{00}$, so it can immediately start doing a local multiplication. However, processor (0, 1) has $A_{01}$ and $B_{01}$, which are not needed together. However, if we rotate the second column of $B$ up by one position, processor (0, 1) will have $A_{01}$, $B_{11}$ and those two do need to be multiplied. Similarly, we rotate the third column of $B$ up by two places so that (0, 2) contains $A_{02}$, $B_{22}$.

How does this story go in the second row? Processor (1, 0) has $A_{10}$, $B_{10}$ which are not needed together. If we rotate the second row of $A$ one position to the left, it contains $A_{11}$, $B_{10}$ which are needed for a partial product. And now processor (1, 1) has $A_{11}$, $B_{11}$.

If we continue this story, we start with a matrix $A$ of which the rows have been rotated left, and $B$ of which the columns have been rotated up. In this setup, processor $(i, j)$ contains $A_{i,i+j}$ and $B_{i+j,j}$, where the addition is done modulo the matrix size. This means

that each processor can start by doing a local product.

Now we observe that $C_{ij} = \sum_k A_{ik}B_{kj}$ implies that the next partial product term comes from increasing $k$ by 1. The corresponding elements of $A, B$ can be moved into the processor by rotating the rows and columns by another location.

<figure><img src="image_90.png" /><figcaption>Figure 6.7: Cannon's algorithm for matrix-matrix multiplication: (a) initial rotation of matrix rows and columns, (b) resulting position in which processor $(i, j)$ can start accumulating $\sum_k A_{ik}B_{kj}$, (c) subsequent rotation of $A, B$ for the next term.</figcaption></figure>
**6.5 Sparse matrix-vector product**

In linear system solving through iterative methods (see section 5.5) the matrix-vector product is computa-
tionally an important kernel, since it is executed in each of potentially hundreds of iterations. In this section
we look at performance aspects of the matrix-vector product on a single processor first; the multi-processor
case will get our attention in section 6.5.

**6.5.1 The single-processor sparse matrix-vector product**

We are not much worried about the dense matrix-vector product in the context of iterative methods, since
one typically does not iterate on dense matrices. In the case that we are dealing with block matrices, refer
to section 1.6.10 for an analysis of the dense product. The sparse product is a lot trickier, since most of that
analysis does not apply.

**Data reuse in the sparse matrix-vector product**

There are some similarities between the dense matrix-vector product, executed by rows, and the CRS sparse product (section 5.4.1.4). In both cases all matrix elements are used sequentially, so any cache line loaded is utilized fully. However, the CRS product is worse at least the following ways:

* The indirect addressing requires loading the elements of an integer vector. This implies that the sparse product has more memory traffic for the same number of operations.

* The elements of the source vector are not loaded sequentially, in fact they can be loaded in effectively a random order. This means that a cacheline contains a source element will likely not be fully utilized. Also, the prefetch logic of the memory subsystem (section 1.3.5) cannot assist here.

For these reasons, an application that is computationally dominated by the sparse matrix-vector product can
very well be run at ≈ 5% of the peak performance of the processor.

It may be possible to improve this performance if the structure of the matrix is regular in some sense. One
such case is where we are dealing with a *block matrix* consisting completely of small dense blocks. This
leads at least to a reduction in the amount of indexing information: if the matrix consists of 2 × 2 blocks
we get a 4× reduction in the amount of integer data transferred.

Exercise 6.9. Give two more reasons why this strategy is likely to improve performance. Hint: cachelines, and reuse.

Such a *matrix tesselation* may give a factor of 2 in performance improvement. Assuming such an improvement, we may adopt this strategy even if the matrix is not a perfect block matrix: if every 2 × 2 block will contain one zero element we may still get a factor of 1.5 performance improvement [161, 20].

**Vectorization in the sparse product**

In other circumstances bandwidth and reuse are not the dominant concerns:

* On old vector computers, such as old Cray machines, memory was fast enough for the processor, but vectorization was paramount. This is a problem for sparse matrices, since the number of zeros in a matrix row, and therefore the vector length, is typically low.
* On GPUs memory bandwidth is fairly high, but it is necessary to find large numbers of identical operations. Matrix can be treated independently, but since the rows are likely of unequal length this is not an appropriate source of parallelism.

For these reasons, a variation on the diagonal storage scheme for sparse matrices has seen a revival recently. The observation here is that if you sort the matrix rows by the number of rows you get a small number of blocks of rows; each block will be fairly large, and in each block the rows have the same number of elements.

A matrix with such a structure is good for vector architectures [31]. In this case the product is computed by diagonals.

**Exercise 6.10.** Write pseudo-code for this case. How did the sorting of the rows improve the situation?

This sorted storage scheme also solves the problem we noted on GPUs [17]. In this case we the traditional CRS product algorithm, and we have an amount of parallelism equal to the number of rows in a block.

Of course there is the complication that we have permuted the matrix: the input and output vectors will need to be permuted accordingly. If the product operation is part of an iterative method, doing this permutation back and forth in each iteration will probably negate any performance gain. Instead we could permute the whole linear system and iterate on the permuted system.

**Exercise 6.11.** Can you think of reasons why this would work? Reasons why it wouldn't?

### 6.5.2 The parallel sparse matrix-vector product

In section 5.4 you saw a first discussion of sparse matrices, limited to use on a single processor. We will now go into parallelism aspects.

The dense matrix-vector product, as you saw above, required each processor to communicate with every other, and to have a local buffer of essentially the size of the global vector. In the sparse case, considerably less buffer space, as well as less communication, is needed. Let us analyze this case. We will assume that the matrix is distributed by block rows, where processor $p$ owns the matrix rows with indices in some set $I_p$.

The line $y_i = y_i + a_{ij}x_j$ now has to take into account that $a_{ij}$ can be zero. In particular, we need to consider that, for some pairs $i \in I_p, j \notin I_p$ no communication will be needed. Declaring for each $i \in I_p$ a sparsity pattern set

$$S_{p;i} = \{j: j \notin I_p, a_{ij} \neq 0\}$$

our multiplication instruction becomes

$$y_i = y_i + a_{ij}x_j \quad \text{if } j \in S_{p;i}.$$

If we want to avoid, as above, a flood of small messages, we combine all communication into a single message per processor. Defining

$$S_p = \cup_{i \in I_p} S_{p;i},$$

the algorithm now becomes:
<figure><img src="image_91.png" /><figcaption>Figure 6.8: A difference stencil applied to a two-dimensional square domain, distributed over processors. A cross-processor connection is indicated.</figcaption></figure>

<figure><img src="image_92.png" /><figcaption>Figure 6.9: The ghost region of a processor, induced by a stencil</figcaption></figure>

* Collect all necessary off-processor elements $x_j$ with $j \in S_p$ into one buffer;

* Perform the matrix-vector product, reading all elements of $x$ from local storage.

This whole analysis of course also applies to dense matrices. This becomes different if we consider where sparse matrices come from. Let us start with a simple case.

Recall figure 4.1, which illustrated a discretized boundary value problem on the simplest domain, a square, and let us now parallelize it. We do this by partitioning the domain; each processor gets the matrix rows corresponding to its subdomain. Figure 6.8 shows how this gives rise to connections between processors: the elements $a_{ij}$ with $i \in I_p, j \notin I_p$ are now the ‘legs’ of the stencil that reach beyond a processor boundary. The set of all such $j$, formally defined as

$$G = \{j \notin I_p : \exists i \in I_p : a_{ij} \neq 0\}$$

is known as the *ghost region* of a processor; see figure 6.9.

**Exercise 6.12.** Show that a one-dimensional partitioning of the domain leads to a partitioning of the matrix into block rows, but a two-dimensional partitioning of the domain does not. You can do this in the abstract, or you can illustrate it: take a 4×4 domain (giving a matrix of size 16), and partition it over 4 processors. The one-dimensional domain partitioning corresponds to giving each processor one line out of the domain, while
the two-dimensional partitioning gives each processor a 2 × 2 subdomain. Draw the matrices for these two cases.

Our crucial observation is now that, for each processor, the number of other processors it is involved with is strictly limited. This has implications for the efficiency of the operation.

### 6.5.3 Parallel efficiency of the sparse matrix-vector product

In the case of the dense matrix-vector product (section 6.2.2), partitioning the matrix over the processors by (block) rows did not lead to a scalable algorithm. Part of the reason was the increase in the number of neighbours that each processors needs to communicate with. Figure 6.8 shows that, for the matrix of a 5-five point stencil, this number is limited to four.

**Exercise 6.13.** Take a square domain and a partitioning of the variables of the processors as in figure 6.8. What is the maximum number of neighbours a processor needs to communication with for the box stencil in figure 4.3? In three space dimensions, what is the number of neighbours if a 7-point central difference stencil is used?

The observation that each processor communicates with just a few neighbours stays intact if we go beyond square domains to more complicated physical objects. If a processor receives a more or less contiguous subdomain, the number of its neighbours will be limited. This implies that even in complicated problems each processor will only communicate with a small number of other processors. Compare this to the dense case where each processor had to receive data from every other processor. It is obvious that the sparse case is far more friendly to the interconnection network. (The fact that it also more common for large systems may influence the choice of network to install if you are about to buy a new parallel computer.)

For square domains, this argument can easily be made formal. Let the unit domain $[0, 1]^2$ be partitioned over $P$ processors in a $\sqrt{P} \times \sqrt{P}$ grid. From figure 6.8 we see that every processor communicates with at most four neighbours. Let the amount of work per processor be $w$ and the communication time with each neighbour $c$. Then the time to perform the total work on a single processor is $T_1 = Pw$, and the parallel time is $T_P = w + 4c$, giving a speed up of

$$S_P = Pw/(w + 4c) = P/(1 + 4c/w) \approx P(1 - 4c/w).$$

**Exercise 6.14.** Express $c$ and $w$ as functions of $N$ and $P$, and show that the speedup is asymptotically optimal, under weak scaling of the problem.

**Exercise 6.15.** In this exercise you will analyze the parallel sparse matrix-vector product for a hypothetical, but realistic, parallel machine. Let the machine parameters be characterized by (see section 1.3.2):

* Latency: $\alpha = 1\mu s = 10^{-6}s$.

* Bandwidth: $1Gb/s$ corresponds to $\beta = 10^{-9}$.

* Computation rate: A per-core flops rate of $1G$ flops means $\gamma = 10^9$. This number may seem low, but note that the matrix-vector product has less reuse than the matrix-matrix product, which can achieve close to peak performance, and that the sparse matrix-vector product is even more bandwidth-bound.
We assume $10^4$ processors, and a five-point stencil matrix of size $N = 25 \cdot 10^{10}$.
This means each processor stores $5 \cdot 8 \cdot N/p = 10^9$ bytes. If the matrix comes
from a problem on a square domain, this means the domain was size $n \times n$ where
$n = \sqrt{N} = 5 \cdot 10^5$.

Case 1. Rather than dividing the matrix, we divide the domain, and we do this first by horizontal slabs of size $n \times (n/p)$. Argue that the communication complexity is $2(\alpha + n\beta)$ and computation complexity is $10 \cdot n \cdot (n/p)$. Show that the resulting computation outweighs the communication by a factor 250.

Case 2. We divide the domain into patches of size $(n/\sqrt{p}) \times (n/\sqrt{p})$. The memory and computation time are the same as before. Derive the communication time and show that it is better by a factor of 50.

Argue that the first case does not weakly scale: under the assumption that $N/p$ is
constant the efficiency will go down. (Show that speedup still goes up asymptotically
as $\sqrt{p}$.) Argue that the second case does scale weakly.

The argument that a processor will only connect with a few neighbours is based on the nature of the scientific computations. It is true for FDM and FEM methods. In the case of the *Boundary Element Method* (BEM), any subdomain needs to communicate with everything in a radius $r$ around it. As the number of processors goes up, the number of neighbours per processor will also go up.

**Exercise 6.16.** Give a formal analysis of the speedup and efficiency of the BEM algorithm.
Assume again a unit amount of work $w$ per processor and a time of communication $c$ per neighbour. Since the notion of neighbour is now based on physical distance, not on graph properties, the number of neighbours will go up. Give $T_1, T_p, S_p, E_p$ for this case.

There are also cases where a sparse matrix needs to be handled similarly to a dense matrix. For instance, *Google's PageRank* algorithm (see section 9.3) has at its heart the repeated operation $x \leftarrow Ax$ where $A$ is a sparse matrix with $A_{ij} \neq 0$ if web page $j$ links to page $i$; see section 9.3. This makes $A$ a very sparse matrix, with no obvious structure, so every processor will most likely communicate with almost every other.

**6.5.4 Memory behaviour of the sparse matrix-vector product**

In section 1.6.10 you saw an analysis of the sparse matrix-vector product in the dense case, on a single
processor. Some of the analysis carries over immediately to the sparse case, such as the fact that each
matrix element is used only once and that the performance is bound by the bandwidth between processor
and memory.

Regarding reuse of the input and the output vector, if the matrix is stored by rows, such as in CRS format (section 5.4.1.3), access to the output vector will be limited to one write per matrix row. On the other hand, the loop unrolling trick for getting reuse of the input vector can not be applied here. Code that combines two iterations is as follows:

```c
for (i=0; i<M; i+=2) {
    s1 = s2 = 0;
    for (j) {
        s1 = s1 + a[i][j] * x[j];
```
```c
s2 = s2 + a[i+1][j] * x[j];
}
y[i] = s1; y[i+1] = s2;
}
```

The problem here is that if $a_{ij}$ is nonzero, it is not guaranteed that $a_{i+1,j}$ is nonzero. The irregularity of the sparsity pattern makes optimizing the matrix-vector product hard. Modest improvements are possible by identifying parts of the matrix that are small dense blocks [20, 36, 160].

On a GPU the sparse matrix-vector product is also limited by memory bandwidth. Programming is now harder because the GPU has to work in data parallel mode, with many active threads.

An interesting optimization becomes possible if we consider the context in which the sparse matrix-vector product typically appears. The most common use of this operation is in iterative solution methods for linear systems (section 5.5), where it is applied with the same matrix in possibly hundreds of iterations. Thus we could consider leaving the matrix stored on the GPU and only copying the input and output vectors for each product operation.

### 6.5.5 The transpose product

In section 5.4.1.3 you saw that the code for both the regular and the transpose matrix-vector product are limited to loop orderings where rows of the matrix are traversed. (In section 1.5.2 you saw a discussion of computational effects of changes in loop order; in this case we are limited to row traversal by the storage format.)

In this section we will briefly look at the parallel transpose product. Equivalently to partitioning the matrix by rows and performing the transpose product, we look at a matrix stored and partitioned by columns and perform the regular product.

The algorithm for the product by columns can be given as:

$$ 
\begin{array}{l}
y \leftarrow 0 \\
\text{for } j: \\
\quad \text{for } i: \\
\quad \quad y_i \leftarrow y_i + a_{ij}x_j
\end{array}
$$

Both in shared and distributed memory we distribute outer iterations over processors. The problem is then that each outer iteration updates the whole output vector. This is a problem: with shared memory it leads to multiple writes to locations in the output and in distributed memory it requires communication that is as yet unclear.

One way to solve this would be to allocate a private output vector $y^{(p)}$ for each process:

$$ 
\begin{array}{l}
y^{(p)} \leftarrow 0 \\
\text{for } j \in \text{the rows of processor } p \\
\quad \text{for all } i: \\
\quad \quad y_i^{(p)} \leftarrow y_i^{(p)} + a_{ji}x_j
\end{array}
$$

after which we sum $y \leftarrow \sum_p y^{(p)}$.
6.5.6 Setup of the sparse matrix-vector product

While the dense matrix-vector product relies on collectives (see section 6.2), the sparse case uses point-to-
point communication^4 Every processor sends to just a few neighbours, and receives from just a few.

However, there is an asymmetry between sending and receiving. It is fairly easy for a processor to find out
what other processors it will be receiving from.

Exercise 6.17. Assume that the matrix is divided over the processors by block rows; see figure 6.2 for an illustration. Also assume that each processor knows which rows each other processor stores. (How would you implement that knowledge?) Sketch the algorithm by which a processor can find out who it will be receiving from; this algorithm should not involve any communication itself.

Discovering who to send *to* is harder.

Exercise 6.18. Argue that this is easy in the case of a structurally symmetric matrix: $a_{ij} \neq 0 \Leftrightarrow a_{ji} \neq 0$.

In the general case, a processor can in principle be asked to send to any other, so the simple algorithm is as
follows:

* Each processor makes an inventory of what non-local indices it needs. Under the above assumption that it knows what range of indices each other processor owns, it then decides which indices to get from what neighbours.

* Each processor sends a list of indices to each of its neighbours; this list will be empty for most of the neighbours, but we can not omit sending it.

* Each processor then receives these lists from all others, and draws up lists of which indices to send.

You will note that, even though the communication during the matrix-vector product involves only a few neighbours for each processor, giving a cost that is $O(1)$ in the number of processors, the setup involves all-to-all communications, which have time complexity $O(\alpha P)$

If a processor has only a few neighbours, the above algorithm is wasteful. Ideally, you would want space and running time proportional to the number of neighbours. We could bring the receive time in the setup if we knew how many messages were to be expected. That number can be found:

* Each processor makes an array *need* of length *P*, where *need[i]* is 1 if the processor needs any data from processor *i*, and zero otherwise.

* A *reduce-scatter* collective on this array, with a sum operator, then leaves on each processor a number indicating how many processors need data from it.

* The processor can execute that many receive calls.

The *reduce-scatter* call has time complexity $\alpha \log P + \beta P$, which is of the same order as the previous algorithm, but probably with a lower proportionality constant.

The time and space needed for the setup can be reduced to $O(\log P)$ with some sophisticated trickery [54, 89].

<footnote n="4">4. This makes sense for the type of sparse matrices that come from PDEs which have a clear structure, as you saw in section 4.2.3. However, there are sparse matrices that are so random that you essentially have to use dense techniques; see section 9.4.</footnote>
## 6.6 Computational aspects of iterative methods

All iterative methods feature the following operations:

* A matrix-vector product; this was discussed for the sequential case in section 5.4 and for the parallel case in section 6.5. In the parallel case, construction of FEM matrices has a complication that we will discuss in section 6.6.2.

* The construction of the preconditioner matrix $K \approx A$, and the solution of systems $Kx = y$. This was discussed in the sequential case in section 5.5.5. Below we will go into parallelism aspects in section 6.7.

* Some vector operations (including inner products, in general). These will be discussed next.

### 6.6.1 Vector operations

There are two types of vector operations in a typical iterative method: vector additions and inner products.

**Exercise 6.19.** Consider the CG method of section 5.5.10, figure 5.8, applied to the matrix from a 2D BVP; equation (4.14), First consider the unpreconditioned case $M = I$. Show that there is a roughly equal number of floating point operations are performed in the matrix-vector product and in the vector operations. Express everything in the matrix size $N$ and ignore lower order terms. How would this balance be if the matrix had 20 nonzeros per row?

Next, investigate this balance between vector and matrix operations for the FOM scheme in section 5.5.8. Since the number of vector operations depends on the iteration, consider the first 50 iterations and count how many floating point operations are done in the vector updates and inner product versus the matrix-vector product. How many nonzeros does the matrix need to have for these quantities to be equal?

**Exercise 6.20.** Flop counting is not the whole truth. What can you say about the efficiency of the vector and matrix operations in an iterative method, executed on a single processor?

#### 6.6.1.1 Vector additions

The vector additions are typically of the form $x \leftarrow x + \alpha y$ or $x \leftarrow \alpha x + y$. If we assume that all vectors are distributed the same way, this operation is fully parallel.

#### 6.6.1.2 Inner products

Inner products are vector operations, but they are computationally more interesting than updates, since they involve communication.

When we compute an inner product, most likely every processor needs to receive the computed value. We
use the following algorithm:

**for** processor *p* **do**

compute $a_p \leftarrow x_p^t y_p$ where $x_p, y_p$ are the part of $x, y$ stored on processor $p$

do a global reduction to compute $a = \sum_p a_p$

broadcast the result

**Algorithm 1:** Compute $a \leftarrow x^t y$ where $x, y$ are distributed vectors

The reduction and broadcast (which can be joined into an Allreduce) combine data over all processors, so they have a communication time that increases with the number of processors. This makes the inner product potentially an expensive operation, and people have suggested a number of ways to reducing their impact on the performance of iterative methods.

**Exercise 6.21.** Iterative methods are typically used for sparse matrices. In that context, you can argue that the communication involved in an inner product can have a larger influence on overall performance than the communication in the matrix-vector product. What is the complexity of the matrix-vector product and the inner product as a function of the number of processors?

Here are some of the approaches that have been taken.

* The CG method has two inner products per iteration that are inter-dependent. It is possible to rewrite the method so that it computes the same iterates (in exact arithmetic, at least) but so that the two inner products per iteration can be combined. See [27, 32, 124, 167].

* It may be possible to overlap the inner product calculation with other, parallel, calculations [34].

* In the GMRES method, use of the classical Gram-Schmidt (GS) method takes far fewer independent inner products than the modified GS method, but it is less stable. People have investigated strategies for deciding when it is allowed to use the classic GS method [109].

Since computer arithmetic is not associative, inner products are a prime source of results that differ when the same calculation is executed of two different processor configurations. In section 3.3.7 we sketched a solution.

## 6.6.2 Finite element matrix construction

The FEM leads to an interesting issue in parallel computing. For this we need to sketch the basic outline of how this method works. The FEM derives its name from the fact that the physical objects modeled are divided into small two or three dimensional shapes, the elements, such as triangles and squares in 2D, or pyramids and bricks in 3D. On each of these, the function we are modeling is then assumed to polynomial, often of a low degree, such as linear or bilinear.

The crucial fact is that a matrix element $a_{ij}$ is then the sum of computations, specifically certain integrals, over all elements that contain both variables $i$ and $j$:

$$a_{ij} = \sum_{e: i,j \in e} a_{ij}^{(e)}.$$

The computations in each element share many common parts, so it is natural to assign each element $e$ uniquely to a processor $P_e$, which then computes all contributions $a_{ij}^{(e)}$. In figure 6.10 element 2 is assigned to processor 0 and element 4 to processor 1.
<figure><img src="image_93.png" /><figcaption>Figure 6.10: A finite element domain, parallelization of the matrix construction, and parallelization of matrix element storage</figcaption></figure>

Now consider variables *i* and *j* and the matrix element $a_{ij}$. It is constructed as the sum of computations over domain elements 2 and 4, which have been assigned to different processors. Therefore, no matter what processor row *i* is assigned to, at least one processor will have to communicate its contribution to matrix element $a_{ij}$.

Clearly it is not possibly to make assignments $P_e$ of elements and $P_i$ of variables such that $P_e$ computes in full the coefficients $a_{ij}$ for all $i \in e$. In other words, if we compute the contributions locally, there needs to be some amount of communication to assemble certain matrix elements. For this reason, modern linear algebra libraries such as PETSc (see tutorial section 24) allow any processor to set any matrix element.

### 6.6.3 A simple model for iterative method performance

Above, we have already remarked that iterative methods have very little opportunity for data reuse, and they are therefore characterized as *bandwidth-bound* algorithms. This allows us to make a simple prediction as to the *floating point performance of iterative methods*. Since the number of iterations of an iterative method is hard to predict, by performance here we mean the performance of a single iteration. Unlike with *direct methods for linear systems* (see for instance section 6.8.1), the number of flops to solution is very hard to establish in advance.

First we argue that we can restrict ourselves to the *performance of the sparse matrix vector product*: the time spent in this operation is considerably more than in the vector operations. Additionally, most preconditioners have a computational structure that is quite similar to the matrix-vector product.

Let us then consider the *performance of the CRS matrix-vector product*.

* First we observe that the arrays of matrix elements and column indices have no reuse, so the performance of loading them is completely determined by the available bandwidth. Caches and prefect streams only hide the latency, but do not improve the bandwidth. For each multiplication of matrix element times input vector element we load one floating point number and one integer. Depending one whether the indices are 32-bit or 64-bit, this means 12 or 16 bytes loaded for each multiplication.
* The demands for storing the result vector are less important: an output vector element is written only once for each matrix row.

* The input vector can also be ignored in the bandwidth calculation. At first sight you might think that indirect indexing into the input vector is more or less random and therefore expensive. However, let's take the operator view of a matrix-vector product and consider the space domain of the PDE from which the matrix stems; see section 4.2.3. We now see that the seemingly random indexing is in fact into vector elements that are grouped closely together. This means that these vector elements will likely reside in L3 cache, and therefore accessible with a higher bandwidth (say, by a factor of at least 5) than data from main memory.

For the parallel performance of the sparse matrix-vector product we consider that in a PDE context each processor communicates only with a few neighbours. Furthermore, a surface-to-volume argument shows that the message volume is of a lower order than the on-node computation.

In sum, we conclude that a very simple model for the sparse matrix vector product, and thereby for the whole iterative solver, consists of measuring the effective bandwidth and computing the performance as one addition and one multiplication operation per 12 or 16 bytes loaded.

## 6.7 Parallel preconditioners

Above (section 5.5.5 and in particular 5.5.5.1) we saw a couple of different choices of $K$. In this section we will begin the discussion of parallelization strategies. The discussion is continued in detail in the next sections.

### 6.7.1 Jacobi preconditioning

The Jacobi method (section 5.5.3) uses the diagonal of $A$ as preconditioner. Applying this is as parallel as is possible: the statement $y \leftarrow K^{-1}x$ scales every element of the input vector independently. Unfortunately the improvement in the number of iterations with a Jacobi preconditioner is rather limited. Therefore we need to consider more sophisticated methods such ILU. Unlike with the Jacobi preconditioner, parallelism is then not trivial.

### 6.7.2 The trouble with parallel ILU

Above we saw that, in a flop counting sense, applying an ILU preconditioner (section 5.5.5.1) is about as expensive as doing a matrix-vector product. This is no longer true if we run our iterative methods on a parallel computer.

At first glance the operations are similar. A matrix-vector product $y = Ax$ looks like

```fortran
for i=1..n
  y[i] = sum over j=1..n a[i,j]*x[j]
```

In parallel this would look like
```
for i=myfirstrow..mylastrow
  y[i] = sum over j=1..n a[i,j]*x[j]
```

Suppose that a processor has local copies of all the elements of A and x that it will need, then this operation is fully parallel: each processor can immediately start working, and if the work load is roughly equal, they will all finish at the same time. The total time for the matrix-vector product is then divided by the number of processors, making the speedup more or less perfect.

Consider now the forward solve $Lx = y$, for instance in the context of an ILU preconditioner:

```
for i=1..n
  x[i] = (y[i] - sum over j=1..i-1 ell[i,j]*x[j]) / a[i,i]
```

We can simply write the parallel code:

```
for i=myfirstrow..mylastrow
  x[i] = (y[i] - sum over j=1..i-1 ell[i,j]*x[j]) / a[i,i]
```

but now there is a problem. We can no longer say 'suppose a processor has local copies of everything in the right hand side', since the vector x appears both in the left and right hand side. While the matrix-vector product is in principle fully parallel over the matrix rows, this triangular solve code is recursive, hence sequential.

In a parallel computing context this means that, for the second processor to start, it needs to wait for certain components of x that the first processor computes. Apparently, the second processor can not start until the first one is finished, the third processor has to wait for the second, and so on. The disappointing conclusion is that in parallel only one processor will be active at any time, and the total time is the same as for the sequential algorithm. This is actually not a big problem in the dense matrix case, since parallelism can be found in the operations for handling a single row (see section 6.12), but in the sparse case it means we can not use incomplete factorizations without some redesign.

In the next few subsections we will see different strategies for finding preconditioners that perform efficiently in parallel.

### 6.7.3 Block Jacobi methods

Various approaches have been suggested to remedy this sequentiality the triangular solve. For instance, we could simply let the processors ignore the components of x that should come from other processors:

```
for i=myfirstrow..mylastrow
  x[i] = (y[i] - sum over j=myfirstrow..i-1 ell[i,j]*x[j]) / a[i,i]
```

This is not mathematically equivalent to the sequential algorithm (technically, it is called a *block Jacobi method* with ILU as the *local solve*), but since we're only looking for an approximation $K \approx A$, this is simply a slightly cruder approximation.
Exercise 6.22. Take the Gauss-Seidel code you wrote above, and simulate a parallel run. What is the effect of increasing the (simulated) number of processors?

The idea behind block methods can easily be appreciated pictorially; see figure 6.11. In effect, we make an

<figure><img src="image_94.png" /><figcaption>Figure 6.11: Sparsity pattern corresponding to a block Jacobi preconditioner</figcaption></figure>

ILU of the matrix that we get by ignoring all connections between processors. Since in a BVP all points influence each other (see section 4.2.1), using a less connected preconditioner will increase the number of iterations if executed on a sequential computer. However, block methods are parallel and, as we observed above, a sequential preconditioner is very inefficient in a parallel context, so we put up with this increase in iterations.

### 6.7.4 Parallel ILU

The Block Jacobi preconditioner operates by decoupling domain parts. While this may give a method that is highly parallel, it may give a higher number of iterations than a true ILU preconditioner. (A theoretical argument can be made that this decoupling decreases the efficiency of the iterative method; see section 4.2.1.) Fortunately it is possible to have a parallel ILU method.

In section 6.8.2 you saw the combination of *graph colouring* and *permutation*. Let $P$ be the permutation that groups like-coloured variables together, then $\tilde{A} = P^t A P$ is a matrix with the following structure:

* $\tilde{A}$ has a block structure with the number of blocks equal to the number of colours in the adjacency graph of $A$; and

* each diagonal block is a diagonal matrix.

Now, if you are performing an iterative system solution and you are looking for a parallel preconditioner you can use this permuted matrix. Consider solving $Ly = x$ with the permuted system. We write the usual algorithm (section 5.3.5) as

for $c$ in the set of colours:

for $i$ in the variables of colour $c$:

$$y_i \leftarrow x_i - \sum_{j<i} l_{ij} y_j$$

Exercise 6.23. Show that the *flop count* of solving a system $LUx = y$ remains the same (in the highest order term) when you from an ILU factorization in the natural ordering to one in the colour-permuted ordering.
<figure><img src="image_95.png" /><figcaption>Figure 6.12: A partitioned domain with coloured nodes</figcaption></figure>

Where does all this colouring get us? Solving is still sequential... Well, it is true that the outer loop over the colours is sequential, but all the points of one colour are independent of each other, so they can be solved at the same time. So if we use an ordinary domain partitioning and combine that with a multi-colouring (see figure 6.12), the processors are all active during all the colour stages; see figure 6.13. Ok, if you take a close look at that figure you'll see that one processor is not active in the last colour. With large numbers of nodes per processor this is unlikely to happen, but there may be some load imbalance.

## 6.8 Ordering strategies and parallelism

In the foregoing we have remarked on the fact that solving a linear system of equations is inherently a recursive activity. For dense systems, the number of operations is large enough compared to the recursion length that finding parallelism is fairly straightforward. Sparse systems, on the other hand, take more sophistication. In this section we will look at a number of strategies for reordering the equations (or, equivalently, permuting the matrix) that will increase the available parallelism.

These strategies can all be considered as variants of Gaussian elimination. By making incomplete variants of them (see section 5.5.5.1), all these strategies also apply to constructing preconditioners for iterative solution methods.

### 6.8.1 Nested dissection

Above, you have seen several examples of ordering the variables in the domain other than with the *lexico-graphic ordering*. In this section you will see the *nested dissection ordering*, which was initially designed as a way to reduce fill-in. However, it is also advantageous in a parallel computing context.

Nested dissection is a recursive process for determining a nontrivial ordering of the unknowns in a domain. In the first step, the computational domain is split in two parts, with a dividing strip between them; see figure 6.14. To be precise, the *separator* is wide enough that there are no connections between the left and
<figure><img src="image_96.png" /><figcaption>Figure 6.13: Solving a parallel multicolour ILU in four steps</figcaption></figure>

right subdomain. The resulting matrix $A^{DD}$ has a 3 × 3 structure, corresponding to the three divisions of the domain. Since the subdomains $\Omega_1$ and $\Omega_2$ are not connected, the submatrices $A_{12}^{DD}$ and $A_{21}^{DD}$ are zero.

$$A^{DD} = \left( \begin{array}{ccc} A_{11} & \emptyset & A_{13} \\ \emptyset & A_{22} & A_{23} \\ A_{31} & A_{32} & A_{33} \end{array} \right) = \left\{ \begin{array}{c|c} \begin{array}{cccc} \star & \star & & \\ \star & \star & \star & \\ \vdots & \vdots & \vdots & \\ \star & \star & \star & \\ \star & \star & \star & \star \end{array} & \emptyset \\ \hline \begin{array}{c} \emptyset \\ \hline \emptyset \end{array} & \begin{array}{c} \begin{array}{cccc} \star & \star & & \\ \star & \star & \star & \\ \vdots & \vdots & \vdots & \\ \star & \star & \star & \\ \star & \star & \star & \star \end{array} \\ \begin{array}{cccc} 0 & & & \\ \vdots & \vdots & \vdots & \\ 0 & & & \\ \star & & & \end{array} \end{array} & \begin{array}{c} (n^2 - n)/2 \\ (n^2 - n)/2 \\ n \end{array} \end{array} \right\}$$

This process of dividing the domain by a separator is also called *domain decomposition or substructuring*, although this name is also associated with the mathematical analysis of the resulting matrices [10]. In this example of a rectangular domain it is of course trivial to find a separator. However, for the type of equations we get from BVPs it is usually feasible to find a separator efficiently for any domain [115]; see also section 16.5.1.

Let us now consider the *LU* factorization of this matrix. If we factor it in terms of the 3 × 3 block structure, we get

$$A^{DD} = LU = \begin{pmatrix} I & & \\ \emptyset & I & \\ A_{31}A_{11}^{-1} & A_{32}A_{22}^{-1} & I \end{pmatrix} \begin{pmatrix} A_{11} & \emptyset & A_{13} \\ A_{22} & A_{23} & \\ S_{33} & & \end{pmatrix}$$
<figure><img src="image_97.png" /><figcaption>Figure 6.14: Domain dissection into two unconnected subdomains and a separator</figcaption></figure>

where

$$S_{33} = A_{33} - A_{31}A_{11}^{-1}A_{13} - A_{32}A_{22}^{-1}A_{23}.$$

The important fact here is that

* the contributions $A_{31}A_{11}^{-1}A_{13}$ and $A_{32}A_{22}^{-1}A_{23}$ can be computed simultaneously, so the factorization is largely parallel; and

* both in the forward and backward solve, components 1 and 2 of the solution can be computed simultaneously, so the solution process is also largely parallel.

The third block can not trivially be handled in parallel, so this introduces a sequential component in the algorithm. We also need to take a closer look at the structure of $S_{33}$.

**Exercise 6.24.** In section 5.4.3.1 you saw the connection between LU factorization and graph theory: eliminating a node leads to a graph with that node removed, but with certain new connections added. Show that, after eliminating the first two sets of variables, the graph of the remaining matrix on the separator will be fully connected.

The upshot is that after eliminating all the variables in blocks 1 and 2 we are left with a matrix $S_{33}$ that is fully dense of size $n \times n$.

The introduction of a separator gave us a factorization that was two-way parallel. Now we iterate this process: we put a separator inside blocks 1 and 2 (see figure 6.15), which gives the following matrix structure:
<figure><img src="image_98.png" /><figcaption>Figure 6.15: A four-way domain decomposition</figcaption></figure>

$$A^{DD} = \begin{pmatrix} A_{11} & & & A_{15} & & A_{17} \\ & A_{22} & & A_{25} & & A_{27} \\ & & A_{33} & & A_{36} & A_{37} \\ & & & A_{44} & & A_{46} & A_{47} \\ & & & & A_{55} & & A_{57} \\ A_{51} & A_{52} & & A_{63} & A_{64} & A_{66} & A_{67} \\ & & & & A_{75} & A_{76} & A_{77} \end{pmatrix}$$

(Note the similarities with the 'arrow' matrix in section 5.4.3.4, and recall the argument that this led to lower fill-in.) The LU factorization of this is:

$$\left( \begin{array}{c|c|c} I & I & I \\ \hline A_{51}A_{11}^{-1} & A_{52}A_{22}^{-1} & A_{63}A_{33}^{-1} & A_{64}A_{44}^{-1} & I & I \\ \hline A_{71}A_{11}^{-1} & A_{72}A_{22}^{-1} & A_{73}A_{33}^{-1} & A_{74}A_{44}^{-1} & A_{75}S_{5}^{-1} & A_{76}S_{6}^{-1} & I \end{array} \right)$$

$$\left( \begin{array}{c|c|c} A_{11} & A_{15} & A_{17} \\ & A_{22} & A_{27} \\ & & A_{33} & A_{37} \\ & & & A_{44} & A_{46} & A_{47} \\ & & & & S_{5} & A_{57} \\ & & & & & S_{6} & A_{67} \\ & & & & & & S_{7} \end{array} \right)$$
where

$$ 
\begin{aligned}
S_5 &= A_{55} - A_{51}A_{11}^{-1}A_{15} - A_{52}A_{22}^{-1}A_{25}, & S_6 &= A_{66} - A_{63}A_{33}^{-1}A_{36} - A_{64}A_{44}^{-1}A_{46} \\
S_7 &= A_{77} - \sum_{i=1,2,3,4} A_{7i}A_{ii}^{-1}A_{i7} - \sum_{i=5,6} A_{7i}S_i^{-1}A_{17}.
\end{aligned}
$$

Constructing the factorization now goes as follows:

* Blocks $A_{ii}$ are factored in parallel for $i = 1, 2, 3, 4$; similarly the contributions $A_{5i}A_{ii}^{-1}A_{i5}$ for $i = 1, 2, A_{6i}A_{ii}^{-1}A_{i6}$ for $i = 3, 4$, and $A_{7i}A_{ii}^{-1}A_{i7}$ for $i = 1, 2, 3, 4$ can be constructed in parallel.

* The Schur complement matrices $S_5, S_6$ are formed and subsequently factored in parallel, and the contributions $A_{7i}S_i^{-1}A_{17}$ for $i = 5, 6$ are constructed in parallel.

* The Schur complement $S_7$ is formed and factored.

Analogous to the above reasoning, we conclude that after eliminating blocks 1,2,3,4 the updated matrices $S_5, S_6$ are dense of size $n/2$, and after eliminating blocks 5,6 the Schur complement $S_7$ is dense of size $n$.

**Exercise 6.25.** Show that solving a system with $A^{DD}$ has a similar parallelism to constructing the factorization as described above.

For future reference, we will call the sets 1 and 2 each others' siblings, and similarly for 3 and 4. The set 5 is the parent of 1 and 2, 6 is the parent of 3 and 4; 5 and 6 are siblings and 7 is the parent of 5 and 6.

### 6.8.1.1 Domain decomposition

In figure 6.15 we divided the domain four ways by a recursive process. This leads up to our discussion of nested dissection. It is also possible to immediately split a domain in any number of strips, or in a grid of subdomains. As long as the separators are wide enough, this will give a matrix structure with many independent subdomains. As in the above discussion, an LU factorization will be characterized by

* parallel processing of the subdomains, both in the factorization and $L, U$ solves, and

* a system to be solved on the separator structure.

**Exercise 6.26.** The matrix from a two-dimensional BVP has a block tridiagonal structure. Divide the domain in four strips, that is, using three separators (see figure 6.16). Note that the separators are uncoupled in the original matrix.

Now sketch the sparsity structure of the resulting system on the separators are elimination of the subdomains. Show that the system is block tridiagonal.

### 6.8.1.2 Complexity

The nested dissection method repeats the above process until the subdomains are very small. For a theoretical analysis, we keep dividing until we have subdomains of size $1 \times 1$, but in practice one could stop at sizes such as 32, and use an efficient dense solver to factor and invert the blocks.

To derive the complexity of the algorithm, we take another look at figure 6.15, and see that complexity argument, the total space a full recursive nested dissection factorization needs is the sum of

* one dense matrix on a separator of size $n$, plus

* two dense matrices on separators of size $n/2$,
<figure><img src="image_99.png" /><figcaption>Figure 6.16: One-way domain decomposition</figcaption></figure>

* taking together $3/2 n^2$ space and $5/12 n^3$ time;

* the two terms above then get repeated on four subdomains of size $(n/2) \times (n/2)$.

With the observation that $n = \sqrt{N}$, this sums to

$$ 
\begin{aligned}
\text{space} &= 3/2n^2 + 4 \cdot 3/2(n/2)^2 + \dots \\
&= N(3/2 + 3/2 + \dots) \quad \log n \text{ terms} \\
&= O(N \log N)
\end{aligned}
$$

$$ 
\begin{aligned}
\text{time} &= 5/12n^3/3 + 4 \cdot 5/12(n/2)^3/3 + \dots \\
&= 5/12N^{3/2}(1 + 1/4 + 1/16 + \dots) \\
&= O(N^{3/2})
\end{aligned}
$$

Apparently, we now have a factorization that is parallel to a large extent, and that is done in $O(N \log N)$ space, rather than $O(N^{3/2})$ (see section 5.4.3.3). The factorization time has also gone down from $O(N^2)$ to $O(N^{3/2})$.

Unfortunately, this space savings only happens in two dimensions: in three dimensions we need

* one separator of size $n \times n$, taking $(n \times n)^2 = N^{4/3}$ space and $1/3 \cdot (n \times n)^3 = 1/3 \cdot N^2$ time,

* two separators of size $n \times n/2$, taking $N^{3/2}/2$ space and $1/3 \cdot N^2/4$ time,

* four separators of size $n/2 \times n/2$, taking $N^{3/2}/4$ space and $1/3 \cdot N^2/16$ time,

* adding up to $7/4N^{3/2}$ space and $21/16N^2/3$ time;
* on the next level there are 8 subdomains that contribute these terms with $n \to n/2$ and therefore $N \to N/8$.

This makes the total space

$$ \frac{7}{4}N^{3/2}(1 + (1/8)^{4/3} + \cdots) = O(N^{3/2}) $$

and the total time

$$ \frac{21}{16}N^2(1 + 1/16 + \cdots)/3 = O(N^2). $$

We no longer have the tremendous savings of the 2D case. A much more complicated analysis shows that the order improvement holds for general problems in 2D, and that 3D in general has a higher complexity [114].

### 6.8.1.3 Parallelism

The nested dissection method clearly introduces a lot of parallelism, and we can characterize it as task parallelism (section 2.6.3): associated with each separator is a task of factoring its matrix, and later one of solving a linear system on its variables. However, the tasks are not independent: in figure 6.15 the factorization on domain 7 has to wait for 5 and 6, and they have to wait for 1,2,3,4. Thus, we have tasks with dependencies in the form of a tree: each separator matrix can be factored only when its children have been factored.

Mapping these tasks to processors is not trivial. First of all, if we are dealing with shared memory we can use a simple task queue:

Queue $\leftarrow \{\}$
**for** all bottom level subdomains $d$ **do**
    add $d$ to the Queue
**while** Queue is not empty **do**
    **if** a processor is idle **then**
        assign a queued task to it
    **if** a task is finished AND its sibling is finished **then**
        add its parent to the queue

The main problem here is that at some point we will have more processors than tasks, thus causing load unbalance. This problem is made more severe by the fact that the last tasks are also the most substantial, since the separators double in size from level to level. (Recall that the work of factoring a dense matrix goes up with the third power of the size!) Thus, for the larger separators we have to switch from task parallelism to medium-grained parallelism, where processors collaborate on factoring a block.

With distributed memory, we can now solve the parallelism problem with a simple task queue, since it would involve moving large amounts of data. (But recall that work is a higher power of the matrix size, which this time works in our favour, making communication relatively cheap.) The solution is then to use some form of domain decomposition. In figure 6.15 we could have four processors, associated with block 1,2,3,4. Processors 1 and 2 would then negotiate which one factors block 5 (and similarly processors 3 and 4 and block 6), or they could both do it redundantly.
<figure><img src="image_100.png" /><figcaption>Figure 6.17: Red-black ordering of the points on a line</figcaption></figure>

### 6.8.1.4 Preconditioning

As with all factorizations, it is possible to turn the nested dissection method into a preconditioner by making the factorization incomplete. (For the basic idea of incomplete factorizations, see section 5.5.5.1). However, here the factorization is formulated completely in terms of *block matrices*, and the division by the pivot element becomes an inversion or system solution with the pivot block matrix. We will not go into this further; for details see the literature [4, 50, 125].

## 6.8.2 Variable reordering and colouring: independent sets

Another permutation of the problem variables is based on graph colouring (section 16.3). The direct goal here is to maximize available parallelism.

Let us take a simple example, where $A$ is a tridiagonal matrix. The equation $Ax = b$ looks like

$$ \begin{pmatrix} a_{11} & a_{12} & & & \\ a_{21} & a_{22} & a_{23} & & \\ & a_{32} & a_{33} & a_{34} & \\ & & \ddots & \ddots & \\ & & & \ddots & \\ & & & & \emptyset \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \\ \vdots \end{pmatrix} = \begin{pmatrix} y_1 \\ y_2 \\ y_3 \\ \vdots \end{pmatrix} $$

We observe that $x_i$ directly depends on $x_{i-1}$ and $x_{i+1}$, but not $x_{i-2}$ or $x_{i+1}$. Thus, let us see what happens if we permute the indices to group every other component together.

Pictorially, we take the points 1, ..., n and colour them red and black (figure 6.17), then we permute them to first take all red points, and subsequently all black ones. The correspondingly permuted matrix looks as follows:

$$ \begin{pmatrix} a_{11} & & & & a_{12} & & \\ & a_{33} & & & a_{32} & a_{34} & \\ & & a_{55} & & & \ddots & \\ & & & \ddots & & & \\ a_{21} & a_{23} & & & a_{22} & & \\ & a_{43} & a_{45} & & a_{44} & & \\ & & \ddots & & & \ddots & \end{pmatrix} \begin{pmatrix} x_1 \\ x_3 \\ x_5 \\ \vdots \\ x_2 \\ x_4 \\ \vdots \end{pmatrix} = \begin{pmatrix} y_1 \\ y_3 \\ y_5 \\ \vdots \\ y_2 \\ y_4 \\ \vdots \end{pmatrix} $$
With this permuted A, the Gauss-Seidel matrix $D_A + L_A$ looks like

<figure><img src="image_101.png" /></figure>

What does this buy us? Well, let's spell out the solution of a system $Lx = y$.

for $i = 1, 3, 5, \dots$ do

$$ \text{solve } x_i \leftarrow y_i / a_{ii} $$

for $i = 2, 4, 6, \dots$ do

$$ \text{compute } t = a_{ii-1}x_{i-1} + a_{ii+1}x_{i+1} $$

$$ \text{solve } x_i \leftarrow (y_i - t) / a_{ii} $$

Apparently the algorithm has three stages that are each parallel over half the domain points. This is illustrated in figure 6.18. Theoretically we could accommodate a number of processors that is half the number of

<figure><img src="image_102.png" /><figcaption>Figure 6.18: Red-black solution on a 1d domain</figcaption></figure>

the domain points, but in practice each processor will have a subdomain. Now you can see in figure 6.19 how this causes a very modest amount of communication: each processor sends at most the data of two red points to its neighbours.

Red-black ordering can be applied to two-dimensional problems too. Let us apply a red-black ordering to the points $(i, j)$ where $1 \le i, j \le n$. Here we first apply a successive numbering to the odd points on the first line (1, 1), (3, 1), (5, 1), ..., then the even points of the second line (2, 2), (4, 2), (6, 2), ..., the odd points on the third line, etcetera. Having thus numbered half the points in the domain, we continue with the even points in the first line, the odd points in the second, etcetera. As you can see in figure 6.20, now the red points are only connected to black points, and the other way around. In graph theoretical terms, you have found a *colouring* (see appendix 16 for the definition of this concept) of the matrix graph with two colours.
<figure><img src="image_103.png" /><figcaption>Figure 6.19: Parallel red-black solution on a 1d domain</figcaption></figure>

<figure><img src="image_104.png" /><figcaption>Figure 6.20: Red-black ordering of the variables of a two-dimensional domain</figcaption></figure>

Exercise 6.27. Apply the red-black ordering to the 2D BVP (4.10). Sketch the resulting matrix structure.

The red-black ordering is a simple example of *graph colouring* (sometimes called *multi-colouring*). In simple cases, such as the unit square domain we considered in section 4.2.3 or its extension to 3D, the *colour number* of the adjacency graph is easily determined.

Exercise 6.28. You saw that a red-black ordering of unknowns coupled with the regular five-point star stencil give two subsets of variables that are not connected among themselves, that is, they form a two-colouring of the matrix graph. Can you find a colouring if nodes are connected by the second stencil in figure 4.3?

There is a simple bound for the number of colours needed for the graph of a sparse matrix: the number of colours is at most $d + 1$ where $d$ is the degree of the graph. To see that we can colour a graph with degree $d$ using $d + 1$ colours, consider a node with degree $d$. No matter how its neighbours are coloured, there is always an unused colour among the $d + 1$ available ones.
Exercise 6.29. Consider a sparse matrix, where the graph can be coloured with *d* colours. Permute the matrix by first enumerating the unknowns of the first colour, then the second colour, et cetera. What can you say about the sparsity pattern of the resulting permuted matrix?

If you are looking for a direct solution of the linear system you can repeat the process of colouring and permuting on the matrix that remains after you have eliminated one colour. In the case of a tridiagonal matrix you saw that this remaining matrix was again tridiagonal, so it is clear how to continue the process. This is called recursive doubling. If the matrix is not tridiagonal but block tridiagonal, this operation can be performed on blocks.

## 6.9 Operator splitting

In some contexts, it is necessary to perform implicit calculations through all directions of a two or three-dimensional array. For example, in section 4.3 you saw how the implicit solution of the heat equation gave rise to repeated systems

$$ (\alpha I + \frac{d^2}{dx^2} + \frac{d^2}{dy^2})u^{(t+1)} = u^{(t)} \qquad (6.2) $$

Without proof, we state that the time-dependent problem can also be solved by

$$ (\beta I + \frac{d^2}{dx^2})(\beta I + \frac{d^2}{dy^2})u^{(t+1)} = u^{(t)} \qquad (6.3) $$

for suitable $\beta$. This scheme will not compute the same values on each individual time step, but it will converge to the same steady state. The scheme can also be used as a preconditioner in the BVP case.

This approach has considerable advantages, mostly in terms of operation counts: the original system has to be solved either making a factorization of the matrix, which incurs fill-in, or by solving it iteratively.

Exercise 6.30. Analyze the relative merits of these approaches, giving rough operation counts. Consider both the case where $\alpha$ has dependence on $t$ and where it does not. Also discuss the expected speed of various operations.

A further advantage appears when we consider the parallel solution of (6.3). Note that we have a two-dimensional set of variables $u_{ij}$, but the operator $I + d^2u/dx^2$ only connects $u_{ij}$, $u_{ij-1}$, $u_{ij+1}$. That is, each line corresponding to an $i$ value can be processed independently. Thus, both operators can be solved fully parallel using a one-dimensional partition on the domain. The solution of a the system in (6.2), on the other hand, has limited parallelism.

Unfortunately, there is a serious complication: the operator in $x$ direction needs a partitioning of the domain in on direction, and the operator in $y$ in the other. The solution usually taken is to transpose the $u_{ij}$ value matrix in between the two solves, so that the same processor decomposition can handle both. This transposition can take a substantial amount of the processing time of each time step.

Exercise 6.31. Discuss the merits of and problems with a two-dimensional decomposition of the domain, using a grid of $P = p \times p$ processors. Can you suggest a way to ameliorate the problems?
One way to speed up these calculations, is to replace the implicit solve, by an explicit operation; see section 6.10.3.

## 6.10 Parallelism and implicit operations

In the discussion of IBVPs (section 4.1.2.2) you saw that implicit operations can have great advantages from the point of numerical stability. However, you also saw that they make the difference between methods based on a simple operation such as the matrix-vector product, and ones based on the more complicated linear system solution. There are further problems with implicit methods when you start computing in parallel.

Exercise 6.32. Let A be the matrix

$$A = \begin{pmatrix} a_{11} & & & \emptyset \\ a_{21} & a_{22} & & \\ & \ddots & \ddots & \\ & & a_{n,n-1} & a_{nn} \end{pmatrix}. \qquad (6.4)$$

Show that the matrix vector product $y \leftarrow Ax$ and the system solution $x \leftarrow A^{-1}y$, obtained by solving the triangular system $Ax = y$, not by inverting A, have the same operation count.

Now consider parallelizing the product $y \leftarrow Ax$. Suppose we have n processors, and each processor i stores $x_i$ and the i-th row of A. Show that the product $Ax$ can be computed without idle time on any processor but the first.

Can the same be done for the solution of the triangular system $Ax = y$? Show that the straightforward implementation has every processor idle for an $(n-1)/n$ fraction of the computation.

We will now see a number of ways of dealing with this inherently sequential component.

### 6.10.1 Wavefronts

Above, you saw that solving a lower triangular system of size N can have sequential time complexity of N steps. In practice, things are often not quite that bad. Implicit algorithms such as solving a triangular system are inherently sequential, but the number of steps can be less than is apparent at first.

Exercise 6.33. Take another look at the matrix from a two-dimensional BVP on the unit square, discretized with central differences. Derive the matrix structure if we order the unknowns by diagonals. What can you say about the sizes of the blocks and the structure of the blocks themselves?

Let us take another look at figure 4.1 that describes the difference stencil of a two-dimensional BVP. The corresponding picture for the lower triangular factor is in figure 6.21. This describes the sequentiality of the lower triangular solve process $x \leftarrow L^{-1}y$:

$$x_k = y_k - \ell_{k,k-1}x_{k-1} - \ell_{k,k-n}x_{k-n}$$
<figure><img src="image_105.png" /><figcaption>Figure 6.21: The difference stencil of the *L* factor of the matrix of a two-dimensional BVP</figcaption></figure>

In other words, the value at point *k* can be found if its neighbours to the left (that is, variable *k* − 1) and below (variable *k* − *n*) are known.

Turning this around, we see that, if we know $x_1$, we can not only find $x_2$, but also $x_{n+1}$. In the next step we can determine $x_3$, $x_{n+2}$, and $x_{2n+1}$. Continuing this way, we can solve *x* by *wavefronts*: the values of *x* on each wavefront are independent, so they can be solved in parallel in the same sequential step.

**Exercise 6.34.** Finish this argument. What is the maximum number of processors we can employ, and what is the number of sequential steps? What is the resulting efficiency?

Of course you don't have to use actual parallel processing to exploit this parallelism. Instead you could use
a *vector processor*, *vector instructions*, or a *GPU* [117].

In section 5.4.3.5 you saw the *Cuthill-McKee ordering* for reducing the fill-in of a matrix. We can modify this algorithm as follows to give wavefronts:

1. Take an arbitrary node, and call that 'level zero'.

2. For level *n* + 1, find points connected to level *n*, that are not themselves connected.

3. For the so-called 'reverse Cuthill-McKee ordering', reverse the numbering of the levels.

**Exercise 6.35.** This algorithm is not entirely correct. What is the problem; how can you correct it? Show that the resulting permutated matrix is no longer tridiagonal, but will likely still have a band structure.

**6.10.2 Recursive doubling**

One strategy for dealing with recurrences is *recursive doubling*, which you already saw in exercise 1.4.
Here we will discuss it in a more systematic manner. First, take the matrix from (6.4) and scale it to be of
the form

$$ \begin{pmatrix} 1 & & & \emptyset \\ b_{21} & 1 & & \\ & \ddots & \ddots & \\ & & b_{n,n-1} & 1 \end{pmatrix} $$

which we write as $A = I + B$.

**Exercise 6.36.** How does solving the system $(I+B)x = y$ help in solving $Ax = y$? What are the operation counts of solving the system in the two different ways?

Now we do something that looks like Gaussian elimination, except that we do not start with the first row, but the second. (What would happen if you did Gaussian elimination or LU decomposition on the matrix $I+B$?) We use the second row to eliminate $b_{32}$:

$$ \begin{pmatrix} 1 & & & \emptyset \\ & b_{21} & 1 & \\ & & \ddots & \\ & & & b_{n,n-1} \end{pmatrix} \times \begin{pmatrix} 1 & & & \emptyset \\ & b_{32} & 1 & \\ & & \ddots & \\ & & & b_{n,n-1} \end{pmatrix} = \begin{pmatrix} 1 & & & \emptyset \\ & b_{21} & 1 & \\ & & \ddots & \\ & & & b_{n,n-1} \end{pmatrix} $$

which we write as $L^{(2)}A = A^{(2)}$. We also compute $L^{(2)}y = y^{(2)}$ so that $A^{(2)}x = y^{(2)}$ has the same solution as $Ax = y$. Solving the transformed system gains us a little: after we compute $x_1, x_2$ and $x_3$ can be computed in parallel.

Now we repeat this elimination process by using the fourth row to eliminate $b_{54}$, the sixth row to eliminate $b_{76}$, et cetera. The final result is, summarizing all $L^{(i)}$ matrices:

$$ \begin{pmatrix} 1 & & & \emptyset \\ & b_{21} & 1 & \\ & & \ddots & \\ & & & b_{n,n-1} \end{pmatrix} \times (I+B) = \begin{pmatrix} 1 & & & \emptyset \\ & b_{32}b_{21} & 1 & \\ & & \ddots & \\ & & & b_{n,n-1} \end{pmatrix} $$

which we write as $L(I + B) = C$, and solving $(I + B)x = y$ now becomes $Cx = L^{-1}y$.

* First of all, computing $y' = L^{-1}y$ is simple. (Work out the details. How much parallelism is available?)

* Solving $Cx = y'$ is still sequential, but it no longer takes $n$ steps: from $x_1$ we can get $x_3$, from that we get $x_5$, et cetera. In other words, there is only a sequential relationship between the odd numbered components of $x$.

* The even numbered components of $x$ do not depend on each other, but only on the odd components: $x_2$ follows from $x_1, x_4$ from $x_3$, et cetera. Once the odd components have been computed, admittedly sequentially, this step is fully parallel.
We can describe the sequential solving of the odd components by itself:

$$ \begin{pmatrix} 1 & & & \emptyset \\ c_{21} & 1 & & \\ & \ddots & \ddots & \\ & & c_{n,n-1} & 1 \end{pmatrix} \begin{pmatrix} x_1 \\ x_3 \\ \vdots \\ x_n \end{pmatrix} = \begin{pmatrix} y'_1 \\ y'_3 \\ \vdots \\ y'_n \end{pmatrix} $$

where $c_{i+1i} = -b_{2n+1,2n}b_{2n,2n-1}$. In other words, we have reduced a size $n$ sequential problem to a sequential problem of the size kind and a parallel problem, both of size $n/2$. Now we can repeat this procedure recursively, reducing the original problem to a sequence of parallel operations, each half the size of the former.

The process of computing all partial sums through recursive doubling is also referred to as a parallel *prefix operation*. Here we use a prefix sum, but in the abstract it can be applied to any associative operator.

### 6.10.3 Approximating implicit by explicit operations, series expansion

There are various reasons why it is sometimes allowed to replace an implicit operation, which, as you saw above, can be problematic in practice, by a different one that is practically more advantageous.

* Using an explicit method for the heat equation (section 4.3) instead of an implicit one is equally legitimate, as long as we observe step size restrictions on the explicit method.

* Tinkering with the preconditioner (section 5.5.7) in an iterative method is allowed, since it will only affect the speed of convergence, not the solution the method converges to. You already saw one example of this general idea in the *block Jacobi method*; section 6.7.3. In the rest of this section you will see how recurrences in the preconditioner, which are implicit operations, can be replaced by explicit operations, giving various computational advantages.

Solving a linear system is a good example of an implicit operation, and since this comes down to solving two triangular systems, let us look at ways of finding a computational alternative to solving a lower triangular system. If $U$ is upper triangular and nonsingular, we let $D$ be the diagonal of $U$, and we write $U = D(I - B)$ where $B$ is an upper triangular matrix with a zero diagonal, also called a *strictly upper triangular matrix*; we say that $I - B$ is a *unit upper triangular matrix*.

**Exercise 6.37.** Let $A = LU$ be an LU factorization where $L$ has ones on the diagonal. Show how solving a system $Ax = b$ can be done, involving only the solution of unit upper and lower triangular systems. Show that no divisions are needed during the system solution.

Our operation of interest is now solving the system $(I - B)x = y$. We observe that

$$ (I - B)^{-1} = I + B + B^2 + \dots \qquad (6.5) $$

and $B^n = 0$ where $n$ is the matrix size (check this!), so we can solve $(I - B)x = y$ exactly by

$$ x = \sum_{k=0}^{n-1} B^k y. $$
Of course, we want to avoid computing the powers $B^k$ explicitly, so we observe that

$$ \sum_{k=0}^{1} B^k y = (I+B)y, \quad \sum_{k=0}^{2} B^k y = (I+B(I+B))y, \quad \sum_{k=0}^{3} B^k y = (I+B(I+B((I+B))))y, \quad (6.6) $$

et cetera. The resulting algorithm for evaluating $\sum_{k=0}^{n-1} B^k y$ is called Horner's rule, and you see that it avoids computing matrix powers $B^k$.

**Exercise 6.38.** Suppose that $I-B$ is bidiagonal. Show that the above calculation takes $n(n+1)$ operations. What is the operation count for computing $(I-B)x = y$ by triangular solution?

We have now turned an implicit operation into an explicit one, but unfortunately one with a high operation count. In practical circumstances, however, we can truncate the sum of matrix powers.

**Exercise 6.39.** Let $A$ be the tridiagonal matrix

$$ A = \begin{pmatrix} 2 & -1 & & & & \emptyset \\ -1 & 2 & -1 & & & \\ & \ddots & \ddots & \ddots & & \\ & & \emptyset & -1 & & \\ & & & -1 & 2 & \\ \emptyset & & & & & \emptyset \end{pmatrix} $$

of the one-dimensional BVP from section 4.2.2.

1. Recall the definition of diagonal dominance in section 5.3.4. Is this matrix diagonally dominant?

2. Show that the pivots in an LU factorization of this matrix (without pivoting) satisfy a recurrence. Hint: show that after $n$ elimination steps ($n \ge 0$) the remaining matrix looks like

$$ A^{(n)} = \begin{pmatrix} d_n & -1 & & & & \emptyset \\ -1 & 2 & -1 & & & \\ & \ddots & \ddots & \ddots & & \\ & & \emptyset & -1 & & \\ & & & -1 & 2 & \\ \emptyset & & & & & \emptyset \end{pmatrix} $$

and show the relation between $d_{n+1}$ and $d_n$.

3. Show that the sequence $n \mapsto d_n$ is descending, and derive the limit value.

4. Write out the $L$ and $U$ factors in terms of the $d_n$ pivots.

5. Are the $L$ and $U$ factors diagonally dominant?

The above exercise implies (note that we did not actually prove it!) that for matrices from BVPs we find that $B^k \downarrow 0$, in element size and in norm. This means that we can approximate the solution of $(I-B)x = y$ by, for instance, $x = (I+B)y$ or $x = (I+B+B^2)y$. Doing this still has a higher operation count than the direct triangular solution, but it is computationally advantageous in at least two ways:

* The explicit algorithm has a better pipeline behaviour.
* The implicit algorithm has problems in parallel, as you have seen; the explicit algorithm is easily parallelized.

Of course, this approximation may have further implications for the stability of the overall numerical algorithm.

Exercise 6.40. Describe the parallelism aspects of Horner's rule; equation (6.6).

## 6.11 Grid updates

One of the conclusions of chapter 4 was that explicit methods for time-dependent problems are computationally easier than implicit ones. For instance, they typically involve a matrix-vector product rather than a system solution, and parallelizing explicit operations is fairly simple: each result value of the matrix-vector product can be computed independently. That does not mean that there are other computational aspects worth remarking on.

Since we are dealing with sparse matrices, stemming from some computational stencil, we take the operator point of view. In figures 6.8 and 6.9 you saw how applying a stencil in each point of the domain induces certain relations between processors: in order to evaluate the matrix-vector product $y \leftarrow Ax$ on a processor, that processor needs to obtain the $x$-values of its ghost region. Under reasonable assumptions on the partitioning of the domain over the processors, the number of messages involved will be fairly small.

Exercise 6.41. Reason that, in a FEM or FDM context, the number of messages is $O(1)$ as $h \downarrow 0$.

In section 1.5.1 you saw that the matrix-vector product has little data reuse, though there is some locality to the computation; in section 5.4.1.4 it was pointed out that the locality of the sparse matrix-vector product is even worse because of indexing schemes that the sparsity necessitates. This means that the sparse product is largely a *bandwidth-bound algorithm*.

Looking at just a single product there is not much we can do about that. However, often we do a number of such products in a row, for instance as the steps in a time-dependent process. In that case there may be rearrangements of the operations that lessen the bandwidth demands. Consider as a simple example

$$ \forall i: x_i^{(n+1)} = f(x_i^{(n)}, x_{i-1}^{(n)}, x_{i+1}^{(n)}) \qquad (6.7) $$

and let's assume that the set $\{x_i^{(n)}\}_i$ is too large to fit in cache. This is a model for, for instance, the explicit scheme for the heat equation in one space dimension; section 4.3.1.1. Schematically:

<figure><img src="image_106.png" /></figure>

In the ordinary computation, where we first compute all $x_i^{(n+1)}$, then all $x_i^{(n+2)}$, the intermediate values at level $n+1$ will be flushed from the cache after they were generated, and then brought back into cache as input for the level $n+2$ quantities.
However, if we compute not one, but two iterations, the intermediate values may stay in cache. Consider $x_0^{(n+2)}$: it requires $x_0^{(n+1)}$, $x_1^{(n+1)}$, which in turn require $x_0^{(n)}$, ..., $x_2^{(n)}$.

Now suppose that we are not interested in the intermediate results, but only the final iteration. Figure 6.22 shows a simple example. The first processor computes 4 points on level $n + 2$. For this it needs 5 points

<figure><img src="image_107.png" /><figcaption>Figure 6.22: Computation of blocks of grid points over multiple iterations</figcaption></figure>

from level $n + 1$, and these need to be computed too, from 6 points on level $n$. We see that a processor apparently needs to collect a *ghost region* of width two, as opposed to just one for the regular single step update. One of the points computed by the first processor is $x_3^{(n+2)}$, which needs $x_4^{(n+1)}$. This point is also needed for the computation of $x_4^{(n+2)}$, which belongs to the second processor.

The easiest solution is to let this sort of point on the intermediate level *redundantly computed*, in the computation of both blocks where it is needed, on two different processors.

**Exercise 6.42.** Can you think of cases where a point would be redundantly computed by more than two processors?

We can give several interpretations to this scheme of computing multiple update steps by blocks.

* First of all, as we motivated above, doing this on a single processor increases locality: if all points in a coloured block (see the figure) fit in cache, we get reuse of the intermediate points.

* Secondly, if we consider this as a scheme for distributed memory computation, it reduces message traffic. Normally, for every update step the processors need to exchange their boundary data. If we accept some redundant duplication of work, we can now eliminate the data exchange for the intermediate levels. The decrease in communication will typically outweigh the increase in work.

**Exercise 6.43.** Discuss the case of using this strategy for multicore computation. What are the savings? What are the potential pitfalls?

### 6.11.1 Analysis

Let's analyze the algorithm we have just sketched. As in equation (6.7) we limit ourselves to a 1D set of points and a function of three points. The parameters describing the problem are these:

* $N$ is the number of points to be updated, and $M$ denotes the number of update steps. Thus, we perform $MN$ function evaluations.
* $\alpha, \beta, \gamma$ are the usual parameters describing latency, transmission time of a single point, and time for an operation (here taken to be an $f$ evaluation).

* $b$ is the number of steps we block together.

Each halo communication consists of $b$ points, and we do this $\sqrt{N}/b$ many times. The work performed consists of the $MN/p$ local updates, plus the redundant work because of the halo. The latter term consists of $b^2/2$ operations, performed both on the left and right side of the processor domain.

Adding all these terms together, we find a cost of

$$ \frac{M}{b}\alpha + M\beta + \left(\frac{MN}{p} + Mb\right)\gamma. $$

We observe that the overhead of $\alpha M/b + \gamma Mb$ is independent of $p$,

**Exercise 6.44.** Compute the optimal value of $b$, and remark that it only depends on the architectural parameters $\alpha, \beta, \gamma$ but not on the problem parameters.

### 6.11.2 Communication and work minimizing strategy

We can make this algorithm more efficient by overlapping computation with communication. As illustrated in figure 6.23, each processor start by communicating its halo, and overlapping this communication with the part of the communication that can be done locally. The values that depend on the halo will then be computed last.

<figure><img src="image_108.png" /><figcaption>Figure 6.23: Computation of blocks of grid points over multiple iterations</figcaption></figure>

**Exercise 6.45.** What is a great practical problem with organizing your code (with the emphasis on ‘code’!) this way?

If the number of points per processor is large enough, the amount of communication is low relative to the computation, and you could take $b$ fairly large. However, these grid updates are mostly used in iterative methods such as the CG method (section 5.5.10), and in that case considerations of roundoff prevent you from taking $b$ too large[27].

**Exercise 6.46.** Go through the complexity analysis for the non-overlapping algorithm in case the points are organized in a 2D grid. Assume that each point update involves four neighbours, two in each coordinate direction.
<figure><img src="image_109.png" /><figcaption>Figure 6.24: Computation of blocks of grid points over multiple iterations</figcaption></figure>

A further refinement of the above algorithm is possible. Figure 6.24 illustrates that it is possible to use a halo region that uses different points from different time steps. This algorithm (see [35]) cuts down on the amount of redundant computation. However, now the halo values that are communicated first need to be computed, so this requires splitting the local communication into two phases.

## 6.12 Block algorithms on multicore architectures

In section 5.3.7 you saw that certain linear algebra algorithms can be formulated in terms of submatrices. This point of view can be beneficial for the efficient execution of linear algebra operations on shared memory architectures such as current multicore processors.

As an example, let us consider the Cholesky factorization, which computes $A = LL^t$ for a symmetric positive definite matrix $A$; see also section 5.3.2. Recursively, we can describe the algorithm as follows:

$$ \mathrm{Chol} \begin{pmatrix} A_{11} & A_{21}^t \\ A_{21} & A_{22} \end{pmatrix} = LL^t \quad \text{where} \quad L = \begin{pmatrix} L_{11} & 0 \\ \tilde{A}_{21} & \mathrm{Chol}(A_{22} - \tilde{A}_{21}\tilde{A}_{21}^t) \end{pmatrix} $$

and where $\tilde{A}_{21} = A_{21}L_{11}^{-t}$, $A_{11} = L_{11}L_{11}^t$.

In practice, the block implementation is applied to a partitioning

$$ \left( \begin{array}{c|c} \text{finished} & \\ \hline A_{kk} & A_{k,>k} \\ \hline A_{>k,k} & A_{>k,>k} \end{array} \right) $$

where $k$ is the index of the current block row, and the factorization is finished for all indices $< k$. The factorization is written as follows, using Blas names for the operations:

for $k = 1, n_{\text{blocks}}$:

*   **Chol**: factor $L_k L_k^t \leftarrow A_{kk}$

*   **Trsm**: solve $\tilde{A}_{>k,k} \leftarrow A_{>k,k} L_k^{-t}$

*   **Gemm**: form the product $\tilde{A}_{>k,k} \tilde{A}_{>k,k}^t$

*   **Syrk**: symmetric rank-k update $A_{>k,k} \leftarrow A_{>k,k} - \tilde{A}_{>k,k} \tilde{A}_{>k,k}^t$
The key to parallel performance is to partition the indices > k and write the algorithm in terms of these
blocks:

$$
\left(
\begin{array}{c|ccc}
 & \multicolumn{3}{c}{\text{finished}} \\
\hline
A_{kk} & A_{k,k+1} & A_{k,k+2} \cdots \\
A_{k+1,k} & A_{k+1,k+1} & A_{k+1,k+2} \cdots \\
A_{k+2,k} & A_{k+2,k+2} & \\
\vdots & \vdots & \\
\end{array}
\right)
$$

The algorithm now gets an extra level of inner loops:

for $k = 1, n$ blocks:

Chol: factor $L_k L_k^t \leftarrow A_{kk}$

for $\ell > k$:

Trsm: solve $\tilde{A}_{\ell,k} \leftarrow A_{\ell,k} L_k^{-t}$

for $\ell_1, \ell_2 > k$:

Gemm: form the product $\tilde{A}_{\ell_1,k} \tilde{A}_{\ell_2,k}^t$

for $\ell_1, \ell_2 > k, \ell_1 \le \ell_2$:

Syrk: symmmetric rank-k update $A_{\ell_1, \ell_2} \leftarrow A_{\ell_1, \ell_2} - \tilde{A}_{\ell_1, k} \tilde{A}_{\ell_2, k}^t$

Now it is clear that the algorithm has a good deal of parallelism: the iterations in every *l*-loop can be
processed independently. However, these loops get shorter in every iteration of the outer *k*-loop, so it is not
immediate how many processors we can accomodate. Moreover, it is not necessary to preserve the order of
operations of the algorithm above. For instance, after

$$
L_1 L_1^t = A_{11}, \quad A_{21} \leftarrow A_{21} L_1^{-t}, \quad A_{22} \leftarrow A_{22} - A_{21} A_{21}^t
$$

the factorization $L_2 L_2^t = A_{22}$ can start, even if the rest of the $k = 1$ iteration is still unfinished. Thus, there
is probably a lot more parallelism than we would get from just parallelizing the inner loops.

The best way to approach parallelism in this case is to shift away from a *control flow* view of the algorithm,
where the sequence of operations is prescribed, to a *data flow* view. In the latter only data dependencies are
indicated, and any ordering of operations that obeys these dependencies is allowed. (Technically, we aban-
don the program order of the tasks and replace it with a *partial ordering*<fnref n="5" />.) The best way of representing
the data flow of an algorithm is by constructing a *Directed Acyclic Graph (DAG)* (see section 16 for a brief
tutorial on graphs) of tasks. We add an edge (*i*, *j*) to the graph if task *j* uses the output of task *i*.

Exercise 6.47. In section 2.7.1.5 you learned the concept of sequential consistency: a threaded
parallel code program should give the same results when executed in parallel as when
it's executed sequentially. We have just stated that DAG-based algorithms are free to
execute tasks in any order that obeys the partial order of the graph nodes. Discuss
whether sequential consistency is a problem in this context.

In our example, we construct a DAG by making a vertex task for every inner iteration. Figure 6.25 shows the
DAG of all tasks of matrix of 4 × 4 blocks. This graph is constructed by simulating the Cholesky algorithm
above,

<footnote n="5">5. Let's write $a \le b$ if $a$ is executed before $b$, then the relation $\le$ is a partial order if $a \le b \land b \le a \Rightarrow a = b$ and $a \le b \land b \le c \Rightarrow a \le c$. The difference with a total ordering, such as program ordering, is that it is not true that $a \le b \lor b \le a$: there can be pairs that are not ordered, meaning that their time ordering is not prescribed.</footnote>
Exercise 6.48. What is the diameter of this graph? Identify the tasks that lie on the path that determines the diameter. What is the meaning of these tasks in the context of the algorithm? This path is called the *critical path*. Its length determines the execution time of the computation in parallel, even if an infinite number of processors is available.

Exercise 6.49. Assume there are *T* tasks that all take a unit time to execute, and assume we have *p* processors. What is the theoretical minimum time to execute the algorithm? Now amend this formula to take into account the critical path; call its length *C*.

In the execution of the tasks a DAG, several observations can be made.

* If more than one update is made to a block, it is probably advantageous to have these updates be computed by the same process. This simplifies maintaining *cache coherence*.

* If data is used and later modified, the use must be finished before the modification can start. This can even be true if the two actions are on different processors, since the memory subsystem typically maintains cache coherence, so the modifications can affect the process that is reading the data. This case can be remedied by having a copy of the data in main memory, giving a reading process data that is reserved (see section 1.4.1).
<figure><img src="image_110.png" /><figcaption>Figure 6.25: Graph of task dependencies in a 4 × 4 Cholesky factorization</figcaption></figure>
6. High performance linear algebra
# PART II

## APPLICATIONS
# Chapter 7

## Molecular dynamics

Molecular dynamics is a technique for simulating the atom-by-atom behavior of molecules and deriving macroscopic properties from these atomistic motions. It has application to biological molecules such as proteins and nucleic acids, as well as natural and synthetic molecules in materials science and nanotechnology. Molecular dynamics falls in the category of particle methods, which includes N-body problems in celestial mechanics and astrophysics, and many of the ideas presented here will carry over to these other fields. In addition, there are special cases of molecular dynamics including ab initio molecular dynamics where electrons are treated quantum mechanically and thus chemical reactions can be modeled. We will not treat these special cases, but will instead concentrate on *classical* molecular dynamics.

The idea behind molecular dynamics is very simple: a set of particles interact according to Newton's law of motion, $F = ma$. Given the initial particle positions and velocities, the particle masses and other parameters, as well as a model of the forces that act between particles, Newton's law of motion can be integrated numerically to give a trajectory for each of the particles for all future (and past) time. Commonly, the particles reside in a computational box with periodic boundary conditions.

A molecular dynamics time step is thus composed of two parts:

1. compute forces on all particles

2. update positions (integration).

The computation of the forces is the expensive part. State-of-the-art molecular dynamics simulations are performed on parallel computers because the force computation is costly and a vast number of time steps are required for reasonable simulation lengths. In many cases, molecular dynamics is applied to simulations on molecules with a very large number of atoms as well, e.g., up to a million for biological molecules and long time scales, and up to billions for other molecules and shorter time scales.

Numerical integration techniques are also of interest in molecular dynamics. For simulations that take a large number of time steps and for which the preservation of quantities such as energy is more important than order of accuracy, the solvers that must be used are different than the traditional ODE solvers presented in Chapter 4.

In the following, we will introduce force fields used for biomolecular simulations and discuss fast methods for computing these forces. Then we devote sections to the parallelization of molecular dynamics for short-range forces and the parallelization of the 3-D FFT used in fast computations of long-range forces. We
end with a section introducing the class of integration techniques that are suitable for molecular dynamics simulations. Our treatment of the subject of molecular dynamics in this chapter is meant to be introductory and practical; for more information, the text [59] is recommended.

# 7.1 Force Computation

## 7.1.1 Force Fields

In classical molecular dynamics, the model of potential energy and of the forces that act between atoms is called a *force field*. The force field is a tractable but approximate model of quantum mechanical effects which are computationally too expensive to determine for large molecules. Different force fields are used for different types of molecules, as well as for the same molecule by different researchers, and none are ideal.

In biochemical systems, commonly-used force fields model the potential energy function as the sum of bonded, van der Waals, and electrostatic (Coulomb) energy:

$$E = E_{\text{bonded}} + E_{\text{Coul}} + E_{\text{vdW}}.$$

The potential is a function of the positions of all the atoms in the simulation. The force on an atom is the negative gradient of this potential at the position of the atom.

The bonded energy is due to covalent bonds in a molecule,

$$E_{\text{bonded}} = \sum_{\text{bonds}} k_i (r_i - r_{i,0})^2 + \sum_{\text{angles}} k_i (\theta_i - \theta_{i,0})^2 + \sum_{\text{torsions}} V_n (1 + \cos(n\omega - \gamma))$$

where the three terms are, respectively, sums over all covalent bonds, sums over all angles formed by two bonds, and sums over all dihedral angles formed by three bonds. The fixed parameters $k_i$, $r_{i,0}$, etc. depend on the types of atoms involved, and may differ for different force fields. Additional terms or terms with different functional forms are also commonly used.

The remaining two terms for the potential energy $E$ are collectively called the nonbonded terms. Computationally, they form the bulk of the force calculation. The electrostatic energy is due to atomic charges and is modeled by the familiar

$$E_{\text{Coul}} = \sum_{i} \sum_{j>i} \frac{q_i q_j}{4\pi \epsilon_0 r_{ij}}$$

where the sum is over all pairs of atoms, $q_i$ and $q_j$ are the charges on atoms $i$ and $j$, and $r_{ij}$ is the distance between atoms $i$ and $j$. Finally, the van der Waals energy approximates the remaining attractive and repulsive effects, and is commonly modeled by the Lennard-Jones function

$$E_{\text{vdW}} = \sum_{i} \sum_{j>i} 4\epsilon_{ij} \left[ \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{12} - \left( \frac{\sigma_{ij}}{r_{ij}} \right)^{6} \right]$$
where $\epsilon_{ij}$ and $\sigma_{ij}$ are force field parameters depending on atom types. At short distances, the repulsive ($r^{12}$) term is in effect, while at long distances, the dispersive (attractive, $-r^6$) term is in effect.

Parallelization of the molecular dynamics force calculation depends on parallelization each of these individual types of force calculations. The bonded forces are local computations in the sense that for a given atom, only nearby atom positions and data are needed. The van der Waals forces are also local and are termed short-range because they are negligible for large atom separations. The electrostatic forces are long-range, and various techniques have been developed to speed up these calculations. In the next two subsections, we separately discuss the computation of short-range and long-range nonbonded forces.

## 7.1.2 Computing Short-Range Nonbonded Forces

The computation of short-range nonbonded forces for a particle can be truncated beyond a cutoff radius, $r_c$, of that particle. The naive approach to perform this computation for a particle $i$ is by examining all other particles and computing their distance to particle $i$. For $n$ particles, the complexity of this approach is $O(n^2)$, which is equivalent to computing forces between all pairs of particles. There are two data structures, *cell lists* and *Verlet neighbor lists*, that can be used independently for speeding up this calculation, as well as an approach that combines the two.

<figure><img src="image_111.png" /><figcaption>Figure 7.1: Computing nonbonded forces within a cutoff, $r_c$. To compute forces involving the highlighted particle, only particles in the shaded regions are considered.</figcaption></figure>

### 7.1.2.1 Cell Lists

The idea of cell lists appears often in problems where a set of points that are nearby a given point is sought. Referring to Fig. 7.1(a), where we illustrate the idea with a 2-D example, a grid is laid over the set of particles. If the grid spacing is no less than $r_c$, then to compute the forces on particle $i$, only the particles in the cell containing $i$ and the 8 adjacent cells need to be considered. One sweep through all the particles is used to construct a list of particles for each cell. These cell lists are used to compute the forces for all particles. At the next time step, since the particles have moved, the cell lists must be regenerated or updated. The complexity of this approach is $O(n)$ for computing the data structure and $O(n \times n_c)$ for the
force computation, where $n_c$ is the average number of particles in 9 cells (27 cells in 3-D). The storage required for the cell list data structure is $O(n)$.

### 7.1.2.2 Verlet Neighbor Lists

The cell list structure is somewhat inefficient because, for each particle $i$, $n_c$ particles are considered, but this is much more than the number of particles within the cutoff $r_c$. A Verlet neighbor list is a list of particles within the cutoff for a particle $i$. Each particle has its own list, and thus the storage required is $O(n \times n_v)$ where $n_v$ is the average number of particles within the cutoff. Once these lists are constructed, computing the forces is then very fast, requiring the minimal complexity $O(n \times n_v)$. Constructing the list is more expensive, requiring examining all the particles for each particle, i.e., no less than the original complexity of $O(n^2)$. The advantage, however, is that the neighbor lists can be reused for many time steps if an expanded cutoff, $r_v$, is used. Referring to a 2-D example in Fig. 7.1(b), the neighbor list can be reused as long as no particle from outside the two circles moves inside the inner circle. If the maximum speed of the particles can be estimated or bounded, then one can determine a number of time steps for which it is safe to reuse the neighbor lists. (Alternatively, it may be possible to signal when any particle crosses to a position within the cutoff.) Technically, the Verlet neighbor list is the list of particles within the expanded cutoff, $r_v$.

### 7.1.2.3 Using Cell and Neighbor Lists Together

The hybrid approach is simply to use Verlet neighbor lists but to use cell lists to construct the neighbor lists. This reduces the high cost when neighbor lists need to be regenerated. This hybrid approach is very effective and is often the approach used in state-of-the-art molecular dynamics software.

Both cell lists and Verlet neighbor lists can be modified to exploit the fact that the force $f_{ij}$ on particle $i$ due to particle $j$ is equal to $-f_{ji}$ (Newton's third law) and only needs to be computed once. For example, for cell lists, only 4 of the 8 cells (in 2-D) need to be considered.

## 7.1.3 Computing Long-Range Forces

Electrostatic forces are challenging to compute because they are long-range: each particle feels a non-negligible electrostatic force from all other particles in the simulation. An approximation that is sometimes used is to truncate the force calculation for a particle after a certain cutoff radius (as is done for short-range van der Waals forces). This generally produces unacceptable artifacts in the results, however.

There are several more accurate methods for speeding up the computation of electrostatic forces, avoiding the $O(n^2)$ sum over all pairs of $n$ particles. We briefly outline some of these methods here.

### 7.1.3.1 Hierarchical N-body Methods

Hierarchical N-body methods, including the Barnes-Hut method and the fast multipole method, are popular for astrophysical particle simulations, but are typically too costly for the accuracy required in biomolecular simulations. In the Barnes-Hut method, space is recursively divided into 8 equal cells (in 3-D) until each cell contains zero or one particles. Forces between nearby particles are computed individually, as normal, but
for distant particles, forces are computed between one particle and a set of distant particles within a cell. An accuracy measure is used to determine if the force can be computed using a distant cell or must be computed by individually considering its children cells. The Barnes-Hut method has complexity $O(n \log n)$. The fast multipole method has complexity $O(n)$; this method calculates the potential and does not calculate forces directly.

### 7.1.3.2 Particle-Mesh Methods

In particle-mesh methods, we exploit the Poisson equation

$$ \nabla^2 \phi = -\frac{1}{\epsilon} \rho $$

which relates the potential $\phi$ to the charge density $\rho$, where $1/\epsilon$ is a constant of proportionality. To utilize this equation, we discretize space using a mesh, assign charges to the mesh points, solve Poisson's equation on the mesh to arrive at the potential on the mesh. The force is the negative gradient of the potential (for conservative forces such as electrostatic forces). A number of techniques have been developed for distributing point charges in space to a set of mesh points and also for numerically interpolating the force on the point charges due to the potentials at the mesh points. Many fast methods are available for solving the Poisson equation, including multigrid methods and fast Fourier transforms. With respect to terminology, particle-mesh methods are in contrast to the naive particle-particle method where forces are computed between all pairs of particles.

It turns out that particle-mesh methods are not very accurate, and a more accurate alternative is to split each force into a short-range, rapidly-varying part and a long-range, slowly-varying part:

$$ f_{ij} = f_{ij}^{sr} + f_{ij}^{lr} $$

One way to accomplish this easily is to weigh $f$ by a function $h(r)$, which emphasizes the short-range part (small $r$) and by $1 - h(r)$ which emphasizes the long-range part (large $r$). The short-range part is computed by computing the interaction of all pairs of particles within a cutoff (a particle-particle method) and the long-range part is computed using the particle-mesh method. The resulting method, called particle-particle-particle-mesh (PPPM, or P³M) is due to Hockney and Eastwood, in a series of papers beginning in 1973.

### 7.1.3.3 Ewald Method

The Ewald method is the most popular of the methods described so far for electrostatic forces in biomolecular simulations and was developed for the case of periodic boundary conditions. The structure of the method is similar to PPPM in that the force is split between short-range and long-range parts. Again, the short-range part is computed using particle-particle methods, and the long-range part is computed using Fourier transforms. Variants of the Ewald method are very similar to PPPM in that the long-range part uses a mesh, and fast Fourier transforms are used to solve the Poisson equation on the mesh. For additional details, see, for example [59]. In Section 7.3, we describe the parallelization of the 3-D FFT to solve the 3-D Poisson equation.
## 7.2 Parallel Decompositions

We now discuss the parallel computation of forces. Plimpton [136] created a very useful categorization of molecular dynamics parallelization methods, identifying *atom*, *force*, and *spatial* decomposition methods. Here, we closely follow his description of these methods. We also add a fourth category which has come to be recognized as differing from the earlier categories, called *neutral territory* methods, a name coined by Shaw [146]. Neutral territory methods are currently used by many state-of-the-art molecular dynamics codes. Spatial decompositions and neutral territory methods are particularly advantageous for parallelizing cutoff-based calculations.

### 7.2.1 Atom Decompositions

In an atom decomposition, each particle is assigned to one processor, and that processor is responsible for computing the particle's forces and updating its position for the entire simulation. For the computation to be roughly balanced, each processor is assigned approximately the same number of particles (a random distribution works well). An important point of atom decompositions is that each processor generally needs to communicate with all other processors to share updated particle positions.

<figure><img src="image_112.png" /><figcaption>Figure 7.2: Atom decomposition, showing a force matrix of 16 particles distributed among 8 processors. A dot represents a nonzero entry in the force matrix. On the left, the matrix is symmetric; on the right, only one element of a pair of skew-symmetric elements is computed, to take advantage of Newton's third law.</figcaption></figure>

An atom decomposition is illustrated by the *force matrix* in Fig. 7.2(a). For *n* particles, the force matrix is an *n*-by-*n* matrix; the rows and columns are numbered by particle indices. A nonzero entry $f_{ij}$ in the matrix denotes a nonzero force on particle *i* due to particle *j* which must be computed. This force may be a nonbonded and/or a bonded force. When cutoffs are used, the matrix is sparse, as in this example. The matrix is dense if forces are computed between all pairs of particles. The matrix is skew-symmetric because of Newton's third law, $f_{ij} = -f_{ji}$. The lines in Fig. 7.2(a) show how the particles are partitioned. In the figure, 16 particles are partitioned among 8 processors.
Algorithm 1 shows one time step from the point of view of one processor. At the beginning of the time step, each processor holds the positions of particles assigned to it.

**Algorithm 1 Atom decomposition time step**

1: send/receive particle positions to/from all other processors

2: (if nonbonded cutoffs are used) determine which nonbonded forces need to be computed

3: compute forces for particles assigned to this processor

4: update positions (integration) for particles assigned to this processor

An optimization is to halve the amount of computation, which is possible because the force matrix is skew-symmetric. To do this, we choose exactly one of $f_{ij}$ or $f_{ji}$ for all skew-symmetric pairs such that each processor is responsible for computing approximately the same number of forces. Choosing the upper or lower triangular part of the force matrix is a bad choice because the computational load is unbalanced. A better choice is to compute $f_{ij}$ if $i + j$ is even in the upper triangle, or if $i + j$ is odd in the lower triangle, as shown in Fig. 7.2(b). There are many other options.

When taking advantage of skew-symmetry in the force matrix, all the forces on a particle owned by a processor are no longer computed by that processor. For example, in Fig. 7.2(b), the forces on particle 1 are no longer computed only by the first processor. To complete the force calculation, processors must communicate to send forces that are needed by other processors and receive forces that are computed by other processors. The above algorithm must now be modified by adding a communication step (step 4) as shown in Algorithm 2.

**Algorithm 2 Atom decomposition time step, without redundant calculations**

1: send/receive particle positions to/from all other processors

2: (if nonbonded cutoffs are used) determine which nonbonded forces need to be computed

3: compute *partial* forces for particles assigned to this processor

4: send particle forces needed by other processors and receive particle forces needed by this processor

5: update positions (integration) for particles assigned to this processor

This algorithm is advantageous if the extra communication is outweighed by the savings in computation. Note that the amount of communication doubles in general.

### 7.2.2 Force Decompositions

In a force decomposition, the forces are distributed among the processors for computation. A straightforward way to do this is to partition the force matrix into blocks and to assign each block to a processor. Fig. 7.3(a) illustrates this for the case of 16 particles and 16 processors. Particles also need to be assigned to processors (as in atom decompositions) for the purpose of having processors assigned to update particle positions. In the example of the Figure, processor *i* is assigned to update the positions of particle *i*; in practical problems, a processor would be assigned to update the positions of many particles. Note that, again, we first consider the case of a skew-symmetric force matrix.
<figure><img src="image_113.png" /><figcaption>Figure 7.3: Force decomposition, showing a force matrix of 16 particles and forces partitioned among 16 processors.</figcaption></figure>

We now examine the communication required in a time step for a force decomposition. Consider processor 3, which computes partial forces for particles 0, 1, 2, 3, and needs positions from particles 0, 1, 2, 3, and also 12, 13, 14, 15. Thus processor 3 needs to perform communication with processors 0, 1, 2, 3, and processors 12, 13, 14, 15. After forces have been computed by all processors, processor 3 needs to collect forces on particle 3 computed by other processors. Thus processor 2 needs to perform communication again with processors 0, 1, 2, 3.

Algorithm 3 shows what is performed in one time step, from the point-of-view of one processor. At the beginning of the time step, each processor holds the positions of all the particles assigned to it.

**Algorithm 3** Force decomposition time step

1: send positions of my assigned particles which are needed by other processors; receive *row* particle positions needed by my processor (this communication is between processors in the same processor row, e.g., processor 3 communicates with processors 0, 1, 2, 3)

2: receive *column* particle positions needed by my processor (this communication is generally with processors in another processor row, e.g., processor 3 communicates with processors 12, 13, 14, 15)

3: (if nonbonded cutoffs are used) determine which nonbonded forces need to be computed

4: compute forces for my assigned particles

5: send forces needed by other processors; receive forces needed for my assigned particles (this communication is between processors in the same processor row, e.g., processor 3 communicates with processors 0, 1, 2, 3)

6: update positions (integration) for my assigned particles

In general, if there are *p* processors (and *p* is square, for simplicity), then the force matrix is partitioned
into $\sqrt{p}$ by $\sqrt{p}$ blocks. The force decomposition just described requires a processor to communicate in three steps, with $\sqrt{p}$ processors in each step. This is much more efficient than atom decompositions which require communications among all $p$ processors.

We can also exploit Newton’s third law in force decompositions. Like for atom decompositions, we first choose a modified force matrix where only one of $f_{ij}$ and $f_{ji}$ is computed. The forces on particle $i$ are computed by a row of processors and now also by a column of processors. Thus an extra step of communication is needed by each processor to collect forces from a column of processors for particles assigned to it. Whereas there were three communication steps, there are now four communication steps when Newton’s third law is exploited (the communication is not doubled in this case as in atom decompositions).

A modification to the force decomposition saves some communication. In Fig. 7.4, the columns are re-ordered using a *block-cyclic* ordering. Consider again processor 3, which computes partial forces for particles 0, 1, 2, 3. It needs positions from particles 0, 1, 2, 3, as before, but now also with processors 3, 7, 11, 15. The latter are processors in the same processor column as processor 3. Thus all communications are within the same processor row or processor column, which may be advantageous on mesh-based network architectures. The modified method is shown as Algorithm 4.

<figure><img src="image_114.png" /><figcaption>Figure 7.4: Force decomposition, with permuted columns in the force matrix. Note that columns 3, 7, 11, 15 are now in the block column corresponding to processors 3, 7, 11, 15 (the same indices), etc.</figcaption></figure>

### 7.2.3 Spatial Decompositions

In a spatial decomposition, space is decomposed into cells. Each cell is assigned to a processor which is responsible for computing the forces on particles that lie inside the cell. Fig. 7.5(a) illustrates a spatial decomposition into 64 cells for the case of a 2-D simulation. (This is a decomposition of space and is not to be confused with a force matrix.) Typically, the number of cells is chosen to be equal to the number of processors. Since particles move during the simulation, the assignment of particles to cells changes as well. This is in contrast to atom and force decompositions.
**Algorithm 4** Force decomposition time step, with permuted columns of force matrix

1: send positions of my assigned particles which are needed by other processors; receive *row* particle positions needed by my processor (this communication is between processors in the same processor row, e.g., processor 3 communicates with processors 0, 1, 2, 3)

2: receive *column* particle positions needed by my processor (this communication is generally with processors the same processor column, e.g., processor 3 communicates with processors 3, 7, 11, 15)

3: (if nonbonded cutoffs are used) determine which nonbonded forces need to be computed

4: compute forces for my assigned particles

5: send forces needed by other processors; receive forces needed for my assigned particles (this communication is between processors in the same processor row, e.g., processor 3 communicates with processors 0, 1, 2, 3)

6: update positions (integration) for my assigned particles

Fig. 7.5(b) shows one cell (center square) and the region of space (shaded) that contains particles that are potentially within the cutoff radius, $r_c$, with particles in the given cell. The shaded region is often called the *import region*, since the given cell must import positions of particles lying in this region to perform its force calculation. Note that not all particles in the given cell must interact with all particles in the *import* region, especially if the *import* region is large compared to the cutoff radius.

<figure><img src="image_115.png" /><figcaption>Figure 7.5: Spatial decomposition, showing particles in a 2-D computational box, (a) partitioned into 64 cells, (b) *import* region for one cell.</figcaption></figure>

**Algorithm 5** shows what each processor performs in one time step. We assume that at the beginning of the time step, each processor holds the positions of the particles in its cell.

To exploit Newton’s third law, the shape of the *import* region can be halved. Now each processor only computes a partial force on particles in its cell, and needs to receive forces from other processors to compute the total force on these particles. Thus an extra step of communication is involved. We leave it as an exercise to the reader to work out the details of the modified *import* region and the pseudocode for this case.
**Algorithm 5 Spatial decomposition time step**

1: send positions needed by other processors for particles in their import regions; receive positions for particles in my import region

2: compute forces for my assigned particles

3: update positions (integration) for my assigned particles

In the implementation of a spatial decomposition method, each cell is associated with a list of particles in its import region, similar to a Verlet neighbor list. Like a Verlet neighbor list, it is not necessary to update this list at every time step, if the import region is expanded slightly. This allows the import region list to be reused for several time steps, corresponding to the amount of time it takes a particle to traverse the width of the expanded region. This is exactly analogous to Verlet neighbor lists.

In summary, the main advantage of spatial decomposition methods is that they only require communication between processors corresponding to nearby particles. A disadvantage of spatial decomposition methods is that, for very large numbers of processors, the import region is large compared to the number of particles contained inside each cell.

### 7.2.4 Neutral Territory Methods

Our description of neutral territory methods follows closely that of Shaw [146]. A neutral territory method can be viewed as combining aspects of spatial decompositions and force decompositions. To parallelize the integration step, particles are assigned to processors according to a partitioning of space. To parallelize the force computation, each processor computes the forces between two sets of particles, but these particles may be unrelated to the particles that have been assigned to the processor for integration. As a result of this additional flexibility, neutral territory methods may require much less communication than spatial decomposition methods.

An example of a neutral territory method is shown in Fig. 7.6 for the case of a 2-D simulation. In the method shown in the Figure, the given processor is assigned the computation of forces between particles lying in the horizontal bar with particles lying in the vertical bar. These two regions thus form the import region for this method. By comparing to Fig. 7.6(b), the import region for this neutral territory method is much smaller than that for the corresponding spatial decomposition method. The advantage is greater when the size of the cells corresponding to each processor is small compared to the cutoff radius.

After the forces are computed, the given processor sends the forces it has computed to the processors that need these forces for integration. We thus have Algorithm 6.

**Algorithm 6 Neutral territory method time step**

1: send and receive particle positions corresponding to import regions

2: compute forces assigned to this processor

3: send and receive forces required for integration

4: update positions (integration) for particles assigned to this processor

Like other methods, the import region of the neutral territory method can be modified to take advantage of
<figure><img src="image_116.png" /><figcaption>Figure 7.6: Neutral territory method, showing particles in a 2-D computational box and the import region (shaded) for one cell (center square). This Figure can be compared directly to the spatial decomposition case of Fig. 7.5(b). See Shaw [146] for additional details.</figcaption></figure>

Newton’s third law. We refer to Shaw [146] for additional details and for illustrations of neutral territory methods in 3-D simulations.

7.3 Parallel Fast Fourier Transform

A common component of many methods for computing long-range forces is the 3-D FFT for solving the Poisson equation on a 3-D mesh. The Fourier transform diagonalizes the Poisson operator (called the Laplacian) and one forward and one inverse FFT transformation are required in the solution. Consider the discrete Laplacian operator L (with periodic boundary conditions) and the solution of φ in −Lφ = ρ . Let F denote the Fourier transform. The original problem is equivalent to

$$
-(FLF^{-1})F\phi = F\rho \\
\phi = -F^{-1}(FLF^{-1})^{-1}F\rho.
$$

The matrix FLF − 1 is diagonal. The forward Fourier transform F is applied to ρ , then the Fourier-space components are scaled by the inverse of the diagonal matrix, and finally, the inverse Fourier transform F − 1 is applied to obtain the solution φ .

For realistic protein sizes, a mesh spacing of approximately 1 Ångstrom is typically used, leading to a 3-D mesh that is quite small by many standards: 64 × 64 × 64, or 128 × 128 × 128. Parallel computation would often not be applied to a problem of this size, but parallel computation must be used because the data ρ is already distributed among the parallel processors (assuming a spatial decomposition is used).

A 3-D FFT is computed by computing 1-D FFTs in sequence along each of the three dimensions. For the 64 × 64 × 64 mesh size, this is 4096 1-D FFTs of dimension 64. The parallel FFT calculation is typically bound by communication. The best parallelization of the FFT depends on the size of the transforms and
the architecture of the computer network. Below, we first describe some concepts for parallel 1-D FFTs and then describe some concepts for parallel 3-D FFTs. For current software and research dedicated to the parallelization and efficient computation (using SIMD operations) of large 1-D transforms, we refer to the SPIRAL and FFTW packages. These packages use autotuning to generate FFT codes that are efficient for the user's computer architecture.

## 7.3.1 Parallel 1-D FFT

### 7.3.1.1 1-D FFT without Transpose

Fig. 7.7 shows the data dependencies (data flow diagram) between the inputs (left) and outputs (right) for the 16-point radix-2 decimation-in-frequency FFT algorithm. (Not shown is a bit-reversal permutation that may be necessary in the computation.) The Figure also shows a partitioning of the computation among four processors. In this parallelization, the initial data is not moved among processors, but communication occurs during the computation. In the example shown in the Figure, communication occurs in the first two FFT stages; the final two stages do not involve communication. When communication does occur, every processor communicates with exactly one other processor.

### 7.3.1.2 1-D FFT with Transpose

Use of transposes is common to parallelize FFT computations. Fig. 7.8(a) shows the same data flow diagram as in Fig. 7.7, but horizontal lines have been removed and additional index labels have been added for clarity. As before, the first two FFT stages are performed without communication. The data is then transposed among the processors. With this transposed data layout, the last two FFT stages can be performed without communication. The final data is not in the original order; an additional transpose may be needed, or the data may be used in this transposed order. Fig. 7.8(b) shows how the indices are partitioned among four processors before and after the transpose. From these two Figures, notice that the first two stages have data dependencies that only involve indices in the same partition. The same is true for the second two stages for the partitioning after the transpose. Observe also that the structure of the computations before and after the transpose are identical.

## 7.3.2 Parallel 3-D FFT

### 7.3.2.1 3-D FFT with Block Decomposition

Fig. 7.9(a) shows a block decomposition of the FFT input data when a spatial decomposition is used for a mesh of size 8 × 8 × 8 distributed across 64 processors arranged in a 4 × 4 × 4 topology. The parallel 1-D FFT algorithms can be applied in each of the dimensions. For the example shown in the Figure, each 1-D FFT computation involves 4 processors. Each processor performs multiple 1-D FFTs simultaneously (four in this example). Within each processor, data is ordered contiguously if traversing one of the dimensions, and thus data access is strided for computation in the other two dimensions. Strided data access can be slow, and thus it may be worthwhile to reorder the data within each processor when computing the FFT for each of the dimensions.
<figure><img src="image_117.png" /><figcaption>Figure 7.7: Data flow diagram for 1-D FFT for 16 points. The shaded regions show a decomposition for 4 processors (one processor per region). In this parallelization, the first two FFT stages have no communication; the last two FFT stages do have communication.</figcaption></figure>

### 7.3.2.2 3-D FFT with Slab Decomposition

The slab decomposition is shown in Fig. 7.9(b) for the case of 4 processors. Each processor holds one or more planes of the input data. This decomposition is used if the input data is already distributed in slabs, or if it can be easily redistributed this way. The two 1-D FFTs in the plane of the slabs require no communication. The remaining 1-D FFTs require communication and could use one of the two approaches for parallel 1-D FFTs described above. A disadvantage of the slab decomposition is that for large numbers of processors, the number of processors may exceed the number of points in the 3-D FFT along any one dimension. An alternative is the pencil decomposition below.

### 7.3.2.3 3-D FFT with Pencil Decomposition

The pencil decomposition is shown in Fig. 7.9(c) for the case of 16 processors. Each processor holds one or more pencils of the input data. If the original input data is distributed in blocks as in Fig. 7.9(a), then communication among a row of processors (in a 3-D processor mesh) can distribute the data into the pencil decomposition. The 1-D FFTs can then be performed with no communication. To perform the 1-D FFT in another dimension, the data needs to be redistributed into pencils in another dimension. In total, four communication stages are needed for the entire 3-D FFT computation.
<figure><img src="image_118.png" /><figcaption>(a) Data flow diagram (shown without horizontal lines for clarity) for 1-D FFT for 16 points.</figcaption></figure>

<caption>(b) Partitioning of the indices before (left) and after (right) the transpose.</caption>

Figure 7.8: 1-D FFT with transpose. The first two stages do not involve communication. The data is then transposed among the processors. As a result, the second two stages also do not involve communication.

## 7.4 Integration for Molecular Dynamics

To numerically integrate the system of ordinary differential equations in molecular dynamics, special methods are required, different than the traditional ODE solvers that were studied in Chapter 4. These special methods, called symplectic methods, are better than other methods at producing solutions that have constant energy, for example, for systems that are called Hamiltonian (which include systems from molecular dynamics). When Hamiltonian systems are integrated with many time steps over a long time interval, preservation of structure such as total energy is often more important than the order of accuracy of the method. In this section, we motivate some ideas and give some details of the Störmer-Verlet method, which is sufficient for simple molecular dynamics simulations.

Hamiltonian systems are a class of dynamical systems which conserve energy and which can be written in a form called Hamilton's equations. Consider, for simplicity, the simple harmonic oscillator

$$u'' = -u$$

where $u$ is the displacement of a single particle from an equilibrium point. This equation could model a particle with unit mass on a spring with unit spring constant. The force on a particle at position $u$ is $-u$. This system does not look like a molecular dynamics system but is useful for illustrating several ideas.

The above second order equation can be written as a system of first order equations

$$q' = p$$

$$p' = -q$$
<figure><img src="image_119.png" /><figcaption>Figure 7.9: Three data decompositions for 3-D FFTs.</figcaption></figure>

where $q = u$ and $p = u'$ which is common notation used in classical mechanics. The general solution is

$$ \begin{pmatrix} q \\ p \end{pmatrix} = \begin{pmatrix} \cos t & \sin t \\ -\sin t & \cos t \end{pmatrix} \begin{pmatrix} q \\ p \end{pmatrix}. $$

The kinetic energy of the simple harmonic oscillator is $p^2/2$ and the potential energy is $q^2/2$ (the negative gradient of potential energy is the force, $-q$). Thus the total energy is proportional to $q^2 + p^2$.

Now consider the solution of the system of first order equations by three methods, explicit Euler, implicit Euler, and a method called the Störmer-Verlet method. The initial condition is $(q, p) = (1, 0)$. We use a time step of $h = 0.05$ and take 500 steps. We plot $q$ and $p$ on the horizontal and vertical axes, respectively (called a *phase plot*). The exact solution, as given above, is a unit circle centered at the origin.

Figure 7.10 shows the solutions. For explicit Euler, the solution spirals outward, meaning the displacement and momentum of the solution increases over time. The opposite is true for the implicit Euler method. A plot of the total energy would show the energy increasing and decreasing for the two cases, respectively. The solutions are better when smaller time steps are taken or when higher order methods are used, but these methods are not at all appropriate for integration of symplectic systems over long periods of time. Figure 7.10(c) shows the solution using a symplectic method called the Störmer-Verlet method. The solution shows that $q^2 + p^2$ is preserved much better than in the other two methods.

The Störmer-Verlet method is derived very easily. We derive it for the second order equation

$$ u'' = f(t, u) $$

by simply replacing the left-hand side with a finite difference approximation

$$ \frac{u_{k+1} - 2u_k + u_{k-1}}{h^2} = f(t_k, u_k) $$

which can be rearranged to obtain the method

$$ u_{k+1} = 2u_k - u_{k-1} + h^2 f(t_k, u_k). $$

The formula can equivalently be derived from Taylor series. The method is similar to linear multistep methods in that some other technique is needed to supply the initial step of the method. The method is also
<figure><img src="image_120.png" /><figcaption>Figure 7.10: Phase plot of the solution of the simple harmonic oscillator for three methods with initial value (1,0), time step 0.05, and 500 steps. For explicit Euler, the solution spirals outward; for implicit Euler, the solution spirals inward; the total energy is best preserved with the Störmer-Verlet method.</figcaption></figure>

time-reversible, because the formula is the same if $k+1$ and $k-1$ are swapped. To explain why this method is symplectic, unfortunately, is beyond the scope of this introduction.

The method as written above has a number of disadvantages, the most severe being that the addition of the small $h^2$ term is subject to catastrophic cancellation. Thus this formula should not be used in this form, and a number of mathematically equivalent formulas (which can be derived from the formula above) have been developed.

One alternative formula is the *leap-frog* method:

$$ 
\begin{aligned}
u_{k+1} &= u_k + h v_{k+1/2} \\
v_{k+1/2} &= v_{k-1/2} + h f(t_k, u_k)
\end{aligned}
 $$

where $v$ is the first derivative (velocity) which is offset from the displacement $u$ by a half step. This formula does not suffer from the same roundoff problems and also makes available the velocities, although they need to be re-centered with displacements in order to calculate total energy at a given step. The second of this pair of equations is basically a finite difference formula.

A third form of the Störmer-Verlet method is the *velocity Verlet* variant:

$$ 
\begin{aligned}
u_{k+1} &= u_k + h v_k + \frac{h^2}{2} f(t_k, u_k) \\
v_{k+1} &= v_k + \frac{h^2}{2} (f(t_k, u_k) + f(t_{k+1}, u_{k+1}))
\end{aligned}
 $$

where now the velocities are computed at the same points as the displacements. Each of these algorithms can be implemented such that only two sets of quantities need to be stored (two previous positions, or a position and a velocity). These variants of the Störmer-Verlet method are popular because of their simplicity, requiring only one costly force evaluation per step. Higher-order methods have generally not been practical.
The velocity Verlet scheme is also the basis of multiple time step algorithms for molecular dynamics. In these algorithms, the slowly-varying (typically long-range) forces are evaluated less frequently and update the positions less frequently than the quickly-varying (typically short-range) forces. Finally, many state-of-the-art molecular dynamics integrate a Hamiltonian system that has been modified in order to control the simulation temperature and pressure. Much more complicated symplectic methods have been developed for these systems.
# Chapter 8

## Combinatorial algorithms

In this chapter we will briefly consider a few combinatorial algorithms: sorting, and prime number finding with the Sieve of Eratosthenes.

### 8.1 Sorting

Sorting is not a common operation in scientific computing: one expects it to be more important in databases, whether these be financial or biological (for instance in sequence alignment). However, it sometimes comes up, for instance in *Adaptive Mesh Refinement (AMR)* and other applications where significant manipulations of data structures occurs.

In this section we will briefly look at some basic algorithms and how they can be done in parallel. For more details, see [104] and the references therein.

#### 8.1.1 Brief introduction to sorting

##### 8.1.1.1 Complexity

There are many sorting algorithms. Traditionally, they have been distinguished by their computational complexity, that is, given an array of $n$ elements, how many operations does it take to sort them, as a function of $n$.

Theoretically one can show that a sorting algorithm has to have at least complexity $O(n \log n)^1$. There are indeed several algorithms that are guaranteed to attain this complexity, but a very popular algorithm, called *Quicksort* has only an ‘expected’ complexity of $O(n \log n)$, and a worst case complexity of $O(n^2)$.

<footnote n="1">1. One can consider a sorting algorithm as a decision tree: a first comparison is made, depending on it two other comparisons are made, et cetera. Thus, an actual sorting becomes a path through this decision tree. If every path has running time $h$, the tree has $2^h$ nodes. Since a sequence of $n$ elements can be ordered in $n!$ ways, the tree needs to have enough paths to accomodate all of these; in other words, $2^h \ge n!$. Using Stirling’s formula, this means that $n \ge O(n \log n)$</footnote>
### 8.1.1.2 Sorting networks

There is another way of looking at sorting algorithms, distinguishing them by whether they work independent of the state of the input data, or not. As an example of an algorithm that depends on its input data, imagine a sorting algorithm that makes a pass over the data to check if it was already sorted. This algorithm has linear complexity (that is, $n$ operations for $n$ elements) for a sorted list, and something higher for an unsorted list.

The popular Quicksort algorithm has such data dependencies, since it involves finding 'pivot' elements.

On the other hand, the so-called *bubble sort* algorithm can be described independently of the data to be sorted. This algorithm is given by:

```
for t from n - 1 down to 1 do
    for e from 1 to t do
        if elements e and e + 1 are ordered the wrong way then
            exchange them
```

**Algorithm 2:** The bubble sort algorithm

It is easy to see that this algorithm has a complexity of $O(n^2)$: the inner loop does $t$ comparisons and up to $t$ exchanges. Summing this from 1 to $n-1$ gives approximately $n^2/2$ comparisons and at most the same number of exchanges.

Such a data-independent sorting algorithm is sometimes called *sorting network*. It can be considered as custom hardware that implements just one algorithm. The basic hardware element is the *compare-and-swap* element, which has two inputs and two outputs. For two inputs $x, y$ the outputs are $\max(x, y)$ and $\min(x, y)$.

Below we will consider the *Bitonic sort* algorithm as a prime example of a sorting network.

### 8.1.1.3 Parallel complexity

Above we remarked that sorting sequentially takes at least $O(N \log N)$ time. If we can have perfect speedup, using for simplicity $P = N$ processors, we would have parallel time $O(\log N)$. If the parallel time is more than that, we define the *sequential complexity* as the total number of operations. For instance, below you will see *odd-even transposition sort*, which takes $N$ steps, giving a sequential complexity of $N^2$, and bitonic sort which has a parallel time of $(\log N)^2$ and a sequential complexity of $N(\log N)^2$.

## 8.1.2 Odd-even transposition sort

The *odd-even transposition sort* algorithm is a simple algorithm, with as main virtue that it is easily implemented on a linear area of processors. A single step of the algorithm consists of two substeps:

* Every even-numbered processor does a compare-and-swap with its right neighbour; then

* Every odd-numbered processor does a compare-and-swap with its right neighbour.

**Theorem 2** After $N$ steps, each consisting of the two substeps just given, a sequence is sorted.

With a parallel time of $2N$, this gives a sequential complexity $O(N^2)$.
### 8.1.3 Quicksort

Quicksort is a recursive algorithm, that, unlike bubble sort, is not deterministic. It is a two step procedure, based on a reordering of the sequence<fnref n="2" />:

**Algorithm:** Dutch National Flag ordering of an array

**Input:** An array of elements, and a 'pivot' value

**Output:** The input array with elements ordered as red-white-blue, where red elements are larger than the pivot, white elements are equal to the pivot, and blue elements are less than the pivot

We state without proof that this can be done in $O(n)$ operations. With this, quicksort becomes:

**Algorithm:** Quicksort

**Input:** An array of elements

**Output:** The input array, sorted

**while** The array is longer than one element do

pick an arbitrary value as pivot

apply the Dutch National Flag reordering to this array

Quicksort( the blue elements )

Quicksort( the red elements )

The indeterminacy of this algorithm, and the variance in its complexity, stems from the pivot choice. In the worst case, the pivot is always the (unique) smallest element of the array. There will then be no blue elements, the only white element is the pivot, and the recursive call will be on the array of $n - 1$ red elements. It is easy to see that the running time will then be $O(n^2)$. On the other hand, if the pivot is always (close to) the median, that is, the element that is intermediate in size, then the recursive calls will have an about equal running time, and we get a recursive formula for the running time:

$$T_n = 2T_{n/2} + O(n)$$

which is (again without proof) $O(n \log n)$.

We will now consider parallel implementations of quicksort.

#### 8.1.3.1 Quicksort in shared memory

A simple parallelization of the quicksort algorithm can be achieved by executing the two recursive calls in parallel. This is easiest realized with a shared memory model, and threads (section 2.7.1) for the recursive calls. However, this implementation is not efficient.

On an array of length $n$, and with perfect pivot choice, there will be $n$ threads active in the final stage of the algorithm. Optimally, we would want a parallel algorithm to run in $O(\log n)$ time, but here the time is dominated by the initial reordering of the array by the first thread.

<footnote n="2">2. The name is explained by its origin with the Dutch computer scientist Edsger Dijkstra; see http://en.wikipedia.org/wiki/Dutch_national_flag Problem.</footnote>
**Exercise 8.1.** Make this argument precise. What is the total running time, the speedup, and the efficiency of parallelizing the quicksort algorithm this way?

Is there a way to make splitting the array more efficient? As it turns out, yes, and the key is to use a parallel prefix operation; see appendix 18. If the array of values is $x_1, \dots, x_n$, we use a parallel prefix to compute how many elements are less than the pivot $\pi$:

$$X_i = \#\{x_j: j < i \land x_j < \pi\}.$$

With this, if a processor looks at $x_i$, and $x_i$ is less than the pivot, it needs to be moved to location $X_i + 1$ in the array where elements are split according to the pivot. Similarly, one would could how many elements there are over the pivot, and move those accordingly.

This shows that each pivoting step can be done in $O(\log n)$ time, and since there $\log n$ steps to the sorting algorithm, the total algorithm runs in $O((\log n)^2)$ time.

### 8.1.3.2 Quicksort on a hypercube

As was apparent from the previous section, for an efficient parallelization of the quicksort algorithm, we need to make the Dutch National Flag reordering parallel too. Let us then assume that the array has been partitioned over the $p$ processors of a hypercube of dimension $d$ (meaning that $p = 2^d$).

In the first step of the parallel algorithm, we choose a pivot, and broadcast it to all processors. All processors will then apply the reordering independently on their local data.

In order to bring together the red and blue elements in this first level, every processor is now paired up with one that has a binary address that is the same in every bit but the most significant one. In each pair, the blue elements are sent to the processor that has a 1 value in that bit; the red elements go to the processor that has a 0 value in that bit.

After this exchange (which is local, and therefore fully parallel), the processors with an address $1xxxxx$ have all the red elements, and the processors with an address $0xxxxx$ have all the blue elements. The previous steps can now be repeated on the subcubes.

This algorithm keeps all processors working in every step; however, it is susceptible to load imbalance if the chosen pivots are far from the median. Moreover, this load imbalance is not lessened during the sort process.

### 8.1.3.3 Quicksort on a general parallel processor

Quicksort can also be done on any parallel machine that has a linear ordering of the processors. We assume at first that every processor holds exactly one array element, and, because of the flag reordering, sorting will always involve a consecutive set of processors.

Parallel quicksort of an array (or subarray in a recursive call) starts by constructing a binary tree on the processors storing the array. A pivot value is chosen and broadcast through the tree. The tree structure is then used to count on each processor how many elements in the left and right subtree are less than, equal to, or more than the pivot value.
With this information, the root processor can compute where the red/white/blue regions are going to be stored. This information is sent down the tree, and every subtree computes the target locations for the elements in its subtree.

If we ignore network contention, the reordering can now be done in unit time, since each processor sends at most one element. This means that each stage only takes time in summing the number of blue and red elements in the subtrees, which is $O(\log n)$ on the top level, $O(\log n/2)$ on the next, etcetera. This makes for almost perfect speedup.

### 8.1.4 Bitonic sort

To motivate bitonic sorting, suppose a sequence $x = \langle x_0, \dots, x_n - 1 \rangle$ consists of an ascending followed by a descending part. Now split this sequence into two subsequences of equal length defined by:

<figure><img src="image_121.png" /><figcaption>Figure 8.1: Illustration of a splitting a bitonic sequence</figcaption></figure>

$$ 
\begin{aligned}
s_1 &= \langle \min\{x_0, x_{n/2}\}, \dots, \min\{x_{n/2-1}, x_{n-1}\} \rangle \\
s_2 &= \langle \max\{x_0, x_{n/2}\}, \dots, \max\{x_{n/2-1}, x_{n-1}\} \rangle
\end{aligned} 
\quad (8.1) 
$$

From the picture it is easy to see that $s_1, s_2$ are again sequences with an ascending and descending part. Moreover, all the elements in $s_1$ are less than all the elements in $s_2$.

We call (8.1) an ascending bitonic sorter, since the second subsequence contains elements larger than in the first. Likewise we can construct a descending sorter by reversing the roles of maximum and minimum.

It's not hard to imagine that this is a step in a sorting algorithm: starting out with a sequence on this form, recursive application of formula (8.1) gives a sorted sequence. Figure 8.2 shows how 4 bitonic sorters, over distances 8,4,2,1 respectively, will sort a sequence of length 16.

The actual definition of a bitonic sequence is slightly more complicated. A sequence is bitonic if it consists of an ascending part followed by a descending part, or is a cyclic permutation of such a sequence.

**Exercise 8.2.** Prove that splitting a bitonic sequence according to formula (8.1) gives two bitonic sequences.

So the question is how to get a bitonic sequence. The answer is to use larger and larger bitonic networks.

* A bitonic sort of two elements gives you a sorted sequence.

* If you have two sequences of length two, one sorted up, the other sorted down, that is a bitonic sequence.
<figure><img src="image_122.png" /><figcaption>Figure 8.2: Illustration of a bitonic network that sorts a bitonic sequence of length 16</figcaption></figure>

* So this sequence of length four can be sorted in two bitonic steps.

* And two sorted sequences of length four form a bitonic sequence of length;

* which can be sorted in three bitonic steps; et cetera.

<figure><img src="image_123.png" /><figcaption>Figure 8.3: Full bitonic sort for 16 elements</figcaption></figure>

From this description you see that you $log_2 N$ stages to sort $N$ elements, where the $i$-th stage is of length $log_2 i$. This makes the total sequential complexity of bitonic sort $(log_2 N)^2$.

The sequence of operations in figure 8.3 is called a *sorting network*, built up out of simple *compare-and-swap* elements. There is no dependence on the value of the data elements, as is the case with quicksort.

## 8.2 Prime number finding

The sieve of Eratosthenes is a very old method for finding prime numbers. In a way, it is still the basis of many more modern methods.
Write down all natural numbers from 2 to some upper bound $N$, and we are going to mark these numbers as prime or definitely not-prime. All numbers are initially unmarked.

* The first unmarked number is 2: mark it as prime, and mark all its multiples as not-prime.

* The first unmarked number is 3: mark it as prime, and mark all its multiples as not-prime.

* The next number is 4, but it has been marked already, so mark 5 and its multiples.

* The next number is 6, but it has been marked already, so mark 7 and its multiples.

* Numbers 8,9,10 have been marked, so continue with 11.

* Et cetera.
# Chapter 9

## Graph analytics

Various problems in scientific computing can be formulated as graph problems (for an introduction to graph theory see Appendix 16); for instance, you have encountered the problem of load balancing (section 2.11.1.3) and finding independent sets (section 6.8.2).

Many traditional graph algorithms are not immediately, or at least not efficiently, applicable, since the graphs are often distributed, and traditional graph theory assumes global knowledge of the whole graph. Moreover, graph theory is often concerned with finding the optimal algorithm, which is usually not a parallel one. Therefore, parallel graph algorithms are a field of study by themselves.

Recently, new types of graph computations have arisen in scientific computing. Here the graph are no longer tools, but objects of study themselves. Examples are the *World Wide Web* or the *social graph* of *Facebook*, or the graph of all possible *protein interactions* in a living organism.

For this reason, combinatorial computational science is becoming a discipline in its own right. In this section we look at *graph analytics*: computations on large graphs. We start by discussing some classic algorithms, but we give them in an algebraic framework that will make parallel implementation much easier.

### 9.1 Traditional graph algorithms

We start by looking at a few ‘classic’ graph algorithms, and we discuss how they can be implemented in parallel. The connection between graphs and sparse matrices (see Appendix 16) will be crucial here: many graph algorithms have the structure of a *sparse matrix-vector multiplication*.

#### 9.1.1 Shortest path algorithms

There are several types of shortest path algorithms. For instance, in the *single source shortest path* algorithm one wants the shortest path from a given node to any other node. In *all-pairs shortest path* algorithms one wants to know the distance between any two nodes. Computing the actual paths is not part of these algorithms; however, it is usually easy to include some information by which the path can be reconstructed later.
We start with a simple algorithm: finding the single source shortest paths in an unweighted graph. This is simple to do with a Breadth-First Search (BFS):

**Input** : A graph, and a starting node $s$

**Output**: A function $d(v)$ that measures the distance from $s$ to $v$

Let $s$ be given, and set $d(s) = 0$

Initialize the finished set as $U = \{s\}$

Set $c = 1$

**while not finished do**

Let $V$ the neighbours of $U$ that are not themselves in $U$

if $V = \emptyset$ then

We're done

else

Set $d(v) = c + 1$ for all $v \in V$.

$U \leftarrow U \cup V$

Increase $c \leftarrow c + 1$

This way of formulating an algorithm is useful for theoretical purposes: typically you can formulate a *predicate* that is true for every iteration of the while loop. This then allows you to prove that the algorithm terminates and that it computes what you intended it to compute. And on a traditional processor this would indeed be how you program a graph algorithm. However, these days graphs such as from Facebook can be enormous, and you want to program your graph algorithm in parallel.

In that case, the traditional formulations fall short:

* They are often based on queues to which nodes are added or subtracted; this means that there is some form of shared memory.

* Statements about nodes and neighbours are made without wondering if these satisfy any sort of spatial locality; if a node is touched more than once there is no guarantee of temporal locality.

We will now show that graph algorithms can often be considered as *sparse-matrix algorithms*, which means that we can apply all the concepts and analysis we have developed for these.

If $G$ is the *adjacency matrix* of the graph, we can also formulate the shortest path algorithm analogous to a series of matrix-vector multiplications (see appendix section 16.4.4). Let $x$ be the vector tabulating distances from the source, that is, $x_i$ is the distance of node $i$ from the source. For any neighbour $j$ of $i$, the distance to the source is then $x_i + G_{ij}$, unless a shorter distance was already known. In other words, we can define a product

$$y^t = x^t G \equiv \forall i: y_j = \min\{x_j, \min_{i: G_{ij} \neq 0} \{x_i + 1\}\},$$

and the iterations of the above while-loop correspond to subsequent matrix-vector products under this definition.

This algorithm works because we can set $d(v)$ to its final value the first time we visit it: this happens precisely after a number of outer iterations equal to the path length. The total number of inner loop executions equals the number of edges in the graph. A weighted graph is somewhat trickier, since a path with more
stages can actually be shorter as measured in the sum of the weights of the stages. Here is the Bellman-Ford algorithm:

Let $s$ be given, and set $d(s) = 0$
Set $d(v) = \infty$ for all other nodes $v$

for $|E| - 1$ times do

for all edges $e = (u, v)$ do

Relax: **if** $d(u) + w_{uv} < d(v)$ **then**

Set $d(v) \leftarrow d(u) + w_{uv}$

This algorithm is correct since, for a given node $u$, after $k$ steps of the outer iteration it has considered all path $s \to u$ of $k$ stages.

**Exercise 9.1.** What is the complexity of this algorithm? Can the length of the outer loop be reduced if you have some knowledge of the graph?

We can again write this as a series of matrix-vector products, if we define the product as

$$y^t = x^t G \equiv \forall i: y_j = \min\{x_j, \min_{i: G_{ij} \neq 0} \{x_i + g_{ij}\}\},$$

This has essentially the same basis as above: the minimum distance to $j$ is minimum of an already computed distance, or the minimum distance to any node $i$ plus the transition $g_{ij}$.

## 9.1.2 Parallelization

Graph algorithms in their traditional formulation can be parallelized. Typically there are macro steps that are sequential, but in which a number of variables are considered independently. Thus there is indeed parallelism to be exploited.

However, these traditional formulations often feature sets of variables that are gradually built up or depleted. This is implementable by using a *task queue*, but this limits the implementation to some form of shared memory. Also, assignment of tasks to processing elements is completely arbitrary, so the algorithm will have little spatial or temporal locality.

For these reasons, a linear algebra formulation can be preferable. We definitely need this approach once distributed memory is considered, but even on multicore architectures it can pay off to encourage locality.

In section 6.5 we discussed the parallel evaluation of the sparse matrix-vector product. Because of the sparsity, only a partitioning by block rows or columns made sense. In effect, we let the partitioning be determined by one of the problem variables. This is also the only strategy that makes sense for single-source shortest path algorithms.

**Exercise 9.2.** Can you make an *a priori* argument for basing the parallelization on a distribution of the vector? How much data is involved in this operation, how much work, and how many sequential steps?

## 9.1.3 Floyd-Warshall all-pairs shortest path

The *Floyd-Warshall algorithm* is an example of an all-pairs shortest paths algorithm. It is a *dynamic programming* algorithm that is based on gradually increasing the set of intermediate nodes for the paths.
Specifically, in step k all paths u ↪ → v are considered that have intermediate nodes in the set k = {0, ..., k − 1}, and Δ$_{k}$(u, v) is the defined as the path length from u to v where all intermediate nodes are in k . Initially, this means that only graph edges are considered, and when k ≡ | V | we have considered all possible paths and we are done.

The computational step is

$$
\Delta_{k+1}(u, v) = \min\{\Delta_k(u, v), \Delta_k(u, k) + \Delta_k(k, v)\}. \tag{9.1}
$$

that is, the k-th estimate for the distance Δ(u, v) is the minimum of the old one, and a new path that has become feasible now that we are considering node k . This latter path is found by concatenating paths u ↪ → k and k ↪ → v .

Writing it algorithmically:

$$
\begin{array}{ll}
\textbf{for } k \textbf{ from zero to } |V| & \textbf{do} \\
& \textbf{for all nodes } u, v \textbf{ do} \\
& \qquad \Delta_{uv} \leftarrow f(\Delta_{uv}, \Delta_{uk}, \Delta_{kv})
\end{array}
$$

we see that this algorithm has a similar structure to Gaussian elimination, except that there the inner loop
would be ‘for all u, v > k’.

Algebraically:

$$
\begin{array}{ll}
\textbf{for } k \textbf{ from zero to } |V| & \textbf{do} \\
& D \leftarrow D_{\min} [D(:, k) \min \cdot_{+} D(k, :)]
\end{array}
$$

The Floyd-Warshall algorithm does not tell you the actual path. Storing those paths during the distance calculation above is costly, both in time and memory. A simpler solution is possible: we store a second matrix n(i, j) that has the highest node number of the path between i and j.

**Exercise 9.3.** Include the calculation of *n*(i, j) in the Floyd-Warshall algorithm, and describe
how to use this to find the shortest path between *i* and *j* once *d*(·, ·) and *n*(·, ·) are
known.

**Parallelization** In the single-source shortest path algorithm above we didn’t have much choice but to parallelize by distributing the vector rather than the matrix. This type of distribution is possible here too, and it corresponds to a one-dimensional distribution of the *D*(·, ·) quantities.

**Exercise 9.4.** Sketch the parallel implementation of this variant of the algorithm. Show that each *k*-th iteration involves a broadcast with processor *k* as root.

However, this approach runs into the same scaling problems as the matrix-vector product using a one-
dimensional distribution of the matrix; see section 6.2.2. Therefore we need to use a two-dimensional
distribution.

**Exercise 9.5.** Do the scaling analysis. In a weak scaling scenario with constant memory, what is the asymptotic efficiency?

**Exercise 9.6.** Sketch the Floyd-Warshall algorithm using a two-dimensional distribution of the *D*(·, ·) quantities.
Exercise 9.7. The parallel Floyd-Warshall algorithm performs quite some operations on zeros, certainly in the early stages. Can you design an algorithm that avoids those?

### 9.1.4 Spanning trees

In an undirected graph $G = (V, E)$, we call $T \subset E$ a 'tree' if it is connected and acyclic. It is called a *spanning tree* if additionally its edges contain all vertices. If the graph has edge weights $w_i: i \in E$, the tree has weight $\sum_{e \in T} w_e$, and we call a tree a *minimum spanning tree* if it has minimum weight. A minimum spanning tree need not be unique.

*Prim's algorithm*, a minimal variant of *Dijkstra's shortest path algorithm*, computes a spanning tree starting at a root. The root has path length zero, and all other nodes infinity. In each step, all nodes connected to the known tree nodes are considered, and their best known path length updated.

**for all vertices v do**

$$ \ell(v) \leftarrow \infty $$

$$ \ell(s) \leftarrow 0 $$

$$ Q \leftarrow V - \{s\} \text{ and } T \leftarrow \{s\} $$

**while** $Q \neq \emptyset$ **do**

let $u$ be the element in $Q$ with minimal $\ell(u)$ value

remove $u$ from $Q$, and add it to $T$

**for** $v \in Q$ **with** $(u, v) \in E$ **do**

**if** $\ell(u) + w_{uv} < \ell(v)$ **then**

$$ \text{Set } \ell(v) \leftarrow \ell(u) + w_{uv} $$

**Theorem 3** The above algorithm computes the shortest distance from each node to the root node.

*Proof.* The key to the correctness of this algorithm is the fact that we choose $u$ to have minimal $\ell(u)$ value. Call the true shortest path length to a vector $L(v)$. Since we start with an $\ell$ value of infinity and only ever decrease it, we always have $L(v) \le \ell(v)$.

Our induction hypothesis will be that, at any stage in the algorithm, for the nodes in the current $T$ the path length is determined correctly:

$$ u \in T \Rightarrow L(u) = \ell(u). $$

This is certainly true when the tree consists only of the root $s$. Now we need to prove the induction step: if for all nodes in the current tree the path length is correct, then we will also have $L(u) = \ell(u)$.

Suppose this is not true, and there is another path that is shorter. This path needs to go through some node $y$ that is not currently in $T$; this illustrated in figure 9.1. Let $x$ be the node in $T$ on the purported shortest path right before $y$. Now we have $\ell(u) > L(u)$ since we do not have the right pathlength yet, and $L(u) > L(x) + w_{xy}$ since there is at least one edge (which has positive weight) between $y$ and $u$. But $x \in T$, so $L(x) = \ell(x)$ and $L(x) + w_{xy} = \ell(x) + w_{xy}$. Now we observe that when $x$ was added to $T$ its neighbours were updated, so $\ell(y)$ is $\ell_x + w_{xy}$ or less. Stringing together these inequalities we find

$$ \ell(y) < \ell(x) + w_{xy} = L(x) + w_{xy} < L(u) < \ell(u) $$

which contradicts the fact that we choose $u$ to have minimal $\ell$ value.
<figure><img src="image_124.png" /><figcaption>Figure 9.1: Illustration of the correctness of Dijkstra’s algorithm</figcaption></figure>

To parallelize this algorithm we observe that the inner loop is independent and therefore easily parallelized. However, the outer loop has a choice that minimizes a function value. Computing this choice is a reduction operator, and subsequently it needs to be broadcast. This strategy makes the sequential time equal to $d \log P$ where $d$ is the depth of the spanning tree.

On a single processor, finding the minimum value in an array is naively an $O(N)$ operation, but through the use of a *priority queue* this can be reduced to $O(\log N)$. For the parallel version of the spanning tree algorithm the corresponding term is $O(\log(N/P))$, not counting the $O(\log P)$ cost of the reduction.

## 9.2 ‘Real world’ graphs

In discussions such as in section 4.2.3 you have seen how the discretization of PDEs leads to computational problems that have a graph aspect to them. Such graphs have properties that make them amenable to certain kinds of problems. For instance, using FDMs or FEMs to model two or three-dimensional objects leads graphs where each node is connected to just a few neighbours. This makes it easy to find *separators*, which in turn allows such solution methods as *nested dissection*; see section 6.8.1.

There are however applications with computationally intensive graph problems that do not look like FEM graphs. We will briefly look at the example of the world-wide web, and algorithms such *Google’s PageRank* which try to find authoritative nodes.

For now, we will call such graphs *random graphs*, although this term has a technical meaning too [52].
### 9.2.1 Properties of random graphs

The graphs we have seen in most of this course have properties that stem from the fact that they model objects in our three-dimensional world. Thus, the typical distance between two nodes is typically $O(N^{1/3})$ where $N$ is the number of nodes. Random graphs do not behave like this: they often have a *small world* property where the typical distance is $O(\log N)$. A famous example is the graph of film actors and their connection by having appeared in the same movie: according to ‘Six degrees of separation’, no two actors have a distance more than six in this graph. In graph terms this means that the diameter of the graph is six.

Small-world graphs have other properties, such as the existence of cliques (although these feature too in higher order FEM problems) and hubs: nodes of a high degree. This leads to implications such as the following: deleting a random node in such a graph does not have a large effect on shortest paths.

**Exercise 9.8.** Considering the graph of airports and the routes that exist between them. If there are only hubs and non-hubs, argue that deleting a non-hub has no effect on shortest paths between other airports. On the other hand, consider the nodes ordered in a two-dimensional grid, and delete an arbitrary node. How many shortest paths are affected?

## 9.3 Hypertext algorithms

There are several algorithms based on linear algebra for measuring the importance of web sites [110]. We will briefly define a few and discuss computational implications.

### 9.3.1 HITS

In the HITS (Hypertext-Induced Text Search) algorithm, sites have a *hub* score that measures how many other sites it points to, and an *authority* score that measures how many sites point to it. To calculate such scores we define an *incidence matrix* $L$, where

$$L_{ij} = \begin{cases} 1 & \text{document } i \text{ points to document } j \\ 0 & \text{otherwise} \end{cases}$$

The authority scores $x_i$ are defined as the sum of the hub scores $y_j$ of everything that points to $i$, and the other way around. Thus

$$x = L^t y$$

$$y = Lx$$

or $x = LL^t x$ and $y = L^t Ly$, showing that this is an eigenvalue problem. The eigenvector we need has only nonnegative entries; this is known as the *Perron vector* for a *nonnegative matrix*, see appendix 12.4. The Perron vector is computed by a *power method*; see section 12.3.

A practical search strategy is:

* Find all documents that contain the search terms;
* Build the subgraph of these documents, and possible one or two levels of documents related to them;

* Compute authority and hub scores on these documents, and present them to the user as an ordered list.

### 9.3.2 PageRank

The PageRank [133] basic idea is similar to HITS. Again we define a connectivity matrix

$$M_{ij} = \begin{cases} 1 & \text{if page } j \text{ links to } i \\ 0 & \text{otherwise} \end{cases}$$

With $e = (1, \dots, 1)$, the vector $d^t = e^t M$ counts how many links there are on a page: $d_i$ is the number of links on page $i$. We construct a diagonal matrix $D = \text{diag}(d_1, \dots)$ we normalize $M$ to $T = MD^{-1}$.

Now the columns sums (that is, the sum of the elements in any column) of $T$ are all 1, which we can express as $e^t T = e^t$ where $e^t = (1, \dots, 1)$. Such a matrix is called stochastic matrix. It has the following interpretation:

If $p$ is a probability vector, that is, $p_i$ is the probability that the user is looking at page $i$, then $Tp$ is the probability vector after the user has clicked on a random link.

**Exercise 9.9.** Mathematically, a probability vector is characterized by the fact that the sum of its elements is 1. Show that the product of a stochastic matrix and a probability vector is indeed again a probability vector.

The PageRank algorithm now does a mathematical investigation of the question ‘if the user keeps clicking on links that are somehow the most desirable on a page, what will overall be the set of the most desirable links’. In principle this would correspond to taking an arbitrary stochastic vector $p$, computing the power method $Tp, T^2p, T^3p, \dots$ and seeing if that sequence converges to something.

There is, however, a small problem with this algorithm. If a page has no outgoing links, a user that would wind up there would never leave. Mathematically this problem also comes up. Consider for instance an web with only 2 pages and the following adjacency matrix:

$$A = \begin{pmatrix} 1/2 & 0 \\ 1/2 & 1 \end{pmatrix}.$$

Check for yourself that this corresponds to the second page having no outgoing links. Now let $p$ be the starting vector $p^t = (1, 1)$, and compute a few iterations of the power method. Do you see that the probability of the user being on the second page goes up to 1? The problem here is that we are dealing with a reducible matrix.

To prevent this problem, PageRank introduces another element: sometimes the user will get bored from clicking, and will go to an arbitrary page (there are also provisions for pages with no outgoing links). If we call $s$ the chance that the user will click on a link, then the chance of going to an arbitrary page is $1 - s$. Together, we now have the process

$$p' \leftarrow sTp + (1-s)e,$$
that is, if $p$ is a vector of probabilities then $p'$ is a vector of probabilities that describes where the user is after making one page transition, either by clicking on a link or by 'teleporting'.

The PageRank vector is the stationary point of this process; you can think of it as the probability distribution after the user has made infinitely many transitions. The PageRank vector satisfies

$$p = sTp + (1-s)e \Leftrightarrow (I - sT)p = (1-s)e.$$

Thus, we now have to wonder whether $I - sT$ has an inverse. If the inverse exists it satisfies

$$(I - sT)^{-1} = I + sT + s^2T^2 + \dots$$

It is not hard to see that the inverse exists: with the Gershgorin theorem (appendix 12.5) you can see that the eigenvalues of $T$ satisfy $|\lambda| \le 1$. Now use that $s < 1$, so the series of partial sums converges.

The above formula for the inverse also indicates a way to compute the PageRank vector $p$ by using a series of matrix-vector multiplications.

**Exercise 9.10.** Write pseudo-code for computing the PageRank vector, given the matrix $T$. Show that you never need to compute the powers of $T$ explicitly. (This is an instance of Horner's rule).

In the case that $s = 1$, meaning that we rule out teleportation, the PageRank vector satisfies $p = Tp$, which is again the problem of finding the Perron vector; see appendix 12.4.

We find the Perron vector by a power iteration (section 12.3)

$$p^{(i+1)} = Tp^{(i)}.$$

This is a sparse matrix vector product, but unlike in the BVP case the sparsity is unlikely to have a structure such as bandedness. Computationally, one probably has to use the same parallelism arguments as for a dense matrix: the matrix has to be distributed two-dimensionally [130].

## 9.4 Large-scale computational graph theory

In the preceding sections you have seen that many graph algorithms have a computational structure that makes the matrix-vector product their most important kernel. Since most graphs are of low degree relative to the number of nodes, the product is a sparse matrix-vector product.

In many cases we can then, as in section 6.5, make a one-dimensional distribution of the matrix, induced by a distribution of the graph nodes: if a processor owns a graph node $i$, it owns all the edges $i, j$.

However, often the computation is very unbalanced. For instance, in the single-source shortest path algorithm only the vertices along the front are active. For this reason, sometimes a distribution of edges rather than vertices makes sense. For an even balancing of the load even random distributions can be used.

The parallelization of the Floyd-Warshall algorithm (section 9.1.3) proceeds along different lines. Here we don't compute a quantity per node, but a quantity that is a function of pairs of nodes, that is, a matrix-like quantity. Thus, instead of distributing the nodes, we distribute the pair distances.
Sometimes we need to go beyond a one-dimensional decomposition in dealing with sparse graph matrices. Let's assume that we're dealing with a graph that has a structure that is more or less random, for instance in the sense that the chance of there being an edge is the same for any pair of vertices. Also assuming that we have a large number of vertices and edges, every processor will then store a certain number of vertices. The conclusion is then that the chance that any pair of processors needs to exchange a message is the same, so the number of messages is $O(P)$. (Another way of visualizing this is to see that nonzeros are randomly distributed through the matrix.) This does not give a scalable algorithm.

The way out is to treat this sparse matrix as a dense one, and invoke the arguments from section 6.2.2 to decide on a two-dimensional distribution of the matrix. (See [169] for an application to the BFS problem; they formulate their algorithm in graph terms, but the structure of the 2D matrix-vector product is clearly recognizable.) The two-dimensional product algorithm only needs collectives in processor rows and columns, so the number of processors involved is $O(\sqrt{P})$.
# Chapter 10

## N-body problems

In chapter 4 we looked at continuous phenomena, such as the behaviour of a heated rod in the entire interval [0, 1] over a certain time period. There are also applications where you may be interested in a finite number of points. One such application is the study of collections of particles, possibly very big particles such as planets or stars, under the influence of a force such as gravity or the electrical force. (There can also be external forces, which we will ignore; also we assume there are no collisions, otherwise we need to incorporate nearest-neighbour interactions.) This type of problems is known as N-body problems; for an introduction see http://www.scholarpedia.org/article/N-body_simulations_(gravitational).

A basic algorithm for this problem is easy enough:

* choose some small time interval,

* calculate the forces on each particle, given the locations of all the particles,

* move the position of each particle as if the force on it stays constant throughout that interval.

For a small enough time interval this algorithm gives a reasonable approximation to the truth.

The last step, updating the particle positions, is easy and completely parallel: the problem is in evaluating the forces. In a naive way this calculation is simple enough, and even completely parallel:

for each particle i

for each particle j

let $\bar{r}_{ij}$ be the vector between i and j;

then the force on i because of j is

$$f_{ij} = -\bar{r}_{ij} \frac{m_i m_j}{|r_{ij}|}$$

(where $m_i, m_j$ are the masses or charges) and

$$f_{ji} = -f_{ij}.$$

The main objection to this algorithm is that it has quadratic computational complexity: for N particles, the number of operations is $O(N^2)$.

<figure><img src="image_125.png" /><figcaption>Figure 10.1: Summing all forces on a particle</figcaption></figure>
**Exercise 10.1.** If we had *N* processors, the computations for one update step would take time *O*(*N*). What is the communication complexity? Hint: is there a collective operations you can use?

Several algorithms have been invented to get the sequential complexity down to $O(N \log N)$ or even $O(N)$. As might be expected, these are harder to implement than the naive algorithm. We will discuss a popular method: the *Barnes-Hut algorithm* [5], which has $O(N \log N)$ complexity.

## 10.1 The Barnes-Hut algorithm

The basic observation that leads to a reduction in complexity is the following. If you are calculating the forces on two particles $i_1, i_2$ that are close together, coming from two particles $j_1, j_2$ that are also close together, you can clump $j_1, j_2$ together into one particle, and use that for both $i_1, i_2$.

Next, the algorithm uses a recursive division of space, in two dimensions in quadrants and in three dimensions in octants; see figure 10.2.

<figure><img src="image_126.png" /><figcaption>Figure 10.2: Recursive subdivision of a domain in quadrants with levels indicated (left); actual subdivision with one particle per box (right)</figcaption></figure>

The algorithm is then as follows. First total mass and center of mass are computed for all cells on all levels:

* for each level $l$, from fine to coarse:

* for each cell $c$ on level $l$:
  * compute the total mass and center of mass for cell $c$
  * for cell $c$ by considering its children
* if there are no particles in this cell,
  * set its mass to zero

Then the levels are used to compute the interaction with each particle:
for each particle $p$:

for each cell $c$ on the top level

if $c$ is far enough away from $p$:

use the total mass and center of mass of $c$;

otherwise consider the children of $c$

The test on whether a cell is far enough away is typically implemented as the ratio of its diameter to its distance being small enough. This is sometimes referred to as the ‘cell opening criterion’. In this manner, each particle interacts with a number of concentric rings of cells, each next ring of double width; see figure 10.3.

This algorithm is easy to realize if the cells are organized in a tree. In the three-dimensional case, each cell has eight children, so this is known as an octtree.

The computation of centres of masses has to be done each time after the particles move. Updating can be less expensive than computing from scratch. Also, it can happen that a particle crosses a cell border, in which case the data structure needs to be updated. In the worst case, a particle moves into a cell that used to be empty.

<figure><img src="image_127.png" /><figcaption>Figure 10.3: Boxes with constant distance/diameter ratio</figcaption></figure>

## 10.2 The Fast Multipole Method

The *Fast Multipole Method (FMM)* computes an expression for the potential at every point, not the force as does Barnes-Hut. FMM uses more information than the mass and center of the particles in a box. This more complicated expansion is more accurate, but also more expensive. In compensation, the FMM uses a fixed set of boxes to compute the potential, rather than a set varying with the accuracy parameter $\theta$, and location of the center of mass.

However, computationally the FMM is much like the Barnes-Hut method so we will discuss their implementation jointly.

## 10.3 Full computation

Despite the above methods for judicious approximation, there are also efforts at full calculation of the $N^2$ interactions; see for instance the NBODY6 code of Sverre Aarseth; see http://www.ast.cam.ac.uk/~sverre/web/pages/home.htm. Such codes use high order integrators and adaptive time steps. Fast implementation on the *Grape computer* exist; general parallelization is typically hard because of the need for regular load balancing.
## 10.4 Implementation

Octtree methods offer some challenges on high performance architectures. First of all, the problem is irregular, and secondly, the irregularity dynamically changes. The second aspect is mostly a problem in distributed memory, and it needs load rebalancing; see section 2.11.1. In this section we concentrated on the force calculation in a single step.

### 10.4.1 Vectorization

The structure of a problem as in figure 10.2 is quite irregular. This is a problem for vectorization on the small scale of SSE/AVX instructions and on the large scale of vector pipeline processors (see section 2.4.1 for an explanation of both). Program steps ‘for all children of a certain box do something’ will be of irregular length, and data will possibly be not stored in a regular manner.

This problem can be alleviated by subdividing the grid even if this means having empty boxes. If the bottom level is fully divided, there will always be eight (in three dimension) particles to operate on. Higher levels can also be filled in, but this means an increasing number of empty boxes on the lower levels, so there is a trade-off between increased work and increasing efficiency.

### 10.4.2 Shared memory implementation

Executed on a sequential architecture, this algorithm has complexity $O(N \log N)$. It is clear that this algorithm will also work on shared memory if each particle is turned into a task. Since not all cells contain particles, tasks will have a different running time.

### 10.4.3 Distributed memory implementation

The above shared-memory version of the Barnes-Hut algorithm can not immediately be used in a distributed memory context, since each particle can in principle access information from any part of the total data. It is possible to realize an implementation along these lines using a hashed octtree, but we will not pursue this.

We observe data access is more structured than it seems at first. Consider a particle $p$ and the cells on level $\ell$ that it interacts with. Particles located close to $p$ will interact with the same cells, so we can rearrange interaction by looking at cells on level $\ell$ and the other cells on the same level that they interact with.

This gives us the following algorithm [97]: the calculation of centres of mass become a calculation of the force $g_p^{(\ell)}$ exerted by a particle $p$ on level $\ell$:

* for level $\ell$ from one above the finest to the coarsest:

* for each cell $c$ on level $\ell$

* let $g_c^{(\ell)}$ be the combination of the $g_i^{(\ell+1)}$ for all children $i$ of $c$

With this we compute the force on a cell:

* for level $\ell$ from one below the coarsest to the finest:

* for each cell $c$ on level $\ell$:

* let $f_c^{(\ell)}$ be the sum of
1. the force $f_p^{(\ell-1)}$ on the parent $p$ of $c$, and

2. the sums $g_i^{(\ell)}$ for all $i$ on level $\ell$ that satisfy the cell opening criterium

We see that on each level, each cell now only interacts with a small number of neighbours on that level. In the first half of the algorithm we go up the tree using only parent-child relations between cells. Presumably this is fairly easy.

The second half of the algorithm uses more complicated data access. The cells $i$ in the second term are all at some distance from the cell $c$ on which we are computing the force. In graph terms these cells can be described as cousins: children of a sibling of $c$'s parent. If the opening criterium is made sharper, we use second cousins: grandchildren of the sibling of $c$'s grandparent, et cetera.

**Exercise 10.2.** Argue that this force calculation operation has much in common, structurally, with the sparse matrix-vector product.

In the shared memory case we already remarked that different subtrees take different time to process, but, since we are likely to have more tasks than processor cores, this will all even out. With distributed memory we lack the possibility to assign work to arbitrary processors, so we need to assign load carefully. Space-Filling Curves (SFCs) can be used here to good effect (see section 2.11.1.5).
# Chapter 11

## Monte Carlo Methods

Monte Carlo simulation is a broad term for methods that use random numbers and statistical sampling to solve problems, rather than exact modeling. From the nature of this sampling, the result will have some uncertainty, but the statistical ‘law of large numbers’ will ensure that the uncertainty goes down as the number of samples grows.

### 11.1 Motivation

Let's start with a simple example: measuring an area, for instance, $\pi$ is the area of a circle inscribed in a square with sides 2. If you picked a random point in the square, the chance of it falling in the circle is $\pi/4$, so you could estimate this ratio by taking many random points $(x, y)$ and seeing in what proportion their length $\sqrt{x^2 + y^2}$ is less than 1.

You could even do this as a physical experiment: suppose you have a pond of an irregular shape in your backyard, and that the yard itself is rectangular with known dimensions. If you would now throw pebbles into your yard so that they are equally likely to land at any given spot, then the ratio of pebbles falling in the pond to those falling outside equals the ratio of the areas.

Less fanciful and more mathematically, we need to formalize the idea of falling inside or outside the shape you are measuring. Therefore, let $\Omega \in [0, 1]^2$ be the shape, and let a function $f(\bar{x})$ describe the boundary of $\Omega$, that is

$$ \begin{cases} f(\bar{x}) < 0 & x \notin \Omega \\ f(\bar{x}) > 0 & x \in \Omega \end{cases} $$

Now take random points $\bar{x}_0, \bar{x}_1, \bar{x}_2 \in [0, 1]^2$, then we can estimate the area of $\Omega$ by counting how often $f(\bar{x}_i)$ is positive or negative.

We can extend this idea to integration. The average value of a function on an interval $(a, b)$ is defined as

$$ \langle f \rangle = \frac{1}{b-a} \int_a^b f(x)dx $$
On the other hand, we can also estimate the average as

$$ \langle f \rangle \approx \frac{1}{N} \sum_{i=1}^{n} f(x_i) $$

if the points $x_i$ are reasonably distributed and the function $f$ is not too wild. This leads us to

$$ \int_{a}^{b} f(x) dx \approx (b - a) \frac{1}{N} \sum_{i=1}^{n} f(x_i) $$

Statistical theory, that we will not go into, tells us that the uncertainty $\sigma_I$ in the integral is related to the standard deviation $\sigma_f$ by

$$ \sigma_I \sim \frac{1}{\sqrt{N}} \sigma_f $$

for normal distributions.

### 11.1.1 What is the attraction?

So far, Monte Carlo integration does not look much different from classical integration by *Riemann sums*. The difference appears when we go to higher dimensions. In that case, for classical integration we would need $N$ points in each dimension, leading to $N^d$ points in $d$ dimensions. In the Monte Carlo method, on the other hand, the points are taken at random from the $d$-dimensional space, and a much lower number of points suffices.

Computationally, Monte Carlo methods are attractive since all function evaluations can be performed in parallel.

The statistical law that underlies this is as follows: if $N$ independent observations are made of a quantity with standard deviation $\sigma$, then the standard deviation of the mean is $\sigma/\sqrt{N}$. This means that more observations will lead to more accuracy; what makes Monte Carlo methods interesting is that this gain in accuracy is not related to dimensionality of the original problem.

Monte Carlo techniques are of course natural candidatates for simulating phenomena that are statistical in nature, such as radioactive decay, or Brownian motion. Other problems where Monte Carlo simulation is attractive are outside the realm of scientific computing. For instance, the *Black-Scholes model* for stock option pricing [11] uses Monte Carlo simulation.

Some problems that you have seen before, such as solving a linear system of equations, can be tackled with Monte Carlo techniques. However, this is not a typical application. Below we will discuss two applications where exact methods would take far too much time to compute and where statistical sampling can quickly give a reasonably accurate answer.

An important tool for statistical sampling is a random number generator. We start by briefly discussing the problems in generating random numbers, especially in parallel.