<table><thead><tr><th>computation</th><th>correct result</th></tr></thead><tbody><tr><td>y<sub>0</sub> = ln 6 - ln 5 = .182|322 × 10<sup>1</sup>...</td><td>1.82</td></tr><tr><td>y<sub>1</sub> = .900 × 10<sup>-1</sup></td><td>.884</td></tr><tr><td>y<sub>2</sub> = .500 × 10<sup>-1</sup></td><td>.0580</td></tr><tr><td>y<sub>3</sub> = .830 × 10<sup>-1</sup></td><td>going up? .0431</td></tr><tr><td>y<sub>4</sub> = -.165</td><td>negative? .0343</td></tr></tbody></table>

We see that the computed results are quickly not just inaccurate, but actually nonsensical. We can analyze why this is the case.

If we define the error $\epsilon_n$ in the n-th step as

$$ \tilde{y}_n - y_n = \epsilon_n, $$

then

$$ \tilde{y}_n = 1/n - 5\tilde{y}_{n-1} = 1/n + 5n_{n-1} + 5\epsilon_{n-1} = y_n + 5\epsilon_{n-1} $$

so $\epsilon_n \ge 5\epsilon_{n-1}$. The error made by this computation shows exponential growth.

#### 3.3.6.4 Linear system solving

Sometimes we can make statements about the numerical precision of a problem even without specifying what algorithm we use. Suppose we want to solve a linear system, that is, we have an $n \times n$ matrix $A$ and a vector $b$ of size $n$, and we want to compute the vector $x$ such that $Ax = b$. (We will actually considering algorithms for this in chapter 5.) Since the vector $b$ will the result of some computation or measurement, we are actually dealing with a vector $\tilde{b}$, which is some perturbation of the ideal $b$:

$$ \tilde{b} = b + \Delta b. $$

The perturbation vector $\Delta b$ can be of the order of the machine precision if it only arises from representation error, or it can be larger, depending on the calculations that produced $\tilde{b}$.

We now ask what the relation is between the exact value of $x$, which we would have obtained from doing an exact calculation with $A$ and $b$, which is clearly impossible, and the computed value $\tilde{x}$, which we get from computing with $A$ and $\tilde{b}$. (In this discussion we will assume that $A$ itself is exact, but this is a simplification.)

Writing $\tilde{x} = x + \Delta x$, the result of our computation is now

$$ A\tilde{x} = \tilde{b} $$

or

$$ A(x + \Delta x) = b + \Delta b. $$

Since $Ax = b$, we get $A\Delta x = \Delta b$. From this, we get (see appendix 12 for details)

$$ \left\{ \begin{array}{l} \Delta x = A^{-1}\Delta b \\ Ax = b \end{array} \right\} \Rightarrow \left\{ \begin{array}{l} \|A\|\|x\| \ge \|b\| \\ \|\Delta x\| \le \|A^{-1}\|\|\Delta b\| \end{array} \right\} \Rightarrow \frac{\|\Delta x\|}{\|x\|} \le \|A\|\|A^{-1}\| \frac{\|\Delta b\|}{\|b\|} \quad (3.2) $$
The quantity $||A||||A^{-1}||$ is called the condition number of a matrix. The bound (3.2) then says that any perturbation in the right hand side can lead to a perturbation in the solution that is at most larger by the condition number of the matrix $A$. Note that it does not say that the perturbation in $x$ needs to be anywhere close to that size, but we can not rule it out, and in some cases it indeed happens that this bound is attained.

Suppose that $b$ is exact up to machine precision, and the condition number of $A$ is $10^4$. The bound (3.2) is often interpreted as saying that the last 4 digits of $x$ are unreliable, or that the computation ‘loses 4 digits of accuracy’.

Equation (3.2) can also be interpreted as follows: when we solve a linear system $Ax = b$ we get an approximate solution $x + \Delta x$ which is the exact solution of a perturbed system $A(x + \Delta x) = b + \Delta b$. The fact that the perturbation in the solution can be related to the perturbation in the system, is expressed by saying that the algorithm exhibits backwards stability.

The analysis of the accuracy of linear algebra algorithms is a field of study in itself; see for instance the book by Higham [86].

### 3.3.7 Roundoff error in parallel computations

As we discussed in section 3.3.5, and as you saw in the above example of summing a series, addition in computer arithmetic is not associative. A similar fact holds for multiplication. This has an interesting consequence for parallel computations: the way a computation is spread over parallel processors influences the result.

As a very simple example, consider computing the sum of four numbers $a+b+c+d$. On a single processor, ordinary execution corresponds to the following associativity:

$$((a + b) + c) + d.$$

On the other hand, spreading this computation over two processors, where processor 0 has $a, b$ and processor 1 has $c, d$, corresponds to

$$((a + b) + (c + d)).$$

Generalizing this, we see that reduction operations will most likely give a different result on different numbers of processors. (The MPI standard declares that two program runs on the same set of processors should give the same result.) It is possible to circumvent this problem by replace a reduction operation by a gather operation to all processors, which subsequently do a local reduction. However, this increases the memory requirements for the processors.

There is an intriguing other solution to the parallel summing problem. If we use a mantissa of 4000 bits to store a floating point number, we do not need an exponent, and all calculations with numbers thus stored are exact since they are a form of fixed-point calculation [102, 101]. While doing a whole application with such numbers would be very wasteful, reserving this solution only for an occasional inner product calculation may be the solution to the reproducibility problem.
## 3.4 Compilers and round-off

From the above discussion it should be clear that some simple statements that hold for mathematical real numbers do not hold for floating-point numbers. For instance, in floating-point arithmetic

$$ (a + b) + c \neq a + (b + c). $$

This implies that a compiler can not perform certain optimizations without it having an effect on round-off behaviour<fnref n="7" />. In some codes such slight differences can be tolerated, for instance because the method has built-in safeguards. For instance, the stationary iterative methods of section 5.5 damp out any error that is introduced.

On the other hand, if the programmer has written code to account for round-off behaviour, the compiler has no such liberties. This was hinted at in exercise 3.5 above. We use the concept of *value safety* to describe how a compiler is allowed to change the interpretation of a computation. At its strictest, the compiler is not allowed to make any changes that affect the result of a computation.

These matters are also of importance if you care about *reproducibility* of results. If a code is compiled with two different compilers, should runs with the same input give the same output? If a code is run in parallel on two different processor configurations? These questions are very subtle. In the first case, people sometimes insist on *bitwise reproducibility*, whereas in the second case some differences are allowed, as long as the result stays ‘scientifically’ equivalent. Of course, that concept is hard to make rigorous.

Here are some issues that are relevant when considering the influence of the compiler on code behaviour and reproducibility.

**Re-association** Foremost among changes that a compiler can make to a computation is *re-association*, the technical term for grouping $a + b + c$ as $a + (b + c)$. The C language standard and the C++ language standard prescribe strict left-to-right evaluation of expressions without parentheses, so re-association is in fact not allowed by the standard. The Fortran language standard has no such prescription.

A common source of re-association is *loop unrolling*; see section 1.6.2. Under strict value safety, a compiler is limited in how it can unroll a loop, which has implications for performance. The amount of loop unrolling, and whether it's performed at all, depends on the compiler optimization level, the choice of compiler, and the target platform.

A more subtle source of re-association is *parallel execution*; see section 3.3.7. This implies that the output of a code need not be strictly reproducible between two runs on different parallel configurations.

**Expression evaluation** In evaluating the expression $a + (b + c)$, a processor will generate an intermediate result for $b + c$ which is not assigned to any variable. Many processors are able to assign a higher *precision* of the intermediate result. A compiler can have a flag to dictate whether to use this facility.

<footnote n="7">7. This section borrows from documents by Microsoft http://msdn.microsoft.com/en-us/library/aa289157(vs.71).aspx and Intel http://software.intel.com/sites/default/files/article/164389/fp-consistency-122712_1.pdf; for detailed discussion the reader is referred to these.</footnote>
**Behaviour of the floating point unit** Rounding behaviour and treatment of gradual underflow may be controlled by library functions.

**Library functions** The IEEE 754 standard only prescribes simple operations; there is as yet no standard that treats sine or log functions. Therefore, their implementation may be a source of variability.

## 3.5 More about floating point arithmetic

### 3.5.1 Programming languages

Different languages have different approaches to declaring integers and floating point numbers.

*   **Fortran** In Fortran there are various ways of specifying the storage format for integer and real variables. For instance, it is possible to declare the number of bytes that it takes to store a variable: `INTEGER*2`, `REAL*8`. One advantage of this approach is the easy interoperability with other languages, or the MPI library.

Often it is possible to write a code using only `INTEGER`, `REAL`, and use *compiler flags* to indicate the size of an integer and real number in bytes.

More sophisticated, modern versions of Fortran can indicate the number of digits of precision a floating point number needs to have:

```
integer, parameter :: k9 = selected_real9(9)
real(kind=k9)  :: r
r = 2._k9; print *, sqrt(r) ! prints 1.4142135623730
```

The 'kind' values will usually be 4,8,16 but this is compiler dependent.

*   **C** In C, the type identifiers have no standard length. For integers there is `short int`, `int`, `long int`, and for floating point `float`, `double`. The `sizeof()` operator gives the number of bytes used to store a datatype.

*   **C99, Fortran2003** Recent standards of the C and Fortran languages incorporate the C/Fortran interoperability standard, which can be used to declare a type in one language so that it is compatible with a certain type in the other language.

### 3.5.2 Other computer arithmetic systems

Other systems have been proposed to dealing with the problems of inexact arithmetic on computers. One solution is *extended precision* arithmetic, where numbers are stored in more bits than usual. A common use of this is in the calculation of inner products of vectors: the accumulation is internally performed in extended precision, but returned as a regular floating point number. Alternatively, there are libraries such as `GMPlib` [65] that allow for any calculation to be performed in higher precision.

Another solution to the imprecisions of computer arithmetic is 'interval arithmetic' [91], where for each calculation interval bounds are maintained. While this has been researched for considerable time, it is not practically used other than through specialized libraries [18].
### 3.5.3 Extended precision

When the IEEE 754 standard was drawn up, it was envisioned that processors could have a whole range of precisions. In practice, only single and double precision as defined have been used. However, one instance of extended precision still survives: processors sometimes have 80-bit registers for storing intermediate results. (This goes back to the Intel 80287 co-processor.) This strategy makes sense in FMA instructions, and in the accumulation of inner products.

### 3.5.4 Fixed-point arithmetic

A fixed-point number (for a more thorough discussion than found here, see [168]) can be represented as $⟨N, F⟩$ where $N \ge \beta^0$ is the integer part and $F < 1$ is the fractional part. Another way of looking at this, is that a fixed-point number is an integer stored in $N + F$ digits, with an implied decimal point after the first $N$ digits.

Fixed-point calculations can overflow, with no possibility to adjust an exponent. Consider the multiplication $⟨N_1, F_1⟩ \times ⟨N_2, F_2⟩$, where $N_1 \ge \beta^{n_1}$ and $N_2 \ge \beta^{n_2}$. This overflows if $n_1 + n_2$ is more than the number of positions available for the integer part. (Informally, the number of digits of the product is the sum of the digits of the operands.) This means that, in a program that uses fixed-point, numbers will need to have a number of zero digits, if you are ever going to multiply them, which lowers the numerical accuracy. It also means that the programmer has to think harder about calculations, arranging them in such a way that overflow will not occur, and that numerical accuracy is still preserved to a reasonable extent.

So why would people use fixed-point numbers? One important application is in embedded low-power devices, think a battery-powered digital thermometer. Since fixed-point calculations are essentially identical to integer calculations, they do not require a floating-point unit, thereby lowering chip size and lessening power demands. Also, many early video game systems had a processor that either had no floating-point unit, or where the integer unit was considerably faster than the floating-point unit. In both cases, implementing non-integer calculations as fixed-point, using the integer unit, was the key to high throughput.

Another area where fixed point arithmetic is still used is in signal processing. In modern CPUs, integer and floating point operations are of essentially the same speed, but converting between them is relatively slow. Now, if the sine function is implemented through table lookup, this means that in $sin(\sin x)$ the output of a function is used to index the next function application. Obviously, outputting the sine function in fixed point obviates the need for conversion between real and integer quantities, which simplifies the chip logic needed, and speeds up calculations.

### 3.5.5 Complex numbers

Some programming languages have complex numbers as a native data type, others not, and others are in between. For instance, in Fortran you can declare

```fortran
COMPLEX z1,z2, z(32)
COMPLEX*16 zz1, zz2, zz(36)
```
A complex number is a pair of real numbers, the real and imaginary part, allocated adjacent in memory. The first declaration then uses 8 bytes to store to `REAL*4` numbers, the second one has `REAL*8`s for the real and imaginary part. (Alternatively, use `DOUBLE COMPLEX` or in Fortran90 `COMPLEX(KIND=2)` for the second line.)

By contrast, the C language does not natively have complex numbers, but both C99 and C++ have a `complex.h` header file<fnref n="8" />. This defines as complex number as in Fortran, as two real numbers.

Storing a complex number like this is easy, but sometimes it is computationally not the best solution. This becomes apparent when we look at arrays of complex numbers. If a computation often relies on access to the real (or imaginary) parts of complex numbers exclusively, striding through an array of complex numbers, has a stride two, which is disadvantageous (see section 1.3.4.6). In this case, it is better to allocate one array for the real parts, and another for the imaginary parts.

**Exercise 3.13.** Suppose arrays of complex numbers are stored the Fortran way. Analyze the memory access pattern of pairwise multiplying the arrays, that is, $\forall i: c_i \leftarrow a_i \cdot b_i$, where $a()$, $b()$, $c()$ are arrays of complex numbers.

**Exercise 3.14.** Show that an $n \times n$ linear system $Ax = b$ over the complex numbers can be written as a $2n \times 2n$ system over the real numbers. Hint: split the matrix and the vectors in their real and imaginary parts. Argue for the efficiency of storing arrays of complex numbers as separate arrays for the real and imaginary parts.

## 3.6 Conclusions

Computations done on a computer are invariably beset with numerical error. In a way, the reason for the error is the imperfection of computer arithmetic: if we could calculate with actual real numbers there would be no problem. (There would still be the matter of measurement error in data, and approximations made in numerical methods; see the next chapter.) However, if we accept roundoff as a fact of life, then various observations hold:

* Mathematically equivalent operations need not behave identically from a point of stability; see the ‘abc-formula’ example.

* Even rearrangements of the same computations do not behave identically; see the summing example.

Thus it becomes imperative to analyze computer algorithms with regard to their roundoff behaviour: does roundoff increase as a slowly growing function of problem parameters, such as the number of terms evaluated, or is worse behaviour possible? We will not address such questions in further detail in this book.

<footnote n="8">8. These two header files are not identical, and in fact not compatible. Beware, if you compile C code with a C++ compiler [47].</footnote>
# Chapter 4

## Numerical treatment of differential equations

In this chapter we will look at the numerical solution of Ordinary Differential Equations (ODEs) and Partial Differential Equations (PDEs). These equations are commonly used in physics to describe phenomena such as the flow of air around an aircraft, or the bending of a bridge under various stresses. While these equations are often fairly simple, getting specific numbers out of them ('how much does this bridge sag if there are a hundred cars on it') is more complicated, often taking large computers to produce the desired results. Here we will describe the techniques that turn ODEs and PDEs into computable problems.

First of all, we will look at Initial Value Problems (IVPs), which describes processes that develop in time. Here we only consider ODEs: scalar functions that are only dependent on time. The name derives from the fact that typically the function is specified at an initial time point.

Next, we will look at Boundary Value Problems (BVPs), describing processes in space. In realistic situations, this will concern multiple space variables, so we have a PDE. The name BVP is explained by the fact that the solution is specified on the boundary of the domain of definition.

Finally, we will consider the 'heat equation', an Initial Boundary Value Problem (IBVP) which has aspects of both IVPs and BVPs: it describes heat spreading through a physical object such as a rod. The initial value describes the initial temperature, and the boundary values give prescribed temperatures at the ends of the rod.

Our aim in this chapter is to show the origins of an important class of computational problems. Therefore we will not go into theoretical matters of existence, uniqueness, or conditioning of solutions. For this, see [81] or any book that is specifically dedicated to ODEs or PDEs. For ease of analysis we will also assume that all functions involved have sufficiently many higher derivatives, and that each derivative is sufficiently smooth.

### 4.1 Initial value problems

Many physical phenomena change over time, and typically the laws of physics give a description of the change, rather than of the quantity of interest itself. For instance, Newton's second law

$$F = ma$$
is a statement about the change in position of a point mass: expressed as

$$a = \frac{d^2}{dt^2}x = F/m$$

it states that acceleration depends linearly on the force exerted on the mass. A closed form description $x(t) = ...$ for the location of the mass can sometimes be derived analytically, but in many cases some form of approximation or numerical computation is needed.

Newton's equation is an ODE since it describes a function of one variable, time. It is an IVP since it describes the development in time, starting with some initial conditions. As an ODE, it is 'of second order' since it involves a second derivative, We can reduce this to first order, involving only first derivatives, if we allow vector quantities.

We introduce a two-component vector $u$ that combines the location $x$ and velocity $x'^1$:

$$u(t) = (x(t), x'(t))^t,$$

Newton's equation expressed in $u$ then becomes:

$$u' = Au + B, \quad A = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}, \quad B = \begin{pmatrix} 0 \\ F/a \end{pmatrix}$$

For simplicity, in this course we will only consider scalar equations; our reference equation is then

$$u'(t) = f(t, u(t)), \quad u(0) = u_0, \quad t > 0. \qquad (4.1)$$

Equation (4.1) allows for an explicit time dependence of the process, but in general we only consider equations that do not have this explicit dependence, the so-called 'autonomous' ODEs of the form

$$u'(t) = f(u(t)) \qquad (4.2)$$

in which the right hand side does not explicitly depend on $t^2$.

Typically, the initial value in some starting point (often chosen as $t = 0$) is given: $u(0) = u_0$ for some value $u_0$, and we are interested in the behaviour of $u$ as $t \to \infty$. As an example, $f(x) = x$ gives the equation $u'(t) = u(t)$. This is a simple model for population growth: the equation states that the rate of growth is equal to the size of the population. The equation (4.1) can be solved analytically for some choices of $f$, but we will not consider this. Instead, we only consider the numerical solution and the accuracy of this process.

We will now consider numerical methods for the solution In a numerical method, we consider discrete size time steps to approximate the solution of the continuous time-dependent process. Since this introduces a certain amount of error, we will analyze the error introduced in each time step, and how this adds up to a global error. In some cases, the need to limit the global error will impose restrictions on the numerical scheme.

<footnote n="1">1. We use the prime symbol to indicate differentiation in case of functions of a single variable.</footnote>

<footnote n="2">2. Non-autonomous ODEs can be transformed to autonomous ones, so this is no limitation. If $u = u(t)$ is a scalar function and $f = f(t, u)$, we define $u_2(t) = t$ and consider the equivalent autonomous system ($\begin{pmatrix} u' \\ u_2' \end{pmatrix} = \begin{pmatrix} f(u_2, u) \\ f_1(u_2, u) \end{pmatrix}$)</footnote>
### 4.1.1 Error and stability

Since numerical computation will always involve the inaccuracies stemming from the use of machine arithmetic, we want to avoid the situation where a small perturbation in the initial value leads to large perturbations in the solution. Therefore, we will call a differential equation ‘stable’ if solutions corresponding to different initial values $u_0$ converge to one another as $t \to \infty$.

A sufficient criterium for stability is:

$$ \frac{\partial}{\partial u} f(u) = \begin{cases} > 0 & \text{unstable} \\ = 0 & \text{neutrally stable} \\ < 0 & \text{stable} \end{cases} $$

Proof. If $u^*$ is a zero of $f$, meaning $f(u^*) = 0$, then the constant function $u(t) \equiv u^*$ is a solution of $u' = f(u)$, a so-called ‘equilibrium’ solution. We will now consider how small perturbations from the equilibrium behave. Let $u$ be a solution of the PDE, and write $u(t) = u^* + \eta(t)$, then we have

$$ \begin{aligned} \eta' &= u' = f(u) = f(u^* + \eta) = f(u^*) + \eta f'(u^*) + O(\eta^2) \\ &= \eta f'(u^*) + O(\eta^2) \end{aligned} $$

Ignoring the second order terms, this has the solution

$$ \eta(t) = e^{f'(x^*)t} $$

which means that the perturbation will damp out if $f'(x^*) < 0$.

We will often refer to the simple example $f(u) = -\lambda u$, with solution $u(t) = u_0 e^{-\lambda t}$. This problem is stable if $\lambda > 0$.

### 4.1.2 Finite difference approximation: Euler explicit and implicit methods

In order to solve the problem numerically, we turn the continuous problem into a discrete one by looking at finite time/space steps. Assuming all functions are sufficiently smooth, a straightforward Taylor series expansion<fnref n="3" /> gives:

$$ u(t + \Delta t) = u(t) + u'(t)\Delta t + u''(t)\frac{\Delta t^2}{2!} + u'''(t)\frac{\Delta t^3}{3!} + \dots $$

This gives for $u'$:

$$ \begin{aligned} u'(t) &= \frac{u(t+\Delta t)-u(t)}{\Delta t} + \frac{1}{\Delta t} \left( u''(t)\frac{\Delta t^2}{2!} + u'''(t)\frac{\Delta t^3}{3!} + \dots \right) \\ &= \frac{u(t+\Delta t)-u(t)}{\Delta t} + \frac{1}{\Delta t} O(\Delta t^2) \\ &= \frac{u(t+\Delta t)-u(t)}{\Delta t} + O(\Delta t) \end{aligned} \quad (4.3) $$

We can approximate the infinite sum of higher derivatives by a single $O(\Delta t^2)$ if all derivatives are bounded; alternatively, appendix 15 shows that this sum is equal to $\Delta t^2 u''(t + \alpha \Delta t)$ with $0 < \alpha < 1$.

<footnote n="3">3. See appendix 15 if you are unfamiliar with this.</footnote>
We see that we can approximate a differential operator by a finite difference, with an error that is known in its order of magnitude as a function of the time step.

Substituting this in $u' = f(t, u)$ gives<fnref n="4" />

$$ \frac{u(t + \Delta t) - u(t)}{\Delta t} = f(t, u(t)) + O(\Delta t) $$

or

$$ u(t + \Delta t) = u(t) + \Delta t f(t, u(t)) + O(\Delta t^2) $$

We use this equation to derive a numerical scheme: with $t_0 = 0$, $t_{k+1} = t_k + \Delta t = \cdots = (k + 1)\Delta t$, we get a difference equation

$$ u_{k+1} = u_k + \Delta t f(t_k, u_k) $$

for $u_k$ quantities, and we hope that $u_k$ will be a good approximation to $u(t_k)$. This is known as the ‘Explicit Euler’ or ‘Euler forward’ method.

The process of going from a differential equation to a difference equation is often referred to as *discretization*, since we compute function values only in a discrete set of points. The values computed themselves are still real valued. Another way of phrasing this: the numerical solution is found in the finite dimensional space $\mathbb{R}^k$ if we compute $k$ time steps. The solution to the original problem is found in the space of functions $\mathbb{R} \to \mathbb{R}$.

In (4.3) we approximated one operator by another, and in doing so made a truncation error of order $O(\Delta t)$ as $\Delta t \downarrow 0$ (see appendix 13 for a more formal introduction to this notation for orders of magnitude.). This does *not* immediately imply that the difference equation computes a solution that is close to the true solution. For that some more analysis is needed.

We start by analyzing the ‘local error’: if we assume the computed solution is exact at step $k$, that is, $u_k = u(t_k)$, how wrong will we be at step $k + 1$? We have

$$ \begin{aligned} u(t_{k+1}) &= u(t_k) + u'(t_k)\Delta t + u''(t_k)\frac{\Delta t^2}{2!} + \cdots \\ &= u(t_k) + f(t_k, u(t_k))\Delta t + u''(t_k)\frac{\Delta t^2}{2!} + \cdots \end{aligned} $$

and

$$ u_{k+1} = u_k + f(t_k u_k) \Delta t $$

So

$$ \begin{aligned} L_{k+1} &= u_{k+1} - u(t_{k+1}) = u_k - u(t_k) + f(t_k, u_k) - f(t_k, u(t_k)) - u''(t_k)\frac{\Delta t^2}{2!} + \cdots \\ &= -u''(t_k)\frac{\Delta t^2}{2!} + \cdots \end{aligned} $$

<footnote n="4">4. The following equation is a mathematical equality, and should not be interpreted as a way of computing $u'$ for a given function $u$. Recalling the discussion in section 3.3.4 you can see that this formula would quickly lead to cancellation for small $\Delta t$. For a discussion of numerical differentiation, see a numerical analysis textbook.</footnote>
This shows that in each step we make an error of $O(\Delta t^2)$. If we assume that these errors can be added, we find a global error of

$$E_k \approx \Sigma_k L_k = k \Delta t \frac{\Delta t^2}{2!} = O(\Delta t)$$

Since the global error is of first order in $\Delta t$, we call this a ‘first order method’. Note that this error, which measures the distance between the true and computed solutions, is of the same order $O(\Delta t)$ as the truncation error, which is the error in approximating the operator.

### 4.1.2.1 Stability of the Euler explicit method

Consider the IVP $u' = f(t, u)$ for $t \ge 0$, where $f(t, u) = -\lambda u$ and an initial value $u(0) = u_0$ is given. This has an exact solution of $u(t) = u_0 e^{-\lambda t}$. From the above discussion, we conclude that this problem is stable, meaning that small perturbations in the solution ultimately damp out, if $\lambda > 0$. We will now investigate the question of whether the numerical solution behaves the same way as the exact solution, that is, whether numerical solutions also converge to zero.

The Euler forward, or explicit Euler, scheme for this problem is

$$u_{k+1} = u_k - \Delta t \lambda u_k = (1 - \lambda \Delta t) u_k \Rightarrow u_k = (1 - \lambda \Delta t)^k u_0.$$

For stability, we require that $u_k \to 0$ as $k \to \infty$. This is equivalent to

$$\begin{align*}
u_k \downarrow 0 &\Leftrightarrow |1 - \lambda \Delta t| < 1 \\
&\Leftrightarrow -1 < 1 - \lambda \Delta t < 1 \\
&\Leftrightarrow -2 < -\lambda \Delta t < 0 \\
&\Leftrightarrow 0 < \lambda \Delta t < 2 \\
&\Leftrightarrow \Delta t < 2/\lambda
\end{align*}$$

We see that the stability of the numerical solution scheme depends on the value of $\Delta t$: the scheme is only stable if $\Delta t$ is small enough. For this reason, we call the explicit Euler method conditionally stable. Note that the stability of the differential equation and the stability of the numerical scheme are two different questions. The continuous problem is stable if $\lambda > 0$; the numerical problem has an additional condition that depends on the discretization scheme used.

Note that the stability analysis we just performed was specific to the differential equation $u' = -\lambda u$. If you are dealing with a different IVP you have to perform a separate analysis. However, you will find that explicit methods typically give conditional stability.

### 4.1.2.2 The Euler implicit method

The explicit method you just saw was easy to compute, but the conditional stability is a potential problem. For instance, it could imply that the number of time steps would be a limiting factor. There is an alternative to the explicit method that does not suffer from the same objection.
Instead of expanding $u(t + \Delta t)$, consider the following expansion of $u(t - \Delta t)$:

$$u(t - \Delta t) = u(t) - u'(t)\Delta t + u''(t)\frac{\Delta t^2}{2!} + \dots$$

which implies

$$u'(t) = \frac{u(t) - u(t - \Delta t)}{\Delta t} + u''(t)\Delta t/2 + \dots$$

As before, we take the equation $u'(t) = f(t, u(t))$ and approximate $u'(t)$ by a difference formula:

$$\frac{u(t) - u(t - \Delta t)}{\Delta t} = f(t, u(t)) + O(\Delta t) \Rightarrow u(t) = u(t - \Delta t) + \Delta t f(t, u(t)) + O(\Delta t^2)$$

Again we define fixed points $t_k = kt$, and we define a numerical scheme:

$$u_{k+1} = u_k + \Delta t f(t_{k+1}, u_{k+1})$$

where $u_k$ is an approximation of $u(t_k)$.

An important difference with the explicit scheme is that $u_{k+1}$ now also appears on the right hand side of the equation. That is, computation of $u_{k+1}$ is now implicit. For example, let $f(t, u) = -u^3$, then $u_{k+1} = u_k - \Delta t u_{k+1}^3$. In other words, $u_{k+1Z}$ is the solution for $x$ of the equation $\Delta t x^3 + x = u_k$. This is a nonlinear equation, which typically can be solved using the Newton method.

### 4.1.2.3 Stability of the implicit Euler method

Let us take another look at the example $f(t, u) = -\lambda u$. Formulating the implicit method gives

$$u_{k+1} = u_k - \lambda \Delta t u_{k+1} \Leftrightarrow (1 + \lambda \Delta t) u_{k+1} = u_k$$

so

$$u_{k+1} = \left( \frac{1}{1 + \lambda \Delta t} \right) u_k \Rightarrow u_k = \left( \frac{1}{1 + \lambda \Delta t} \right)^k u_0.$$

If $\lambda > 0$, which is the condition for a stable equation, we find that $u_k \to 0$ for all values of $\lambda$ and $\Delta t$. This method is called *unconditionally stable*. One advantage of an implicit method over an explicit one is clearly the stability: it is possible to take larger time steps without worrying about unphysical behaviour. Of course, large time steps can make convergence to the *steady state* (see Appendix 14.4) slower, but at least there will be no divergence.

On the other hand, implicit methods are more complicated. As you saw above, they can involve nonlinear systems to be solved in every time step. In cases where $u$ is vector-valued, such as in the heat equation, discussed below, you will see that the implicit method requires the solution of a system of equations.
Exercise 4.1. Analyse the accuracy and computational aspects of the following scheme for the IVP $u'(x) = f(x)$:

$$u_{i+1} = u_i + h(f(x_i) + f(x_{i+1}))/2$$

which corresponds to adding the Euler explicit and implicit schemes together. You do not have to analyze the stability of this scheme.

Exercise 4.2. Consider the initial value problem $y'(t) = y(t)(1 - y(t))$. Observe that $y \equiv 0$ and $y \equiv 1$ are solutions. These are called ‘equilibrium solutions’.

1. A solution is stable, if perturbations ‘converge back to the solution’, meaning that for $\epsilon$ small enough,

$$\text{if } y(t) = \epsilon \text{ for some } t, \text{ then } \lim_{t \to \infty} y(t) = 0$$

and

$$\text{if } y(t) = 1 + \epsilon \text{ for some } t, \text{ then } \lim_{t \to \infty} y(t) = 1$$

This requires for instance that

$$y(t) = \epsilon \Rightarrow y'(t) < 0.$$

Investigate this behaviour. Is zero a stable solution? Is one?

2. Consider the explicit method

$$y_{k+1} = y_k + \Delta t y_k (1 - y_k)$$

for computing a numerical solution to the differential equation. Show that

$$y_k \in (0, 1) \Rightarrow y_{k+1} > y_k, \quad y_k > 1 \Rightarrow y_{k+1} < y_k$$

3. Write a small program to investigate the behaviour of the numerical solution under various choices for $\Delta t$. Include program listing and a couple of runs in your homework submission.

4. You see from running your program that the numerical solution can oscillate. Derive a condition on $\Delta t$ that makes the numerical solution monotone. It is enough to show that $y_k < 1 \Rightarrow y_{k+1} < 1$, and $y_k > 1 \Rightarrow y_{k+1} > 1$.

5. Now consider the implicit method

$$y_{k+1} - \Delta t y_{k+1} (1 - y_{k+1}) = y_k$$

and show that $y_{k+1}$ is easily computed from $y_k$. Write a program, and investigate the behaviour of the numerical solution under various choices for $\Delta t$.

6. Show that the numerical solution of the implicit scheme is monotone for all choices of $\Delta t$.
## 4.2 Boundary value problems

In the previous section you saw initial value problems, which model phenomena that evolve over time. We will now move on to ‘boundary value problems’, which are in general stationary in time, but which describe a phenomenon that is location dependent. Examples would be the shape of a bridge under a load, or the heat distribution in a window pane, as the temperature outside differs from the one inside.

The general form of a (second order, one-dimensional) BVP is<fnref n="5" />

$$u''(x) = f(x, u, u') \text{ for } x \in [a, b] \text{ where } u(a) = u_a, u(b) = u_b$$

but here we will only consider the simple form

$$-u''(x) = f(x) \text{ for } x \in [0, 1] \text{ with } u(0) = u_0, u(1) = u_1. \qquad (4.4)$$

in one space dimension, or

$$-u_{xx}(\bar{x}) - u_{yy}(\bar{x}) = f(\bar{x}) \text{ for } \bar{x} \in \Omega = [0, 1]^2 \text{ with } u(\bar{x}) = u_0 \text{ on } \delta\Omega. \qquad (4.5)$$

in two space dimensions. Here, $\delta\Omega$ is the boundary of the domain $\Omega$. Since we prescribe the value of $u$ on the boundary, such a problem is called a Boundary Value Problem (BVP).

### 4.2.1 General PDE theory

There are several types of PDE, each with distinct mathematical properties. The most important property is that of region of influence: if we tinker with the problem so that the solution changes in one point, what other points will be affected.

* Hyperbolic PDEs are of the form

$$Au_{xx} + Bu_{yy} + \text{lower order terms} = 0$$

with $A, B$ of opposite sign. Such equations describe waves, or more general convective phenomena, that are conservative, and do not tend to a steady state.

Intuitively, changing the solution of a wave equation at any point will only change certain future points, since waves have a propagation speed that makes it impossible for a point to influence points in the near future that are too far away in space. This type of PDE will not be discussed in this book.

* Parabolic PDEs are of the form

$$Au_x + Bu_{yy} + \text{no higher order terms in } x = 0$$

and they describe diffusion-like phenomena; these often tend to a steady state. The best way to characterize them is to consider that the solution in each point in space and time is influenced by a certain finite region at each previous point in space<fnref n="6" />. The heat equation (section 4.3) is the standard example of this type.

<footnote n="5">5. Actually, the boundary conditions are can be more general, involving derivatives on the interval end points. Here we only look at Dirichlet boundary conditions which prescribe function values on the boundary of the domain.</footnote>

<footnote n="6">6. This leads to a condition limiting the time step in IBVP, known as the Courant-Friedrichs-Lewy condition http://en.wikipedia.org/wiki/CourantFriedrichsLewy_condition. It describes the notion that in the exact problem $u(x, t)$ depends on a range of $u(x', t - \Delta t)$ values; the time step of the numerical method has to be small enough that the numerical solution takes all these points into account.</footnote>
* Elliptic PDEs have the form

$$Au_{xx} + Bu_{yy} + \text{lower order terms} = 0$$

where $A, B > 0$; they typically describe processes that have reached a steady state, for instance as $t \to \infty$ in a parabolic problem. They are characterized by the fact that all points influence each other. These equations often describe phenomena in structural mechanics, such as a beam or a membrane. It is intuitively clear that pressing down on any point of a membrane will change the elevation of every other point, no matter how little. The Poisson equation (section 4.2.2) is the standard example of this type.

### 4.2.2 The Poisson equation in one space dimension

We call the operator $\Delta$, defined by

$$\Delta u = u_{xx} + u_{yy},$$

a second order differential operator, and equation (4.5) a second-order PDE. Specifically, the problem

$$-\Delta u = -u_{xx}(\bar{x}) - u_{yy}(\bar{x}) = f(\bar{x}) \text{ for } \bar{x} \in \Omega = [0, 1]^2 \text{ with } u(\bar{x}) = u_0 \text{ on } \delta\Omega. \quad (4.6)$$

is called the Poisson equation, in this case defined on the unit square. Second order PDEs are quite common, describing many phenomena in fluid and heat flow and structural mechanics.

At first, for simplicity, we consider the one-dimensional Poisson equation

$$-u_{xx} = f(x).$$

and we consider the two-dimensional case below; the extension to three dimensions will then be clear.

In order to find a numerical scheme we use Taylor series as before, expressing $u(x + h)$ and $u(x - h)$ in terms of $u$ and its derivatives at $x$. Let $h > 0$, then

$$u(x+h) = u(x) + u'(x)h + u''(x)\frac{h^2}{2!} + u'''(x)\frac{h^3}{3!} + u^{(4)}(x)\frac{h^4}{4!} + u^{(5)}(x)\frac{h^5}{5!} + \dots$$

and

$$u(x-h) = u(x) - u'(x)h + u''(x)\frac{h^2}{2!} - u'''(x)\frac{h^3}{3!} + u^{(4)}(x)\frac{h^4}{4!} - u^{(5)}(x)\frac{h^5}{5!} + \dots$$

Our aim is now to approximate $u''(x)$. We see that the $u'$ terms in these equations would cancel out under addition, leaving $2u(x)$:

$$u(x+h) + u(x-h) = 2u(x) + u''(x)h^2 + u^{(4)}(x)\frac{h^4}{12} + \dots$$

so

$$-u''(x) = \frac{2u(x) - u(x+h) - u(x-h)}{h^2} + u^{(4)}(x)\frac{h^2}{12} + \dots \quad (4.7)$$
The basis for a numerical scheme for (4.4) is then the observation

$$ \frac{2u(x) - u(x+h) - u(x-h)}{h^2} = f(x, u(x), u'(x)) + O(h^2), $$

which shows that we can approximate the differential operator by a difference operator, with an $O(h^2)$ truncation error as $h \downarrow 0$.

To derive a numerical method, we divide the interval $[0, 1]$ into equally spaced points: $x_k = kh$ where $h = 1/(n + 1)$ and $k = 0...n + 1$. With these, the Finite Difference (FD) formula (4.7) leads to a numerical scheme that forms a system of equations:

$$ -u_{k+1} + 2u_k - u_{k-1} = h^2 f(x_k) \quad \text{for } k = 1, \dots, n \qquad (4.8) $$

This process of using the FD formula (4.7) for the approximate solution of a PDE is known as the Finite Difference Method (FDM).

For most values of $k$ this equation relates $u_k$ unknown to the unknowns $u_{k-1}$ and $u_{k+1}$. The exceptions are $k = 1$ and $k = n$. In that case we recall that $u_0$ and $u_{n+1}$ are known boundary conditions, and we write the equations with unknowns on the left and known quantities on the right as

$$ \begin{cases} 2u_1 - u_2 = h^2 f(x_1) + u_0 \\ 2u_n - u_{n-1} = h^2 f(x_n) + u_{n+1}. \end{cases} $$

We can now summarize these equations for $u_k$, $k = 1...n-1$ as a matrix equation:

$$ \begin{pmatrix} 2 & -1 & & \\ -1 & 2 & -1 & \\ & \ddots & \ddots & \ddots \\ & & & h \end{pmatrix} \begin{pmatrix} u_1 \\ u_2 \\ \vdots \\ u_n \end{pmatrix} = \begin{pmatrix} h^2 f_1 + u_0 \\ h^2 f_2 \\ \vdots \\ u_n \end{pmatrix} \qquad (4.9) $$

This has the form $Au = f$ with $A$ a fully known matrix, $f$ a fully known vector, and $u$ a vector of unknowns. Note that the right hand side vector has the boundary values of the problem in the first and last locations. This means that, if you want to solve the same differential equation with different boundary conditions, only the vector $f$ changes.

**Exercise 4.3.** A condition of the type $u(0) = u_0$ is called a *Dirichlet boundary condition*. Physically, this corresponds to knowing the temperature at the end point of a rod. Other boundary conditions exist. Specifying a value for the derivative, $u'(0) = u'_0$, rather than for the function value, would be appropriate if we are modeling fluid flow and the outflow rate at $x = 0$ is known. This is known as a *Neumann boundary condition*.

A Neumann boundary condition $u'(0) = u'_0$ can be modeled by stating

$$ \frac{u_0 - u_1}{h} = u'_0. $$

Show that, unlike in the case of the Dirichlet boundary condition, this affects the matrix of the linear system.
Show that having a Neumann boundary condition at both ends gives a singular matrix, and therefore no unique solution to the linear system. (Hint: guess the vector that has eigenvalue zero.)

Physically this makes sense. For instance, in an elasticity problem, Dirichlet boundary conditions state that the rod is clamped at a certain height; a Neumann boundary condition only states its angle at the end points, which leaves its height undetermined.

Let us list some properties of A that you will see later are relevant for solving such systems of equations:

* The matrix is very sparse: the percentage of elements that is nonzero is low. The nonzero elements are not randomly distributed but located in a band around the main diagonal. We call this a banded matrix in general, and a tridiagonal matrix in this specific case.

* The matrix is symmetric. This property does not hold for all matrices that come from discretizing BVPs, but it is true if there are no odd order (meaning first, third, fifth, ...) derivatives, such as $u_x, u_{xxx}, u_{xy}$.

* Matrix elements are constant in each diagonal, that is, in each set of points $\{(i, j) : i - j = c\}$ for some c. This is only true for very simple problems. It is no longer true if the differential equation has location dependent terms such as $\frac{d}{dx}(a(x)\frac{d}{dx}u(x))$. It is also no longer true if we make h variable through the interval, for instance because we want to model behaviour around the left end point in more detail.

* Matrix elements conform to the following sign pattern: the diagonal elements are positive, and the off-diagonal elements are nonpositive. This property depends on the numerical scheme used, but it is often true. Together with the following property of definiteness, this is called an M-matrix. There is a whole mathematical theory of these matrices [9].

* The matrix is positive definite: $x^tAx > 0$ for all nonzero vectors $x$. This property is inherited from the original continuous problem, if the numerical scheme is carefully chosen. While the use of this may not seem clear at the moment, later you will see methods for solving the linear system that depend on it.

Strictly speaking the solution of equation (4.9) is simple: $u = A^{-1}f$. However, computing $A^{-1}$ is not the best way of finding $u$. As observed just now, the matrix A has only $3N$ nonzero elements to store. Its inverse, on the other hand, does not have a single nonzero. Although we will not prove it, this sort of statement holds for most sparse matrices. Therefore, we want to solve $Au = f$ in a way that does not require $O(n^2)$ storage.

### 4.2.3 The Poisson equation in two space dimensions

The one-dimensional BVP above was atypical in a number of ways, especially related to the resulting linear algebra problem. In this section we will investigate the two-dimensional Poisson problem. You'll see that it constitutes a non-trivial generalization of the one-dimensional problem. The three-dimensional case is very much like the two-dimensional case, so we will not discuss it. (For one essential difference, see the discussion in section 6.8.1.)

The one-dimensional problem above had a function $u = u(x)$, which now becomes $u = u(x, y)$ in two dimensions. The two-dimensional problem we are interested is then

$$ -u_{xx} - u_{yy} = f, \quad (x, y) \in [0, 1]^2, \qquad (4.10) $$
where the values on the boundaries are given. We get our discrete equation by applying equation (4.7) in x and y directions:

$$ 
\begin{aligned}
-u_{xx}(x, y) &= \frac{2u(x,y)-u(x+h,y)-u(x-h,y)}{h^2} + u^{(4)}(x,y)\frac{h^2}{12} + \dots \\
-u_{yy}(x, y) &= \frac{2u(x,y)-u(x,y+h)-u(x,y-h)}{h^2} + u^{(4)}(x,y)\frac{h^2}{12} + \dots
\end{aligned}
$$

or, taken together:

$$ 4u(x, y) - u(x + h, y) - u(x - h, y) - u(x, y + h) - u(x, y - h) = 1/h^2 f(x, y) + O(h^2) \quad (4.11) $$

Let again $h = 1/(n + 1)$ and define $x_i = ih$ and $y_j = jh$; let $u_{ij}$ be the approximation to $u(x_i, y_j)$, then our discrete equation becomes

$$ 4u_{ij} - u_{i+1,j} - u_{i-1,j} - u_{i,j+1} - u_{i,j-1} = h^{-2}f_{ij}. \quad (4.12) $$

We now have $n \times n$ unknowns $u_{ij}$. To render this in a linear system as before we need to put them in a linear ordering, which we do by defining $I = I_{ij} = j + i \times n$. This is called the lexicographic ordering since it sorts the coordinates $(i, j)$ as if they are strings.

Using this ordering gives us $N = n^2$ equations

$$ 4u_I - u_{I+1} - u_{I-1} - u_{I+n} - u_{I-n} = h^{-2}f_I, \quad I = 1, \dots, N \quad (4.13) $$

and the linear system looks like

$$ A = \begin{pmatrix}
4 & -1 & & & & \emptyset & & -1 & & & & & & \emptyset \\
-1 & 4 & -1 & & & & & & & & & & & \\
& \ddots & \ddots & \ddots & & & & & & & & & & & \\
& & \ddots & \ddots & -1 & & & & & & & & & & \\
& & & \emptyset & -1 & 4 & & \emptyset & & & & & & \\
-1 & & & & \emptyset & & 4 & -1 & & & & & & \\
& -1 & & & & -1 & 4 & -1 & & & & & & \\
\uparrow & & \ddots & & & & \uparrow & \uparrow & \uparrow & & & & & \\
k-n & & & -1 & & & k-1 & k & k+1 & -1 & & & k+n \\
& & & & & & & & -1 & 4 & & & \\
& & & & & & & & & \ddots & & & \ddots
\end{pmatrix} \quad (4.14) $$

As in the one-dimensional case, we see that a BVP gives rise to a sparse matrix.

It may seem the natural thing to consider this system of linear equations in its matrix form. However, it can be more insightful to render these equations in a way that makes clear the two-dimensional connections of the unknowns. For this, figure 4.1 pictures the variables in the domain, and how equation (4.12) relates them. From now on, when making such a picture of the domain, we will just use the indices of the variables, and omit the 'u' identifier.
<figure><img src="image_73.png" /><figcaption>Figure 4.1: A difference stencil applied to a two-dimensional square domain</figcaption></figure>

The matrix of equation 4.14 is banded as in the one-dimensional case, but unlike in the one-dimensional case, there are zeros inside the band. (This has some important consequences when we try to solve the linear system; see section 5.4.3.) Because the matrix has five nonzero diagonals, it is said to be of *penta-diagonal* structure.

You can also put a block structure on the matrix, by grouping the unknowns together that are in one row of the domain. This is called a *block matrix*, and, on the block level, it has a *tridiagonal* structure, so we call this a *block tridiagonal* matrix. Note that the diagonal blocks themselves are tridiagonal; the off-diagonal blocks are minus the identity matrix.

This matrix, like the one-dimensional example above, has constant diagonals, but this is again due to the simple nature of the problem. In practical problems it will not be true. That said, such ‘constant coefficient’ problems occur, and when they are on rectangular domains, there are very efficient methods for solving the linear system with $N \log N$ time complexity.

**Exercise 4.4.** The block structure of the matrix, with all diagonal blocks having the same size, is due to the fact that we defined our BVP on a square domain. Sketch the matrix structure that arises from discretizing equation (4.10), again with central differences, but this time defined on a triangular domain; see figure 4.2. Show that, again, there is a block tridiagonal matrix structure, but that the blocks are now of varying sizes. Hint: start by sketching a small example. For $n = 4$ you should get a $10 \times 10$ matrix with a $4 \times 4$ block structure.

For domains that are even more irregular, the matrix structure will also be irregular.

The regular block structure is also caused by our decision to order the unknowns by rows and columns. This known as the *natural ordering* or *lexicographic ordering*; various other orderings are possible. One common way of ordering the unknowns is the *red-black ordering* or *checkerboard ordering* which has advantages for parallel computation. This will be discussed in section 6.7.
<figure><img src="image_74.png" /><figcaption>Figure 4.2: A triangular domain of definition for the Poisson equation</figcaption></figure>

There is more to say about analytical aspects of the BVP (for instance, how smooth is the solution and how does that depend on the boundary conditions?) but those questions are outside the scope of this course. Here we only focus on the numerical aspects of the matrices. In the chapter on linear algebra, we will come back to the BVP, since solving the linear system is mathematically interesting.

### 4.2.4 Difference stencils

The discretization (4.11) is often phrased as applying the *difference stencil*

$$ \begin{array}{ccc} \cdot & -1 & \cdot \\ -1 & 4 & -1 \\ \cdot & -1 & \cdot \end{array} $$

to the function $u$. Given a physical domain, we apply the stencil to each point in that domain to derive the equation for that point. Figure 4.1 illustrates that for a square domain of $n \times n$ points. Connecting this figure with equation (4.14), you see that the connections in the same line give rise to the main diagonal and first upper and lower offdiagonal; the connections to the next and previous lines become the nonzeros in the off-diagonal blocks.

This particular stencil is often referred to as the ‘5-point star’. There are other difference stencils; the structure of some of them are depicted in figure 4.3. A stencil with only connections in horizontal or vertical direction is called a ‘star stencil’, while one that has cross connections (such as the second in figure 4.3) is called a ‘box stencil’.

**Exercise 4.5.** Consider the third stencil in figure 4.3, used for a BVP on a square domain. What does the sparsity structure of the resulting matrix look like, if we again order the variables by rows and columns?

Other stencils than the 5-point star can be used to attain higher accuracy, for instance giving a truncation error of $O(h^4)$. They can also be used for other differential equations than the one discussed above. For
<figure><img src="image_75.png" /><figcaption>Figure 4.3: The structure of some difference stencils in two dimensions</figcaption></figure>

instance, it is not hard to show that for the equation $u_{xxxx} + u_{yyyy} = f$ we need a stencil that contains both $x, y \pm h$ and $x, y \pm 2h$ connections, such as the third stencil in the figure. Conversely, using the 5-point stencil no values of the coefficients give a discretization of the fourth order problem with less than $O(1)$ truncation error.

While the discussion so far has been about two-dimensional problems, it is easily generalized to higher dimensions for such equations as $-u_{xx} - u_{yy} - u_{zz} = f$. The straightforward generalization of the 5-point stencil, for instance, becomes a 7-point stencil in three dimensions.

### 4.2.5 Other discretization techniques

In the above, we used the FDM to find a numerical solution to a differential equation. There are various other techniques, and in fact, in the case of boundary value problems, they are usually preferred over finite differences. The most popular methods are the *FEM* and the *finite volume method*. Especially the finite element method is attractive, since it can handle irregular shapes more easily, and it is more amenable to approximation error analysis. However, on the simple problems discussed here it gives similar or even the same linear systems as FD methods, so we limit the discussion to Finite Differences, since we are mostly concerned with the computational aspects of the linear systems.

There will be a brief discussion of finite element matrices in section 6.6.2.

## 4.3 Initial boundary value problem

We will now go on to discuss an Initial Boundary Value Problem (IBVP), which, as you may deduce from the name, combines aspects of IVPs and BVPs. Here we will limit ourselves to one space dimension.

The problem we are considering is that of heat conduction in a rod, where $T(x, t)$ describes the temperature in location $x$ at time $t$, for $x \in [a, b]$, $t > 0$. The so-called *heat equation* (see Appendix 14 for a quick introduction to PDEs in general and the heat equation in particular) is:

$$ \frac{\partial T(x, t)}{\partial t} - \alpha \frac{\partial^2 T(x, t)}{\partial x^2} = q(x, t) $$

where

* The initial condition $T(x, 0) = T_0(x)$ describes the initial temperature distribution.
* The boundary conditions $T(a, t) = T_a(t)$, $T(b, t) = T_b(t)$ describe the ends of the rod, which can for instance be fixed to an object of a known temperature.

* The material the rod is made of is modeled by a single parameter $\alpha > 0$, the thermal diffusivity, which describes how fast heat diffuses through the material.

* The forcing function $q(x, t)$ describes externally applied heating, as a function of both time and place.

There is a simple connection between the IBVP and the BVP: if the boundary functions $T_a$ and $T_b$ are constant, and $q$ does not depend on time, only on location, then intuitively $T$ will converge to a steady state. The equation for this is $-\alpha u''(x) = q$.

### 4.3.1 Discretization

We now discretize both space and time, by $x_{j+1} = x_j + \Delta x$ and $t_{k+1} = t_k + \Delta t$, with boundary conditions $x_0 = a$, $x_n = b$, and $t_0 = 0$. We write $T_j^k$ for the numerical solution at $x = x_j, t = t_k$; with a little luck, this will approximate the exact solution $T(x_j, t_k)$.

For the space discretization we use the central difference formula (4.8):

$$ \left. \frac{\partial^2 T(x, t_k)}{\partial x^2} \right|_{x=x_j} \Rightarrow \frac{T_{j-1}^k - 2T_j^k + T_{j+1}^k}{\Delta x^2}. $$

For the time discretization we can use any of the schemes in section 4.1.2. We will investigate again the explicit and implicit schemes, with similar conclusions about the resulting stability.

#### 4.3.1.1 Explicit scheme

With explicit time stepping we approximate the time derivative as

$$ \left. \frac{\partial T(x_j, t)}{\partial t} \right|_{t=t_k} \Rightarrow \frac{T_j^{k+1} - T_j^k}{\Delta t}. \qquad (4.15) $$

Taking this together with the central differences in space, we now have

$$ \frac{T_j^{k+1} - T_j^k}{\Delta t} - \alpha \frac{T_{j-1}^k - 2T_j^k + T_{j+1}^k}{\Delta x^2} = q_j^k $$

which we rewrite as

$$ T_j^{k+1} = T_j^k + \frac{\alpha \Delta t}{\Delta x^2} (T_{j-1}^k - 2T_j^k + T_{j+1}^k) + \Delta t q_j^k \quad (4.16) $$

<figure><img src="image_76.png" /><figcaption>Figure 4.4: The difference stencil of the Euler forward method for the heat equation</figcaption></figure>

Pictorially, we render this as a difference stencil in figure 4.4. This expresses that the function value in each point is determined by a combination of points on the previous time level.
It is convenient to summarize the set of equations (4.16) for a given k and all values of j in vector form as

$$T^{k+1} = \left( I - \frac{\alpha \Delta t}{\Delta x^2} K \right) T^k + \Delta t q^k \quad (4.17)$$

where

$$K = \begin{pmatrix} 2 & -1 & & \\ -1 & 2 & -1 & \\ & \ddots & \ddots & \ddots \\ & & & T_n^k \end{pmatrix}, \quad T^k = \begin{pmatrix} T_1^k \\ \vdots \\ T_n^k \end{pmatrix}.$$

The important observation here is that the dominant computation for deriving the vector $T^{k+1}$ from $T^k$ is a simple matrix-vector multiplication:

$$T^{k+1} \leftarrow AT^k + \Delta t q^k$$

where $A = I - \frac{\alpha \Delta t}{\Delta x^2} K$. This is a first indication that the sparse matrix-vector product is an important operation; see sections 5.4 and 6.5. Actual computer programs using an explicit method often do not form the matrix, but evaluate the equation (4.16). However, the linear algebra formulation (4.17) is more insightful for purposes of analysis.

In later chapters we will consider parallel execution of operations. For now, we note that the explicit scheme is trivially parallel: each point can be updated with just the information of a few surrounding points.

### 4.3.1.2 Implicit scheme

In equation (4.15) we let $T^{k+1}$ be defined from $T^k$. We can turn this around by defining $T^k$ from $T^{k-1}$, as we did for the IVP in section 4.1.2.2. For the time discretization this gives

$$\left. \frac{\partial}{\partial t} T(x, t) \right|_{t=t_k} \Rightarrow \frac{T_j^k - T_j^{k-1}}{\Delta t} \quad \text{or} \quad \left. \frac{\partial}{\partial t} T(x, t) \right|_{t=t_{k+1}} \Rightarrow \frac{T_j^{k+1} - T_j^k}{\Delta t}. \quad (4.18)$$

The implicit time step discretization of the whole heat equation, evaluated in $t_{k+1}$, now becomes:

$$\frac{T_j^{k+1} - T_j^k}{\Delta t} - \alpha \frac{T_{j-1}^{k+1} - 2T_j^{k+1} + T_{j+1}^{k+1}}{\Delta x^2} = q_j^{k+1}$$

or

$$T_j^{k+1} - \frac{\alpha \Delta t}{\Delta x^2} (T_{j-1}^{k+1} - 2T_j^{k+1} + T_{j+1}^{k+1}) = T_j^k + \Delta t q_j^{k+1} \quad (4.19)$$

Figure 4.5 renders this as a stencil; this expresses that each point on the current time level influences a combination of points on the next level. Again we write this in vector form:

<figure><img src="image_77.png" /><figcaption>Figure 4.5: The difference stencil of the Euler backward method for the heat equation</figcaption></figure>

$$\left( I + \frac{\alpha \Delta t}{\Delta x^2} K \right) T^{k+1} = T^k + \Delta t q^{k+1} \quad (4.20)$$
As opposed to the explicit method, where a matrix-vector multiplication sufficed, the derivation of the vector $T^{k+1}$ from $T^k$ now involves solving a linear system

$$T^{k+1} \leftarrow A^{-1}(T^k + \Delta t q^{k+1})$$

where $A = I + \frac{\alpha \Delta t}{\Delta x^2} K$, a harder operation than the matrix-vector multiplication. In this case, it is not possible, as above, to evaluate the equation (4.19) directly. Codes using an implicit method actually form the coefficient matrix, and solve the system (4.20) as such. Solving linear systems will be the focus of much of chapters 5 and 6.

In contrast to the explicit scheme, we now have no obvious parallelization strategy. The parallel solution of linear systems will occupy us in sections 6.6 and on.

**Exercise 4.6.** Show that the flop count for a time step of the implicit method is of the same order as of a time step of the explicit method. (This only holds for a problem with one space dimension.) Give at least one argument why we consider the implicit method as computationally 'harder'.

The numerical scheme that we used here is of first order in time and second order in space: the truncation error (section 4.1.2) is $O(\Delta t + \Delta x^2)$. It would be possible to use a scheme that is second order in time by using central differences in time too. Alternatively, see exercise 4.7.

### 4.3.2 Stability analysis

We now analyse the stability of the explicit and implicit schemes in a simple case. Let $q \equiv 0$, and assume $T_j^k = \beta^k e^{i \ell x_j}$ for some $\ell^7$. This assumption is intuitively defensible: since the differential equation does not 'mix' the $x$ and $t$ coordinates, we surmise that the solution will be a product $T(x, t) = v(x) \cdot w(t)$ of the separate solutions of

$$\begin{cases} v_{xx} = c_1 v, & v(0) = 0, v(1) = 0 \\ w_t = c_2 w, & w(0) = 1 \end{cases}$$

The only meaningful solution occurs with $c_1, c_2 < 0$, in which case we find:

$$v_{xx} = -c^2 v \Rightarrow v(x) = e^{-icx} = e^{-i\ell\pi x}$$

where we substitute $c = \ell\pi$ to take boundary conditions into account, and

$$w(t) = e^{-ct} = e^{-ck\Delta t} = \beta^k, \quad \beta = e^{-ck}.$$

If the assumption on this form of the solution holds up, we need $|\beta| < 1$ for stability.

<footnote n="7">7. Actually, $\beta$ is also dependent on $\ell$, but we will save ourselves a bunch of subscripts, since different $\beta$ values never appear together in one formula.</footnote>
Substituting the surmised form for $T_j^k$ into the explicit scheme gives

$$
\begin{align*}
T_j^{k+1} &= T_j^k + \frac{\alpha \Delta t}{\Delta x^2} (T_{j-1}^k - 2T_j^k + T_{j+1}^k) \\
\Rightarrow \beta^{k+1} e^{i\ell x_j} &= \beta^k e^{i\ell x_j} + \frac{\alpha \Delta t}{\Delta x^2} (\beta^k e^{i\ell x_{j-1}} - 2\beta^k e^{i\ell x_j} + \beta^k e^{i\ell x_{j+1}}) \\
&= \beta^k e^{i\ell x_j} \left[ 1 + \frac{\alpha \Delta t}{\Delta x^2} \left[ e^{-i\ell \Delta x} - 2 + e^{i\ell \Delta x} \right] \right] \\
\Rightarrow \beta &= 1 + 2 \frac{\alpha \Delta t}{\Delta x^2} \left[ \frac{1}{2}(e^{i\ell \Delta x} + e^{-i\ell \Delta x}) - 1 \right] \\
&= 1 + 2 \frac{\alpha \Delta t}{\Delta x^2} (\cos(\ell \Delta x) - 1)
\end{align*}
$$

For stability we need $|\beta| < 1$:

* $\beta < 1 \Leftrightarrow 2\frac{\alpha\Delta t}{\Delta x^2}(\cos(\ell\Delta x) - 1) < 0$: this is true for any $\ell$ and any choice of $\Delta x, \Delta t$.

* $\beta > -1 \Leftrightarrow 2\frac{\alpha\Delta t}{\Delta x^2}(\cos(\ell\Delta x) - 1) > -2$: this is true for all $\ell$ only if $2\frac{\alpha\Delta t}{\Delta x^2} < 1$, that is

$$ \Delta t < \frac{\Delta x^2}{2\alpha} $$

The latter condition poses a big restriction on the allowable size of the time steps: time steps have to be small enough for the method to be stable. This is similar to the stability analysis of the explicit method for the IVP; however, now the time step is also related to the space discretization. This implies that, if we decide we need more accuracy in space and we halve the space discretization $\Delta x$, the number of time steps will be multiplied by four.

Let us now consider the stability of the implicit scheme. Substituting the form of the solution $T_j^k = \beta^k e^{i\ell x_j}$ into the numerical scheme gives

$$
\begin{align*}
T_j^{k+1} - T_j^k &= \frac{\alpha \Delta t}{\Delta x^2} (T_{j1}^{k+1} - 2T_j^{k+1} + T_{j+1}^{k+1}) \\
\Rightarrow \beta^{k+1} e^{i\ell \Delta x} - \beta^k e^{i\ell x_j} &= \frac{\alpha \Delta t}{\Delta x^2} (\beta^{k+1} e^{i\ell x_{j-1}} - 2\beta^{k+1} e^{i\ell x_j} + \beta^{k+1} e^{i\ell x_{j+1}})
\end{align*}
$$

Dividing out $e^{i\ell x_j} \beta^{k+1}$ gives

$$
\begin{align*}
1 &= \beta^{-1} + 2\alpha \frac{\Delta t}{\Delta x^2} (\cos \ell \Delta x - 1) \\
\beta &= \frac{1}{1 + 2\alpha \frac{\Delta t}{\Delta x^2} (1 - \cos \ell \Delta x)}
\end{align*}
$$

Since $1 - \cos \ell \Delta x \in (0, 2)$, the denominator is strictly $> 1$. Therefore the condition $|\beta| < 1$ is always satisfied, regardless the value of $\ell$ and the choices of $\Delta x$ and $\Delta t$: the method is always stable.

**Exercise 4.7.** The schemes we considered here are of first order in time and second order in space: their discretization order are $O(\Delta t) + O(\Delta x^2)$. Derive the Crank-Nicolson method that is obtained by averaging the explicit and implicit schemes, show that it is unconditionally stable, and of second order in time.
# Chapter 5

## Numerical linear algebra

In chapter 4 you saw how the numerical solution of partial differential equations can lead to linear algebra problems. Sometimes this is a simple problem – a matrix-vector multiplication in the case of the Euler forward method – but sometimes it is more complicated, such as the solution of a system of linear equations in the case of Euler backward methods. Solving linear systems will be the focus of this chapter; in other applications, which we will not discuss here, eigenvalue problems need to be solved.

You may have learned a simple algorithm for solving system of linear equations: elimination of unknowns, also called Gaussian elimination. This method can still be used, but we need some careful discussion of its efficiency. There are also other algorithms, the so-called iterative solution methods, which proceed by gradually approximating the solution of the linear system. They warrant some discussion of their own.

Because of the PDE background, we only consider linear systems that are square and nonsingular. Rectangular, in particular overdetermined, systems have important applications too in a corner of numerical analysis known as optimization theory. However, we will not cover that in this book.

The standard work on numerical linear algebra is Golub and Van Loan's *Matrix Computations* [71]. It covers algorithms, error analysis, and computational details. Heath's *Scientific Computing* covers the most common types of computations that arise in scientific computing; this book has many excellent exercises and practical projects.

### 5.1 Elimination of unknowns

In this section we are going to take a closer look at Gaussian elimination, or elimination of unknowns. You may have seen this method before (and if not, it will be explained below), but we will be a bit more systematic here so that we can analyze various aspects of it.

One general thread of this chapter will be the discussion of the efficiency of the various algorithms. When you learned how to solve a system of unknowns by gradually eliminating unknowns, you most likely never applied that method to a matrix larger than 4 × 4. The linear systems that occur in PDE solving can be thousands of times larger, and counting how many operations, as well as how much memory, they require becomes important.
Let us consider an example of the importance of efficiency in choosing the right algorithm. The solution of a linear system can be written with a fairly simple explicit formula, using determinants. This is called ‘Cramer’s rule’. It is mathematically elegant, but completely impractical for our purposes.

If a matrix A and a vector b are given, and a vector x satisfying $Ax = b$ is wanted, then, writing $|A|$ for the determinant,

$$x_i = \begin{vmatrix} a_{11} & a_{12} & \dots & a_{1i-1} & b_1 & a_{1i+1} & \dots & a_{1n} \\ a_{21} & \dots & & b_2 & \dots & a_{2n} & & \vdots \\ \vdots & & & \vdots & & \vdots & & \vdots \\ a_{n1} & \dots & & b_n & \dots & a_{nn} & & |A| \end{vmatrix}$$

For any matrix M the determinant is defined recursively as

$$|M| = \sum_{i} (-1)^{i} m_{1i} |M^{[1,i]}|$$

where $M^{[1,i]}$ denotes the matrix obtained by deleting row 1 and column $i$ from M. This means that computing the determinant of a matrix of dimension $n$ means $n$ times computing a size $n-1$ determinant. Each of these requires $n-1$ determinants of size $n-2$, so you see that the number of operations required to compute the determinant is factorial in the matrix size. This quickly becomes prohibitive, even ignoring any issues of numerical stability. Later in this chapter you will see complexity estimates for other methods of solving systems of linear equations that are considerably more reasonable.

Let us now look at a simple example of solving linear equations with elimination of unknowns. Consider the system

$$\begin{aligned} 6x_1 & -2x_2 +2x_3 = 16 \\ 12x_1 & -8x_2 +6x_3 = 26 \\ 3x_1 & -13x_2 +3x_3 = -19 \end{aligned} \qquad (5.1)$$

We eliminate $x_1$ from the second and third equation by

* multiplying the first equation $\times 2$ and subtracting the result from the second equation, and

* multiplying the first equation $\times 1/2$ and subtracting the result from the third equation.

The linear system then becomes

$$\begin{aligned} 6x_1 & -2x_2 +2x_3 = 16 \\ 0x_1 & -4x_2 +2x_3 = -6 \\ 0x_1 & -12x_2 +2x_3 = -27 \end{aligned}$$

Finally, we eliminate $x_2$ from the third equation by multiplying the second equation by 3, and subtracting the result from the third equation:

$$\begin{aligned} 6x_1 & -2x_2 +2x_3 = 16 \\ 0x_1 & -4x_2 +2x_3 = -6 \\ 0x_1 & +0x_2 -4x_3 = -9 \end{aligned}$$
We can now solve $x_3 = 9/4$ from the last equation. Substituting that in the second equation, we get $-4x_2 = -6 - 2x_2 = -21/2$ so $x_2 = 21/8$. Finally, from the first equation $6x_1 = 16 + 2x_2 - 2x_3 = 16 + 21/4 - 9/2 = 76/4$ so $x_1 = 19/6$.

We can write this more compactly by omitting the $x_i$ coefficients. Write

$$ \begin{pmatrix} 6 & -2 & 2 \\ 12 & -8 & 6 \\ 3 & -13 & 3 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 16 \\ 26 \\ -19 \end{pmatrix} $$

as

$$ \begin{bmatrix} 6 & -2 & 2 & 16 \\ 12 & -8 & 6 & 26 \\ 3 & -13 & 3 & -19 \end{bmatrix} \qquad (5.2) $$

then the elimination process is

$$ \begin{bmatrix} 6 & -2 & 2 & 16 \\ 12 & -8 & 6 & 26 \\ 3 & -13 & 3 & -19 \end{bmatrix} \rightarrow \begin{bmatrix} 6 & -2 & 2 & 16 \\ 0 & -4 & 2 & -6 \\ 0 & -12 & 2 & -27 \end{bmatrix} \rightarrow \begin{bmatrix} 6 & -2 & 2 & 16 \\ 0 & -4 & 2 & -6 \\ 0 & 0 & -4 & -9 \end{bmatrix} $$

In the above example, the matrix coefficients could have been any real (or, for that matter, complex) coefficients, and you could follow the elimination procedure mechanically. There is the following exception. At some point in the computation, we divided by the numbers 6, -4, -4 which are found on the diagonal of the matrix in the last elimination step. These quantities are called the *pivots*, and clearly they are required to be nonzero.

**Exercise 5.1.** The system

$$ \begin{array}{lclcl} 6x_1 & -2x_2 & +2x_3 & = & 16 \\ 12x_1 & -4x_2 & +6x_3 & = & 26 \\ 3x_1 & -13x_2 & +3x_3 & = & -19 \end{array} $$

is the same as the one we just investigated in equation (5.1), except for the (2, 2) element. Confirm that you get a zero pivot in the second step.

The first pivot is an element of the original matrix; as you saw in the preceding exercise, the other pivots can not easily be found without doing the actual elimination. In particular, there is no easy way of predicting zero pivots from looking at the system of equations.

If a pivot turns out to be zero, all is not lost for the computation: we can always exchange two matrix rows; this is known as *pivoting*. It is not hard to show<fnref n="1" /> that with a nonsingular matrix there is always a row exchange possible that puts a nonzero element in the pivot location.

**Exercise 5.2.** Suppose you want to exchange matrix rows 2 and 3 of the system of equations in equation (5.2). What other adjustments would you have to make to make sure you still compute the correct solution? Continue the system solution of the previous exercise by exchanging rows 2 and 3, and check that you get the correct answer.

<footnote n="1">1. And you can find this in any elementary linear algebra textbook.</footnote>
Exercise 5.3. Take another look at exercise 5.1. Instead of exchanging rows 2 and 3, exchange columns 2 and 3. What does this mean in terms of the linear system? Continue the process of solving the system; check that you get the same solution as before.

In general, with floating point numbers and round-off, it is very unlikely that a matrix element will become exactly zero during a computation. Also, in a PDE context, the diagonal is usually nonzero. Does that mean that pivoting is in practice almost never necessary? The answer is no: pivoting is desirable from a point of view of numerical stability. In the next section you will see an example that illustrates this fact.

## 5.2 Linear algebra in computer arithmetic

In most of this chapter, we will act as if all operations can be done in exact arithmetic. However, it is good to become aware of some of the potential problems due to our finite precision computer arithmetic. This allows us to design algorithms that minimize the effect of roundoff. A more rigorous approach to the topic of numerical linear algebra includes a full-fledged error analysis of the algorithms we discuss; however, that is beyond the scope of this course. Error analysis of computations in computer arithmetic is the focus of Wilkinson's classic *Rounding errors in Algebraic Processes* [164] and Higham's more recent *Accuracy and Stability of Numerical Algorithms* [86].

Here, we will only note a few paradigmatic examples of the sort of problems that can come up in computer arithmetic: we will show why pivoting during LU factorization is more than a theoretical device, and we will give two examples of problems in eigenvalue calculations due to the finite precision of computer arithmetic.

### 5.2.1 Roundoff control during elimination

Above, you saw that row interchanging ('pivoting') is necessary if a zero element appears on the diagonal during elimination of that row and column. Let us now look at what happens if the pivot element is not zero, but close to zero.

Consider the linear system

$$ \begin{pmatrix} \epsilon & 1 \\ 1 & 1 \end{pmatrix} x = \begin{pmatrix} 1 + \epsilon \\ 2 \end{pmatrix} $$

which has the solution $x = (1, 1)^t$. Using the (1, 1) element to clear the remainder of the first column gives:

$$ \begin{pmatrix} \epsilon & 1 \\ 0 & 1 - \frac{1}{\epsilon} \end{pmatrix} x = \begin{pmatrix} 1 + \epsilon \\ 2 - \frac{1+\epsilon}{\epsilon} \end{pmatrix}. $$

We can now solve $x_2$ and from it $x_1$.

If $\epsilon$ is small, say $\epsilon < \epsilon_{\text{mach}}$, the $1+\epsilon$ term in the right hand side will be simply 1: our linear system will be

$$ \begin{pmatrix} \epsilon & 1 \\ 1 & 1 \end{pmatrix} x = \begin{pmatrix} 1 \\ 2 \end{pmatrix} $$
but the solution $(1, 1)^t$ will still satisfy the system in machine arithmetic.

Next, $1/\epsilon$ will be very large, so the second component of the right hand side after elimination will be $2 - \frac{1}{\epsilon} = -1/\epsilon$. Also, the $(2, 2)$ element of the matrix is then $-1/\epsilon$ instead of $1 - 1/\epsilon$:

$$ \begin{pmatrix} \epsilon & 1 \\ 0 & -1/\epsilon \end{pmatrix} x = \begin{pmatrix} 1 \\ -1/\epsilon \end{pmatrix} $$

We get the correct value $x_2 = 1$, but

$$ \epsilon x_1 + x_2 = 1 \Rightarrow \epsilon x_1 = 0 \Rightarrow x_1 = 0, $$

which is 100% wrong.

What would have happened if we had pivoted as described above? We exchange the matrix rows, giving

$$ \begin{pmatrix} 1 & 1 \\ \epsilon & 1 \end{pmatrix} x = \begin{pmatrix} 2 \\ 1 + \epsilon \end{pmatrix} \Rightarrow \begin{pmatrix} 1 & 1 \\ 0 & 1 - \epsilon \end{pmatrix} x = \begin{pmatrix} 2 \\ 1 - \epsilon \end{pmatrix} $$

Now we get, regardless the size of epsilon:

$$ x_2 = \frac{1 - \epsilon}{1 - \epsilon} = 1, \quad x_1 = 2 - x_2 = 1 $$

In this example we used a very small value of $\epsilon$; a much more refined analysis shows that even with $\epsilon$ greater than the machine precision pivoting still makes sense. The general rule of thumb is: Always do row exchanges to get the largest remaining element in the current column into the pivot position. In chapter 4 you saw matrices that arise in certain practical applications; it can be shown that for them pivoting is never necessary; see exercise 5.13.

The pivoting that was discussed above is also known as partial pivoting, since it is based on row exchanges only. Another option would be full pivoting, where row and column exchanges are combined to find the largest element in the remaining subblock, to be used as pivot. Finally, diagonal pivoting applies the same exchange to rows and columns. (This is equivalent to renumbering the unknowns of the problem, a strategy which we will consider in section 6.8 for increasing the parallelism of the problem.) This means that pivots are only searched on the diagonal. From now on we will only consider partial pivoting.

## 5.2.2 Influence of roundoff on eigenvalue computations

Consider the matrix

$$ A = \begin{pmatrix} 1 & \epsilon \\ \epsilon & 1 \end{pmatrix} $$

where $\epsilon_{mach} < |\epsilon| < \sqrt{\epsilon_{mach}}$, which has eigenvalues $1 + \epsilon$ and $1 - \epsilon$. If we calculate its characteristic polynomial in computer arithmetic

$$ \begin{vmatrix} 1 - \lambda & \epsilon \\ \epsilon & 1 - \lambda \end{vmatrix} = \lambda^2 - 2\lambda + (1 - \epsilon^2) = \lambda^2 - 2\lambda + 1. $$
we find a double eigenvalue 1. Note that the exact eigenvalues are expressible in working precision; it is the algorithm that causes the error. Clearly, using the characteristic polynomial is not the right way to compute eigenvalues, even in well-behaved, symmetric positive definite, matrices.

An unsymmetric example: let A be the matrix of size 20

$$A = \begin{pmatrix} 20 & 20 & & & & \emptyset \\ & 19 & 20 & & & \\ & & \ddots & \ddots & & \\ & & & 2 & 20 & \\ & \emptyset & & & 1 & \\ & & & & & \emptyset \end{pmatrix} .$$

Since this is a triangular matrix, its eigenvalues are the diagonal elements. If we perturb this matrix by setting $A_{20,1} = 10^{-6}$ we find a perturbation in the eigenvalues that is much larger than in the elements:

$$\lambda = 20.6 \pm 1.9i, 20.0 \pm 3.8i, 21.2, 16.6 \pm 5.4i, \dots$$

Also, several of the computed eigenvalues have an imaginary component, which the exact eigenvalues do not have.

## 5.3 LU factorization

So far, we have looked at eliminating unknowns in the context of solving a single system of linear equations. Suppose you need to solve more than one system with the same matrix, but with different right hand sides. This happens for instance if you take multiple time steps in an implicit method. Can you use any of the work you did in the first system to make solving subsequent ones easier?

The answer is yes. You can split the solution process in a part that only concerns the matrix, and a part that is specific to the right hand side. If you have a series of systems to solve, you have to do the first part only once, and, luckily, that even turns out to be the larger part of the work.

Let us take a look at the same example of section 5.1 again.

$$A = \begin{pmatrix} 6 & -2 & 2 \\ 12 & -8 & 6 \\ 3 & -13 & 3 \end{pmatrix}$$

In the elimination process, we took the 2nd row minus 2× the first and the 3rd row minus 1/2× the first. Convince yourself that this combining of rows can be done by multiplying A from the left by

$$L_1 = \begin{pmatrix} 1 & 0 & 0 \\ -2 & 1 & 0 \\ -1/2 & 0 & 1 \end{pmatrix},$$

which is the identity with the elimination coefficients in the first column, below the diagonal. The first step in elimination of variables is equivalent to transforming the system $Ax = b$ to $L_1Ax = L_1b$.
In the next step, you subtracted 3× the second row from the third. Convince yourself that this corresponds to left-multiplying the current matrix $L_1A$ by

$$L_2 = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -3 & 1 \end{pmatrix}$$

We have now transformed our system $Ax = b$ into $L_2L_1Ax = L_2L_1b$, and $L_2L_1A$ is of 'upper triangular' form. If we define $U = L_2L_1A$, then $A = L_1^{-1}L_2^{-1}U$. How hard is it to compute matrices such as $L_2^{-1}$? Remarkably easy, it turns out to be.

We make the following observations:

$$L_1 = \begin{pmatrix} 1 & 0 & 0 \\ -2 & 1 & 0 \\ -1/2 & 0 & 1 \end{pmatrix} \qquad L_1^{-1} = \begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 1/2 & 0 & 1 \end{pmatrix}$$

and likewise

$$L_2 = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -3 & 1 \end{pmatrix} \qquad L_2^{-1} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 3 & 1 \end{pmatrix}$$

and even more remarkable:

$$L_1^{-1}L_2^{-1} = \begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 1/2 & 3 & 1 \end{pmatrix},$$

that is, $L_1^{-1}L_2^{-1}$ contains the off-diagonal elements of $L_1^{-1}$, $L_2^{-1}$ unchanged, and they in turn contain the elimination coefficients. (This is a special case of *Householder reflectors*; see 12.6.)

**Exercise 5.4.** Show that a similar statement holds, even if there are elements above the diagonal.

If we define $L = L_1^{-1}L_2^{-1}$, we now have $A = LU$; this is called an *LU factorization*. We see that the coefficients of $L$ below the diagonal are the negative of the coefficients used during elimination. Even better, the first column of $L$ can be written while the first column of $A$ is being eliminated, so the computation of $L$ and $U$ can be done without extra storage, at least if we can afford to lose $A$.

### 5.3.1 The algorithm

Let us write out the LU factorization algorithm in more or less formal code.

```
(LU factorization):
    for k = 1, n-1:
        <eliminate values in column k>
    <eliminate values in column k>:
        for i = k + 1 to n:
```
$\langle \text{compute multiplier for row } i \rangle$

$\langle \text{update row } i \rangle$

$\langle \text{compute multiplier for row } i \rangle$

$$a_{ik} \leftarrow a_{ik}/a_{kk}$$

$\langle \text{update row } i \rangle$:

for $j = k + 1$ to $n$:

$$a_{ij} \leftarrow a_{ij} - a_{ik} * a_{kj}$$

Or, putting everything together:

$\langle \text{LU factorization} \rangle$:

for $k = 1, n - 1$:

$$\text{for } i = k + 1 \text{ to } n: \tag{5.3}$$

$$a_{ik} \leftarrow a_{ik}/a_{kk}$$

for $j = k + 1$ to $n$:

$$a_{ij} \leftarrow a_{ij} - a_{ik} * a_{kj}$$

This is the most common way of presenting the LU factorization. However, other ways of computing the same result exist. Algorithms such as the LU factorization can be coded in several ways that are mathematically equivalent, but that have different computational behaviour. This issue, in the context of dense matrices, is the focus of van de Geijn and Quintana's *The Science of Programming Matrix Computations* [157].

### 5.3.2 The Cholesky factorization

The LU factorization of a symmetric matrix does not give an L and U that are each other's transpose: L has ones on the diagonal and U has the pivots. However it is possible to make a factorization of a symmetric matrix A of the form $A = LL^t$. This has the advantage that the factorization takes the same amount of space as the original matrix, namely $n(n + 1)/2$ elements. We a little luck we can, just as in the LU case, overwrite the matrix with the factorization.

We derive this algorithm by reasoning inductively. Let us write $A = LL^t$ on block form:

$$A = \begin{pmatrix} A_{11} & A_{21}^t \\ A_{21} & A_{22} \end{pmatrix} = LL^t = \begin{pmatrix} \ell_{11} & 0 \\ \ell_{21} & L_{22} \end{pmatrix} \begin{pmatrix} \ell_{11} & \ell_{21}^t \\ 0 & L_{22}^t \end{pmatrix}$$

then $\ell_{11}^2 = a_{11}$, from which we get $\ell_{11}$. We also find $\ell_{11}(L^t)_{1j} = \ell_{j1} = a_{1j}$, so we can compute the whole first column of L. Finally, $A_{22} = L_{22}L_{22}^t + \ell_{12}\ell_{12}^t$, so

$$L_{22}L_{22}^t = A_{22} - \ell_{12}\ell_{12}^t,$$

which shows that $L_{22}$ is the Cholesky factor of the updated $A_{22}$ block. Recursively, the algorithm is now defined.

### 5.3.3 Uniqueness

It is always a good idea, when discussing numerical algorithms, to wonder if different ways of computing lead to the same result. This is referred to as the 'uniqueness' of the result, and it is of practical use: if
the computed result is unique, swapping one software library for another will not change anything in the computation.

Let us consider the uniqueness of LU factorization. The definition of an LU factorization algorithm (without pivoting) is that, given a nonsingular matrix A, it will give a lower triangular matrix L and upper triangular matrix U such that $A = LU$. The above algorithm for computing an LU factorization is deterministic (it does not contain instructions 'take any row that satisfies...'), so given the same input, it will always compute the same output. However, other algorithms are possible, so we need worry whether they give the same result.

Let us then assume that $A = L_1U_1 = L_2U_2$ where $L_1, L_2$ are lower triangular and $U_1, U_2$ are upper triangular. Then, $L_2^{-1}L_1 = U_2U_1^{-1}$. In that equation, the left hand side is the product of lower triangular matrices, while the right hand side contains only upper triangular matrices.

**Exercise 5.5.** Prove that the product of lower triangular matrices is lower triangular, and the product of upper triangular matrices is upper triangular. Is a similar statement true for inverses of nonsingular triangular matrices?

The product $L_2^{-1}L_1$ is apparently both lower triangular and upper triangular, so it must be diagonal. Let us call it $D$, then $L_1 = L_2D$ and $U_2 = DU_1$. The conclusion is that LU factorization is not unique, but it is unique 'up to diagonal scaling'.

**Exercise 5.6.** The algorithm in section 5.3.1 resulted in a lower triangular factor $L$ that had ones on the diagonal. Show that this extra condition makes the factorization unique.

**Exercise 5.7.** Show that an alternative condition of having ones on the diagonal of $U$ is also sufficient for the uniqueness of the factorization.

Since we can demand a unit diagonal in $L$ or in $U$, you may wonder if it is possible to have both. (Give a simple argument why this is not strictly possible.) We can do the following: suppose that $A = LU$ where $L$ and $U$ are nonsingular lower and upper triangular, but not normalized in any way. Write

$$L = (I + L')D_L, \quad U = D_U(I + U'), \quad D = D_L D_U.$$

After some renaming we now have a factorization

$$A = (I + L)D(I + U) \qquad (5.4)$$

where $D$ is a diagonal matrix containing the pivots.

**Exercise 5.8.** Show that you can also normalize the factorization on the form

$$A = (D + L)D^{-1}(D + U).$$

How does this $D$ relate to the previous?

**Exercise 5.9.** Consider the factorization of a tridiagonal matrix this way. How do $L$ and $U$ relate to the triangular parts of $A$? Derive a relation between $D$ and $D_A$ and show that this is the equation that generates the pivots.
### 5.3.4 Pivoting

Above, you saw examples where pivoting, that is, exchanging rows, was necessary during the factorization process, either to guarantee the existence of a nonzero pivot, or for numerical stability. We will now integrate pivoting into the LU factorization.

Let us first observe that row exchanges can be described by a matrix multiplication. Let

$$P^{(i,j)} = \begin{array}{c} i \quad j \\ \begin{array}{cc} 1 & 0 \\ 0 & \ddots \\ \vdots & \ddots \\ \end{array} & \begin{array}{cc} 0 & 1 \\ I & \end{array} \\ j & \begin{array}{cc} 1 & 0 \\ \multicolumn{2}{c}{I} \\ \vdots & \ddots \\ \end{array} \end{array}$$

then $P^{(i,j)}A$ is the matrix $A$ with rows $i$ and $j$ exchanged. Since we may have to pivot in every iteration of the factorization process, we introduce a sequence $p(\cdot)$ where $p(i)$ is the $j$ values of the row that row $i$ is switched with. We write $P^{(i)} \equiv P^{(i,p(i))}$ for short.

**Exercise 5.10.** Show that $P^{(i)}$ is its own inverse.

The process of factoring with partial pivoting can now be described as:

* Let $A^{(i)}$ be the matrix with columns $1 \dots i-1$ eliminated, and partial pivoting applied to get the desired element in the $(i, i)$ location.

* Let $\ell^{(i)}$ be the vector of multipliers in the $i$-th elimination step. (That is, the elimination matrix $L_i$ in this step is the identity plus $\ell^{(i)}$ in the $i$-th column.)

* Let $P^{(i+1)}$ (with $j \ge i+1$) be the matrix that does the partial pivoting for the next elimination step as described above.

* Then $A^{(i+1)} = P^{(i+1)}L_i A^{(i)}$.

In this way we get a factorization of the form

$$L_{n-1}P^{(n-2)}L_{n-2}\cdots L_1P^{(1)}A = U.$$

Suddenly it has become impossible to write $A = LU$: instead we write

$$A = P^{(1)}L_1^{-1}\cdots P^{(n-2)}L_{n-1}^{-1}U. \qquad (5.5)$$

**Exercise 5.11.** Recall from sections 1.6.10 and 1.6.8 that blocked algorithms are often desirable from a performance point of view. Why is the 'LU' factorization with interleaved pivoting matrices' in equation (5.5) bad news for performance?
Fortunately, equation (5.5) can be simplified: the $P$ and $L$ matrices 'almost commute'. We show this by looking at an example: $P^{(2)}L_1 = \tilde{L}_1P^{(2)}$ where $\tilde{L}_1$ is very close to $L_1$.

$$ \begin{pmatrix} 1 & & & \\ & 0 & 1 & \\ & & I & \\ & & & I \end{pmatrix} \begin{pmatrix} 1 & & \emptyset & \\ & \vdots & 1 & \\ \ell^{(1)} & & \ddots & \\ & \vdots & & 1 \end{pmatrix} = \begin{pmatrix} 1 & & \emptyset & \\ & \vdots & 0 & 1 \\ \tilde{\ell}^{(1)} & & & \\ & \vdots & 1 & 0 \\ & & \vdots & I \end{pmatrix} = \begin{pmatrix} 1 & & \emptyset & \\ & \vdots & 1 & \\ \tilde{\ell}^{(1)} & & \ddots & \\ & \vdots & & 1 \end{pmatrix} \begin{pmatrix} 1 & & & \\ & 0 & & 1 \\ & & I & \\ & 1 & & 0 \end{pmatrix} $$

where $\tilde{\ell}^{(1)}$ is the same as $\ell^{(1)}$, except that elements $i$ and $p(i)$ have been swapped. You can now easily convince yourself that similarly $P^{(2)}$ et cetera can be 'pulled through' $L_1$.

As a result we get

$$ P^{(n-2)} \cdots P^{(1)} A = \tilde{L}_1^{-1} \cdots L_{n-1}^{-1} U = \tilde{L}U. \qquad (5.6) $$

This means that we can again form a matrix $L$ just as before, except that every time we pivot, we need to update the columns of $L$ that have already been computed.

**Exercise 5.12.** If we write equation (5.6) as $PA = LU$, we get $A = P^{-1}LU$. Can you come up with a simple formula for $P^{-1}$ in terms of just $P$? Hint: each $P^{(i)}$ is symmetric and its own inverse; see the exercise above.

**Exercise 5.13.** Earlier, you saw that 2D BVP (section 4.2.3) give rise to a certain kind of matrix. We stated, without proof, that for these matrices pivoting is not needed. We can now formally prove this, focusing on the crucial property of *diagonal dominance*:

$$ \forall_i a_{ii} \geq \sum_{j \neq i} |a_{ij}|. $$

Assume that a matrix $A$ satisfies $\forall_{j \neq i}: a_{ij} \leq 0$. Show that the matrix is diagonally dominant iff there are vectors $u, v \geq 0$ (meaning that each component is nonnegative) such that $Au = v$.

Show that, after eliminating a variable, for the remaining matrix $\tilde{A}$ there are again vectors $\tilde{u}, \tilde{v} \geq 0$ such that $\tilde{A}\tilde{u} = \tilde{v}$.

Now finish the argument that (partial) pivoting is not necessary if $A$ is symmetric and diagonally dominant. (One can actually prove that pivoting is not necessary for any *symmetric positive definite (SPD)* matrix, and diagonal dominance is a stronger condition than SPD-ness.)

### 5.3.5 Solving the system

Now that we have a factorization $A = LU$, we can use this to solve the linear system $Ax = LUX = b$. If we introduce a temporary vector $y = Ux$, then we see this takes two steps:

$$ Ly = b, \quad Ux = z. $$
The first part, $Ly = b$ is called the ‘lower triangular solve’, since it involves the lower triangular matrix $L$.

$$ \begin{pmatrix} 1 & & & \\ & \ell_{21} & 1 & \\ & \ell_{31} & \ell_{32} & 1 \\ & & \ddots & \\ & & & \ell_{n1} & \ell_{n2} \\ & & & & \ddots \\ & & & & & \ddots \\ & & & & & & \ell_{nn} \end{pmatrix} \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix} = \begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{pmatrix} $$

In the first row, you see that $y_1 = b_1$. Then, in the second row $\ell_{21}y_1 + y_2 = b_2$, so $y_2 = b_2 - \ell_{21}y_1$. You can imagine how this continues: in every i-th row you can compute $y_i$ from the previous y-values:

$$ y_i = b_i - \sum_{j<i} \ell_{ij} y_j. $$

Since we compute $y_i$ in increasing order, this is also known as the forward substitution, forward solve, or forward sweep.

The second half of the solution process, the upper triangular solve, backward substitution, or backward sweep, computes $x$ from $Ux = y$:

$$ \begin{pmatrix} u_{11} & u_{12} & \cdots & u_{1n} \\ & u_{22} & \cdots & u_{2n} \\ & & \ddots & \\ & & & u_{nn} \\ \emptyset & & & & \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{pmatrix} = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix} $$

Now we look at the last line, which immediately tells us that $x_n = u_{nn}^{-1}y_n$. From this, the line before the last states $u_{n-1n-1}x_{n-1} + u_{n-1n}x_n = y_{n-1}$, which gives $x_{n-1} = u_{n-1n-1}^{-1}(y_{n-1} - u_{n-1n}x_n)$. In general, we can compute

$$ x_i = u_{ii}^{-1} (y_i - \sum_{j>i} u_{ij} y_j) $$

for decreasing values of $i$.

Exercise 5.14. In the backward sweep you had to divide by the numbers $u_{ii}$. That is not possible if any of them are zero. Relate this problem back to the above discussion.

### 5.3.6 Complexity

In the beginning of this chapter, we indicated that not every method for solving a linear system takes the same number of operations. Let us therefore take a closer look at the complexity<fnref n="2" />, that is, the number of operations as function of the problem size, of the use of an LU factorization in solving the linear system.

First we look at the computation of $x$ from $LUx = b$ (‘solving the linear system’; section 5.3.5), given that we already have the factorization $A = LU$. Looking at the lower and upper triangular part together, you

<footnote n="2">2. See appendix 13 for an explanation of complexity.</footnote>
see that you perform a multiplication with all off-diagonal elements (that is, elements $l_{ij}$ or $u_{ij}$ with $i \neq j$). Furthermore, the upper triangular solve involves divisions by the $u_{ii}$ elements. Now, division operations are in general much more expensive than multiplications, so in this case you would compute the values $1/u_{ii}$, and store them instead.

**Exercise 5.15.** Take a look at the factorization algorithm, and argue that storing the reciprocals of the pivots does not add to the computational complexity.

Summing up, you see that, on a system of size $n \times n$, you perform $n^2$ multiplications and roughly the same number of additions. This shows that, given a factorization, solving a linear system has the same complexity as a simple matrix-vector multiplication, that is, of computing $Ax$ given $A$ and $x$.

The complexity of constructing the LU factorization is a bit more involved to compute. Refer to the algorithm in section 5.3.1. You see that in the $k$-th step two things happen: the computation of the multipliers, and the updating of the rows. There are $n-k$ multipliers to be computed, each of which involve a division. After that, the update takes $(n-k)^2$ additions and multiplications. If we ignore the divisions for now, because there are fewer of them, we find that the LU factorization takes $\sum_{k=1}^{n-1} 2(n-k)^2$ operations. If we number the terms in this sum in the reverse order, we find

$$ \#ops = \sum_{k=1}^{n-1} 2k^2. $$

Since we can approximate a sum by an integral, we find that this is $2/3n^3$ plus some lower order terms. This is of a higher order than solving the linear system: as the system size grows, the cost of constructing the LU factorization completely dominates.

Of course there is more to algorithm analysis than operation counting. While solving the linear system has the same complexity as a matrix-vector multiplication, the two operations are of a very different nature. One big difference is that both the factorization and forward/backward solution involve recursion, so they are not simple to parallelize. We will say more about that later on.

### 5.3.7 Block algorithms

Often, matrices have a natural block structure, such as in the case of two-dimensional BVPs; section 4.2.3. Many linear algebra operations can be formulated in terms of these blocks. This can have several advantages over the traditional scalar view of matrices. For instance, it improves *cache blocking* (section 1.6.7); it also facilitates scheduling linear algebra algorithms on *multicore* architectures (section 6.12).

For block algorithms we write a matrix as

$$ A = \begin{pmatrix} A_{11} & \dots & A_{1N} \\ \vdots & \ddots & \vdots \\ A_{M1} & \dots & A_{MN} \end{pmatrix} $$

where $M, N$ are the block dimensions, that is, the dimension expressed in terms of the subblocks. Usually, we choose the blocks such that $M = N$ and the diagonal blocks are square.
As a simple example, consider the matrix-vector product $y = Ax$, expressed in block terms.

$$ \begin{pmatrix} Y_1 \\ \vdots \\ Y_M \end{pmatrix} = \begin{pmatrix} A_{11} & \dots & A_{1M} \\ \vdots & \ddots & \vdots \\ A_{M1} & \dots & A_{MM} \end{pmatrix} \begin{pmatrix} X_1 \\ \vdots \\ X_M \end{pmatrix} $$

To see that the block algorithm computes the same result as the old scalar algorithm, we look at a component $Y_{ik}$, that is the $k$-th scalar component of the $i$-th block. First,

$$ Y_i = \sum_j A_{ij} X_j $$

so

$$ Y_{ik} = \left( \sum_j A_{ij} X_j \right)_k = \sum_j (A_{ij} X_j)_k = \sum_j \sum_\ell A_{ij\ell} X_{j\ell} $$

which is the product of the $k$-th row of the $i$-th blockrow of $A$ with the whole of $X$.

A more interesting algorithm is the block version of the LU factorization. The algorithm (5.3) then becomes

$$ \begin{array}{l} \langle LU \text{ factorization} \rangle: \\ \quad \text{for } k = 1, n-1: \\ \quad \quad \text{for } i = k+1 \text{ to } n: \\ \quad \quad \quad A_{ik} \leftarrow A_{ik} A_{kk}^{-1} \\ \quad \quad \text{for } j = k+1 \text{ to } n: \\ \quad \quad \quad A_{ij} \leftarrow A_{ij} - A_{ik} \cdot A_{kj} \end{array} \tag{5.7} $$

which mostly differs from the earlier algorithm in that the division by $a_{kk}$ has been replaced by a multiplication by $A_{kk}^{-1}$. Also, the $U$ factor will now have pivot blocks, rather than pivot elements, on the diagonal, so $U$ is only ‘block upper triangular’, and not strictly upper triangular.

**Exercise 5.16.** We would like to show that the block algorithm here again computes the same result as the scalar one. Doing so by looking explicitly at the computed elements is cumbersome, so we take another approach. First, recall from section 5.3.3 that $LU$ factorizations are unique: if $A = L_1 U_1 = L_2 U_2$ and $L_1, L_2$ have unit diagonal, then $L_1 = L_2, U_1 = U_2$.

Next, consider the computation of $A_{kk}^{-1}$. Show that this can be done easily by first computing an LU factorization of $A_{kk}$. Now use this to show that the block LU factorization can give $L$ and $U$ factors that are strictly triangular. The uniqueness of $LU$ factorizations then proves that the block algorithm computes the scalar result.

## 5.4 Sparse matrices

In section 4.2.3 you saw that the discretization of BVPs (and IBVPs) may give rise to sparse matrices. Since such a matrix has $N^2$ elements but only $O(N)$ nonzeros, it would be a big waste of space to store this as a two-dimensional array. Additionally, we want to avoid operating on zero elements.

In this section we will explore efficient storage schemes for sparse matrices, and the form that familiar linear algebra operations take when using sparse storage.
### 5.4.1 Storage of sparse matrices

It is pointless to look for an exact definition of *sparse matrix*, but an operational definition is that a matrix is called 'sparse' if there are enough zeros to make specialized storage feasible. We will discuss here briefly the most popular storage schemes for sparse matrices. Since a matrix is no longer stored as a simple 2-dimensional array, algorithms using such storage schemes need to be rewritten too.

#### 5.4.1.1 Diagonal storage

In section 4.2.2 you have seen examples of sparse matrices that were banded. In fact, their nonzero elements are located precisely on a number of subdiagonals. For such a matrix, a specialized storage scheme is possible.

Let us take as an example the matrix of the one-dimensional BVP (section 4.2.2). Its elements are located on three subdiagonals: the main diagonal and the first super and subdiagonal. The idea of *storage by diagonals* or *diagonal storage* is to store the diagonals consecutively in memory. The most economical storage scheme for such a matrix would store the $2n - 2$ elements consecutively. However, for various reasons it is more convenient to waste a few storage locations, as shown in figure 5.1.

<figure><img src="image_78.png" /><figcaption>Figure 5.1: Diagonal storage of a banded matrix</figcaption></figure>

Thus, for a matrix with size $n \times n$ and a *bandwidth* $p$, we need a rectangular array of size $n \times p$ to store the matrix. The matrix of equation (4.9) would then be stored as

$$ \begin{pmatrix} \star & 2 & -1 \\ -1 & 2 & -1 \\ \vdots & \vdots & \vdots \\ -1 & 2 & \star \end{pmatrix} \tag{5.8} $$

where the stars indicate array elements that do not correspond to matrix elements: they are the triangles in the top left and bottom right in figure 5.1.
Of course, now we have to wonder about the conversion between array elements A(i, j) and matrix elements A<sub>ij</sub>. This is easiest done in the Fortran language. If we allocate the array with

dimension A(n, -1:1)

then the main diagonal A<sub>ii</sub> is stored in A(*, 0). For instance, A(1, 0) ~ A<sub>11</sub>. The next location in the same row of the matrix A, A(1, 1) ~ A<sub>12</sub>. It is easy to see that together we have the conversion

$$
A(i, j) \sim A_{i,i+j}. \tag{5.9}
$$

Exercise 5.17. What is the reverse conversion, that is, what array location A(?, ?) does the matrix element A<sub>ij</sub> correspond to?

Exercise 5.18. If you are a C programmer, derive the conversion between matrix elements A<sub>ij</sub>
and array elements A[i][j].

If we apply this scheme to the matrix of the two-dimensional BVP (section 4.2.3), it becomes wasteful,
since we would be storing many zeros that exist inside the band. Therefore, we refine this scheme by
storing only the nonzero diagonals: if the matrix has *p* nonzero diagonals, we need an *n* × *p* array. For the
matrix of equation (4.14) this means:

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

Of course, we need an additional integer array telling us the locations of these nonzero diagonals.

In the preceding examples, the matrices had an equal number of nonzero diagonals above and below the
main diagonal. In general this need not be true. For this we introduce the concepts of

* left halfbandwidth: if A has a left halfbandwidth of *p* then *A*<sub>*ij*</sub> = 0 for *i* > *j* + *p*, and

* right halfbandwidth: if A has a right halfbandwidth of *p* then *A*<sub>*ij*</sub> = 0 for *j* > *i* + *p*.

5.4.1.2 Operations on diagonal storage

The most important operation on sparse matrices is the matrix-vector product. With a matrix stored by
diagonals, as described above, it is still possible to perform the ordinary rowwise or columnwise product
using the conversion formula (5.9)<fnref n="3" />. However, with a small bandwidth, this gives short vector lengths and
relatively high loop overhead, so it will not be efficient. It is possible do to much better than that.

If we look at how the matrix elements are used in the matrix-vector product, we see that the main diagonal
is used as

$$
y_i \leftarrow y_i + A_{ii} x_i,
$$

<footnote n="3">3. In fact, this is how Lapack banded routines work.</footnote>
the first superdiagonal is used as

$$y_i \leftarrow y_i + A_{ii+1}x_{i+1} \quad \text{for } i < n,$$

and the first subdiagonal as

$$y_i \leftarrow y_i + A_{ii-1}x_{i-1} \quad \text{for } i > 1.$$

In other words, the whole matrix-vector product can be executed in just three vector operations of length $n$ (or $n-1$), instead of $n$ inner products of length 3 (or 2).

```fortran
for diag = -diag_left, diag_right
    for loc = max(1,1-diag), min(n,n-diag)
        y(loc) = y(loc) + val(loc,diag) * x(loc+diag)
    end
end
```

**Exercise 5.19.** Write a routine that computes $y \leftarrow A^t x$ by diagonals. Implement it in your favourite language and test it on a random matrix.

**Exercise 5.20.** The above code fragment is efficient if the matrix is dense inside the band. This is not the case for, for instance, the matrix of two-dimensional BVPs; see section 4.2.3 and in particular equation (4.14). Write code for the matrix-vector product by diagonals that only uses the nonzero diagonals.

**Exercise 5.21.** Multiplying matrices is harder than multiplying a matrix times a vector. If matrix $A$ has left and halfbandwidth $p_A$, $q_A$, and matrix $B$ has $p_B$, $q_B$, what are the left and right halfbandwidth of $C = AB$? Assuming that an array of sufficient size has been allocated for $C$, write a routine that computes $C \leftarrow AB$.

### 5.4.1.3 Compressed row storage

If we have a sparse matrix that does not have a simple band structure, or where the number of nonzero diagonals becomes impractically large, we use the more general Compressed Row Storage (CRS) scheme. As the name indicates, this scheme is based on compressing all rows, eliminating the zeros; see figure 5.2. Since this loses the information what columns the nonzeros originally came from, we have to store this explicitly. Consider an example of a sparse matrix:

$$A = \begin{pmatrix}
10 & 0 & 0 & 0 & -2 & 0 \\
3 & 9 & 0 & 0 & 0 & 3 \\
0 & 7 & 8 & 7 & 0 & 0 \\
3 & 0 & 8 & 7 & 5 & 0 \\
0 & 8 & 0 & 9 & 9 & 13 \\
0 & 4 & 0 & 0 & 2 & -1
\end{pmatrix}. \qquad (5.10)$$

After compressing all rows, we store all nonzeros in a single real array. The column indices are similarly stored in an integer array, and we store pointers to where the columns start. Using 0-based indexing this gives:
<figure><img src="image_80.png" /><figcaption>Figure 5.2: Compressing a row of a sparse matrix in the CRS format</figcaption></figure>

<table><tr><td>val</td><td>10</td><td>-2</td><td>3</td><td>9</td><td>3</td><td>7</td><td>8</td><td>7</td><td>3 ... 9</td><td>13</td><td>4</td><td>2</td><td>-1</td></tr><tr><td>col_ind</td><td>0</td><td>4</td><td>0</td><td>1</td><td>5</td><td>1</td><td>2</td><td>3</td><td>0 ... 4</td><td>5</td><td>1</td><td>4</td><td>5</td></tr><tr><td>row_ptr</td><td>0</td><td>2</td><td>5</td><td>8</td><td>12</td><td>16</td><td>19</td><td></td><td></td><td></td><td></td><td></td><td></td></tr></table>

A simple variant of CRS is *Compressed Column Storage (CCS)* where the elements in columns are stored contiguously. This is also known as the *Harwell-Boeing matrix format* [48]. Another storage scheme you may come across is *coordinate storage*, where the matrix is stored as a list of triplets $\langle i, j, a_{ij} \rangle$. The popular *Matrix Market* website [129] uses a variant of this scheme.

### 5.4.1.4 Algorithms on compressed row storage

In this section we will look at the form some algorithms take in CRS.

First we consider the implementation of the sparse matrix-vector product.

```c
for (row=0; row<nrows; row++) {
    s = 0;
    for (icol=ptr資[row]; icol<ptr資[row+1]; icol++) {
        int col = ind[icol];
        s += a[icol] * x[col];
    }
    y資[row] = s;
}
```

You recognize the standard matrix-vector product algorithm for $y = Ax$, where the inner product is taken of each row $A_{i*}$ and the input vector $x$. However, note that the inner loop no longer has the column number as index, but rather the location where that number is to be found. This extra step is known as *indirect addressing*.
**Exercise 5.22.** Compare the data locality of the dense matrix-vector product, executed by rows, with the sparse product given just now. Show that, for general sparse matrices, the spatial locality in addressing the input vector $x$ has now disappeared. Are there matrix structures for which you can still expect some spatial locality?

Now, how about if you wanted to compute the product $y = A^t x$? In that case you need rows of $A^t$, or, equivalently, columns of $A$. Finding arbitrary columns of $A$ is hard, requiring lots of searching, so you may think that this algorithm is correspondingly hard to compute. Fortunately, that is not true.

If we exchange the $i$ and $j$ loop in the standard algorithm for $y = Ax$, we get

$$y \leftarrow 0$$

for $i$:

⇒

for $j$:

$$y \leftarrow 0$$

for $j$:

for $i$:

$$y_i \leftarrow y_i + a_{ij}x_j$$

$$y_i \leftarrow y_i + a_{ji}x_j$$

We see that in the second variant, columns of $A$ are accessed, rather than rows. This means that we can use the second algorithm for computing the $A^t x$ product by rows.

**Exercise 5.23.** Write out the code for the transpose product $y = A^t x$ where $A$ is stored in CRS format. Write a simple test program and confirm that your code computes the right thing.

**Exercise 5.24.** What if you need access to both rows and columns at the same time? Implement an algorithm that tests whether a matrix stored in CRS format is symmetric. Hint: keep an array of pointers, one for each row, that keeps track of how far you have progressed in that row.

**Exercise 5.25.** The operations described so far are fairly simple, in that they never make changes to the sparsity structure of the matrix. The CRS format, as described above, does not allow you to add new nonzeros to the matrix, but it is not hard to make an extension that does allow it.

Let numbers $p_i$, $i = 1...n$, describing the number of nonzeros in the $i$-th row, be given. Design an extension to CRS that gives each row space for $q$ extra elements. Implement this scheme and test it: construct a matrix with $p_i$ nonzeros in the $i$-th row, and check the correctness of the matrix-vector product before and after adding new elements, up to $q$ elements per row.

Now assume that the matrix will never have more than a total of $qn$ nonzeros. Alter your code so that it can deal with starting with an empty matrix, and gradually adding nonzeros in random places. Again, check the correctness.

We will revisit the transpose product algorithm in section 6.5.5 in the context of shared memory parallelism.

## 5.4.2 Sparse matrices and graph theory

Many arguments regarding sparse matrices can be formulated in terms of graph theory. To see why this can be done, consider a matrix $A$ of size $n$ and observe that we can define a graph $\langle E, V \rangle$ by $V = \{1, ..., n\}$, $E = \{(i, j): a_{ij} \neq 0\}$. This is called the adjacency graph of the matrix. For simplicity, we assume that $A$
has a nonzero diagonal. If necessary, we can attach weights to this graph, defined by $w_{ij} = a_{ij}$. The graph is then denoted $\langle E, V, W \rangle$. (If you are not familiar with the basics of graph theory, see appendix 16.)

Graph properties now correspond to matrix properties; for instance, the degree of the graph is the maximum number of nonzeros per row, not counting the diagonal element. As another example, if the graph of the matrix is an *undirected graph*, this means that $a_{ij} \neq 0 \Leftrightarrow a_{ji} \neq 0$. We call such a matrix *structurally symmetric*: it is not truly symmetric in the sense that $\forall_{ij}: a_{ij} = a_{ji}$, but every nonzero in the upper triangle corresponds to one in the lower triangle and vice versa.

One advantage of considering the graph of a matrix is that graph properties do not depend on how we order the nodes, that is, they are invariant under permutation of the matrix.

**Exercise 5.26.** Let us take a look at what happens with a matrix $A$ when the nodes of its graph $G = \langle V, E, W \rangle$ are renumbered. As a simple example, we number the nodes backwards; that is, with $n$ the number of nodes, we map node $i$ to $n + 1 - i$. Correspondingly, we find a new graph $G' = \langle V, E', W' \rangle$ where

$$ (i, j) \in E' \Leftrightarrow (n + 1 - i, n + 1 - j) \in E, \quad w'_{ij} = w_{n+1-i,n+1-j}. $$

What does this renumbering imply for the matrix $A'$ that corresponds to $G'$? If you exchange the labels $i, j$ on two nodes, what is the effect on the matrix $A$?

Some graph properties can be hard to see from the sparsity pattern of a matrix, but are easier deduced from the graph.

**Exercise 5.27.** Let $A$ be the tridiagonal matrix of the one-dimensional BVP (see section 4.2.2) of size $n$ with $n$ odd. What does the graph of $A$ look like? Consider the permutation that results from putting the nodes in the following sequence: 1, 3, 5, ..., $n$, 2, 4, 6, ..., $n-1$. What does the sparsity pattern of the permuted matrix look like? Renumbering strategies such as this will be discussed in more detail in section 6.8.2.

**Exercise 5.28.** Take again the matrix from the previous exercise, and zero the offdiagonal elements closest to the 'middle' of the matrix: let $a_{(n+1)/2, (n+1)/2+1} = a_{(n+1)/2+1, (n+1)/2} = 0$. Describe what that does to the graph of $A$. Such a graph is called *reducible*. Now apply the permutation of the previous exercise and sketch the resulting sparsity pattern. Note that the reducibility of the graph is now harder to read from the sparsity pattern.

### 5.4.3 LU factorizations of sparse matrices

In section 4.2.2 the one-dimensional BVP led to a linear system with a tridiagonal coefficient matrix. If we do one step of Gaussian elimination, the only element that needs to be eliminated is in the second row:

$$ \begin{pmatrix} 2 & -1 & 0 & \dots \\ -1 & 2 & -1 & \\ 0 & -1 & 2 & -1 \\ \vdots & \vdots & \vdots & \ddots \end{pmatrix} \Rightarrow \begin{pmatrix} 2 & -1 & 0 & \dots \\ 0 & 2 - \frac{1}{2} & -1 & \\ 0 & -1 & 2 & -1 \\ \vdots & \vdots & \vdots & \ddots \end{pmatrix} $$
There are two important observations to be made: one is that this elimination step does not change any zero elements to nonzero. The other observation is that the part of the matrix that is left to be eliminated is again tridiagonal. Inductively, during the elimination no zero elements change to nonzero: the sparsity pattern of $L+U$ is the same as of $A$, and so the factorization takes the same amount of space to store as the matrix.

The case of tridiagonal matrices is unfortunately not typical, as we will shortly see in the case of two-dimensional problems. But first we will extend the discussion on graph theory of section 5.4.2 to factorizations.

### 5.4.3.1 Graph theory of sparse LU factorization

Graph theory is often useful when discussing the factorization of a sparse matrix. Let us investigate what eliminating the first unknown (or sweeping the first column) means in graph theoretic terms. We are assuming a structurally symmetric matrix.

We consider eliminating an unknown as a process that takes a graph $G = \langle V, E \rangle$ and turns it into a graph $G' = \langle V', E' \rangle$. The relation between these graphs is first that a vertex, say $k$, has been removed from the vertices: $k \notin V'$, $V' \cup \{k\} = V$.

The relationship between $E$ and $E'$ is more complicated. In the Gaussian elimination algorithm the result of eliminating variable $k$ is that the statement

$$a_{ij} \leftarrow a_{ij} - a_{ik}a_{kk}^{-1}a_{kj}$$

is executed for all $i, j \neq k$. If $a_{ij} \neq 0$ originally, then its value is merely altered. In case $a_{ij} = 0$ in the original matrix, there will be a nonzero element, termed a fill-in element, after the $k$ unknown is eliminated:

<figure><img src="image_81.png" /><figcaption>Figure 5.3: Eliminating a vertex introduces a new edge in the quotient graph</figcaption></figure>

in $E$ there was no edge $(i, j)$, and this edge is present in $E'$. This is illustrated in figure 5.3.

Summarizing, eliminating an unknown gives a graph that has one vertex less, and that has edges for all $i, j$ such that there were edges between $i$ or $j$ and the eliminated variable $k$.

**Exercise 5.29.** Go back to exercise 5.27. Use a graph argument to determine the sparsity pattern after the odd variables have been eliminated.

**Exercise 5.30.** Prove the generalization of the above argument about eliminating a single vertex. Let $I \subset V$ be any set of vertices, and let $J$ be the vertices connected to $I$:

$$J \cap I = \emptyset, \quad \forall i \in I \exists j \in J: (i, j) \in E.$$
Now show that eliminating the variables in *I* leads to a graph $\langle V', E' \rangle$ where all nodes in *J* are connected in the remaining graph, if there was a path between them through *I*:

$$ \forall_{j_1, j_2 \in J} : \text{ there is a path } j_i \to j_2 \text{ through } I \text{ in } E \Rightarrow (j_1, j_2) \in E'. $$

### 5.4.3.2 Fill-in

We now return to the factorization of the matrix from two-dimensional problems. We write such matrices of size $N \times N$ as block matrices with block dimension $n$, each block being of size $n$. (Refresher question: where do these blocks come from?) Now, in the first elimination step we need to zero two elements, $a_{21}$ and $a_{n+1,1}$.

$$ \left( \begin{array}{cccc|c} 4 & -1 & 0 & \cdots & -1 \\ -1 & 4 & -1 & 0 & \cdots & 0 & -1 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\ -1 & 0 & \cdots & \cdots & -1 & \cdots & -1 \\ 0 & -1 & 0 & \cdots & -1 & 4 & -1 \end{array} \right) \quad \Rightarrow \quad \left( \begin{array}{c|cccc|c} 4 & -1 & 0 & \cdots & -1 \\ \hline 4 - \frac{1}{4} & -1 & 0 & \cdots & -1/4 & -1 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \ddots \\ -1/4 & -1 & 0 & \cdots & 4 - \frac{1}{4} & -1 \\ -1 & 0 & \cdots & \cdots & -1 & 4 & -1 \end{array} \right) $$

You see that eliminating $a_{21}$ and $a_{n+1,1}$ causes two fill elements to appear: in the original matrix $a_{2,n+1}$ and $a_{n+1,2}$ are zero, but in the modified matrix these locations are nonzero. We define fill locations as locations $(i, j)$ where $a_{ij} = 0$, but $(L+U)_{ij} \neq 0$.

Clearly the matrix fills in during factorization. With a little imagination you can also see that every element in the band outside the first diagonal block will fill in. However, using the graph approach of section 5.4.3.1 it becomes easy to visualize the fill-in connections that are created.

In figure 5.4 this is illustrated for the graph of the 2d BVP example. (The edges corresponding to diagonal elements have not been pictured here.) Each variable in the first row that is eliminated creates connections between the next variable and the second row, and between variables in the second row. Inductively you see that after the first row is eliminated the second row is fully connected. (Connect this to exercise 5.30.)

**Exercise 5.31.** Finish the argument. What does the fact that variables in the second row are fully connected imply for the matrix structure? Sketch in a figure what happens after the first variable in the second row is eliminated.

**Exercise 5.32.** The LAPACK software for dense linear algebra has an LU factorization routine that overwrites the input matrix with the factors. Above you saw that is possible since the columns of $L$ are generated precisely as the columns of $A$ are eliminated. Why is such an algorithm not possible if the matrix is stored in sparse format?

### 5.4.3.3 Fill-in estimates

In the above example you saw that the factorization of a sparse matrix can take much more space than the matrix itself, but still less than storing an entire square array of size the matrix dimension. We will now give some bounds for the space complexity of the factorization, that is, the amount of space needed to execute the factorization algorithm.
<figure><img src="image_82.png" /><figcaption>Figure 5.4: Creation of fill-in connection in the matrix graph</figcaption></figure>

**Exercise 5.33.** Prove the following statements.

1. Assume that the matrix $A$ has a halfbandwidth $p$, that is, $a_{ij} = 0$ if $|i-j| > p$. Show that, after a factorization without pivoting, $L+U$ has the same halfbandwidth.

2. Show that, after a factorization with partial pivoting, $L$ has a left halfbandwidth of $p$, whereas $U$ has a right halfbandwidth of $2p$.

3. Assuming no pivoting, show that the fill-in can be characterized as follows:

Consider row $i$. Let $j_{\min}$ be the leftmost nonzero in row $i$, that is $a_{ij} = 0$ for $j < j_{\min}$. Then there will be no fill-in in row $i$ to the left of column $j_{\min}$. Likewise, if $i_{\min}$ is the topmost nonzero in column $j$, there will be no fill-in in column $j$ above row $i_{\min}$.

As a result, $L$ and $U$ have a ‘skyline’ profile. Given a sparse matrix, it is now easy to allocate enough storage to fit a factorization without pivoting: this is known as *skyline storage*.

This exercise shows that we can allocate enough storage for the factorization of a banded matrix:

* for the factorization without pivoting of a matrix with bandwidth $p$, an array of size $N \times p$ suffices;

* the factorization with partial pivoting of a matrix left halfbandwidth $p$ and right halfbandwidth $q$ can be stored in $N \times (p + 2q + 1)$.

* A skyline profile, sufficient for storing the factorization, can be constructed based on the specific matrix.

We can apply this estimate to the matrix from the two-dimensional BVP, section 4.2.3.

**Exercise 5.34.** Show that in equation (4.14) the original matrix has $O(N) = O(n^2)$ nonzero elements, $O(N^2) = O(n^4)$ elements in total, and the factorization has $O(nN) = O(n^3) = O(N^{3/2})$ nonzeros.

These estimates show that the storage required for an *LU* factorization can be more than what is required for *A*, and the difference is not a constant factor, but related to the matrix size. Without proof we state that
the inverses of the kind of sparse matrices you have seen so far are fully dense, so storing them takes even more. This is an important reason that solving linear systems $Ax = y$ is not done in practice by computing $A^{-1}$ and multiplying $x = A^{-1}y$. (Numerical stability is another reason that this is not done.) The fact that even a factorization can take a lot of space is one reason for considering iterative methods, as we will do in section 5.5.

Above, you saw that the factorization of a dense matrix of size $n \times n$ takes $O(n^3)$ operations. How is this for a sparse matrix? Let us consider the case of a matrix with halfbandwidth $p$, and assume that the original matrix is dense in that band. The pivot element $a_{11}$ is used to zero $p$ elements in the first column, and for each the first row is added to that row, involving $p$ multiplications and additions. In sum, we find that the number of operations is roughly

$$ \sum_{i=1}^{n} p^2 = p^2 \cdot n $$

plus or minus lower order terms.

**Exercise 5.35.** The assumption of a band that is initially dense is not true for the matrix of a two-dimensional BVP. Why does the above estimate still hold, up to some lower order terms?

In exercise 5.33 above you derived an estimate for the amount of fill-in that is easy to apply. However, it can be a considerable overestimate. It is desirable to compute or estimate the amount of fill-in with less work than doing the actual factorization. We will now sketch an algorithm for finding the exact number of nonzeros in $L+U$, with a cost that is linear in this number. We will do this in the (structurally) symmetric case. The crucial observation is the following. Suppose column $i$ has more than one nonzero below the diagonal:

$$ \left( \begin{array}{ccc} \cdots & & \\ & a_{ii} & a_{ij} & a_{ik} \\ & \ddots & & \\ & & a_{ji} & a_{jj} \\ & & & \ddots \\ & & & & a_{ki} & ?a_{kj}? & a_{kk} \\ & & & & & & \end{array} \right) $$

Eliminating $a_{ki}$ in the $i$-th step causes an update of $a_{kj}$, or a fill-in element if originally $a_{kj} = 0$. However, we can infer the existence of this nonzero value: eliminating $a_{ji}$ causes a fill-in element in location $(j, k)$, and we know that structural symmetry is preserved. In other words, if we are only counting nonzeros, it is enough to look at the effects of eliminating the $(j, i)$ location, or in general the first nonzero below the pivot. Following this argument through, we only need to record the nonzeros in one row per pivot, and the entire process has a complexity linear in the number of nonzeros in the factorization.

### 5.4.3.4 Fill-in reduction

Graph properties of a matrix, such as degree and diameter, are invariant under renumbering the variables. Other properties, such as fill-in during a factorization, are affected by renumbering. In fact, it is worthwhile
investigating whether it is possible to reduce the amount of fill-in by renumbering the nodes of the matrix graph, or equivalently, by applying a *permutation* to the linear system.

**Exercise 5.36.** Consider the ‘arrow’ matrix with nonzeroes only in the first row and column and on the diagonal:

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

What is the number of nonzeroes in the matrix, and in the factorization, assuming that no addition ever results in zero? Can you find a symmetric permutation of the variables of the problem such that the new matrix has no fill-in?

This example is not typical, but it is true that fill-in estimates can sometimes be improved upon by clever permuting of the matrix (see for instance section 6.8.1). Even with this, as a rule the statement holds that an *LU* factorization of a sparse matrix takes considerably more space than the matrix itself. This is one of the motivating factors for the iterative methods in the next section.

### 5.4.3.5 Fill-in reducing orderings

Some matrix properties are invariant under symmetric permutations.

**Exercise 5.37.** In linear algebra classes, you typically look at matrix properties and whether they are invariant under a change of basis, in particular under *unitary basis transformations*:

$$B = VAV^t, \quad \text{where } VV^t = I.$$

Show that a symmetric permutation a particular change of basis is. Name some matrix properties that do not change under unitary transformations.

Other properties are not: in the previous section you saw that the amount of fill-in is one of those. Thus, you may wonder what the best ordering is to reduce the fill-in of factoring a given matrix. This problem is intractable in practice, but various heuristics exist. Some of these heuristics can also be justified from a point of view of parallelism; in fact, the *nested dissection* ordering will only be discussed in the section on parallelism 6.8.1. Here we briefly show two other heuristics that predate the need for parallelism.

The *Cuthill-McKee ordering* [30] orders the variables in *level sets*. It considers the *adjacency graph* of the matrix, and proceeds as follows:

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

2. Given a level $n$, assign all nodes that are connecting to level $n$, and that are not yet in a level, to level $n+1$.

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

**Exercise 5.38.** Show that permuting a matrix according to the Cuthill-McKee ordering has a *block tridiagonal* structure.
We will revisit this algorithm in section 6.10.1 when we consider parallelism.

Another ordering is motivated by the observation that the amount of fill-in is related to the degree of nodes.

Exercise 5.39. Show that eliminating a node with degree *d* leads to at most 2*d* fill elements
The so-called *minimum degree ordering* proceeds as follows:

* Find the node with lowest degree;

* eliminate that node and update the degree information for the remaining nodes;

* repeat from the first step, with the updated matrix graph.

Exercise 5.40. Indicate a difference between the two above-mentioned methods. Both are based on inspection of the matrix graph; however, the minimum degree method requires much more flexibility in the data structures used. Explain why and discuss two aspects in detail.

**5.5 Iterative methods**

Gaussian elimination, the use of an *LU* factorization, is a simple way to find the solution of a linear system, but as we saw above, in the sort of problems that come from discretized PDEs, it can create a lot of fill-in. In this section we will look at a completely different approach, where the solution of the system is found by a sequence of approximations.

The computational scheme looks, very roughly, like:

$$
\left\{
\begin{array}{l}
\text{Choose any starting vector } x_0 \text{ and repeat for } i \ge 0: \\
x_{i+1} = Bx_i + c \\
\quad \text{until some stopping test is satisfied.}
\end{array}
\right.
$$

The important feature here is that no systems are solved with the original coefficient matrix; instead, every
iteration involves a matrix-vector multiplication or a solution of a much simpler system. Thus we have
replaced a complicated operation, constructing an *LU* factorization and solving a system with it, by a
repeated simpler and cheaper operation. This makes iterative methods easier to code, and potentially more
efficient.

Let us consider a simple example to motivate the precise definition of the iterative methods. Suppose we
want to solve the system

$$
\begin{pmatrix} 10 & 0 & 1 \\ \frac{1}{2} & 7 & 1 \\ 1 & 0 & 6 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 21 \\ 9 \\ 8 \end{pmatrix}
$$

which has the solution (2, 1, 1). Suppose you know (for example, from physical considerations) that so-
lution components are roughly the same size. Observe the dominant size of the diagonal, then, to decide
that

$$
\begin{pmatrix} 10 & & \\ & 7 & \\ & & 6 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 21 \\ 9 \\ 8 \end{pmatrix}
$$
might be a good approximation. This has the solution (2.1, 9/7, 8/6). Clearly, solving a system that only involves the diagonal of the original system is both easy to do, and, at least in this case, fairly accurate.

Another approximation to the original system would be to use the lower triangle. The system

$$ \begin{pmatrix} 10 \\ \frac{1}{2} & 7 \\ 1 & 0 & 6 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix} = \begin{pmatrix} 21 \\ 9 \\ 8 \end{pmatrix} $$

has the solution (2.1, 7.95/7, 5.9/6). Solving triangular systems is a bit more work than diagonal systems, but still a lot easier than computing an LU factorization. Also, we have not generated any fill-in in the process of finding this approximate solution.

Thus we see that there are easy to compute ways of getting reasonably close to the solution. Can we somehow repeat this trick?

Formulated a bit more abstractly, what we did was instead of solving $Ax = b$ we solved $L\tilde{x} = b$. Now define $\Delta x$ as the distance from the true solution: $\tilde{x} = x + \Delta x$. This gives $A\Delta x = A\tilde{x} - b \equiv r$, where $r$ is the residual. Next we solve again $L\tilde{\Delta}x = r$ and update $\tilde{x} = \tilde{x} - \tilde{\Delta}x$.

<table><thead><tr><th>iteration</th><th>1</th><th>2</th><th>3</th></tr></thead><tbody><tr><th>$x_1$</th><td>2.1000</td><td>2.0017</td><td>2.000028</td></tr><tr><th>$x_2$</th><td>1.1357</td><td>1.0023</td><td>1.000038</td></tr><tr><th>$x_3$</th><td>0.9833</td><td>0.9997</td><td>0.999995</td></tr></tbody></table>

In this case we get two decimals per iteration, which is not typical.

It is now clear why iterative methods can be attractive. Solving a system by Gaussian elimination takes $O(n^3)$ operations, as shown above. A single iteration in a scheme as the above takes $O(n^2)$ operations if the matrix is dense, and possibly as low as $O(n)$ for a sparse matrix. If the number of iterations is low, this makes iterative methods competitive.

**Exercise 5.41.** When comparing iterative and direct methods, the flop count is not the only relevant measure. Outline some issues relating to the efficiency of the code in both cases. Also compare the cases of solving a single linear system and solving multiple.

### 5.5.1 Abstract presentation

It is time to do a formal presentation of the iterative scheme of the above example. Suppose we want to solve $Ax = b$, and a direct solution is too expensive, but multiplying by $A$ is feasible. Suppose furthermore that we have a matrix $K \approx A$ such that solving $Kx = b$ can be done cheaply.

Instead of solving $Ax = b$ we solve $Kx = b$, and define $x_0$ as the solution: $Kx_0 = b$. This leaves us with an error $e_0 = x_0 - x$, for which we have the equation $A(x_0 - e_0) = b$ or $Ae_0 = Ax_0 - b$. We call $r_0 \equiv Ax_0 - b$ the residual; the error then satisfies $Ae_0 = r_0$.

If we could solve the error from the equation $Ae_0 = r_0$, we would be done: the true solution is then found as $x = x_0 - e_0$. However, since solving with $A$ was too expensive the last time, we can not do so this time either, so we determine the error correction approximately. We solve $K\tilde{e}_0 = r_0$ and set $x_1 := x_0 - \tilde{e}_0$; the story can now continue with $e_1 = x_1 - x$, $r_1 = Ax_1 - b$, $K\tilde{e}_1 = r_1$, $x_2 = x_1 - \tilde{e}_1$, et cetera.
The iteration scheme is then:

Let $x_0$ be given

For $i \ge 0$:

let $r_i = Ax_i - b$

compute $e_i$ from $Ke_i = r_i$

update $x_{i+1} = x_i - e_i$

We call the basic scheme

$$x_{i+1} = x_i - K^{-1}r_i \qquad (5.11)$$

a stationary iteration. It is stationary because every update is performed the same way, without any dependence on the iteration number. This scheme has a simple analysis, but unfortunately limited applicability.

There are several questions we need to answer about iterative schemes:

* Does this scheme always take us to the solution?

* If the scheme converges, how quickly?

* When do we stop iterating?

* How do we choose $K$?

We will now devote some attention to these matters, though a full discussion is beyond the scope of this book.

### 5.5.2 Convergence and error analysis

We start with the question of whether the iterative scheme converges, and how quickly. Consider one iteration step:

$$ \begin{aligned} r_1 &= Ax_1 - b = A(x_0 - \tilde{e}_0) - b \\ &= r_0 - AK^{-1}r_0 \\ &= (I - AK^{-1})r_0 \end{aligned} \qquad (5.12) $$

Inductively we find $r_n = (I - AK^{-1})^n r_0$, so $r_n \downarrow 0$ if all eigenvalues satisfy $|\lambda(I - AK^{-1})| < 1^4$.

This last statement gives us both a condition for convergence, by relating $K$ to $A$, and a geometric convergence rate, if $K$ is close enough.

**Exercise 5.42.** Derive a similar inductive relation for $e_n$.

It is hard to determine if the condition $|\lambda(I - AK^{-1})| < 1$ is satisfied by computing the actual eigenvalues. However, sometimes the Gershgorin theorem (appendix 12.5) gives us enough information.

**Exercise 5.43.** Consider the matrix $A$ of equation (4.14) that we obtained from discretization of a two-dimensional BVP. Let $K$ be matrix containing the diagonal of $A$, that is $k_{ii} = a_{ii}$ and $k_{ij} = 0$ for $i \ne j$. Use the Gershgorin theorem to show that $|\lambda(I - AK^{-1})| < 1$.

<footnote n="4">4. This is fairly easy to see in the case where the matrix is diagonalizable and has a full basis of eigenvectors. However, it is true in the general case too.</footnote>
The argument in this exercise is hard to generalize for more complicated choices of $K$, such as you will see below. Here we only remark that for certain matrices $A$, these choices of $K$ will always lead to convergence, with a speed that decreases as the matrix size increases. We will not go into the details beyond stating that for $M$-matrices (see section 4.2.2) these iterative methods converge. For more details on the convergence theory of stationary iterative methods, see [159]

### 5.5.3 Computational form

Above, in section 5.5.1, we derived stationary iteration as a process that involves multiplying by $A$ and solving with $K$. However, in some cases a simpler implementation is possible. Consider the case where $A = K - N$, and we know both $K$ and $N$. Then we write $Ax = b$ as

$$ Kx = Nx + b \tag{5.13} $$

and we observe that an $x$ satisfying (5.13) is a fixed point of the iteration

$$ Kx^{(n+1)} = Nx^{(n)} + b. $$

It is easy to see that this is a stationary iteration:

$$ \begin{aligned} Kx^{(n+1)} &= Nx^{(n)} + b \\ &= Kx^{(n)} - Ax^{(n)} + b \\ &= Kx^{(n)} - r^{(n)} \\ \Rightarrow x^{(n+1)} &= x^{(n)} - K^{-1}r^{(n)}. \end{aligned} $$

which is the basic form of equation (5.11). The convergence criterium $|\lambda(I - AK^{-1})| < 1$ (see above) now simplifies to $|\lambda(NK^{-1})| < 1$.

Let us consider some special cases. First of all, let $K = D_A$, that is, the matrix containing the diagonal part of $A$: $k_{ii} = a_{ii}$ and $k_{ij} = 0$ for all $i \neq j$. Likewise, $n_{ii} = 0$ and $n_{ij} = -a_{ij}$ for all $i \neq j$.

This is known as the *Jacobi method*. The iteration scheme $Kx^{(n+1)} = Nx^{(n)} + b$ now becomes

for $t = 1, \dots$ until convergence, do:

for $i = 1 \dots n$:

$$ \begin{aligned} // a_{ii}x_i^{(t+1)} &= \sum_{j \neq i} a_{ij}x_j^{(t)} + b_i \text{ becomes:} \\ x_i^{(t+1)} &= a_{ii}^{-1} \left( \sum_{j \neq i} a_{ij}x_j^{(t)} + b_i \right) \end{aligned} $$

(Bearing in mind that divisions are relatively costly, section 1.2, we would actually store the $a_{ii}^{-1}$ quantities explicitly, and replace the division by a multiplication.)

This requires us to have one vector $x$ for the current iterate $x^{(t)}$, and one temporary $u$ for the next vector $x^{(t+1)}$. The easiest way to write this is probably:

for $t = 1, \dots$ until convergence, do:

for $i = 1 \dots n$:

$$ u_i = a_{ii}^{-1} \left( - \sum_{j \neq i} a_{ij} x_j^{(t)} + b_i \right) $$

copy $x \leftarrow u$
<figure><img src="image_84.png" /><figcaption>Figure 5.5: Data movement pattern in the Jacobi iteration on a one-dimensional problem</figcaption></figure>

For the simple case of a one-dimensional problem this is illustrated in figure 5.5: in each $x_i$ point the values of the two neighbours are combined with the current value to generate a new value. Since the computations in all the $x_i$ points are independent, this can be done in parallel on a parallel computer.

But, you might think, in the sum $\sum_{j \neq i} a_{ij}x_j$ why not use the $x^{(t+1)}$ values for as far as already computed? In terms of the vectors $x^{(t)}$ this means

for $k = 1, \dots$ until convergence, do:

for $i = 1, \dots, n$:

$$x_i^{(t+1)} = a_{ii}^{-1} \left( -\sum_{j<i} a_{ij}x_j^{(t+1)} - \sum_{j>i} a_{ij}x_j^{(t)} + b_i \right)$$

Surprisingly, the implementation is simpler than the Jacobi method:

for $t = 1, \dots$ until convergence, do:

for $i = 1, \dots, n$:

$$x_i = a_{ii}^{-1} \left( -\sum_{j \neq i} a_{ij}x_j + b_i \right)$$

If you write this out as a matrix equation, you see that the newly computed elements $x_i^{(t+1)}$ are multiplied with elements of $D_A + L_A$, and the old elements $x_j^{(t)}$ by $U_A$, giving

$$(D_A + L_A)x^{(k+1)} = -U_Ax^{(k)} + b$$

which is called the Gauss-Seidel method.

For the one-dimensional case, the Gauss-Seidel method is illustrated in figure 5.6; every $x_i$ point still

<figure><img src="image_85.png" /><figcaption>Figure 5.6: Data movement pattern in the Gauss-Seidel iteration on a one-dimensional problem</figcaption></figure>

combines its neighbours' values, but now the left value is actually from the next outer iteration. Since the computations in the $x_i$ points are now dependent, this type of iteration is not simple to implement on a parallel computer. Some strategies for dealing with this are discussed in section 6.7.
Finally, we can insert a damping parameter into the Gauss-Seidel scheme, giving the Successive Over-Relaxation (SOR) method:

for $t = 1, \dots$ until convergence, do:

for $i = 1 \dots n$:

$$x_i^{(t+1)} = \omega a_{ii}^{-1} \left( -\sum_{j<i} a_{ij} x_j^{(t+1)} - \sum_{j>i} a_{ij} x_j^{(t)} + b_i \right) + (1-\omega)x^{(t)}$$

Surprisingly for something that looks like an interpolation, the method actually works with value for $\omega$ in the range $\omega \in (0, 2)$, the optimal value being larger than 1 [79]. Computing the optimal $\omega$ is not simple.

### 5.5.4 Convergence of the method

We are interested in two questions: firstly whether the iterative method converges at all, and if so, with what speed. The theory behind these questions goes far beyond this book. Above we remarked that convergence can often be guaranteed for M-matrices; with regard to the convergence speed a full analysis is usually only possible in model cases. For the matrices from BVPs, as described in section 4.2.3, we state without proof that the smallest eigenvalue of the coefficient matrix is $O(h^2)$. The geometric convergence ratio $|\lambda(I - AK^{-1})|$ derived above can then be shown to be as follows:

* For the Jacobi method, the ratio is $1 - O(h^2)$;

* For the Gauss-Seidel iteration it also is $1 - O(h^2)$, but the method converges twice as fast;

* For the SOR method, the optimal omega can improve the convergence factor to $1 - O(h)$.

### 5.5.5 Choice of K

The convergence and error analysis above showed that the closer $K$ is to $A$, the faster the convergence will be. In the initial examples we already saw the diagonal and lower triangular choice for $K$. We can describe these formally by letting $A = D_A + L_A + U_A$ be a splitting into diagonal, lower triangular, upper triangular part of $A$. Here are some methods with their traditional names:

* Richardson iteration: $K = \alpha I$.

* Jacobi method: $K = D_A$ (diagonal part),

* Gauss-Seidel method: $K = D_A + L_A$ (lower triangle, including diagonal)

* The SOR method: $K = \omega^{-1}D_A + L_A$

* Symmetric SOR (SSOR) method: $K = (D_A + L_A)D_A^{-1}(D_A + U_A)$.

* Iterative refinement: $K = LU$ where $LU$ is a true factorization of $A$. In exact arithmetic, solving a system $LUx = y$ gives you the exact solution, so using $K = LU$ in an iterative method would give convergence after one step. In practice, roundoff error will make the solution be inexact, so people will sometimes iterate a few steps to get higher accuracy.

**Exercise 5.44.** What is the extra cost of a few steps of iterative refinement over a single system solution, assuming a dense system?

**Exercise 5.45.** The Jacobi iteration for the linear system $Ax = b$ is defined as

$$x_{i+1} = x_i - K^{-1}(Ax_i - b)$$
where $K$ is the diagonal of $A$. Show that you can transform the linear system (that is, find a different coefficient matrix and right hand side vector that will still have the same solution) so that you can compute the same $x_i$ vectors but with $K = I$, the identity matrix.

What are the implications of this strategy, in terms of storage and operation counts? Are there special implications if $A$ is a sparse matrix?

Suppose $A$ is symmetric. Give a simple example to show that $K^{-1}A$ does not have to be symmetric. Can you come up with a different transformation of the system so that symmetry of the coefficient matrix is preserved and that has the same advantages as the transformation above? You can assume that the matrix has positive diagonal elements.

**Exercise 5.46.** Show that the transformation of the previous exercise can also be done for the Gauss-Seidel method. Give several reasons why this is not a good idea.

There are many different ways of choosing the preconditioner matrix $K$. Some of them are defined algebraically, such as the incomplete factorization discussed below. Other choices are inspired by the differential equation. For instance, if the operator is

$$ \frac{\delta}{\delta x}(a(x, y) \frac{\delta}{\delta x}u(x, y)) + \frac{\delta}{\delta y}(b(x, y) \frac{\delta}{\delta y}u(x, y)) = f(x, y) $$

then the matrix $K$ could be derived from the operator

$$ \frac{\delta}{\delta x}(\tilde{a}(x) \frac{\delta}{\delta x}u(x, y)) + \frac{\delta}{\delta y}(\tilde{b}(y) \frac{\delta}{\delta y}u(x, y)) = f(x, y) $$

for some choices of $\tilde{a}, \tilde{b}$. The second set of equations is called a separable problem, and there are fast solvers for them, meaning that they have $O(N \log N)$ time complexity; see [163].

### 5.5.5.1 Constructing $K$ as an incomplete LU factorization

We briefly mention one other choice of $K$, which is inspired by Gaussian elimination. As in Gaussian elimination, we let $K = LU$, but now we use an *Incomplete LU (ILU)* factorization. Remember that a regular $LU$ factorization is expensive because of the fill-in phenomenon. In an incomplete factorization, we limit the fill-in artificially.

If we write Gauss elimination as

$$ \begin{array}{l} \text{for } k, i, j: \\ \quad a[i, j] = a[i, j] - a[i, k] * a[k, j] / a[k, k] \end{array} $$

we define an incomplete variant by

$$ \begin{array}{l} \text{for } k, i, j: \\ \quad \text{if } a[i, j] \text{ not zero:} \\ \quad \quad a[i, j] = a[i, j] - a[i, k] * a[k, j] / a[k, k] \end{array} $$
* The resulting factorization is no longer exact: $LU \approx A$, so it is called an Incomplete LU (ILU) factorization.

* An ILU factorization takes much less space than a full factorization: the sparsity of $L+U$ is the same as of $A$.

The algorithm above is called 'ILU(0)', where the zero refers to the fact that absolutely no fill-in is allowed during the incomplete factorization. Other schemes that allow a limited amount of fill-in exist. Much more can be said about this method; we will only remark that for *M*-matrices this scheme typically gives a converging method [122].

**Exercise 5.47.** How do operation counts of the matrix-vector product and solving a system with an ILU factorization compare?

You have seen that a full factorization of sparse matrix can need a higher order storage ($N^{3/2}$ for the factorization versus $N$ for the matrix), but that an incomplete factorization takes $O(N)$, just like the matrix. It may therefore come as a surprise that the error matrix $R = A - LU$ is not dense, but itself sparse.

**Exercise 5.48.** Let $A$ be the matrix of the Poisson equation, $LU$ an incomplete factorization, and $R = A - LU$. Show that $R$ is a bi-diagonal matrix:

* Consider that $R$ consists of those elements that are discarded during the factorization. Where are they located in the matrix?

* Alternatively, write out the sparsity pattern of the product $LU$ and compare that to the sparsity pattern of $A$.

### 5.5.5.2 Cost of constructing a preconditioner

In the example of the heat equation (section 4.3) you saw that each time step involves solving a linear system. As an important practical consequence, any setup cost for solving the linear system, such as constructing the preconditioner, will be amortized over the sequence of systems that is to be solved. A similar argument holds in the context of nonlinear equations, a topic that we will not discuss as such. Nonlinear equations are solved by an iterative process such as the *Newton method*, which in its multidimensional form leads to a sequence of linear systems. Although these have different coefficient matrices, it is again possible to amortize setup costs by reusing the preconditioner for a number of Newton steps.

### 5.5.5.3 Parallel preconditioners

Constructing and using a preconditioner is a balancing act of many considerations: a more accurate preconditioner leads to convergence in fewer iterations, but these iterations can be more expensive; in addition, a more accurate preconditioner can be more costly to construct. In parallel, this story is even more complicated, because certain preconditioners are not very parallel to begin with. Therefore, we may accept a preconditioner that is parallel, but that has a worse iteration count than a serial preconditioner. For more discussion, see section 6.7.

## 5.5.6 Stopping tests

The next question we need to tackle is when to stop iterating. Above we saw that the error decreases geometrically, so clearly we will never reach the solution exactly, even if that were possible in computer
arithmetic. Since we only have this relative convergence behaviour, how do we know when we are close enough?

We would like the error $e_i = x - x_i$ to be small, but measuring this is impossible. Above we observed that $Ae_i = r_i$, so

$$ \|e_i\| \le \|A^{-1}\|\|r_i\| \le \lambda_{\max}(A^{-1})\|r_i\| $$

If we know anything about the eigenvalues of $A$, this gives us a bound on the error. (The norm of $A$ is only the largest eigenvalue for symmetric $A$. In general, we need singular values here.)

Another possibility is to monitor changes in the computed solution. If these are small:

$$ \|x_{n+1} - x_n\| / \|x_n\| < \epsilon $$

we can also conclude that we are close to the solution.

**Exercise 5.49.** Prove an analytic relationship between the distance between iterates and the distance to the true solution. If your equation contains constants, can they be determined theoretically or in practice?

**Exercise 5.50.** Write a simple program to experiment with linear system solving. Take the matrix from the 1D BVP (use an efficient storage scheme) and program an iterative method using the choice $K = D_A$. Experiment with stopping tests on the residual and the distance between iterates. How does the number of iterations depend on the size of the matrix?

Change the matrix construction so that a certain quantity is added the diagonal, that is, add $\alpha I$ to the original matrix. What happens when $\alpha > 0$? What happens when $\alpha < 0$? Can you find the value where the behaviour changes? Does that value depend on the matrix size?

### 5.5.7 Theory of general iterative methods

Above, you saw iterative methods of the form $x_{i+1} = x_i - K^{-1}r_i$, and we will now see iterative methods of the more general form

$$ x_{i+1} = x_i + \sum_{j \le i} K^{-1} r_j \alpha_{ji}, \qquad (5.14) $$

that is, using all previous residuals to update the iterate. One might ask, ‘why not introduce an extra parameter and write $x_{i+1} = \alpha_i x_i + \cdots$?’ Here we give a short argument that the former scheme describes a large class of methods. Indeed, the current author is not aware of methods that fall outside this scheme.

We defined the residual, given an approximate solution $\tilde{x}$, as $\tilde{r} = A\tilde{x} - b$. For this general discussion we precondition the system as $K^{-1}Ax = K^{-1}b$. (See section 5.5.5 where we discussed transforming the linear system.) The corresponding residual for the initial guess $\tilde{x}$ is

$$ \tilde{r} = K^{-1} A \tilde{x} - K^{-1} b. $$
We now find that

$$x = A^{-1}b = \tilde{x} - A^{-1}K\tilde{r} = \tilde{x} - (K^{-1}A)^{-1}\tilde{r}.$$

Now, the *Cayley-Hamilton theorem* states that for every $A$ there exists a polynomial $\phi(x)$ (the characteristic polynomial) such that

$$\phi(A) = 0.$$

We observe that we can write this polynomial $\phi$ as

$$\phi(x) = 1 + x\pi(x)$$

where $\pi$ is another polynomial. Applying this to $K^{-1}A$, we have

$$0 = \phi(K^{-1}A) = I + K^{-1}A\pi(K^{-1}A) \Rightarrow (K^{-1}A)^{-1} = -\pi(K^{-1}A)$$

so that $x = \tilde{x} + \pi(K^{-1}A)\tilde{r}$. Now, if we let $x_0 = \tilde{x}$, then $\tilde{r} = K^{-1}r_0$, giving the equation

$$x = x_0 + \pi(K^{-1}A)K^{-1}r_0.$$

This equation suggests an iterative scheme: if we can find a series of polynomials $\pi^{(i)}$ of degree $i$ to approximate $\pi$, it will give us a sequence of iterates

$$x_{i+1} = x_0 + \pi^{(i)}(K^{-1}A)K^{-1}r_0 = x_0 + K^{-1}\pi^{(i)}(AK^{-1})r_0 \quad (5.15)$$

that ultimately reaches the true solution. Multiplying this equation by $A$ and subtracting $b$ on both sides gives

$$r_{i+1} = r_0 + \tilde{\pi}^{(i)}(AK^{-1})r_0$$

where $\tilde{\pi}^{(i)}(x) = x\pi^{(i)}(x)$. This immediately gives us

$$r_i = \hat{\pi}^{(i)}(AK^{-1})r_0 \quad (5.16)$$

where $\hat{\pi}^{(i)}$ is a polynomial of degree $i$ with $\hat{\pi}^{(i)}(0) = 1$. This statement can be used as the basis of a convergence theory of iterative methods. However, this goes beyond the scope of this book.

Let us look at a couple of instances of equation (5.16). For $i = 1$ we have

$$r_1 = (\alpha_1 AK^{-1} + \alpha_2 I)r_0 \Rightarrow AK^{-1}r_0 = \beta_1 r_1 + \beta_0 r_0$$

for some values $\alpha_i, \beta_i$. For $i = 2$

$$r_2 = (\alpha_2(AK^{-1})^2 + \alpha_1 AK^{-1} + \alpha_0)r_0$$

for different values $\alpha_i$. But we had already established that $AK_0^{-1}$ is a combination of $r_1, r_0$, so now we have that

$$(AK^{-1})^2 r_0 \in [r_2, r_1, r_0],$$
and it is clear how to show inductively that

$$ (AK^{-1})^i r_0 \in [r_i, \dots, r_0]. \qquad (5.17) $$

Substituting this in (5.15) we finally get

$$ x_{i+1} = x_0 + \sum_{j \le i} K^{-1} r_j \alpha_{ji}. \qquad (5.18) $$

It is easy to see that the scheme (5.14) is of the form (5.18) and that the reverse implication also holds.

Summarizing, the basis of iterative methods is a scheme where iterates get updated by all residuals computed so far:

$$ x_{i+1} = x_i + \sum_{j \le i} K^{-1} r_j \alpha_{ji}. \qquad (5.19) $$

Compare that to the stationary iteration (section 5.5.1) where the iterates get updated from just the last residual, and with a coefficient that stays constant.

We can say more about the $\alpha_{ij}$ coefficients. If we multiply equation (5.19) by A, and subtract b from both sides, we find

$$ r_{i+1} = r_i + \sum_{j \le i} AK^{-1}r_j \alpha_{ji}. \qquad (5.20) $$

Let us consider this equation for a moment. If we have a starting residual $r_0$, the next residual is computed as

$$ r_1 = r_0 + AK^{-1}r_0 \alpha_{00}. $$

From this we get that $AK^{-1}r_0 = \alpha_{00}^{-1}(r_1 - r_0)$, so for the next residual,

$$ \begin{aligned} r_2 &= r_1 + AK^{-1}r_1 \alpha_{11} + AK^{-1}r_0 \alpha_{01} \\ &= r_1 + AK^{-1}r_1 \alpha_{11} + \alpha_{00}^{-1} \alpha_{01} (r_1 - r_0) \\ \Rightarrow AK^{-1}r_1 &= \alpha_{11}^{-1} (r_2 - (1 + \alpha_{00}^{-1} \alpha_{01})r_1 + \alpha_{00}^{-1} \alpha_{01} r_0) \end{aligned} $$

We see that we can express $AK^{-1}r_1$ as a sum $r_2\beta_2 + r_1\beta_1 + r_0\beta_0$, and that $\sum_i \beta_i = 0$.

Generalizing this, we find (with different $\alpha_{ij}$ than above)

$$ \begin{aligned} r_{i+1} &= r_i + AK^{-1}r_i \delta_i + \sum_{j \le i+1} r_j \alpha_{ji} \\ r_{i+1}(1 - \alpha_{i+1,i}) &= AK^{-1}r_i \delta_i + r_i(1 + \alpha_{ii}) + \sum_{j<i} r_j \alpha_{ji} \end{aligned} $$

$$ r_{i+1} \alpha_{i+1, i} = AK^{-1} r_i \delta_i + \sum_{j \le i} r_j \alpha_{ji} $$

$$ \begin{array}{l} \text{substituting} \quad \alpha_{ii} := 1 + \alpha_{ii} \\ \alpha_{i+1,i} := 1 - \alpha_{i+1,i} \end{array} $$

note that $\alpha_{i+1,i} = \sum_{j \le i} \alpha_{ji}$

$$ r_{i+1} \alpha_{i+1, i} \delta_i^{-1} = AK^{-1} r_i + \sum_{j \le i} r_j \alpha_{ji} \delta_i^{-1} $$

$$ r_{i+1} \alpha_{i+1, i} \delta_i^{-1} = AK^{-1} r_i + \sum_{j \le i} r_j \alpha_{ji} \delta_i^{-1} $$

$$ r_{i+1} \gamma_{i+1, i} = AK^{-1} r_i + \sum_{j \le i} r_j \gamma_{ji} $$

$$ \text{substituting } \gamma_{ij} = \alpha_{ij} \delta_j^{-1} $$
and we have that $\gamma_{i+1,i} = \sum_{j \le i} \gamma_{ji}$.

We can take this last equation and write it as $AK^{-1}R = RH$ where

$$H = \begin{pmatrix} -\gamma_{11} & -\gamma_{12} & \cdots & & \\ \gamma_{21} & -\gamma_{22} & -\gamma_{23} & \cdots & \\ 0 & \gamma_{32} & -\gamma_{33} & -\gamma_{34} & \\ \vdots & \vdots & \ddots & \ddots & \vdots \end{pmatrix}$$

In this, $H$ is a so-called Hessenberg matrix: it is upper triangular plus a single lower subdiagonal. Also we note that the elements of $H$ in each column sum to zero.

Because of the identity $\gamma_{i+1,i} = \sum_{j \le i} \gamma_{ji}$ we can subtract $b$ from both sides of the equation for $r_{i+1}$ and ‘divide out $A$, giving

$$x_{i+1}\gamma_{i+1,i} = K^{-1}r_i + \sum_{j \le i} x_j\gamma_{ji}.$$

This gives us the general form for iterative methods:

$$\begin{cases} r_i = Ax_i - b \\ x_{i+1}\gamma_{i+1,i} = K^{-1}r_i + \sum_{j \le i} x_j\gamma_{ji} \\ r_{i+1}\gamma_{i+1,i} = AK^{-1}r_i + \sum_{j \le i} r_j\gamma_{ji} \end{cases} \quad \text{where } \gamma_{i+1,i} = \sum_{j \le i} \gamma_{ji}. \qquad (5.21)$$

This form holds for many iterative methods, including the stationary iterative methods you have seen above. In the next sections you will see how the $\gamma_{ij}$ coefficients follow from orthogonality conditions on the residuals.

### 5.5.8 Iterating by orthogonalization

The stationary methods described above (section 5.5.1) have been around in some form or another for a long time: Gauss described some variant in a letter to a student. They were perfected in the thesis of Young [170] in 1950; the final reference is probably the book by Varga [159]. These methods are little used these days, except in the specialized context of multigrid smoothers, a topic not discussed in this course.

At almost the same time, the field of methods based on orthogonalization was kicked off by two papers [107, 85], though it took a few decades for them to find wide applicability. (For further history, see [69].)

The basic idea is as follows:

If you can make all your residuals orthogonal to each other, and the matrix is of dimension $n$, then after $n$ iterations you have to have converged: it is not possible to have an $n+1$-st residual that is orthogonal to all previous and nonzero. Since a zero residual means that the corresponding iterate is the solution, we conclude that after $n$ iterations we have the true solution in hand.

With the size of matrices that contemporary applications generate this reasoning is no longer relevant: it is not computationally realistic to iterate for $n$ iterations. Moreover, roundoff will probably destroy any accuracy of the solution. However, it was later realized [141] that such methods are a realistic option in the case of symmetric positive definite (SPD) matrices. The reasoning is then:
The sequence of residuals spans a series of subspaces of increasing dimension, and by orthogonalizing, the new residuals are projected on these spaces. This means that they will have decreasing sizes.

<figure><img src="image_86.png" /><figcaption>Figure 5.7: The optimal update $u_m$ make the new residual orthogonal to the $AK_m$ subspace</figcaption></figure>

This is illustrated in figure 5.7.

In this section you will see the basic idea of iterating by orthogonalization. The method presented here is only of theoretical interest; next you will see the Conjugate Gradients (CG) and Generalized Minimum Residual (GMRES) methods that are the basis of many real-life applications.

Let us now take the basic scheme (5.21) and orthogonalize the residuals. Instead of the normal inner product we use the $K^{-1}$-inner product:

$$ (x, y)_{K^{-1}} = x^T K^{-1} y $$

and we will force residuals to be $K^{-1}$-orthogonal:

$$ \forall i \neq j: r_i \perp_{K^{-1}} r_j \Leftrightarrow \forall i \neq j: r_i K^{-1} r_j = 0 $$

This is known as the Full Orthogonalization Method (FOM) scheme:

Let $r_0$ be given

For $i \ge 0$:

let $s \leftarrow K^{-1}r_i$

let $t \leftarrow AK^{-1}r_i$

for $j \le i$:

let $\gamma_j$ be the coefficient so that $t - \gamma_j r_j \perp r_j$

for $j \le i$:

form $s \leftarrow s - \gamma_j x_j$

and $t \leftarrow t - \gamma_j r_j$

let $x_{i+1} = (\sum_j \gamma_j)^{-1}s, r_{i+1} = (\sum_j \gamma_j)^{-1}t$.

You may recognize the Gram-Schmidt orthogonalization in this (see appendix 12.2 for an explanation): in each iteration $r_{i+1}$ is initially set to $AK^{-1}r_i$, and orthogonalized against $r_j$ with $j \le i$.

We can use *modified Gram-Schmidt* by rewriting the algorithm as:
Let $r_0$ be given

For $i \ge 0$:

let $s \leftarrow K^{-1}r_i$

let $t \leftarrow AK^{-1}r_i$

for $j \le i$:

let $\gamma_j$ be the coefficient so that $t - \gamma_j r_j \perp r_j$

form $s \leftarrow s - \gamma_j x_j$

and $t \leftarrow t - \gamma_j r_j$

let $x_{i+1} = (\sum_j \gamma_j)^{-1}s, r_{i+1} = (\sum_j \gamma_j)^{-1}t$.

These two versions of the FOM algorithm are equivalent in exact arithmetic, but differ in practical circumstances in two ways:

* The modified Gram-Schmidt method is more numerically stable;

* The unmodified method allows you to compute all inner products simultaneously. We discuss this below in section 6.6.

Even though the FOM algorithm is not used in practice, these computational considerations carry over to the GMRES method below.

### 5.5.9 Coupled recurrences form of iterative methods

Above, you saw the general equation (5.21) for generating iterates and search directions. This equation is often split as

* An update of the $x_i$ iterate from a single search direction:

$$x_{i+1} = x_i - \delta_i p_i,$$

and

* A construction of the search direction from the residuals known so far:

$$p_i = K^{-1}r_i + \sum_{j<i} \beta_{ij} K^{-1}r_j.$$

It is not hard to see inductively that we can also define

$$p_i = K^{-1}r_i + \sum_{j<i} \gamma_{ji} p_j,$$

and this last form is the one that is used in practice.

The iteration dependent coefficients are typically chosen to let the residuals satisfy various orthogonality conditions. For instance, one can choose to let the method be defined by letting the residuals be orthogonal $(r_i^t r_j = 0 \text{ if } i \neq j)$, or A-orthogonal $(r_i^t A r_j = 0 \text{ if } i \neq j)$. Many more schemes exist. Such methods can converge much faster than stationary iteration, or converge for a wider range of matrix and preconditioner types. Below we will see two such methods; their analysis, however, is beyond the scope of this course.
5.5.10 The method of Conjugate Gradients

In this section, we will derive the Conjugate Gradients (CG) method, which is a specific implementation of the FOM algorithm. In particular, it has pleasant computational properties in the case of an SPD matrix A.

The CG method takes as its basic form the coupled recurrences formulation described above, and the coefficients are defined by demanding that the sequence of residuals $r_0, r_1, r_2, \dots$ satisfy

$$r_i^t K^{-1} r_j = 0 \quad \text{if } i \neq j.$$

We start by deriving the CG method for nonsymmetric systems, and then show how it simplifies in the symmetric case. (The approach here is taken from [51]).

The basic equations are

$$\begin{cases}
x_{i+1} = x_i - \delta_i p_i \\
r_{i+1} = r_i - \delta_i A p_i \\
p_{i+1} = K^{-1} r_{i+1} + \sum_{j \le i} \gamma_{ji+1} p_j,
\end{cases} \qquad (5.22)$$

where the first and third equation were introduced above, and the second can be found by multiplying the first by A (check this!).

We will now derive the coefficients in this method by induction. In essence, we assume that we have current residual $r_{cur}$, a residuals to be computed $r_{new}$, and a collection of known residuals $R_{old}$. Rather than using subscripts ‘old, cur, new’, we use the following convention:

* $x_1, r_1, p_1$ are the current iterate, residual, and search direction. Note that the subscript 1 does not denote the iteration number here.

* $x_2, r_2, p_2$ are the iterate, residual, and search direction that we are about to compute. Again, the subscript does not equal the iteration number.

* $X_0, R_0, P_0$ are all previous iterates, residuals, and search directions bundled together in a block of vectors.

In terms of these quantities, the update equations are then

$$\begin{cases}
x_2 = x_1 - \delta_1 p_1 \\
r_2 = r_1 - \delta_1 A p_1 \\
p_2 = K^{-1} r_2 + v_{12} p_1 + P_0 u_{02}
\end{cases} \qquad (5.23)$$

where $\delta_1, v_{12}$ are scalars, and $u_{02}$ is a vector with length the number of iterations before the current. We now derive $\delta_1, v_{12}, u_{02}$ from the orthogonality of the residuals. To be specific, the residuals have to be orthogonal under the $K^{-1}$ inner product: we want to have

$$r_2^t K^{-1} r_1 = 0, \quad r_2^t K^{-1} R_0 = 0.$$

Combining these relations gives us, for instance,

$$\left. \begin{array}{l} r_1^t K^{-1} r_2 = 0 \\ r_2 = r_1 - \delta_i A K^{-1} p_1 \end{array} \right\} \Rightarrow \delta_1 = \frac{r_1^t K^{-1} r_1}{r_1^t K^{-1} A p_1}.$$
Finding $v_{12}, u_{02}$ is a little harder. For this, we start by summarizing the relations for the residuals and search directions in equation (5.22) in block form as

$$ (R_0, r_1, r_2) \begin{pmatrix} 1 & & & \\ & -1 & 1 & \\ & & \ddots & \ddots \\ & & & -1 \end{pmatrix} \begin{pmatrix} & & & \\ & 1 & & \\ & & \ddots & \ddots \\ & & & -1 \end{pmatrix} = A(P_0, p_1, p_2) \operatorname{diag}(D_0, d_1, d_2) $$

$$ (P_0, p_1, p_2) \begin{pmatrix} I - U_{00} & -u_{01} & -u_{02} \\ & 1 & -v_{12} \\ & & 1 \end{pmatrix} = K^{-1} (R_0, r_1, r_2) $$

or abbreviated $RJ = APD$, $P(I - U) = R$ where $J$ is the matrix with identity diagonal and minus identity subdiagonal. We then observe that

* $R^t K^{-t} R$ is diagonal, expressing the orthogonality of the residuals.

* Combining that $R^t K^{-t} R$ is diagonal and $P(I - U) = R$ gives that $R^t P = R^t K^{-t} R(I - U)^{-1}$. We now reason that $(I - U)^{-1}$ is upper diagonal, so $R^t P$ is upper triangular. This tells us quantities such as $r_2^t p_1$ are zero.

* Combining the relations for $R$ and $P$, we get first that

$$ R^t K^{-t} A P = R^t K^{-t} R J D^{-1} $$

which tells us that $R^t K^{-t} A P$ is lower bidiagonal. Expanding $R$ in this equation gives

$$ P^t A P = (I - U)^{-t} R^t R J D^{-1}. $$

Here $D$ and $R^t K^{-t} R$ are diagonal, and $(I - U)^{-t}$ and $J$ are lower triangular, so $P^t A P$ is lower triangular.

* This tells us that $P_0^t A p_2 = 0$ and $p_1^t A p_2 = 0$.

* Taking the product of $P_0^t A$, $p_1^t A$ with the definition of $p_2$ in equation (5.23) gives

$$ u_{02} = -(P_0^t A P_0)^{-1} P_0^t A K^{-1} r_2, \quad v_{12} = -(p_1^t A p_1)^{-1} p_1^t A K^{-1} r_2. $$

* If $A$ is symmetric, $P^t A P$ is lower triangular (see above) and symmetric, so it is in fact diagonal. Also, $R^t K^{-t} A P$ is lower bidiagonal, so, using $A = A^t$, $P^t A K^{-1} R$ is upper bidiagonal. Since $P^t A K^{-1} R = P^t A P(I - U)$, we conclude that $I - U$ is upper bidiagonal, so, only in the symmetric case, $u_{02} = 0$.

Some observations about this derivation.

* Strictly speaking we are only proving necessary relations here. It can be shown that these are sufficient too.

* There are different formulas that wind up computing the same vectors, in exact arithmetic. For instance, it is easy to derive that $p_1^t r_1 = r_1^t r_1$, so this can be substituted in the formulas just derived. The implementation of the CG method as it is typically implemented, is given in figure 5.8.
Compute $r^{(0)} = Ax^{(0)} - b$ for some initial guess $x^{(0)}$

for $i = 1, 2, \dots$

solve $Kz^{(i-1)} = r^{(i-1)}$

$\rho_{i-1} = r^{(i-1)T} z^{(i-1)}$

if $i = 1$

$p^{(1)} = z^{(0)}$

else

$\beta_{i-1} = \rho_{i-1}/\rho_{i-2}$

$p^{(i)} = z^{(i-1)} + \beta_{i-1}p^{(i-1)}$

endif

$q^{(i)} = Ap^{(i)}$

$\delta_i = \rho_{i-1}/p^{(i)T} q^{(i)}$

$x^{(i)} = x^{(i-1)} - \delta_i p^{(i)}$

$r^{(i)} = r^{(i-1)} - \delta_i q^{(i)}$

check convergence; continue if necessary

end

<caption>Figure 5.8: The Preconditioned Conjugate Gradient Method</caption>

* In the *k*-th iteration, computing $P_0^t A r_2$ (which is needed for $u_{02}$) takes *k* inner products. First of all, inner products are disadvantageous in a parallel context. Secondly, this requires us to store all search directions indefinitely. This second point implies that both work and storage go up with the number of iterations. Contrast this with the stationary iteration scheme, where storage was limited to the matrix and a few vectors, and work in each iteration was the same.

* The objections just raised disappear in the symmetric case. Since $u_{02}$ is zero, the dependence on $P_0$ disappears, and only the dependence on $p_1$ remains. Thus, storage is constant, and the amount of work per iteration is constant. The number of inner products per iteration can be shown to be just two.

**Exercise 5.51.** Do a flop count of the various operations in one iteration of the CG method. Assume that $A$ is the matrix of a five-point stencil and that the preconditioner $M$ is an incomplete factorization of $A$ (section 5.5.5.1). Let $N$ be the matrix size.

### 5.5.11 Derivation from minimization

The above derivation of the CG method is not often found in the literature. The typical derivation starts with a minimization problem with a *symmetric positive definite (SPD)* matrix $A$:

$$ \text{For which vector } x \text{ with } \|x\| = 1 \text{ is } f(x) = \frac{1}{2} x^t A x - b^t x \text{ minimal?} \quad (5.24) $$

If we accept the fact that the function $f$ has a minimum, which follows from the positive definiteness, we find the minimum by computing the derivative

$$ f'(x) = Ax - b. $$

and asking where $f'(x) = 0$. And, presto, there we have the original linear system.
**Exercise 5.52.** Derive the derivative formula above. (Hint: write out the definition of derivative as $\lim_{h\downarrow 0} \dots$.) Note that this requires A to be symmetric.

For the derivation of the iterative method, we state that the iterate $x_i$ is updated with a certain stepsize $\delta_i$ along a search direction $p_i$:

$$x_{i+1} = x_i + p_i \delta_i$$

The optimal stepsize

$$\delta_i = \frac{r_i^t p_i}{p_1^t A p_i}$$

is then derived as the one that minimizes the function $f$ along the line $x_i + \delta \delta p_i$:

$$\delta_i = \underset{\delta}{\operatorname{argmin}} \| f(x_i + p_i \delta) \|$$

The construction of the search direction from the residuals follows by induction proof from the requirement that the residuals be orthogonal. For a typical proof, see [3].

### 5.5.12 GMRES

In the discussion of the CG method above, it was pointed out that orthogonality of the residuals requires storage of all residuals, and $k$ inner products in the $k$'th iteration. Unfortunately, it can be proved that the work savings of the CG method can, for all practical purposes, not be found outside of SPD matrices [53].

The GMRES method is a popular implementation of such full orthogonalization schemes. In order to keep the computational costs within bounds, it is usually implemented as a restarted method. That is, only a certain number (say $k = 5$ or 20) of residuals is retained, and every $k$ iterations the method is restarted.

Other methods exist that do not have the growing storage demands of GMRES, for instance QMR [60] and BiCGstab [158]. Even though by the remark above these can not orthogonalize the residuals, they are still attractive in practice.

### 5.5.13 Complexity

The efficiency of Gaussian elimination was fairly easy to assess: factoring and solving a system takes, deterministically, $\frac{1}{3}n^3$ operations. For an iterative method, the operation count is the product of the number of operations per iteration times the number of iterations. While each individual iteration is easy to analyze, there is no good theory to predict the number of iterations. (In fact, an iterative method may not even converge to begin with.) Added to this is the fact that Gaussian elimination can be coded in such a way that there is considerable cache reuse, making the algorithm run at a fair percentage of the computer's peak speed. Iterative methods, on the other hand, are much slower on a flops per second basis.

All these considerations make the application of iterative methods to linear system solving somewhere in between a craft and a black art. In practice, people do considerable experimentation to decide whether an iterative method will pay off, and if so, which method is preferable.
**5.6 Further Reading**

Iterative methods is a very deep field. As a practical introduction to the issues involved, you can read the ‘Templates book’ [6], online at http://netlib.org/templates/. For a deeper treatment of the theory, see the book by Saad [142] of which the first edition can be downloaded at http://www-users.cs.umn.edu/~saad/books.html.
# Chapter 6

## High performance linear algebra

In this section we will discuss a number of issues pertaining to linear algebra on parallel computers. We will take a realistic view of this topic, assuming that the number of processors is finite, and that the problem data is always large, relative to the number of processors. We will also pay attention to the physical aspects of the communication network between the processors.

We will analyze various linear algebra operations, including iterative methods, and their behaviour in the presence of a network with finite bandwidth and finite connectivity. This chapter will conclude with various short remarks regarding complications in algorithms that arise due to parallel execution.

### 6.1 Collective operations

Collective operations play an important part in linear algebra operations. In fact, the scalability of the operations can depend on the cost of these collectives as you will see below. Here we give a short discussion of the essential ideas; see [22] for details.

In computing the cost of a collective operation, three architectural constants are enough to give lower bounds: $\alpha$, the of sending a single message, $\beta$, the inverse of the for sending data (see section 1.3.2), and $\gamma$, the inverse of the , the time for performing an arithmetic operation. Sending $n$ data items then takes time $\alpha + \beta n$. We further assume that a processor can only send one message at a time. We make no assumptions about the connectivity of the processors; thus, the lower bounds derived here will hold for a wide range of architectures.

The main implication of the architectural model above is that the number of active processors can only double in each step of an algorithm. For instance, to do a broadcast, first processor 0 sends to 1, then 0 and 1 can send to 2 and 3, then 0–3 send to 4–7, et cetera. This cascade of messages is called a *minimum spanning tree* of the processor network, and it follows that any collective algorithm has at least $\alpha \log_2 p$ cost associated with the accumulated latencies.

#### 6.1.1 Broadcast

In a *broadcast* operation, a single processor has $n$ data elements that is needs to send to all others: the other processors need a full copy of all $n$ elements. By the above doubling argument, we conclude that
a broadcast to *p* processors takes time at least ⌈log₂*p*⌋ steps with a total latency of ⌈log₂*p*⌋α. Since *n* elements are sent, this adds at least a time *n*β for all elements to leave the sending processor, giving a total cost lower bound of

$$ \lceil \log_2 p \rceil \alpha + n\beta. $$

We can illustrate the spanning tree method as follows:

<table><thead><tr><th></th><th>t = 1</th><th>t = 2</th><th>t = 3</th></tr></thead><tbody><tr><th>p<sub>0</sub></th><td>x<sub>0</sub> ↓, x<sub>1</sub> ↓, x<sub>2</sub> ↓, x<sub>3</sub> ↓</td><td>x<sub>0</sub> ↓, x<sub>1</sub> ↓, x<sub>2</sub> ↓, x<sub>3</sub> ↓</td><td>x<sub>0</sub>, x<sub>1</sub>, x<sub>2</sub>, x<sub>3</sub></td></tr><tr><th>p<sub>1</sub></th><td></td><td>x<sub>0</sub> ↓, x<sub>1</sub> ↓, x<sub>2</sub> ↓, x<sub>3</sub> ↓</td><td>x<sub>0</sub>, x<sub>1</sub>, x<sub>2</sub>, x<sub>3</sub></td></tr><tr><th>p<sub>2</sub></th><td></td><td></td><td>x<sub>0</sub>, x<sub>1</sub>, x<sub>2</sub>, x<sub>3</sub></td></tr><tr><th>p<sub>3</sub></th><td></td><td></td><td>x<sub>0</sub>, x<sub>1</sub>, x<sub>2</sub>, x<sub>3</sub></td></tr></tbody></table>

(On *t* = 1, *p*₀ sends to *p*₁; on *t* = 2 *p*₀, *p*₁ send to *p*₂, *p*₃.) This algorithm has the correct log₂*p* · α term, but processor 0 repeatedly sends the whole vector, so the bandwidth cost is log₂*p* · nβ. If *n* is small, the latency cost dominates, so we may characterize this as a short vector collective operation

The following algorithm implements the broadcast as a combination of a scatter and a bucket brigade algorithm:

<table><thead><tr><th></th><th>t = 0</th><th>t = 1</th><th>etcetera</th></tr></thead><tbody><tr><th>p<sub>0</sub></th><td>x<sub>0</sub> ↓</td><td>x<sub>0</sub></td><td>x<sub>0</sub>, x<sub>2</sub>, x<sub>3</sub></td></tr><tr><th>p<sub>1</sub></th><td>x<sub>1</sub> ↓</td><td>x<sub>0</sub> ↓, x<sub>1</sub></td><td>x<sub>0</sub>, x<sub>1</sub>, x<sub>3</sub></td></tr><tr><th>p<sub>2</sub></th><td>x<sub>2</sub> ↓</td><td>x<sub>1</sub> ↓, x<sub>2</sub></td><td>x<sub>0</sub>, x<sub>1</sub>, x<sub>2</sub></td></tr><tr><th>p<sub>3</sub></th><td>x<sub>3</sub> ↓</td><td>x<sub>2</sub> ↓, x<sub>3</sub></td><td>x<sub>1</sub>, x<sub>2</sub>, x<sub>3</sub></td></tr></tbody></table>

The complexity now becomes

$$ p\alpha + \beta n(p - 1)/p $$

which is not optimal in latency, but is a better algorithm if *n* is large, making this a long vector collective operation.

### 6.1.2 Reduction

In a reduction operation, each processor has *n* data elements, and one processor needs to combine them elementwise, for instance computing *n* sums or products.

By running the broadcast backwards in time, we see that a reduction operation has the same lower bound on the communication of ⌈log₂*p*⌋α + nβ. A reduction operation also involves computation, which would take a total time of (*p* − 1)γ*n* sequentially: each of *n* items gets reduced over *p* processors. Since these operations can potentially be parallelized, the lower bound on the computation is *p*−1/pγ*n*, giving a total of

$$ \lceil \log_2 p \rceil \alpha + n\beta + \frac{p-1}{p}\gamma n. $$
We illustrate the spanning tree algorithm, using the notation $x_i^{(j)}$ for the data item $i$ that was originally on processor $j$, and $x_i^{(j:k)}$ for the sum of the items $i$ of processors $j...k$.

<table><thead><tr><th></th><th>t = 1</th><th>t = 2</th><th>t = 3</th></tr></thead><tbody><tr><td>p<sub>0</sub></td><td>x<sub>0</sub><sup>(0)</sup>, x<sub>1</sub><sup>(0)</sup>, x<sub>2</sub><sup>(0)</sup>, x<sub>3</sub><sup>(0)</sup></td><td>x<sub>0</sub><sup>(0:1)</sup>, x<sub>1</sub><sup>(0:1)</sup>, x<sub>2</sub><sup>(0:1)</sup>, x<sub>3</sub><sup>(0:1)</sup></td><td>x<sub>0</sub><sup>(0:3)</sup>, x<sub>1</sub><sup>(0:3)</sup>, x<sub>2</sub><sup>(0:3)</sup>, x<sub>3</sub><sup>(0:3)</sup></td></tr><tr><td>p<sub>1</sub></td><td>x<sub>0</sub><sup>(1)</sup> ↑, x<sub>1</sub><sup>(1)</sup> ↑, x<sub>2</sub><sup>(1)</sup> ↑, x<sub>3</sub><sup>(1)</sup> ↑</td><td></td><td></td></tr><tr><td>p<sub>2</sub></td><td>x<sub>0</sub><sup>(2)</sup>, x<sub>1</sub><sup>(2)</sup>, x<sub>2</sub><sup>(2)</sup>, x<sub>3</sub><sup>(2)</sup></td><td>x<sub>0</sub><sup>(2:3)</sup> ↑, x<sub>1</sub><sup>(2:3)</sup> ↑, x<sub>2</sub><sup>(2:3)</sup> ↑, x<sub>3</sub><sup>(2:3)</sup> ↑</td><td></td></tr><tr><td>p<sub>3</sub></td><td>x<sub>0</sub><sup>(3)</sup> ↑, x<sub>1</sub><sup>(3)</sup> ↑, x<sub>2</sub><sup>(3)</sup> ↑, x<sub>3</sub><sup>(3)</sup> ↑</td><td></td><td></td></tr></tbody></table>

On time $t = 1$ processors $p_0, p_2$ receive from $p_1, p_3$, and on $t = 2$ $p_0$ receives from $p_2$.

As with the broadcast above, this algorithm does not achieve the lower bound; instead it has a complexity

$$ \lceil \log_2 p \rceil (\alpha + n\beta + \frac{p-1}{p}\gamma n). $$

For short vectors the $\alpha$ term dominates, so this algorithm is sufficient. For long vectors one can, as above, use other algorithms [22].

### 6.1.3 Allreduce

An allreduce operation computes the same elementwise reduction of $n$ elements on each processor, but leaves the result on each processor, rather than just on the root of the spanning tree. This could be implemented as a reduction followed by a broadcast, but more clever algorithms exist.

The lower bound on the cost of an allreduce is, somewhat remarkably, almost the same as of a simple reduction: since in a reduction not all processors are active at the same time, we assume that the extra work can be spread out perfectly. This means that the lower bound on the latency and computation stays the same. For the bandwidth we reason as follows: in order for the communication to be perfectly parallelized, $\frac{p-1}{p}n$ items have to arrive at, and leave each processor. Thus we have a total time of

$$ \lceil \log_2 p \rceil \alpha + 2 \frac{p-1}{p} n\beta + \frac{p-1}{p} n\gamma. $$

### 6.1.4 Allgather

In a gather operation on $n$ elements, each processor has $n/p$ elements, and one processor collects them all, without combining them as in a reduction. The allgather computes the same gather, but leaves the result on all processors.

Again we assume that gathers with multiple targets are active simultaneously. Since every processor originates a minimum spanning tree, we have $\log_2 p\alpha$ latency; since each processor receives $n/p$ elements from $p-1$ processors, there is $(p-1) \times (n/p)\beta$ bandwidth cost. The total cost for constructing a length $n$ vector by allgather is then

$$ \lceil \log_2 p \rceil \alpha + \frac{p-1}{p} n\beta. $$
We illustrate this:

<table><thead><tr><th></th><th>t = 1</th><th>t = 2</th><th>t = 3</th></tr></thead><tbody><tr><th>p<sub>0</sub></th><td>x<sub>0</sub> ↓</td><td>x<sub>0</sub>x<sub>1</sub> ↓</td><td>x<sub>0</sub>x<sub>1</sub>x<sub>2</sub>x<sub>3</sub></td></tr><tr><th>p<sub>1</sub></th><td>x<sub>1</sub> ↑</td><td>x<sub>0</sub>x<sub>1</sub> ↓</td><td>x<sub>0</sub>x<sub>1</sub>x<sub>2</sub>x<sub>3</sub></td></tr><tr><th>p<sub>2</sub></th><td>x<sub>2</sub> ↓</td><td>x<sub>2</sub>x<sub>3</sub> ↑</td><td>x<sub>0</sub>x<sub>1</sub>x<sub>2</sub>x<sub>3</sub></td></tr><tr><th>p<sub>3</sub></th><td>x<sub>3</sub> ↑</td><td>x<sub>2</sub>x<sub>3</sub> ↑</td><td>x<sub>0</sub>x<sub>1</sub>x<sub>2</sub>x<sub>3</sub></td></tr></tbody></table>

At time $t = 1$, there is an exchange between neighbours $p_0, p_1$ and likewise $p_2, p_3$; at $t = 2$ there is an exchange over distance two between $p_0, p_2$ and likewise $p_1, p_3$.

### 6.1.5 Reduce-scatter

In a reduce-scatter operation, each processor has $n$ elements, and an $n$-way reduction is done on them. Unlike in the reduce or allreduce, the result is then broken up, and distributed as in a scatter operation.

Formally, processor $i$ has an item $x_i^{(i)}$, and it needs $\sum_j x_i^{(j)}$. We could implement this by doing a size $p$ reduction, collecting the vector $(\sum_i x_0^{(i)}, \sum_i x_1^{(i)}, \dots)$ on one processor, and scattering the results. However it is possible to combine these operations in a so-called bidirectional exchange algorithm:

<table><thead><tr><th></th><th>t = 1</th><th>t = 2</th><th>t = 3</th></tr></thead><tbody><tr><th>p<sub>0</sub></th><td>x<sub>0</sub><sup>(0)</sup>, x<sub>1</sub><sup>(0)</sup>, x<sub>2</sub><sup>(0)</sup> ↓, x<sub>3</sub><sup>(0)</sup> ↓</td><td>x<sub>0</sub><sup>(0:2:2)</sup>, x<sub>1</sub><sup>(0:2:2)</sup> ↓</td><td>x<sub>0</sub><sup>(0:3)</sup></td></tr><tr><th>p<sub>1</sub></th><td>x<sub>0</sub><sup>(1)</sup>, x<sub>1</sub><sup>(1)</sup>, x<sub>2</sub><sup>(1)</sup> ↓, x<sub>3</sub><sup>(1)</sup> ↓</td><td>x<sub>0</sub><sup>(1:3:2)</sup> ↑, x<sub>1</sub><sup>(1:3:2)</sup></td><td>x<sub>1</sub><sup>(0:3)</sup></td></tr><tr><th>p<sub>2</sub></th><td>x<sub>0</sub><sup>(2)</sup> ↑, x<sub>1</sub><sup>(2)</sup> ↑, x<sub>2</sub><sup>(2)</sup>, x<sub>3</sub><sup>(2)</sup></td><td>x<sub>2</sub><sup>(0:2:2)</sup>, x<sub>3</sub><sup>(0:2:2)</sup> ↓</td><td>x<sub>2</sub><sup>(0:3)</sup></td></tr><tr><th>p<sub>3</sub></th><td>x<sub>0</sub><sup>(3)</sup> ↑, x<sub>1</sub><sup>(3)</sup> ↑, x<sub>2</sub><sup>(3)</sup>, x<sub>3</sub><sup>(3)</sup></td><td>x<sub>0</sub><sup>(1:3:2)</sup> ↑, x<sub>1</sub><sup>(1:3:2)</sup></td><td>x<sub>3</sub><sup>(0:3)</sup></td></tr></tbody></table>

The reduce-scatter can be considered as a allgather run in reverse, with arithmetic added, so the cost is

$$[\log_2 p]\alpha + \frac{p-1}{p}n(\beta + \gamma).$$

## 6.2 Parallel dense matrix-vector product

In this section we will go into great detail into the performance, and in particular the scalability, of the parallel dense matrix-vector product. First we will consider a simple case, and discuss the parallelism aspects in some amount of detail.
### 6.2.1 Implementing the block-row case

In designing a parallel version of an algorithm, one often proceeds by making a data decomposition of the objects involved. In the case of a matrix-vector operations such as the product $y = Ax$, we have the choice of starting with a vector decomposition, and exploring its ramifications on how the matrix can be decomposed, or rather to start with the matrix, and deriving the vector decomposition from it. In this case, it seems natural to start with decomposing the matrix rather than the vector, since it will be most likely of larger computational significance. We now have two choices:

1. We make a one-dimensional decomposition of the matrix, splitting it in block rows or block columns, and assigning each of these – or groups of them – to a processor.

2. Alternatively, we can make a two-dimensional decomposition, assigning to each processor one or more general submatrices.

We start by considering the decomposition in block rows. Consider a processor $p$ and the set $I_p$ of indices of rows that it owns<fnref n="1" />, and let $i \in I_p$ be a row that is assigned to this processor. The elements in row $i$ are used in the operation

$$y_i = \sum_j a_{ij}x_j$$

We now reason:

* If processor $p$ has all $x_j$ values, the matrix-vector product can trivially be executed, and upon completion, the processor has the correct values $y_j$ for $j \in I_p$.

* This means that every processor needs to have a copy of $x$, which is wasteful. Also it raises the question of data integrity: you need to make sure that each processor has the correct value of $x$.

* In certain practical applications (for instance iterative methods, as you have seen before), the output of the matrix-vector product is, directly or indirectly, the input for a next matrix-vector operation. This is certainly the case for the power method which computes $x, Ax, A^2x, \dots$. Since our operation started with each processor having the whole of $x$, but ended with it owning only the local part of $Ax$, we have a mismatch.

* Maybe it is better to assume that each processor, at the start of the operation, has only the local part of $x$, that is, those $x_i$ where $i \in I_p$, so that the start state and end state of the algorithm are the same. This means we have to change the algorithm to include some communication that allows each processor to obtain those values $x_i$ where $i \notin I_p$.

**Exercise 6.1.** Go through a similar reasoning for the case where the matrix is decomposed in block columns. Describe the parallel algorithm in detail, like above, but without giving pseudo code.

Let us now look at the communication in detail: we will consider a fixed processor $p$ and consider the operations it performs and the communication that necessitates. According to the above analysis, in executing the statement $y_i = \sum_j a_{ij}x_j$ we have to be aware what processor the $j$ values ‘belong to’. To acknowledge this, we write

$$y_i = \sum_{j \in I_p} a_{ij}x_j + \sum_{j \notin I_p} a_{ij}x_j \qquad (6.1)$$

<footnote n="1">1. For ease of exposition we will let $I_p$ be a contiguous range of indices, but any general subset is allowed.</footnote>
**Input:** Processor number *p*; the elements *x*<sub>*i*</sub> with *i* ∈ *I*<sub>*p*</sub>; matrix elements *A*<sub>*ij*</sub> with *i* ∈ *I*<sub>*p*</sub>.

**Output:** The elements $y_i$ with $i \in I_p$

for $i \in I_p$ do

$s \leftarrow 0;$

for $j \in I_p$ do

$s \leftarrow s + a_{ij}x_j$

for $j \notin I_p$ do

send $x_j$ from the processor that owns it to the current one, then;

$s \leftarrow s + a_{ij}x_j$

$y_i \leftarrow s$

**Procedure Naive Parallel MVP(A, x<sub>local</sub>, y<sub>local</sub>, p)**

<caption>Figure 6.1: A naively coded parallel matrix-vector product</caption>

If $j \in I_p$, the instruction $y_i \leftarrow y_i + a_{ij}x_j$ involves only quantities that are already local to the processor. Let us therefore concentrate on the case $j \notin I_p$. It would be nice if we could just write the statement

$y(i) = y(i) + a(i, j) * x(j)$

and some lower layer would automatically transfer $x(j)$ from whatever processor it is stored on to a local register. (The PGAS languages of section 2.7.5 aim to do this, but their efficiency is far from guaranteed.) An implementation, based on this optimistic view of parallelism, is given in figure 6.1.

The immediate problem with such a 'local' approach is that too much communication will take place.

* If the matrix $A$ is dense, the element $x_j$ is necessary once for each row $i \in I_p$, and it will thus be fetched once for every row $i \in I_p$.

* For each processor $q \neq p$, there will be (large) number of elements $x_j$ with $j \in I_q$ that need to be transferred from processor $q$ to $p$. Doing this in separate messages, rather than one bulk transfer, is very wasteful.

With shared memory these issues are not much of a problem, but in the context of distributed memory it is better to take a buffering approach.

Instead of communicating individual elements of $x$, we use a local buffer $B_{pq}$ for each processor $q \neq p$ where we collect the elements from $q$ that are needed to perform the product on $p$. (See figure 6.2 for an illustration.) The parallel algorithm is given in figure 6.3.

In addition to preventing an element from being fetched more than once, this also combines many small messages into one large message, which is usually more efficient; recall our discussion of bandwidth and latency in section 2.8.8.

**Exercise 6.2.** Give pseudocode for the matrix-vector product using nonblocking operations (section 2.7.3.6)
<figure><img src="image_87.png" /><figcaption>Figure 6.2: The parallel matrix-vector product with a blockrow distribution.</figcaption></figure>

**Input:** Processor number $p$; the elements $x_i$ with $i \in I_p$; matrix elements $A_{ij}$ with $i \in I_p$.

**Output:** The elements $y_i$ with $i \in I_p$

for $q \neq p$ do

Send elements of $x$ from processor $q$ to $p$, receive in buffer $B_{pq}$.

$y_{local} \leftarrow Ax_{local}$

for $q \neq p$ do

$y_{local} \leftarrow y_{local} + A_{pq}B_q$

**Procedure** Parallel MVP($A$, $x_{local}$, $y_{local}$, $p$)

<caption>Figure 6.3: A buffered implementation of the parallel matrix-vector product</caption>
Above we said that having a copy of the whole of $x$ on each processor was wasteful in space. The implicit argument here is that, in general, we do not want local storage to be function of the number of processors: ideally it should be only a function of the local data. (This is related to weak scaling; section 2.3.4.)

You see that, because of communication considerations, we have actually decided that it is unavoidable, or at least preferable, for each processor to store the whole input vector. Such trade-offs between space and time efficiency are fairly common in parallel programming. For the dense matrix-vector product we can actually defend this overhead, since the vector storage is of lower order than the matrix storage, so our over-allocation is small by ratio. Below (section 6.5), we will see that for the sparse matrix-vector product the overhead can be much less.

It is easy to see that the parallel dense matrix-vector product, as described above, has perfect speedup if we are allowed to ignore the time for communication. In the next couple of sections you will see that the block row implementation above is not optimal if we take communication into account. For scalability we need a two-dimensional decomposition. We start with a discussion of collectives.

## 6.2.2 Scalability of the dense matrix-vector product

In this section, we will give a full analysis of the parallel computation of $y \leftarrow Ax$, where $x, y \in \mathbb{R}^n$ and $A \in \mathbb{R}^{n \times n}$. We will assume that $p$ nodes will be used, but we make no assumptions on their connectivity. We will see that the way the matrix is distributed makes a big difference for the scaling of the algorithm; for the original research see [83, 145, 150], and see section 2.3.4 for the definitions of the various forms of scaling.

### 6.2.2.1 Matrix-vector product, partitioning by rows

Partition

$$A \rightarrow \begin{pmatrix} A_0 \\ A_1 \\ \vdots \\ A_{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_i \in \mathbb{R}^{m_i \times n}$ and $x_i, y_i \in \mathbb{R}^{m_i}$ with $\sum_{i=0}^{p-1} m_i = n$ and $m_i \approx n/p$. We will start by assuming that $A_i, x_i$, and $y_i$ are originally assigned to $\mathcal{P}_i$.

The computation is characterized by the fact that each processor needs the whole vector $x$, but owns only an $n/p$ fraction of it. Thus, we execute an allgather of $x$. After this, the processor can execute the local product $y_i \leftarrow A_i x$; no further communication is needed after that.

An algorithm with cost computation for $y = Ax$ in parallel is then given by

<table><thead><tr><th>Step</th><th>Cost (lower bound)</th></tr></thead><tbody><tr><td>Allgather x<sub>i</sub> so that x is available on all nodes</td><td>$$ \begin{aligned} & \lceil \log_2(p) \rceil \alpha + \frac{p-1}{p} n \beta \\ & \approx \log_2(p) \alpha + n \beta \end{aligned} $$</td></tr><tr><td>Locally compute y<sub>i</sub> = A<sub>i</sub>x</td><td>$$ \begin{aligned} & \approx 2 \frac{n^2}{p} \gamma \end{aligned} $$</td></tr></tbody></table>
**Cost analysis** The total cost of the algorithm is given by, approximately,

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

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

$$S_p^{\text{1D-row}}(n) = \frac{T_1(n)}{T_p^{\text{1D-row}}(n)} = \frac{2n^2\gamma}{2\frac{n^2}{p}\gamma + \log_2(p)\alpha + n\beta} = \frac{p}{1 + \frac{p\log_2(p)}{2n^2}\frac{\alpha}{\gamma} + \frac{p\beta}{2n}\gamma}$$

and the parallel efficiency by

$$E_p^{\text{1D-row}}(n) = \frac{S_p^{\text{1D-row}}(n)}{p} = \frac{1}{1 + \frac{p \log_2(p) \alpha}{2n^2 \gamma} + \frac{p \beta}{2n \gamma}}$$

**An optimist's view** Now, if one fixes $p$ and lets $n$ get large,

$$\lim_{n \to \infty} E_p(n) = \lim_{n \to \infty} \left[ \frac{1}{1 + \frac{p \log_2(p) \alpha}{2n^2 \gamma} + \frac{p \beta}{2n \gamma}} \right] = 1.$$

Thus, if one can make the problem large enough, eventually the parallel efficiency is nearly perfect. However, this assumes unlimited memory, so this analysis is not practical.

**A pessimist's view** In a strong scalability analysis, one fixes $n$ and lets $p$ get large, to get

$$\lim_{p \to \infty} E_p(n) = \lim_{p \to \infty} \left[ \frac{1}{1 + \frac{p \log_2(p) \alpha}{2n^2 \gamma} + \frac{p \beta}{2n \gamma}} \right] = 0.$$

Thus, eventually the parallel efficiency becomes nearly nonexistent.

**A realist's view** In a more realistic view we increase the number of processors with the amount of data. This is called weak scalability, and it makes the amount of memory that is available to store the problem scale linearly with $p$.

Let $M$ equal the number of floating point numbers that can be stored in a single node's memory. Then the aggregate memory is given by $Mp$. Let $n_{max}(p)$ equal the largest problem size that can be stored in the aggregate memory of $p$ nodes. Then, if all memory can be used for the matrix,

$$(n_{\text{max}}(p))^2 = Mp \quad \text{or} \quad n_{\text{max}}(p) = \sqrt{Mp}.$$
The question now becomes what the parallel efficiency for the largest problem that can be stored on p nodes:

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

Now, if one analyzes what happens when the number of nodes becomes large, one finds that

$$ \lim_{p \to \infty} E_p(n_{\text{max}}(p)) = \lim_{p \to \infty} \left[ \frac{1}{1 + \frac{\log_2(p)}{2M} \frac{\alpha}{\gamma} + \frac{\sqrt{p}}{2\sqrt{M}} \frac{\beta}{\gamma}} \right] = 0. $$

Thus, this parallel algorithm for matrix-vector multiplication does not scale.

If you take a close look at this expression for efficiency, you'll see that the main problem is the $1/\sqrt{p}$ part of the expression. This terms involves a factor $\beta$, and if you follow the derivation backward you see that it comes from the time to send data between the processors. Informally this can be described as saying that the message size is too large to make the problem scalable. In fact, the message size is constant $n$, regardless the number of processors.

Alternatively, a realist realizes that there is a limited amount of time, $T_{\text{max}}$, to get a computation done. Under the best of circumstances, that is, with zero communication overhead, the largest problem that we can solve in time $T_{\text{max}}$ is given by

$$ T_p(n_{\text{max}}(p)) = 2 \frac{(n_{\text{max}}(p))^2}{p} \gamma = T_{\text{max}}. $$

Thus

$$ (n_{\text{max}}(p))^2 = \frac{T_{\text{max}}p}{2\gamma} \quad \text{or} \quad n_{\text{max}}(p) = \frac{\sqrt{T_{\text{max}}} \sqrt{p}}{\sqrt{2\gamma}}. $$

Then the parallel efficiency that is attained by the algorithm for the largest problem that can be solved in time $T_{\text{max}}$ is given by

$$ E_{p,n_{\text{max}}} = \frac{1}{1 + \frac{\log_2 p}{T} \alpha + \sqrt{\frac{p \beta}{T \gamma}}} $$

and the parallel efficiency as the number of nodes becomes large approaches

$$ \lim_{p \to \infty} E_p = \sqrt{\frac{T\gamma}{p\beta}}. $$

Again, efficiency cannot be maintained as the number of processors increases and the execution time is capped.