## 11.2 Parallel Random Number Generation

Random numbers are often used in simulations as some examples below will show. True random numbers are very hard to obtain: they could be generated by measuring quantum processes such as radioactive particles, but this is very cumbersome. Instead, we use pseudo-random numbers. This means that we use a deterministic mathematical process, that is sufficiently irregular that for practical purposes no order can be found in it.

### 11.2.1 Sequential random number generators

An easy way to generate random numbers (we leave off the 'pseudo' qualification) is to use linear congruential generators (for all you ever need to know about random numbers, see Knuth [99]), recurrences of the form

$$x_{k+1} = (ax_k + b) \mod m.$$

This sequence is periodic, since it consists of nonnegative integers at most $m-1$, and with period $m$ under certain conditions. A typical period is $2^{31}$. The starting point $x_0$ of the series is known as the 'seed'. Software for random numbers often lets you specify the seed. To get reproducible results you would run your program with the same seed multiple times; to get random behaviour over multiple runs of your program you could for instance derive the seed from clock and calendar functions.

Linear congruential generators may have some amount of correlation between lower bits. A different principle of generating random numbers is the lagged Fibonacci random number generator

$$X_i = X_{i-p} \otimes X_{i-q}$$

where $p, q$ are the lag parameter, and $\otimes$ is any binary operation, such as addition or multiplication modulo $M$.

The main problems with lagged Fibonacci generators are:

* They require setting $max(p, q)$ initial values, and their randomness is sensitive to these choices;

* They theory is not as developed as for congruential generators, so their is a greater reliance on statistical tests to evaluate their 'randomness'.

### 11.2.2 Parallel random number generation

Random number generation is problematic in parallel. To see this, consider a parallel process that uses a random number generator on each subprocess, and consider a single processor emulating the parallel process. Now this single process in effect has a random number generator that consists of interleaving the parallel generator results. This means that, if we use the same generator in all parallel processes, the effective generator over the whole process will produce stretches of identical values.

There are various ways out.
### 11.2.2.1 Master-worker generator

We can generate the random numbers centrally. In shared memory that could mean making its operation atomic. This may introduce a serious bottleneck.

**Exercise 11.1.** Critical sections are usually justified if the work spent there is of lower order than the parallel work. Why does that argument not apply here.

Another solution would be to have one thread or process that generates the random numbers and distributes them to the other processes. Doing this on a number-by-number basis causes considerable overhead. Instead, it would be possible for the generator process to distribute blocks of numbers. However, the manner in which this is done may again cause correlation between processes.

### 11.2.2.2 Sequence splitting solutions

A better solution is to set up independent generators with parameter choices that guarantee statistical randomness. This is not simple. For instance, if two sequences $x_i^{(1)}, x_i^{(2)}$ have the same values of $a, b, m$, and their starting points are close together, the sequences will be strongly correlated. Less trivial examples of correlation exist.

Various techniques for random number generation exist, such as using two sequences, where one generates the starting points for the other sequence, which is the one actually used in the simulation. Software for parallel random number generator can be found at http://sprng.cs.fsu.edu/ [120].

### 11.2.2.3 Blocked random number generators

Some random number generators (see [112]) allow you to calculate a value that is many iteration away from the seed. You could then take the block of values from the seed to that iteration and give it to one processor. Similarly, each processor would get a contiguous block of iterations of the generator.

## 11.3 Examples

### 11.3.1 Monte Carlo simulation of the Ising model

The Ising model (for an introduction, see [28]) was originally proposed to model ferromagnetism. Magnetism is the result of atoms aligning their ‘spin’ direction: let’s say spin can only be ‘up’ or ‘down’, then a material has magnetism if more atoms have spin up than down, or the other way around. The atoms are said to be in a structure called a ‘lattice’.

Now image heating up a material, which loosens up the atoms. If an external field is applied to the material, the atoms will start aligning with the field, and if the field is removed the magnetism disappears again. However, below a certain critical temperature the material will retain its magnetism. We will use Monte Carlo simulation to find the stable configurations that remain.

Let’s say the lattice $\Lambda$ has $N$ atoms, and we denote a configuration of atoms as $\sigma = (\sigma_1, \dots, \sigma_N)$ where each $\sigma_i = \pm 1$. The energy of a lattice is modeled as

$$H = H(\sigma) = -J \sum_i \sigma_i - E \sum_{ij} \sigma_i \sigma_j.$$
The first term models the interaction of individual spins $\sigma_i$ with an external field of strength $J$. The second term, which sums over nearest neighbour pairs, models alignment of atom pairs: the product $\sigma_i\sigma_j$ is positive if the atoms have identical spin, and negative if opposite.

In statistical mechanics, the probability of a configuration is

$$P(\sigma) = \exp(-H(\sigma))/Z$$

where the 'partitioning function' $Z$ is defined as

$$Z = \sum_{\sigma} \exp(H(\sigma))$$

where the sum runs over all $2^N$ configurations.

A configuration is stable if its energy does not decrease under small perturbations. To explore this, we iterate over the lattice, exploring whether altering the spin of atoms lowers the energy. We introduce an element of chance to prevent artificial solutions. (This is the *Metropolis algorithm* [123].)

**for fixed number of iterations do**

**for each atom i do**

calculate the change $\Delta E$ from changing the sign of $\sigma_i$

if $\Delta E < \exp(-\Delta E)$ greater than some random number then

accept the change

This algorithm can be parallelized, if we notice the similarity with the structure of the sparse matrix-vector product. In that algorithm too we compute a local quantity by combining inputs from a few nearest neighbours. This means we can partitioning the lattice, and compute the local updates after each processor collects a *ghost region*.

Having each processor iterate over local points in the lattice corresponds to a particular global ordering of the lattice; to make the parallel computation equivalent to a sequential one we also need a parallel random generator (section 11.2).
# PART III

## APPENDICES
This course requires no great mathematical sophistication. Mostly it assumes that you know the basics of linear algebra: what are matrices and vectors, and the most common operations on them.

In the following appendices we will cover some less common bits of theory that have been moved out of the main storyline of the preceding chapters.
# Chapter 12

## Linear algebra

In this course it is assumed that you know what a matrix and a vector are, simple algorithms such as how to multiply them, and some properties such as invertibility of a matrix. This appendix introduces some concepts and theorems that are not typically part of a first course in linear algebra.

### 12.1 Norms

A norm is a way to generalize the concept of absolute value to multi-dimensional objects such as vectors and matrices. There are many ways of defining a norm, and there is theory of how different norms relate. Here we only give the basic definitions; for more detail see any linear algebra textbook, for instance [70].

#### 12.1.1 Vector norms

A norm is any function $n(\cdot)$ on a vector space $V$ with the following properties:

* $n(x) \ge 0$ for all $x \in V$ and $n(x) = 0$ only for $x = 0$,

* $n(\lambda x) = |\lambda|n(x)$ for all $x \in V$ and $\lambda \in \mathbb{R}$.

* $n(x + y) \le n(x) + n(y)$

For any $p \ge 1$, the following defines a vector norm:

$$|x|_p = \sqrt[p]{\sum_i |x_i|^p}.$$

Common norms are $\| \cdot \|_1$ ('sum of absolute values') and $\| \cdot \|_2$ ('square root of sum of squares'); the $\| \cdot \|_{\infty}$ norm is defined as $\lim_{p \to \infty} \| \cdot \|_p$, and it is not hard to see that this equals

$$\|x\|_{\infty} = \max_{i} |x_{i}|.$$
### 12.1.2 Matrix norms

By considering a matrix of size $n$ as a vector of length $n^2$, we can define the Frobenius matrix norm:

$$ \|A\|_F = \sqrt{\sum_{i,j} |a_{ij}|^2}. $$

However, we will mostly look at associated matrix norms:

$$ \|A\|_p = \sup_{\|x\|_p=1} \|Ax\|_p = \sup_x \frac{\|Ax\|_p}{\|x\|_p}. $$

From their definition, it easily follows that

$$ \|Ax\| \le \|A\| \|x\| $$

for associated norms.

The following are easy to derive:

* $\|A\|_1 = \max_j \sum_i |a_{ij}|$,

* $\|A\|_\infty = \max_i \sum_j |a_{ij}|$.

By observing that $\|A\|_2 = \sup_{\|x\|_2=1} x^t A^t A x$, it is not hard to derive that $\|A\|_2$ is the maximal singular value of $A$, which is the root of the maximal eigenvalue of $A^t A$.

The matrix condition number is defined as

$$ \kappa(A) = \|A\| \|A^{-1}\|. $$

In the symmetric case, and using the 2-norm, this is the ratio between the largest and smallest eigenvalue.

## 12.2 Gram-Schmidt orthogonalization

The Gram-Schmidt (GS) algorithm takes a series of vectors and inductively orthogonalizes them. This can be used to turn an arbitrary basis of a vector space into an orthogonal basis; it can also be viewed as transforming a matrix $A$ into one with orthogonal columns. If $Q$ has orthogonal columns, $Q^tQ$ is diagonal, which is often a convenient property to have.

The basic principle of the GS algorithm can be demonstrated with two vectors $u, v$. Suppose we want a vector $v'$ so that $u, v$ and $u, v'$ span the same space, but $v' \perp u$. For this we let

$$ v' \leftarrow v - \frac{u^t v}{u^t u} u. $$

It is easy to see that this satisfies the requirements.

Suppose we have a set of vectors $u_1, \dots, u_n$ that we want to orthogonalize. We do this by successive applications of the above transformation:
For $i = 1, \dots, n$:

For $j = 1, \dots, i-1$:

$$ \text{let } c_{ji} = u_j^t u_i / u_j^t u_j $$

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

$$ \text{update } u_i \leftarrow u_i - u_j c_{ji} $$

Often the vector $v$ in the algorithm above is normalized; this adds a line

$$ u_i \leftarrow u_i / \|u_i\| $$

to the algorithm. GS orthogonalization with this normalization, applied to a matrix, is also known as the QR factorization.

**Exercise 12.1.** Suppose that we apply the GS algorithm to the columns of a rectangular matrix $A$, giving a matrix $Q$. Prove that there is an upper triangular matrix $R$ such that $A = QR$. (Hint: look at the $c_{ji}$ coefficients above.) If we normalize the orthogonal vector in the algorithm above, $Q$ has the additional property that $Q^t Q = I$. Prove this too.

The GS algorithm as given above computes the desired result, but only in exact arithmetic. A computer implementation can be quite inaccurate if the angle between $v$ and one of the $u_i$ is small. In that case, the Modified Gram-Schmidt (MGS) algorithm will perform better:

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

For $j = 1, \dots, i-1$:

$$ \text{let } c_{ji} = u_j^t u_i / u_j^t u_j $$

$$ \text{update } u_i \leftarrow u_i - u_j c_{ji} $$

To contrast it with MGS, the original GS algorithm is also known as Classical Gram-Schmidt (CGS).

As an illustration of the difference between the two methods, consider the matrix

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

where $\epsilon$ is of the order of the machine precision, so that $1 + \epsilon^2 = 1$ in machine arithmetic. The CGS method proceeds as follows:

* The first column is of length 1 in machine arithmetic, so

$$ q_1 = a_1 = \begin{pmatrix} 1 \\ \epsilon \\ 0 \\ 0 \end{pmatrix}. $$

* The second column gets orthogonalized as $v \leftarrow a_2 - 1 \cdot q_1$, giving

$$ v = \begin{pmatrix} 0 \\ -\epsilon \\ \epsilon \\ 0 \end{pmatrix}, \quad \text{normalized:} \quad q_2 = \begin{pmatrix} 0 \\ -\frac{\sqrt{2}}{2} \\ \frac{\sqrt{2}}{2} \\ 0 \end{pmatrix} $$
* The third column gets orthogonalized as $v \leftarrow a_3 - c_1q_1 - c_2q_2$, where

$$ \begin{cases} c_1 = q_1^t a_3 = 1 \\ c_2 = q_2^t a_3 = 0 \end{cases} \quad \Rightarrow \quad v = \begin{pmatrix} 0 \\ -\epsilon \\ 0 \\ \epsilon \end{pmatrix}; \quad \text{normalized:} \quad q_3 = \begin{pmatrix} 0 \\ \frac{\sqrt{2}}{2} \\ 0 \\ \frac{\sqrt{2}}{2} \end{pmatrix} $$

It is easy to see that $q_2$ and $q_3$ are not orthogonal at all. By contrast, the MGS method differs in the last step:

* As before, $q_1^t a_3 = 1$, so

$$ v \leftarrow a_3 - q_1 = \begin{pmatrix} 0 \\ -\epsilon \\ \epsilon \\ 0 \end{pmatrix}. $$

Then, $q_2^t v = \frac{\sqrt{2}}{2}\epsilon$ (note that $q_2^t a_3 = 0$ before), so the second update gives

$$ v \leftarrow v - \frac{\sqrt{2}}{2}\epsilon q_2 = \begin{pmatrix} 0 \\ \frac{\epsilon}{2} \\ -\frac{\epsilon}{2} \\ \epsilon \end{pmatrix}, \quad \text{normalized:} \quad \begin{pmatrix} 0 \\ \frac{\sqrt{6}}{6} \\ -\frac{\sqrt{6}}{6} \\ 2\frac{\sqrt{6}}{6} \end{pmatrix} $$

Now all $q_i^t q_j$ are on the order of $\epsilon$ for $i \neq j$.

## 12.3 The power method

The vector sequence

$$ x_i = Ax_{i-1}, $$

where $x_0$ is some starting vector, is called the power method since it computes the product of subsequent matrix powers times a vector:

$$ x_i = A^i x_0. $$

There are cases where the relation between the $x_i$ vectors is simple. For instance, if $x_0$ is an eigenvector of A, we have for some scalar $\lambda$

$$ Ax_0 = \lambda x_0 \quad \text{and} \quad x_i = \lambda^i x_0. $$

However, for an arbitrary vector $x_0$, the sequence $\{x_i\}_i$ is likely to consist of independent vectors. Up to a point.
Exercise 12.2. Let $A$ and $x$ be the $n \times n$ matrix and dimension $n$ vector

$$A = \begin{pmatrix} 1 & 1 & & & \\ & 1 & 1 & & \\ & & \ddots & \ddots & \\ & & & 1 & 1 \\ & & & & 1 \end{pmatrix}, \quad x = (0, \dots, 0, 1)^t.$$

Show that the sequence $[x, Ax, \dots, A^ix]$ is an independent set for $i < n$. Why is this no longer true for $i \ge n$?

Now consider the matrix $B$:

$$B = \begin{pmatrix} 1 & 1 & & & \\ & \ddots & \ddots & & \\ & & 1 & 1 & \\ & & & 1 & \\ \hline & & & & 1 & 1 \\ & & \ddots & \ddots & & 1 & 1 \\ & & & & & 1 & 1 \end{pmatrix}, \quad y = (0, \dots, 0, 1)^t$$

Show that the set $[y, By, \dots, B^iy]$ is an independent set for $i < n/2$, but not for any larger values of $i$.

While in general the vectors $x, Ax, A^2x, \dots$ can be expected to be independent, in computer arithmetic this story is no longer so clear.

Suppose the matrix has eigenvalues $\lambda_0 > \lambda_1 \ge \dots \ge \lambda_{n-1}$ and corresponding eigenvectors $u_i$ so that

$$Au_i = \lambda_i u_i.$$

Let the vector $x$ be written as

$$x = c_0 u_0 + \dots + c_{n-1} u_{n-1},$$

then

$$A^i x = c_0 \lambda_0^i u_i + \dots + c_{n-1} \lambda_{n-1}^i u_{n-1}.$$

If we write this as

$$A^i x = \lambda_0^i \left[ c_0 u_i + c_1 \left( \frac{\lambda_1}{\lambda_0} \right)^i + \dots + c_{n-1} \left( \frac{\lambda_{n-1}}{\lambda_0} \right)^i \right],$$

we see that, numerically, $A^ix$ will get progressively closer to a multiple of $u_0$, the dominant eigenvector. Hence, any calculation that uses independence of the $A^ix$ vectors is likely to be inaccurate.
## 12.4 Nonnegative matrices; Perron vectors

If A is a nonnegative matrix, the maximal eigenvalue has the property that its eigenvector is nonnegative: this is the the Perron-Frobenius theorem.

**Theorem 4** If a nonnegative matrix A is irreducible, its eigenvalues satisfy

* The eigenvalue $\alpha_1$ that is largest in magnitude is real and simple:

$$ \alpha_1 > |\alpha_2| \ge \dots $$

* The corresponding eigenvector is positive.

## 12.5 The Gershgorin theorem

Finding the eigenvalues of a matrix is usually complicated. However, there are some tools to estimate eigenvalues. In this section you will see a theorem that, in some circumstances, can give useful information on eigenvalues.

Let A be a square matrix, and x, λ an eigenpair: $Ax = λx$. Looking at one component, we have

$$ a_{ii}x_i + \sum_{j \neq i} a_{ij}x_j = \lambda x_i $$

Taking norms:

$$ (a_{ii} - \lambda) \le \sum_{j \neq i} |a_{ij}| \frac{|x_j|}{|x_i|} $$

Taking the value of i for which $|x_i|$ is maximal, we find

$$ (a_{ii} - \lambda) \le \sum_{j \neq i} |a_{ij}|. $$

This statement can be interpreted as follows:

The eigenvalue λ is located in the circle around $a_{ii}$ with radius $\sum_{j \neq i} |a_{ij}|$.

Since we do not know for which value of i $|x_i|$ is maximal, we can only say that there is some value of i such that λ lies in such a circle. This is the Gershgorin theorem.

**Theorem 5** Let A be a square matrix, and let $D_i$ be the circle with center $a_{ii}$ and radius $\sum_{j \neq i} |a_{ij}|$, then the eigenvalues are contained in the union of circles $\cup_i D_i$.

We can conclude that the eigenvalues are in the interior of these discs, if the constant vector is not an eigenvector.
## 12.6 Householder reflectors

In some contexts the question comes up how to transform one subspace into another. Householder reflectors are in a sense the minimal solution to this. Consider a unit vector $u$, and let

$$H = I - 2uu^t.$$

You can easily check that $Hu = -u$, and if $u \perp v$, then $Hv = v$. In other words, the subspace of multiples of $u$ is flipped, and the orthogonal subspace stays invariant.

Now for the original problem of mapping one space into another. Let the original space be spanned by a vector $x$ and the resulting by $y$, then note that

$$\begin{cases} x = (x + y)/2 + (x - y)/2 \\ y = (x + y)/2 - (x - y)/2 \end{cases}$$

In other words, we can map $x$ into $y$ with the reflector based on $u = (x - y)/2$.

<figure><img src="image_128.png" /><figcaption>Figure 12.1: Householder reflector</figcaption></figure>

We can generalize Householder reflectors to a form

$$H = I - 2uv^t.$$

The matrices $L_i$ used in LU factorization (see section 5.3) can then be seen to be of the form $L_i = I - \ell_i e_i^t$ where $e_i$ has a single one in the $i$-th location, and $\ell_i$ only has nonzero below that location. That form also makes it easy to see that $L_i^{-1} = I + \ell_i e_i^t$:

$$(I - uv^t)(I + uv^t) = I - uv^tuv^t = 0$$

if $v^t u = 0$.
# Chapter 13

## Complexity

At various places in this book we are interested in how many operations an algorithm takes. It depends on the context what these operations are, but often we count additions (or subtractions) and multiplications. This is called the *arithmetic* or *computational complexity* of an algorithm. For instance, summing *n* numbers takes *n* − 1 additions. Another quantity that we may want to describe is the amount of space (computer memory) that is needed. Sometimes the space to fit the input and output of an algorithm is all that is needed, but some algorithms need temporary space. The total required space is called the *space complexity* of an algorithm.

Both arithmetic and space complexity depend on some description of the input, for instance, for summing an array of numbers, the length *n* of the array is all that is needed. We express this dependency by saying ‘summing an array of numbers has time complexity *n* − 1 additions, where *n* is the length of the array’.

The time (or space) the summing algorithm takes is not dependent on other factors such as the values of the numbers. By contrast, some algorithms such as computing the greatest common divisor of an array of integers *can* be dependent on the actual values.

**Exercise 13.1.** What is the time and space complexity of multiplying two square matrices of size *n* × *n*? Assume that an addition and a multiplication take the same amount of time.

Often we aim to simplify the formulas that describe time or space complexity. For instance, if the complexity of an algorithm is $n^2 + 2n$, we see that for $n > 2$ the complexity is less than $2n^2$, and for $n > 4$ it is less than $(3/2)n^2$. On the other hand, for all values of *n* the complexity is at least $n^2$. Clearly, the quadratic term $n^2$ is the most important, and the linear term *n* becomes less and less important by ratio. We express this informally by saying that the complexity is quadratic in *n* as $n \to \infty$: there are constants *c*, *C* so that for *n* large enough the complexity is at least $cn^2$ and at most $Cn^2$.

This is expressed for short by saying that the complexity is of order $n^2$, written as $O(n^2)$ as $n \to \infty$. In chapter 4 you will see phenomena that we want to describe as orders of a parameter that goes to zero. In that case we write for instance $f(h) = O(h^2)$ as $h \downarrow 0$, meaning that *f* is bounded by $ch^2$ and $Ch^2$ for certain constants *c*, *C* and *h* small enough.
# Chapter 14

## Partial Differential Equations

Partial Differential Equations are the source of a large fraction of HPC problems. Here is a quick derivation of two of the most important ones.

### 14.1 Partial derivatives

Derivatives of a function $u(x)$ are a measure of the rate of change. Partial derivatives to the same, but for a function $u(x, y)$ of two variables. Notated $u_x$ and $u_y$, these partial derivatives indicate the rate of change if only one variable changes and the other stays constant.

Formally, we define $u_x, u_y$ by:

$$u_x(x, y) = \lim_{h \to 0} \frac{u(x + h, y) - u(x, y)}{h}, \quad u_y(x, y) = \lim_{h \to 0} \frac{u(x, y + h) - u(x, y)}{h}$$

### 14.2 Poisson or Laplace Equation

Let $T$ be the temperature of a material, then its heat energy is proportional to it. A segment of length $\Delta x$ has heat energy $Q = c\Delta x \cdot u$. If the heat energy in that segment is constant

$$\frac{\delta Q}{\delta t} = c \Delta x \frac{\delta u}{\delta t} = 0$$

but it is also the difference between inflow and outflow of the segment. Since flow is proportional to temperature differences, that is, to $u_x$, we see that also

$$0 = \left. \frac{\delta u}{\delta x} \right|_{x+\Delta x} - \left. \frac{\delta u}{\delta x} \right|_{x}$$

In the limit of $\Delta x \downarrow 0$ this gives $u_{xx} = 0$, which is called the Laplace equation. If we have a source term, for instance corresponding to externally applied heat, the equation becomes $u_{xx} = f$, which is called the Poisson equation.
## 14.3 Heat Equation

Let $T$ be the temperature of a material, then its heat energy is proportional to it. A segment of length $\Delta x$ has heat energy $Q = c\Delta x \cdot u$. The rate of change in heat energy in that segment is

$$ \frac{\delta Q}{\delta t} = c\Delta x \frac{\delta u}{\delta t} $$

but it is also the difference between inflow and outflow of the segment. Since flow is proportional to temperature differences, that is, to $u_x$, we see that also

$$ \frac{\delta Q}{\delta t} = \left. \frac{\delta u}{\delta x} \right|_{x+\Delta x} - \left. \frac{\delta u}{\delta x} \right|_{x} $$

In the limit of $\Delta x \downarrow 0$ this gives $u_t = \alpha u_{xx}$.

## 14.4 Steady state

The solution of an IBVP is a function $u(x, t)$. In cases where the forcing function and the boundary conditions do not depend on time, the solution will converge in time, to a function called the *steady state* solution:

$$ \lim_{t \to \infty} u(x, t) = u_{\text{steady state}}(x). $$

This solution satisfies a BVP, which can be found by setting $u_t \equiv 0$. For instance, for the heat equation

$$ u_t = u_{xx} + q(x) $$

the steady state solution satisfies $-u_{xx} = q(x)$.
# Chapter 15

## Taylor series

Taylor series expansion is a powerful mathematical tool. In this course it is used several times in proving properties of numerical methods.

The Taylor expansion theorem, in a sense, asks how well functions can be approximated by polynomials, that is, for a given function $f$, can we find coefficients $c_i$ with $i = 1, \dots, n$ so that

$$f(x) \approx c_0 + c_1 x + c_2 x^2 + \dots + c_n x^n.$$

This question obviously needs to be refined. What do we mean by 'approximately equal'? This approximation formula can not hold for all functions $f$ and all $x$: the function $\sin x$ is bounded for all $x$, but any polynomial is unbounded for $x \to \pm\infty$, so any polynomial approximation to the $\sin x$ function is unbounded. Clearly we can only approximate on an interval.

We will show that a function $f$ with sufficiently many derivatives can be approximated as follows: if the $n$-th derivative $f^{(n)}$ is continuous on an interval $I$, then there are coefficients $c_0, \dots, c_{n-1}$ such that

$$\forall x \in I: |f(x) - \sum_{i<n} c_i x^i| \le cM_n \quad \text{where } M_n = \max_{x \in I} |f^{(n)}(x)|$$

It is easy to get inspiration for what these coefficients should be. Suppose

$$f(x) = c_0 + c_1 x + c_2 x^2 + \dots$$

(where we will not worry about matters of convergence and how long the dots go on) then filling in

$$x = 0 \text{ gives } c_0 = f(0).$$

Next, taking the first derivative

$$f'(x) = c_1 + 2c_2 x + 3c_3 x^2 + \dots$$

and filling in

$$x = 0 \text{ gives } c_1 = f'(0).$$
From the second derivative

$$f''(x) = 2c_2 + 6c_3x + \dots$$

so filling in $x = 0$ gives

$$c_2 = f''(0)/2.$$

Similarly, in the third derivative

$$\text{filling in } x = 0 \text{ gives } c_3 = \frac{1}{3!} f^{(3)}(0).$$

Now we need to be a bit more precise. Cauchy's form of Taylor's theorem says that

$$f(x) = f(a) + \frac{1}{1!}f'(a)(x-a) + \dots + \frac{1}{n!}f^{(n)}(a)(x-a)^n + R_n(x)$$

where the 'rest term' $R_n$ is

$$R_n(x) = \frac{1}{(n+1)!} f^{(n+1)}(\xi)(x-a)^{n+1} \quad \text{where } \xi \in (a, x) \text{ or } \xi \in (x, a) \text{ depending.}$$

If $f^{(n+1)}$ is bounded, and $x = a+h$, then the form in which we often use Taylor's theorem is

$$f(x) = \sum_{k=0}^{n} \frac{1}{k!} f^{(k)}(a)h^k + O(h^{n+1}).$$

We have now approximated the function $f$ on a certain interval by a polynomial, with an error that decreases geometrically with the inverse of the degree of the polynomial.

For a proof of Taylor's theorem we use integration by parts. First we write

$$\int_a^x f'(t)dt = f(x) - f(a)$$

as

$$f(x) = f(a) + \int_a^x f'(t)dt$$

Integration by parts then gives

$$\begin{aligned} f(x) &= f(a) + [xf'(x) - af'(a)] - \int_a^x tf''(t)dt \\ &= f(a) + [xf'(x) - xf'(a) + xf'(a) - af'(a)] - \int_a^x tf''(t)dt \\ &= f(a) + x \int_a^x f''(t)dt + (x-a)f'(a) - \int_a^x tf''(t)dt \\ &= f(a) + (x-a)f'(a) + \int_a^x (x-t)f''(t)dt \end{aligned}$$

Another application of integration by parts gives

$$f(x) = f(a) + (x-a)f'(a) + \frac{1}{2}(x-a)^2 f''(a) + \frac{1}{2} \int_a^x (x-t)^2 f'''(t)dt$$
Inductively, this gives us Taylor's theorem with

$$R_{n+1}(x) = \frac{1}{n!} \int_{a}^{x} (x-t)^{n} f^{(n+1)}(t) dt$$

By the mean value theorem this is

$$\begin{align*} R_{n+1}(x) &= \frac{1}{(n+1)!} f^{(n+1)}(\xi) \int_{a}^{x} (x-t)^{n} f^{(n+1)}(t) dt \\ &= \frac{1}{(n+1)!} f^{(n+1)}(\xi) (x-a)^{n+1} \end{align*}$$
# Chapter 16

## Graph theory

Graph theory is the branch of mathematics that studies pairwise relations between objects. Graphs both appear as tools for analyzing issues in HPC, and as objects of study themselves. This appendix introduces the basic concepts and some relevant theory.

### 16.1 Definitions

A graph consists of a set of objects, and set of relations between them. The objects, called the *nodes* or *vertices* of the graph, usually form a finite set, so we usually identify them with consecutive integers 1...n or 0...n - 1. The relation that holds between nodes is described by the *edges* of the graph: if i and j are related, we say that (i, j) is an edge of the graph. This relation does not need to be symmetric, take for instance the 'less than' relation.

Formally, then, a graph is a tuple $G = (V, E)$ where $V = \{1, ..., n\}$ for some $n$, and $E \subset \{(i, j): 1 \le i, j \le n, i \neq j\}$.

<figure><img src="image_129.png" /><figcaption>Figure 16.1: A simple graph</figcaption></figure>

A graph is called an *undirected graph* if $(i, j) \in E \Leftrightarrow (j, i) \in E$. The alternative is a *directed graph*, where we indicate an edge $(i, j)$ with an arrow from $i$ to $j$.

Two concepts that often appear in graph theory are the degree and the diameter of a graph.

**Definition 2** The degree denotes the maximum number of nodes that are connected to any node:

$$d(G) = \max_{i} |\{j: j \neq i \land (i, j) \in E\}|.$$
**Definition 3** The diameter of a graph is the length of the longest shortest path in the graph, where a path is defined as a set of vertices $v_1, \dots, v_{k+1}$ such that $v_i \neq v_j$ for all $i \neq j$ and

$$ \forall_{1 \le i \le k} : (v_i, v_{i+1}) \in E. $$

The length of this path is $k$.

The concept of diameter is illustrated in figure 16.2.

<figure><img src="image_130.png" /><figcaption>Figure 16.2: Shortest paths</figcaption></figure>

A path where all nodes are disjoint except for $v_1 = v_{k+1}$ is called a cycle.

Sometimes we are only interested in the mere existence of an edge $(i, j)$, at other times we attach a value or 'weight' $w_{ij}$ to that edge. A graph with weighted edges is called a *weighted graph*. Such a graph can be represented as a tuple $G = (V, E, W)$ where $E$ and $W$ have the same cardinality.

# 16.2 Common types of graphs

## 16.2.1 Directed Acyclic Graphs

A graph that does not have cycles is called *acyclic*. A special case of this type of graph is the *Directed Acyclic Graph (DAG)*. This type of graph can for instance be used to model dependencies between tasks: if there is an edge between $i, j$, it means that task $i$ has to be done before task $j$.
### 16.2.2 Trees

One special case of DAGs is the *tree graph*: here any node can have multiple incoming edges, but only one outgoing edge. Nodes with no incoming edges are *leaf nodes*; a node with no outgoing edges is called a root (can there be more than one root?), and all other nodes are called *interior nodes*.

## 16.3 Graph colouring and independent sets

We can assign labels to the nodes of a graph, which is equivalent to partitioning the set of nodes into disjoint subsets. One type of labeling that is of interest is *graph colouring*: here the labels (or ‘colours’) are chosen so that, if nodes *i* and *j* have the same colour, there is no edge connecting them: $(i, j) \notin E$.

There is a trivial colouring of a graph, where each node has its own colour. More interestingly, the minimum number of colours with which you can colour a graph is called the *colour number* of the graph.

**Exercise 16.1.** Show that, if a graph has degree *d*, the colour number is at most *d* + 1.

A famous graph colouring problem is the ‘four colour theorem’: if a graph depicts countries on a two-dimensional map (a so-called ‘planar’ graph), then the colour number is at most four. In general, finding the colour number is hard (in fact, NP-hard).

The colour sets of a graph colouring are also called *independent sets*, since within each colour no node is connected to a node of the same colour.

There is a trivial way of finding independent sets: declare each node to have its own unique colour. On the other hand, finding the ‘best’ division in independent sets, for instance through finding the colour number of the graph, is difficult. However, often it is enough to find a reasonable partitioning of the nodes into independent sets, for instance in constructing *paralell ILU preconditioners*; section 6.7.4. The following algorithm does that [93, 118]:

* Give each node a unique random number.

* Now find the set of nodes that have a higher number than all of their neighbours; call this the first independent set.

* Remove this set from the graph, and find again the nodes with a higher number than all their neighbours; this will be the second set.

* Repeat this procedure until all nodes are in an independent set.

**Exercise 16.2.** Convince yourself that the sets found this way are indeed independent.

## 16.4 Graphs and matrices

A graph can be rendered in a number of ways. You could of course just list nodes and edges, but little insight can be derived that way. Simple graphs can be visualized by drawing vertices and edges, but for large graphs this becomes unwieldy. Another option is to construct the *adjacency matrix* of the graph. For
a graph $G = \langle V, E \rangle$, the adjacency matrix $M$ (with a size $n$ equal to the number of vertices $|V|$) is defined by

$$M_{ij} = \begin{cases} 1 & (i, j) \in E \\ 0 & \text{otherwise} \end{cases}$$

Conversely, if you have a matrix, especially a sparse matrix, you can construct its adjacency graph. This is illustrated in figure 16.3 for both a dense and a sparse matrix. In this example, the matrices are structurally

<figure><img src="image_131.png" /><figcaption>Figure 16.3: A dense and a sparse matrix, both with their adjacency graph</figcaption></figure>

symmetric, so we use lines instead of arrows in the graphs. There is an edge on each vertex corresponding to the diagonal element; this edge will often be left out of illustrations.

For graphs with edge weights, we set the elements of the adjacency matrix to the weights:

$$M_{ij} = \begin{cases} w_{ij} & (i, j) \in E \\ 0 & \text{otherwise} \end{cases}$$

If a matrix has no zero elements, its adjacency graph has an edge between each pair of vertices. Such a graph is called a *clique*. If the graph is undirected, the adjacency matrix is symmetric, and conversely, if a matrix is *structurally symmetric*, its adjacency graph is undirected.

### 16.4.1 Permutation

Graphs are often used to indicate relations between objects in the real world. One example would be 'friend-of' relations in Facebook. In such cases, the nodes in a graph do not have a natural numbering: they are identified by a name and any numbering is artificial. Thus, we could wonder which graph properties remain invariant, and which ones change, if we apply a different numbering.

Renumbering a set of objects can be modeled algebraically by applying a *permutation matrix*.

**Definition 4** A permutation matrix is a square matrix where each row and column has exactly one element equal to one; all other elements are zero.
**Exercise 16.3.** Let a set of $N$ objects $x_1, \dots, x_N$ be given. What is the permutation matrix that orders them as $x_1, x_3, \dots, x_2, x_4, \dots$? That is, find the matrix $P$ such that

$$ \begin{pmatrix} x_1 \\ x_3 \\ \vdots \\ x_2 \\ x_4 \\ \vdots \end{pmatrix} = P \begin{pmatrix} x_1 \\ \vdots \\ x_N \end{pmatrix} $$

**Exercise 16.4.** Show that the eigenvalues of a matrix are invariant under permutation.

### 16.4.2 Irreducibility

As an example of graph concepts that has an easy interpretation in the adjacency matrix, consider reducibility.

**Definition 5** A graph is called irreducible if for every pair $i, j$ of nodes there is a path from $i$ to $j$ and from $j$ to $i$. A graph is reducible if it is not irreducible.

**Exercise 16.5.** Let $A$ be a matrix

$$ A = \begin{pmatrix} B & C \\ \emptyset & D \end{pmatrix} $$

where $B$ and $D$ are square matrices. Prove the reducibility of the graph of which this is the adjacency matrix.

If we permute graph, its reducibility or irreducibility is not changed. However, it may now no longer be apparent from looking at the adjacency matrix.

### 16.4.3 Graph closure

Here is another example of how adjacency matrices can simplify reasoning about graphs.

**Exercise 16.6.** Let $G = \langle V, E \rangle$ be an undirected graph, and let $G' = \langle V, E' \rangle$ be the graph with the same vertices, but with vertices defined by

$$ (i, j) \in E' \Leftrightarrow \exists_k : (i, k) \in E \land (k, j) \in E. $$

If $M$ is the adjacency matrix of $G$, show that $M^2$ is the adjacency matrix of $G'$, where we use boolean multiplication on the elements: $1 \cdot 1 = 1$, $1 + 1 = 1$.

### 16.4.4 Graph operations in linear algebra

In most of the above, the adjacency matrix was nothing more than a table. We will now show that we can actually do linear algebra with it, making it deserving of the name 'matrix'. As a simple example of using
linear algebra on an adjacency matrix $G$, let $e$ the vector of all 1s, then $Ge$ is the vector that lists the degrees of the nodes.

We can do many operations this way. Consider a weighted graph $G$. Finding for each node $i$ the largest weight $g_{ij}$ can be described as

$$y = G \otimes e \quad \text{where} \quad y_i = \max_j g_{ij} \cdot 1$$

This looks like the regular matrix-vector product $Ge$, but with the sum replaced by a maximum calculation.

In many cases we actually need the left product, that is, multiplying the adjacency matrix from the left by a row vector. Let for example $e_i$ be the vector with a 1 in location $i$ and zero elsewhere. Then $e_i^t G$ has a one in every $j$ that is a neighbour of $i$ in the adjacency graph of $G$.

### 16.4.4.1 Markov chains

This left matrix-vector product has a neat application: Markov chains. Suppose we have a system (see for instance 17) that can be in any of $n$ states. We can then use the adjacency matrix to model state transitions.

Let $G$ by an adjacency matrix with the property that all elements are non-negative, and the elements in each row sum to 1. We now interpret $g_{ij}$ as the probability of going from state $i$ to state $j$. Let $x$ be a probability vector, that is, $x_i$ is the nonegative probability of being in state $i$, then $y^t = x^t G$ describes these probabilities, after one state transition.

**Exercise 16.7.** For a vector $x$ to be a proper probability vector, its elements need to sum to 1. Show that the elements of $y$ again sum to 1.

### 16.4.4.2 General matrix-vector product

The examples above showed that sometimes we perform an operation on an adjacency matrix that has the structure of a matrix-vector product, but may not necessarily use addition and multiplication. (Many more examples of this can be found in section 9.1.)

Motivated by this, we define a general product

$$y = G \oplus \otimes x$$

where $\otimes$ is a binary operator, and $\oplus$ a reduction operator, as

$$y_i = \bigoplus_j (g_{ij} \otimes x_j).$$

In this notation, finding the largest weight would be

$$w = G \max \cdot \times e.$$

This will be used in several algorithms in chapter 9; for other applications see [105].
## 16.5 Spectral graph theory

With a graph $G^1$ and its adjacency matrix $A_G$, we can define a stochastic matrix or Markov matrix by scaling $A_G$ to have unit row sums:

$$W_G = D_G^{-1} A_G \quad \text{where } (D_G)_{ii} = \deg(i).$$

To see how we interpret this, let's look at a simple example. Let's take an unweighted graph with an adjacency matrix

$$A_G = \begin{pmatrix} 1 & 1 & 1 \\ 1 & 1 & 1 \\ 1 & 1 & 1 \\ 1 & 1 & 1 \end{pmatrix}$$

and look at the second row, which says that there are edges (2, 3) and (2, 4). This means that if you are on node 2, you can go to nodes 3 and 4. Scaling this matrix we get

$$W_G = \begin{pmatrix} 1/3 & 1/3 & 1/3 \\ 1/3 & 1/2 & 1/2 \\ 1/3 & 1/3 & 1/3 \\ 1/2 & 1/2 & 1/2 \end{pmatrix}$$

and now the second row says that from node 2 you can get with equal probability to nodes 3 and 4. You can also derive this statement mathematically:

$$(0 \ 1 \ 0 \ 0) W_G = (0 \ 0 \ 1/2 \ 1/2)$$

It is simple to extrapolate that: if $p$ is a vector where the $i$-th component gives the probability of being in node $i$, then $(p^t W_G)_i$ is the probability of being in node $i$ if you take one more step along a graph edge.

**Exercise 16.8.** Prove that $p^t W_G$ is indeed a vector of probabilities. Hint: you can express that $p$ is a probability vector as $p^t e = e$, where $e$ is the vector of all ones.

Another matrix to associate with a graph is the *graph Laplacian*

$$L_G = D_G - A_G.$$

This matrix has zero rowsums and positive diagonal entries, so by the Gershgorin theorem (section 12.5 all its eigenvalues are in the complex right half plane.

**Exercise 16.9.** Show that the vector of all ones is an eigenvector with eigenvalue 1.

This Laplacian matrix gives us a quadratic form:

$$x^t L_G x = \sum_{(i,j) \in E} (x_i - x_j)^2.$$

<footnote n="1">1. This section owes much to Dan Spielman's course on spectral graph theory http://www.cs.yale.edu/homes/spielman/561/.</footnote>
### 16.5.1 Eigenvalue theorems for graph matrices

There are various interesting theorems connected with the graph adjacency and Laplacian matrix. We start by remarking that the following theorems may not seem intuitive, but they closely resemble some simple facts about PDEs.

Connected with the Laplace equation $-\Delta u = f$ is an operator $\mathcal{L}u = -\Delta u$. On the unit interval $[0, 1]$ the eigenfunctions of this operator, that is, the functions for which $\mathcal{L}u = \lambda u$, are $u_n(x) = \sin n\pi x$ for $n > 0$. These have the property that $u_n(x)$ has $n-1$ zeros in the interior of the interval, and they divide the interval in $n$ connected regions where the function is positive of negative. Thus, if you wanted to divide a domain $\Omega$ over $p$ processors, you could consider the $p$-th eigenfunction of the Laplacian on $\Omega$, and find the connected regions where it is positive or negative.

This statement about PDE has a graph equivalent in two versions of Fiedler's theorem. (We will not give any proofs in this section; see [149].)

**Theorem 6** Let $G$ be a weighted path graph on $n$ vertices, let $L_P$ have eigenvalues $0 = \lambda_1 < \lambda_2 \le \dots \le \lambda_n$, and let $v_k$ be an eigenvector of $\lambda_k$. Then $v_k$ changes sign $k-1$ times.

The second theorem is more useful [55]:

**Theorem 7** Let $G = (V, E, w)$ be a weighted connected graph, and let $L_G$ be its Laplacian matrix. Let $0 = \lambda_1 < \lambda_2 \le \dots \le \lambda_n$ be the eigenvalues of $L_G$ and let $v_1, \dots, v_n$ be the corresponding eigenvectors. For any $k \ge 2$, let $W_k = \{i \in V: v_k(i) \ge 0\}$. Then, the graph induced by $G$ on $W_k$ has at most $k-1$ connected components.

The important consequence of this is that the eigenvector to the first nontrivial eigenvalue can be used to partition the graph in two connected pieces: one of nodes where the eigenvector is positive, and one where the eigenvector is negative. This eigenvector is known as the *Fiedler vector*.

The adjacency matrix is nonnegative, and there is an extensive theory for this type of matrix [9]; see the Perron-Frobenius theorem in section 12.4.

**Lemma 1** The largest eigenvalue of the adjacency matrix can be bounded by the maximum node degree:

$$ \alpha_1 \le d_{\max}. $$

### 16.5.2 Cheeger's inequality

Above we remarked that the first non-trivial eigenvalue of the graph Laplacian has a relation to partitioning a graph in two parts. The Cheeger's constant and Cheeger's inequality relate this eigenvalue to a certain quality measure of partitionings.

Let $V$ be the set of vertices and $S \subset V$, then Cheeger's constant of a graph is defined as

$$ C = \min_S \frac{e(S, V-S)}{\min_S \min \operatorname{vol}(S), \operatorname{vol}(V-S)} $$
16. Graph theory

where $e(S, V - S)$ denotes the number of edges connecting $S$ to $V - S$, and the volume of a set of nodes is defined as

$$ \mathrm{vol}(S) = \sum_{e \in S} d(e). $$

Cheeger's inequality then states

$$ 2C \geq \lambda \geq \frac{C^2}{2} $$

where $\lambda$ is the first nontrivial eigenvalue of the graph Laplacian.
# Chapter 17

## Automata theory

Automata are mathematical abstractions of machines. There is an extensive theory of automata; here we will only touch on the basic concepts. Let us start with a simple example.

### 17.1 Finite State Automata

A Finite State Automaton (FSA) is a very simple machine, something on the order of a vending machine that will dispense a candy bar when a quarter has been inserted. There are four actions possible with a vending machine: insert a quarter, press 'coin return' to ask for any inserted money back, open the window to take the candy bar, and close the window again. Whether an action is possible (especially the third) depends on the state the machine is in. There are three states: the begin state, the state where the quarter has been inserted and the window unlocked (let us call this 'ready to dispense'), and the state where the window is open (which we will call 'dispensing').

In certain states, certain actions are not possible. For instance, in the beginning state the window can not be opened.

The mathematical description of this vending machine consists of 1. the list of states, 2. a table of how the possible actions make the machine go from one state to another. However, rather than writing down the table, a graphical representation is usually more insightful.

### 17.2 General discussion

From the vending machine example, you saw an important characteristic of an automaton: it allows certain actions, but only in some circumstances, and in the end it has a state corresponding to 'success', which is only reached if certain sequences of actions are taken.

The formal description of this process is as follows: we call the individual actions an 'alphabet', and a sequence of actions a 'word' from based on that alphabet. The 'success' result of a machine then corresponds to a verdict that a certain word belong to the 'language' that is accepted by the automaton. Thus there is a correspondence between automata theory and language theory.
<figure><img src="image_132.png" /><figcaption>Figure 17.1: A simple real-life automaton</figcaption></figure>

Exercise 17.1. Consider the alphabet $\{a, b\}$, that is, the alphabet with only the letters $a$ and $b$, and consider the language $\{a^m b^n : m, n > 0\}$, that is the words that consist of one or more $a$s followed by one or more $b$s. Draw the automaton that accepts this language.

What makes the FSA the simplest type is that it has no memory. Most vending machines do not complain if you put in more than one quarter: they have no memory beyond 'a quarter has been inserted'. A more complicated machine would count how many quarters you inserted, and then allow you to open that many windows to different candy bars. In the above formal way of describing, that machine would accept the language $\{q^n w^n : n \ge 0\}$, that is, the sequences where you deposit as many quarters ('q') as you open windows ('w'). This language is an example of a so-called *context-free language*; the language of the original vending machine is a *regular language*.

These two language types belong to the four level *Chomsky hierarchy* of languages. The famous *Turing machine*, which recognizes the *recursively enumerable language* type, is on the top level of the hierarchy. The missing step has the *context-sensitive language* type, which is recognized by a *linear bounded automaton*.
# Chapter 18

## Parallel Prefix

For operations to be executable in parallel they need to be independent. That makes recurrences problematic to evaluate in parallel. Recurrences occur in obvious places such as solving a triangular system of equations (section 5.3.5), but they can also appear in sorting and many other operations.

In this appendix we look at *parallel prefix* operations: the parallel execution of an operation that is defined by a recurrence involving an associative operator. Computing the sum of an array of elements is an example of this type of operation (disregarding the *non-associativity* for the moment). Let $\pi(x, y)$ be the binary sum operator:

$$ \pi(x, y) \equiv x + y, $$

then we define the prefix sum of $n \ge 2$ terms as

$$ \Pi(x_1, \dots, x_n) = \begin{cases} \pi(x_1, x_2) & \text{if } n = 2 \\ \pi(\Pi(x_1, \dots, x_{n-1}), x_n) & \text{otherwise} \end{cases} $$

As a non-obvious of a prefix operation, we could count the number of elements of an array that have a certain property.

**Exercise 18.1.** Let $p(\cdot)$ be a predicate, $p(x) = 1$ if it holds for $x$ and $0$ otherwise. Define a binary operator $\pi(x, y)$ so that its reduction over an array of numbers yields the number of elements for which $p$ is true.

So let us now assume the existence of an associative operator $\oplus$, an array of values $x_1, \dots, x_n$. Then we define the prefix problem as the computation of $X_1, \dots, X_n$, where

$$ \begin{cases} X_1 = x_1 \\ X_k = \oplus_{i \le k} x_i \end{cases} $$

### 18.1 Parallel prefix

The key to parallelizing this is the realization that we can compute partial reductions in parallel:

$$ x_1 \oplus x_2, \quad x_3 \oplus x_4, \dots $$
are all independent. Furthermore, partial reductions of these reductions,

$$ (x_1 \oplus x_2) \oplus (x_3 \oplus x_4), \dots $$

are also independent. We use the notation

$$ X_{i,j} = x_i \oplus \dots \oplus x_j $$

for these partial reductions.

You have seen this before in section 2.1: an array of n numbers can be reduced in ⌈log₂n⌉ steps. What is missing to make this a full prefix operation is computation of all intermediate values.

Observing that, for instance, $X_3 = (x_1 \oplus x_2) \oplus x_3 = X_2 \oplus x_3$, you can now easily imagine the whole process; see figure 18.1 for the case of 8 elements. To compute, say, $X_{13}$, you express $13 = 8 + 4 + 1$ and

<figure><img src="image_133.png" /><figcaption>Figure 18.1: A prefix operation applied to 8 elements</figcaption></figure>

compute $X_{13} = X_8 \oplus X_{9,12} \oplus x_{13}$.

In this figure, operations over the same ‘distance’ have been horizontally aligned corresponding to a SIMD type execution. If the execution proceeds with a task graph, some steps can be executed earlier than the figure suggests; for instance $X_3$ can be computed simultaneously with $X_6$.

Regardless the arrangement of the computational steps, it is not hard to see that the whole prefix calculation can be done in $2 \log_2 n$ steps: $\log_2 n$ steps for computing the final reduction $X_n$, then another $\log_2 n$ steps for filling in the intermediate values.

## 18.2 Sparse matrix vector product as parallel prefix

It has been observed that the sparse matrix vector product can be considered a prefix operation; see [14]. The reasoning here is we first compute all $y_{ij} = a_{ij}x_j$, and subsequently compute the sums $y_i = \sum_j y_{ij}$ with a prefix operation.
A prefix sum as explained above does not compute the right result. The first couple of $y_{ij}$ terms do indeed sum to $y_1$, but then continuing the prefix sum gives $y_1 + y_2$, instead of $y_2$. The trick to making this work is to consider two-component quantities $\langle y_{ij}, s_{ij} \rangle$, where

$$s_{ij} = \begin{cases} 1 & \text{if } j \text{ is the first nonzero index in row } i \\ 0 & \text{otherwise} \end{cases}$$

Now we can define prefix sums that are 'reset' every time $s_{ij} = 1$.
# PART IV

## TUTORIALS
In the theory part of this book you learned mathematical models can be translated to algorithms that can be realized efficiently on modern hardware. You learned how data structures and coding decisions influence the performance of your code. In other words, you should now have all the tools to write programs that solve scientific problems.

This would be all you would need to know, if there was any guarantee that a correctly derived algorithm and well designed data structure could immediately be turned into a correct program. Unfortunately, there is more to programming than that. This collection of tutorials will give you the tools to be an effective scientific programmer.

The vast majority of scientific programming is done on the Unix platform so we start out with a tutorial on Unix in section 19, followed by an explanation of how your code is handled by compilers and linkers and such in section 20.

Next you will learn about some tools that will increase your productivity and effectiveness:

* The *Make* utility is used for managing the building of projects; section 21.

* Source control systems store your code in such a way that you can undo changes, or maintain multiple versions; in section 22 you will see the *subversion* software.

* Storing and exchanging scientific data becomes an important matter once your program starts to produce results; in section 23 you will learn the use of *HDF5*.

* A lot of functionality that you need has been coded already by other people; in section 24 you will learn about some of the scientific libraries that are available to you.

* Visual output of program data is important, but too wide a topic to discuss here in great detail; section 25 teaches you the basics of the *gnuplot* package, which is suitable for simple data plotting.

We also consider the activity of program development itself: section 26 considers how to code to prevent errors, and section 27 teaches you to debug code with *gdb*. Section 28 contains some information on how to write a program that uses more than one programming language.

Finally, section 29 teaches you about the *LATEX* document system, so that you can report on your work in beautifully typeset articles.

Many of the tutorials are very hands-on. Do them while sitting at a computer!
# Chapter 19

## Unix intro

Unix is an *Operating System (OS)*, that is, a layer of software between the user or a user program and the hardware. It takes care of files and screen output, and it makes sure that many processes can exist side by side on one system. However, it is not immediately visible to the user. Most of the time that you use Unix, you are typing commands which are executed by an interpreter called the *shell*. The shell makes the actual OS calls. There are a few possible Unix shells available, but in this tutorial we will assume that you are using the *sh* or *bash* shell, although many commands are common to the various shells in existence.

This short tutorial will get you going; if you want to learn more about Unix and shell scripting, see for instance http://www.tldp.org/guides.html. Most of this tutorial will work on any Unix-like platform, including *Cygwin* on Windows. However, there is not just one Unix:

* Traditionally there are a few major flavours of Unix: ATT and BSD. Apple has Darwin which is close to BSD; IBM and HP have their own versions of Unix, and Linux is yet another variant. The differences between these are deep down and if you are taking this tutorial you probably won't see them for quite a while.

* Within Linux there are various *Linux distributions* such as *Red Hat* or *Ubuntu*. These mainly differ in the organization of system files and again you probably need not worry about them.

* As mentioned just now, there are different shells, and they do differ considerably. Here you will learn the *bash* shell, which for a variety of reasons is to be preferred over the *csh* or *tcsh* shell. Other shells are the *ksh* and *zsh*.

### 19.1 Files and such

**Purpose.** In this section you will learn about the Unix file system, which consists of *directories* that store *files*. You will learn about *executable* files and commands for displaying data files.

#### 19.1.1 Looking at files

**Purpose.** In this section you will learn commands for displaying file contents.
The `ls` command gives you a listing of files that are in your present location.

**Exercise.** Type `ls`. Does anything show up?

*Expected outcome.* If there are files in your directory, they will be listed; if there are none, no output will be given. This is standard Unix behaviour: no output does not mean that something went wrong, it only means that there is nothing to report.

The `cat` command is often used to display files, but it can also be used to create some simple content.

**Exercise.** Type `cat > newfilename` (where you can pick any filename) and type some text. Conclude with `Control-d` on a line by itself<fnref n="1" />. Now use `cat` to view the contents of that file: `cat newfilename`.

*Expected outcome.* In the first use of `cat`, text was concatenated from the terminal to a file; in the second the file was cat'ed to the terminal output. You should see on your screen precisely what you typed into the file.

*Caveats.* Be sure to type `Control-d` as the first thing on the last line of input. If you really get stuck, `Control-c` will usually get you out. Try this: start creating a file with `cat > filename` and hit `Control-c` in the middle of a line. What are the contents of your file?

**Remark 1** Instead of `Control-d` you will often see the notation `^D`. The capital letter is for historic reasons: you use the control key and the lowercase letter.

Above you used `ls` to get a directory listing. You can also use the `ls` command on specific files:

**Exercise.** Do `ls newfilename` with the file that you created above; also do `ls nosuchfile` with a file name that does not exist.

*Expected outcome.* For an existing file you get the file name on your screen; for a non-existing file you get an error message.

The `ls` command can give you all sorts of information.

**Exercise.** Read the man page of the `ls` command: `man ls`. Find out the size and creation date of some files, for instance the file you just created.

*Expected outcome.* Did you find the `ls -s` and `ls -l` options? The first one lists the size of each file, usually in kilobytes, the other gives all sorts of information about a file, including things you will learn about later.

*Caveats.* The `man` command puts you in a mode where you can view long text documents. This viewer is common on Unix systems (it is available as the `more` or `less` system command), so memorize the following ways of navigating: Use the space bar to go forward and the `u` key to go back up. Use `g` to go to the beginning of the text, and `G` for the end. Use `q` to exit the viewer. If you really get stuck, `Control-c` will get you out.

<footnote n="1">1. Press the `Control` and hold it while you press the `d` key.</footnote>
(If you already know what command you're looking for, you can use man to get online information about it. If you forget the name of a command, `man -k keyword` can help you find it.)

The `touch` command creates an empty file, or updates the timestamp of a file if it already exists. Use `ls -l` to confirm this behaviour.

Three more useful commands for files are: `cp` for copying, `mv` (short for 'move') for renaming, and `rm` ('remove') for deleting. Experiment with them.

There are more commands for displaying a file, parts of a file, or information about a file.

**Exercise.** Do `ls /usr/share/words` or `ls /usr/share/dict/words` to confirm that a file with words exists on your system. Now experiment with the commands `head`, `tail`, `more`, and `wc` using that file.

*Expected outcome.* `head` displays the first couple of lines of a file, `tail` the last, and `more` uses the same viewer that is used for man pages. Read the man pages for these commands and experiment with increasing and decreasing the amount of output. The `wc` ('word count') command reports the number of words, characters, and lines in a file.

Another useful command is `which`: it tells you what type of file you are dealing with. See what it tells you about one of the text files you just created.

### 19.1.2 Directories

**Purpose.** Here you will learn about the Unix directory tree, how to manipulate it and how to move around in it.

A unix file system is a tree of directories, where a directory is a container for files or more directories. We will display directories as follows:

<table><tr><td>/</td><td>The root of the directory tree</td></tr><tr><td>bin</td><td>Binary programs</td></tr><tr><td>home</td><td>Location of users directories</td></tr></table>

The root of the Unix directory tree is indicated with a slash. Do `ls /` to see what the files and directories there are in the root. Note that the root is not the location where you start when you reboot your personal machine, or when you log in to a server.

**Exercise.** The command to find out your current working directory is `pwd`. Your home directory is your working directory immediately when you log in. Find out your home directory.

*Expected outcome.* You will typically see something like `/home/yourname` or `/Users/yourname`. This is system dependent.

Do `ls` to see the contents of the working directory. In the displays in this section, directory names will be followed by a slash: `dir/` but this character is not part of their name. You can get this output by using `ls -F`, and you can tell your shell to use this output consistently by stating `alias ls=ls -F` at the start of your session. Example:
<figure><img src="image_134.png" /></figure>

The command for making a new directory is `mkdir`.

**Exercise.** Make a new directory with `mkdir newdir` and view the current directory with `ls`

*Expected outcome.* You should see this structure:

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

The command for going into another directory, that is, making it your working directory, is `cd` ('change directory'). It can be used in the following ways:

* `cd` Without any arguments, `cd` takes you to your home directory.

* `cd <absolute path>` An absolute path starts at the root of the directory tree, that is, starts with `/`. The `cd` command takes you to that location.

* `cd <relative path>` A relative path is one that does not start at the root. This form of the `cd` command takes you to `<yourcurrentdir>/<relative path>`.

**Exercise.** Do `cd newdir` and find out where you are in the directory tree with `pwd`. Confirm with `ls` that the directory is empty. How would you get to this location using an absolute path?

*Expected outcome.* `pwd` should tell you `/home/you/newdir`, and `ls` then has no output, meaning there is nothing to list. The absolute path is `/home/you/newdir`.

**Exercise.** Let's quickly create a file in this directory: `touch onefile`, and another directory: `mkdir otherdir`. Do `ls` and confirm that there are a new file and directory.

*Expected outcome.* You should now have:

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

The `ls` command has a very useful option: with `ls -a` you see your regular files and hidden files, which have a name that starts with a dot. Doing `ls -a` in your new directory should tell you that there are the following files:

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

The single dot is the current directory, and the double dot is the directory one level back.
**Exercise.** Predict where you will be after `cd ./otherdir/..` and check to see if you were right.

**Expected outcome.** The single dot sends you to the current directory, so that does not change anything. The `otherdir` part makes that subdirectory your current working directory. Finally, .. goes one level back. In other words, this command puts your right back where you started.

Since your home directory is a special place, there are shortcuts for cd'ing to it: `cd` without arguments, `cd ~`, and `cd $HOME` all get you back to your home.

Go to your home directory, and from there do `ls newdir` to check the contents of the first directory you created, without having to go there.

**Exercise.** What does `ls ..` do?

**Expected outcome.** Recall that .. denotes the directory one level up in the tree: you should see your own home directory, plus the directories of any other users.

**Exercise.** Can you use `ls` to see the contents of someone else's home directory? In the previous exercise you saw whether other users exist on your system. If so, do `ls ../thatotheruser`.

**Expected outcome.** If this is your private computer, you can probably view the contents of the other user's directory. If this is a university computer or so, the other directory may very well be protected – permissions are discussed in the next section – and you get `ls: ../otheruser: Permission denied`.

Make an attempt to move into someone else's home directory with `cd`. Does it work?

You can make copies of a directory with `cp`, but you need to add a flag to indicate that you recursively copy the contents: `cp -r`. Make another directory `somedir` in your home so that you have

<table><thead><tr><th>/home/you/<br>newdir/</th><th>you have been working in this one</th></tr></thead><tbody><tr><td>somedir/</td><td>you just created this one</td></tr></tbody></table>

What is the difference between `cp -r newdir somedir` and `cp -r newdir thirddir` where `thirddir` is not an existing directory name?

### 19.1.3 Permissions

**Purpose.** In this section you will learn about how to give various users on your system permission to do (or not to do) various things with your files.

Unix files, including directories, have permissions, indicating 'who can do what with this file'. Actions that can be performed on a file fall into three categories:

* `reading r:` any access to a file (displaying, getting information on it) that does not change the file;

* `writing w:` access to a file that changes its content, or even its metadata such as 'date modified';

* `executing x:` if the file is executable, to run it; if it is a directory, to enter it.
The people who can potentially access a file are divided into three classes too:

* the user *u*: the person owning the file;

* the group *g*: a group of users to which the owner belongs;

* other *o*: everyone else.

These nine permissions are rendered in sequence

<table><thead><tr><th><em>user</em></th><th><em>group</em></th><th><em>other</em></th></tr></thead><tbody><tr><td><em>rwx</em></td><td><em>rwx</em></td><td><em>rwx</em></td></tr></tbody></table>

For instance `rw-r--r--` means that the owner can read and write a file, the owner's group and everyone else can only read.

Permissions are also rendered numerically in groups of three bits, by letting $r = 4$, $w = 2$, $x = 1$:

<table><tr><td><em>rwx</em></td></tr><tr><td>421</td></tr></table>

Common codes are 7 = `rwx` and 6 = `rw`. You will find many files that have permissions 755 which stands for an executable that everyone can run, but only the owner can change, or 644 which stands for a data file that everyone can see but again only the owner can alter. You can set permissions by the `chmod` command:

`chmod <permissions> file # just one file`

`chmod -R <permissions> directory # directory, recursively`

Examples:

`chmod 766 file # set to rwxrw-rw-`

`chmod g+w file # give group write permission`

`chmod g=rx file # set group permissions`

`chmod o-w file # take away write permission from others`

`chmod o= file # take away all permissions from others.`

`chmod g+r,o-x file # give group read permission`

# remove other execute permission

The man page gives all options.

**Exercise.** Make a file `foo` and do `chmod u-r foo`. Can you now inspect its contents? Make the file readable again, this time using a numeric code. Now make the file readable to your classmates. Check by having one of them read the contents.

*Expected outcome.* When you've made the file 'unreadable' by yourself, you can still `ls` it, but not `cat` it: that will give a 'permission denied' message.

Make a file `com` with the following contents:

```
#!/bin/sh
```

```
echo "Hello world!"
```
This is a legitimate shell script. What happens when you type `./com`? Can you make the script executable?

In the three permission categories it is who 'you' and 'others' refer to. How about 'group'? The command *groups* tells you all the groups you are in, and `ls -l` tells you what group a file belongs to. Analogous to `chmod`, you can use `chgrp` to change the group to which a file belongs, to share it with a user who is also in that group. Adding a user to a group sometimes needs system privileges.

### 19.1.4 Wildcards

You already saw that `ls filename` gives you information about that one file, and `ls` gives you all files in the current directory. To see files with certain conditions on their names, the *wildcard* mechanism exists. The following wildcards exist:

* any number of characters.

? any character.

Example:

```
%% ls
s      sk     ski    skiing   skill

%% ls ski*
ski    skiing   skill
```

The second option lists all files whose name start with `ski`, followed by any number of other characters'; below you will see that in different contexts `ski*` means 'sk followed by any number of i characters'. Confusing, but that's the way it is.

## 19.2 Text searching and regular expressions

**Purpose.** In this section you will learn how to search for text in files.

For this section you need at least one file that contains some amount of text. You can for instance get random text from http://www.lipsum.com/feed/html.

The `grep` command can be used to search for a text expression in a file.

**Exercise.** Search for the letter `q` in your text file with `grep q yourfile` and search for it in all files in your directory with `grep q *`. Try some other searches.

**Expected outcome.** In the first case, you get a listing of all lines that contain a `q`; in the second case, `grep` also reports what file name the match was found in: `qfile>this line has q in it`.

**Caveats.** If the string you are looking for does not occur, `grep` will simply not output anything. Remember that this is standard behaviour for Unix commands if there is nothing to report.
In addition to searching for literal strings, you can look for more general expressions.

<table><tr><td>^</td><td>the beginning of the line</td></tr><tr><td>$</td><td>the end of the line</td></tr><tr><td>.</td><td>any character</td></tr><tr><td>*</td><td>any number of repetitions</td></tr><tr><td>[xyz]</td><td>any of the characters xyz</td></tr></table>

This looks like the wildcard mechanism you just saw (section 19.1.4) but it's subtly different. Compare the example above with:

```
%% cat s
sk
ski
skill
skiing
%% grep "ski*" s
sk
ski
skill
skiing
```

In the second case you search for a string consisting of `sk` and any number of `i` characters, including zero of them.

Some more examples: you can find

* All lines that contain the letter 'q' with `grep q yourfile`;

* All lines that start with an 'a' with `grep "^a" yourfile` (if your search string contains special characters, it is a good idea to use quote marks to enclose it);

* All lines that end with a digit with `grep "[0-9]$" yourfile`.

**Exercise.** Construct the search strings for finding

* lines that start with an uppercase character, and

* lines that contain exactly one character.

*Expected outcome.* For the first, use the range characters `[]`, for the second use the period to match any character.

**Exercise.** Add a few lines `x = 1, x = 2, x = 3` (that is, have different numbers of spaces between `x` and the equals sign) to your test file, and make `grep` commands to search for all assignments to `x`.

The characters in the table above have special meanings. If you want to search that actual character, you have to escape it.

**Exercise.** Make a test file that has both `abc` and `a.c` in it, on separate lines. Try the commands `grep "a.c" file`, `grep a\.c file`, `grep "a\.c" file`.

*Expected outcome.* You will see that the period needs to be escaped, and the search string needs to be quoted. In the absence of either, you will see that `grep` also finds the `abc` string.
### 19.2.1 Stream editing with `sed`

Unix has various tools for processing text files on a line-by-line basis. The stream editor `sed` is one example. If you have used the `vi` editor, you are probably used to a syntax like `s/foo/bar/` for making changes. With `sed`, you can do this on the commandline. For instance

```
sed 's/foo/bar/' myfile > mynewfile
```

will apply the substitute command `s/foo/bar/` to every line of `myfile`. The output is shown on your screen so you should capture it in a new file; see section 19.3.2 for more on output redirection.

### 19.2.2 Cutting up lines with `cut`

Another tool for editing lines is `cut`, which will cut up a line and display certain parts of it. For instance,

```
cut -c 2-5 myfile
```

will display the characters in position 2-5 of every line of `myfile`. Make a test file and verify this example.

Maybe more useful, you can give `cut` a delimiter character and have it split a line on occurrences of that delimiter. For instance, your system will mostly likely have a file `/etc/passwd` that contains user information<fnref n="2" />, with every line consisting of fields separated by colons. For instance:

```
daemon:*:1:1:System Services:/var/root:/usr/bin/false
```

```
nobody:*:-2:-2:Unprivileged User:/var/empty:/usr/bin/false
```

```
root:*:0:0:System Administrator:/var/root:/bin/sh
```

The seventh and last field is the login shell of the user; `/bin/false` indicates that the user is unable to log in.

You can display users and their login shells with:

```
cut -d ":" -f 1,7 /etc/passwd
```

This tells `cut` to use the colon as delimiter, and to print fields 1 and 7.

## 19.3 Command execution

### 19.3.1 Search paths

**Purpose.** In this section you will learn how Unix determines what to do when you type a command name.

<footnote n="2">2. This is traditionally the case; on Mac OS information about users is kept elsewhere and this file only contains system services.</footnote>
If you type a command such as `ls`, the shell does not just rely on a list of commands: it will actually go searching for a program by the name `ls`. This means that you can have multiple different commands with the same name, and which one gets executed depends on which one is found first.

**Exercise.** What you may think of as 'Unix commands' are often just executable files in a system directory. Do which `ls`, and do an `ls -l` on the result

*Expected outcome.* The location of `ls` is something like `/bin/ls`. If you `ls` that, you will see that it is probably owned by root. Its executable bits are probably set for all users.

The locations where unix searches for commands is the 'search path', which is stored in the *environment variable* (for more details see below) `PATH`.

**Exercise.** Do `echo $PATH`. Can you find the location of `cd`? Are there other commands in the same location? Is the current directory `..` in the path? If not, do `export PATH=".:$PATH"`. Now create an executable file `cd` in the current director (see above for the basics), and do `cd`.

*Expected outcome.* The path will be a list of colon-separated directories, for instance `/usr/bin:/usr/local/bin:/usr/X11R6/bin`. If the working directory is in the path, it will probably be at the end: `/usr/X11R6/bin:`. but most likely it will not be there. If you put `..` at the start of the path, unix will find the local `cd` command before the system one.

Some people consider having the working directory in the path a security risk. If your directory is writable, someone could put a malicious script named `cd` (or any other system command) in your directory, and you would execute it unwittingly.

It is possible to define your own commands as aliases of existing commands.

**Exercise.** Do `alias chdir=cd` and convince yourself that now `chdir` works just like `cd`. Do `alias rm='rm -i'`; look up the meaning of this in the man pages. Some people find this alias a good idea; can you see why?

*Expected outcome.* The `-i` 'interactive' option for `rm` makes the command ask for confirmation before each delete. Since unix does not have a trashcan that needs to be emptied explicitly (as on Windows or the Mac OS), this can be a good idea.

### 19.3.2 Redirection and Pipelines

**Purpose.** In this section you will learn how to feed one command into another, and how to connect commands to input and output files.

So far, the unix commands you have used have taken their input from your keyboard, or from a file named on the command line; their output went to your screen. There are other possibilities for providing input from a file, or for storing the output in a file.
### 19.3.2.1 Input redirection

The `grep` command had two arguments, the second being a file name. You can also write `grep string < yourfile`, where the less-than sign means that the input will come from the named file, `yourfile`. This is known as *input redirection*.

### 19.3.2.2 Output redirection

Conversely, `grep string yourfile > outfile` will take what normally goes to the terminal, and *redirect* the output to *outfile*. The output file is created if it didn't already exist, otherwise it is overwritten. (To append, use `grep text yourfile >> outfile`.)

**Exercise.** Take one of the grep commands from the previous section, and send its output to a file. Check that the contents of the file are identical to what appeared on your screen before. Search for a string that does not appear in the file and send the output to a file. What does this mean for the output file?

*Expected outcome.* Searching for a string that does not occur in a file gives no terminal output. If you redirect the output of this `grep` to a file, it gives a zero size file. Check this with `ls` and `wc`.

### 19.3.2.3 Standard files

Unix has three standard files that handle input and output:

`stdin` is the file that provides input for processes.

`stdout` is the file where the output of a process is written.

`stderr` is the file where error output is written.

In an interactive session, all three files are connected to the user terminal. Using input or output redirection then means that the input is taken or the output sent to a different file than the terminal.

### 19.3.2.4 Command redirection

Instead of taking input from a file, or sending output to a file, it is possible to connect two commands together, so that the second takes the output of the first as input. The syntax for this is `cmdone | cmdtwo`; this is called a *pipeline*. For instance, `grep a yourfile | grep b` finds all lines that contain both an *a* and a *b*.

**Exercise.** Construct a pipeline that counts how many lines there are in your file that contain the string `th`. Use the `wc` command (see above) to do the counting.

There are a few more ways to combine commands. Suppose you want to present the result of `wc` a bit nicely. Type the following command

`echo The line count is wc -l foo`

where `foo` is the name of an existing file. The way to get the actual line count echoed is by the backquote:
```
echo The line count is `wc -l foo`
```

Anything in between backquotes is executed before the rest of the command line is evaluated. The way `wc` is used here, it prints the file name. Can you find a way to prevent that from happening?

### 19.3.3 Processes

The Unix operating system can run many programs at the same time, by rotating through the list and giving each only a fraction of a second to run each time. The command `ps` can tell you everything that is currently running.

**Exercise.** Type `ps`. How many programs are currently running? By default `ps` gives you only programs that you explicitly started. Do `ps guwax` for a detailed list of everything that is running. How many programs are running? How many belong to the root user, how many to you?

*Expected outcome.* To count the programs belonging to a user, pipe the `ps` command through an appropriate `grep`, which can then be piped to `wc`.

In this long listing of `ps`, the second column contains the process numbers. Sometimes it is useful to have those. The `cut` command explained above can cut certain position from a line: type `ps guwax | cut -c 10-14`.

To get dynamic information about all running processes, use the `top` command. Read the man page to find out how to sort the output by CPU usage.

When you type a command and hit return, that command becomes, for the duration of its run, the *foreground process*. Everything else that is running at the same time is a *background process*.

Make an executable file `hello` with the following contents:

```sh
#!/bin/sh
while [ 1 ]; do
    sleep 2
    date
done
```

and type `./hello`.

**Exercise.** Type `Control-z`. This suspends the foreground process. It will give you a number like `[1]` or `[2]` indicating that it is the first or second program that has been suspended or put in the background. Now type `bg` to put this process in the background. Confirm that there is no foreground process by hitting return, and doing an `ls`.

*Expected outcome.* After you put a process in the background, the terminal is available again to accept foreground commands. If you hit return, you should see the command prompt. However, the background process still keeps generating output.
**Exercise.** Type `jobs` to see the processes in the current session. If the process you just put in the background was number 1, type `fg %1`. Confirm that it is a foreground process again.

*Expected outcome.* If a shell is executing a program in the foreground, it will not accept command input, so hitting return should only produce blank lines.

**Exercise.** When you have made the `hello` script a foreground process again, you can kill it with `Control-c`. Try this. Start the script up again, this time as `./hello &` which immediately puts it in the background. You should also get output along the lines of `[1] 12345` which tells you that it is the first job you put in the background, and that `12345` is its process ID. Kill the script with `kill %1`. Start it up again, and kill it by using the process number.

*Expected outcome.* The command `kill 12345` using the process number is usually enough to kill a running program. Sometimes it is necessary to use `kill -9 12345`.

### 19.3.4 Shell customization

Above it was mentioned that `ls -F` is an easy way to see which files are regular, executable, or directories; by typing `alias ls='ls -F'` the `ls` command will automatically expanded to `ls -F` every time it is invoked. If you would like this behaviour in every login session, you can add the `alias` command to your `.profile` file. Other shells than `sh/bash` have other files for such customizations.

## 19.4 Scripting

The unix shells are also programming environments. You will learn more about this aspect of unix in this section.

### 19.4.1 Shell environment variables

Above you encountered `PATH`, which is an example of an shell, or environment, variable. These are variables that are known to the shell and that can be used by all programs run by the shell. You can see the full list of all variables known to the shell by typing `env`.

You can get the value of a shell variable by prefixing it with a dollar sign. Type the following two commands and compare the output:

```
echo PATH
echo $PATH
```

**Exercise.** Check on the value of the `HOME` variable by typing `echo $HOME`. Also find the value of `HOME` by piping `env` through `grep`.

Environment variables can be set in a number of ways. The simplest is by an assignment as in other programming languages.

**Exercise.** Type `a=5` on the commandline. This defines a variable `a`; check on its value by using the `echo` command.
*Expected outcome.* The shell will respond by typing the value 5.

*Caveats.* Beware not to have space around the equals sign; also be sure to use the dollar sign to print the value.

A variable set this way will be known to all subsequent commands you issue in this shell, but not to commands in new shells you start up. For that you need the *export* command. Reproduce the following session (the double percent sign is the command prompt):

```
%% a=20
%% echo $a
20
%% /bin/bash
%% echo $a

%% exit
exit
%% export a=21
%% /bin/bash
%% echo $a
21
%% exit
```

You can also temporarily set a variable. Replay this scenario:

1. Find an environment variable that does not have a value:

```
%% echo $b
```

```
%%
```

2. Write a short shell script to print this variable:

```
%% cat > echob
#!/bin/bash
echo $b
```

and of course make it execuable: `chmod +x echob`.

3. Now call the script, preceeding it with a setting of the variable b:

```
%% b=5 ./echob
5
```

The syntax where you set the value, as a prefix without using a separate command, sets the value just for that one command.

4. Show that the variable is still undefined:

```
%% echo $b
```

```
%%
```
That is, you defined the variable just for the execution of a single command.

In section 19.4.2 you will see that the `for` construct also defines a variable; section 19.4.3 shows some more built-in variables that apply in shell scripts.

## 19.4.2 Control structures

Like any good programming system, the shell has some control structures. Their syntax takes a bit of getting used to. (Different shells have different syntax; in this tutorial we only discuss the `bash` shell.)

In the `bash` shell, control structures can be written over several lines:

```bash
if [ $PATH = "" ] ; then
    echo "Error: path is empty"
fi
```

or on a single line:

```bash
if [ `wc -l file` -gt 100 ] ; then echo "file too long" ; fi
```

There are a number of tests defined, for instance `-f somefile` tests for the existence of a file. Change your script so that it will report `-1` if the file does not exist.

There are also loops. A `for` loop looks like

```perl
for var in listofitems ; do
    something with $var
done
```

This does the following:

* for each item in `listofitems`, the variable `var` is set to the item, and

* the loop body is executed.

As a simple example:

```bash
%% for x in a b c ; do echo $x ; done
a
b
c
```

In a more meaningful example, here is how you would make backups of all your `.c` files:

```bash
for cfile in *.c ; do
    cp $cfile $cfile.bak
done
```

Shell variables can be manipulated in a number of ways. Execute the following commands to see that you can remove trailing characters from a variable:
```bash
%% a=b.c
%% echo ${a%.c}
b
```

With this as a hint, write a loop that renames all your `.c` files to `.x` files.

### 19.4.3 Scripting

It is possible to write programs of unix shell commands. First you need to know how to put a program in a file and have it be executed. Make a file `script1` containing the following two lines:

```bash
#!/bin/bash
echo "hello world"
```

and type `./script1` on the command line. Result? Make the file executable and try again.

You can give your script command line arguments. If you want to be able to call

`./script1 foo bar`

you can use variables `$1`, `$2` et cetera in the script:

```bash
#!/bin/bash
echo "The first argument is $1"
echo "There were $# arguments in all"
```

Write a script that takes as input a file name argument, and reports how many lines are in that file.

Edit your script to test whether the file has less than 10 lines (use the `foo -lt bar` test), and if it does, `cat` the file. Hint: you need to use backquotes inside the test.

The number of command line arguments is available as `$#`. Add a test to your script so that it will give a helpful message if you call it without any arguments.

## 19.5 Expansion

The shell performs various kinds of expansion on a command line, that is, replacing part of the command-line with different text.

Brace expansion:

```bash
%% echo a{b,cc,ddd}e
abe acce addde
```

This can for instance be used to delete all extension of some base file name:
```bash
%% rm tmp.{c,s,o} # delete tmp.c tmp.s tmp.o
```

Tilde expansion gives your own, or someone else's home directory:

```bash
%% echo ~
/share/home/00434/eijkhout
%% echo ~eijkhout
/share/home/00434/eijkhout
```

Parameter expansion gives the value of shell variables:

```bash
%% x=5
%% echo $x
5
```

Undefined variables do not give an error message:

```bash
%% echo $y
```

There are many variations on parameter expansion. Above you already saw that you can strip trailing characters:

```bash
%% a=b.c
%% echo ${a%.c}
b
```

Here is how you can deal with undefined variables:

```bash
%% echo ${y:-0}
0
```

The backquote mechanism (section 19.3.2 above) is known as command substitution. It allows you to evalute part of a command and use it as input for another. For example, if you want to ask what type of file the command `ls` is, do

```bash
%% file 'which ls'
```

This first evalutes `which ls`, giving `/bin/ls`, and then evaluates `file /bin/ls`. As another example, here we backquote a whole pipeline, and do a test on the result:

```bash
%% echo 123 > w
%% cat w
123
%% wc -c w
4 w
%% if [ `cat w | wc -c` -eq 4 ] ; then echo four ; fi
four
```
Unix shell programming is very much oriented towards text manipulation, but it is possible to do arithmetic. Arithmetic substitution tells the shell to treat the expansion of a parameter as a number:

```
%% x=1
%% echo $((x*2))
2
```

Integer ranges can be used as follows:

```
%% for i in {1..10} ; do echo $i ; done
1
2
3
4
5
6
7
8
9
10
```

## 19.6 Startup files

In this tutorial you have seen several mechanisms for customizing the behaviour of your shell. For instance, by setting the `PATH` variable you can extend the locations where the shell looks for executables. Other environment variables (section 19.4.1) you can introduce for your own purposes. Many of these customizations will need to apply to every sessions, so you can have shell startup files that will be read at the start of any session.

Unfortunately, there are several startup files, and which one gets read is a complicated functions of circumstances. Here is a good common sense guideline<fnref n="3" />:

* Have a `.profile` that does nothing but read the `.bashrc`:
  ```
  # ~/.profile
  if [ -f ~/.bashrc ]; then
    source ~/.bashrc
  fi
  ```

* Your `.bashrc` does the actual customizations:
  ```
  # ~/.bashrc
  # make sure your path is updated
  if [ -z ``$MYPATH'' ]; then
    export MYPATH=1
  fi
  ```

<footnote n="3">3. Many thanks for Robert McLay for figuring this out.</footnote>
```
export PATH=$HOME/bin:$PATH
fi
```

## 19.7 Shell interaction

Interactive use of Unix, in contrast to script writing (section 19.4), is a complicated conversation between the user and the shell. You, the user, type a line, hit return, and the shell tries to interpret it. There are several cases.

* Your line contains one full command, such as `ls foo`: the shell will execute this command.

* You can put more than one command on a line, separated by semicolons: `mkdir foo; cd foo`. The shell will execute these commands in sequence.

* Your input line is not a full command, for instance `while [ 1 ]`. The shell will recognize that there is more to come, and use a different prompt to show you that it is waiting for the remainder of the command.

* Your input line would be a legitimate command, but you want to type more on a second line. In that case you can end your input line with a backslash character, and the shell will recognize that it needs to hold off on executing your command. In effect, the backslash will hide (escape) the return.

When the shell has collected a command line to execute, by using one or more of your input line or only part of one, as described just now, it will apply expansion to the command line (section 19.5). It will then interpret the commandline as a command and arguments, and proceed to invoke that command with the arguments as found.

There are some subtleties here. If you type `ls *.c`, then the shell will recognize the wildcard character and expand it to a command line, for instance `ls foo.c bar.c`. Then it will invoke the `ls` command with the argument list `foo.c bar.c`. Note that `ls` does not receive `*.c` as argument! In cases where you do want the unix command to receive an argument with a wildcard, you need to escape it so that the shell will not expand it. For instance, `find . -name \*.c` will make the shell invoke `find` with arguments `.`.

`-name *.c`.

## 19.8 The system and other users

Unix is a multi-user operating system. Thus, even if you use it on your own personal machine, you are a user with an *account* and you may occasionally have to type in your username and password.

If you are on your personal machine, you may be the only user logged in. On university machines or other servers, there will often be other users. Here are some commands relating to them.

* `whoami` show your login name.

* `who` show the other users currently logged in.

* `finger otheruser` get information about another user; you can specify a user's login name here, or their real name, or other identifying information the system knows about.
`top` which processes are running on the system; use `top -u` to get this sorted the amount of cpu
time they are currently taking. (On Linux, try also the `vmstat` command.)

`uptime` how long has it been since your last reboot?

## 19.9 The `sed` and `awk` tools

Apart from fairly small utilities such as `tr` and `cut`, Unix some more powerful ones. In this section you will see two tools for line-by-line transformations on text files. Of course this tutorial merely touches on the depth of these tools; for more information see [1, 46].

### 19.9.1 `sed`

The streaming editor `sed` is like an editor by remote control, doing simple line edits with a commandline interface. Most of the time you will use `sed` as follows:

`cat somefile | sed 's/abc/def/:g' > newfile`

(The use of `cat` here is not strictly necessary.) The `s/abc/def/` part has the effect of replacing `abc` by `def` in every line; the `:g` modifier applies it to every instance in every line rather than just the first.

* If you have more than one edit, you can specify them with

`sed -e 's/one/two/' -e 's/three/four/'`

* If an edit needs to be done only on certain lines, you can specify that by prefixing the edit with the match string. For instance

`sed '/^a/s/b/c/'`

only applies the edit on lines that start with an `a`. (See section 19.2 for regular expressions.)

* Traditionally, `sed` could only function in a stream, so you the output file was always different from the input. The GNU version, which is standard on Linux systems, has a flag `-i` which edits 'in place':

`sed -e 's/ab/cd/' -e 's/ef/gh/' -i thefile`

### 19.9.2 `awk`

The `awk` utility also operates on each line, but it can be described as having a memory. An `awk` program consists of a sequence of pairs, where each pair consists of a match string and an action. The simplest `awk` program is

`cat somefile | awk '{ print }'`

where the match string is omitted, meaning that all lines match, and the action is to print the line. Awk breaks each line into fields separated by whitespace. A common application of `awk` is to print a certain field:
awk '{print $2}' file

prints the second field of each line.

Suppose you want to print all subroutines in a Fortran program; this can be accomplished with

awk '/subroutine/ {print}' yourfile.f

Exercise 19.1. Build a commandpipeline that prints of each subroutine header only the sub-
routine name. For this you first use sed to replace the parentheses by spaces, then
awk to print the subroutine name field.

Awk has variables with which it can remember things. For instance, instead of just printing the second field of every line, you can make a list of them and print that later:

cat myfile | awk 'BEGIN {v="Fields:"} {v=v " " $2} END {print v}'

As another example of the use of variables, here is how you would print all lines in between a `BEGIN` and `END` line:

cat myfile | awk '/END/ {p=0} p==1 {print} /BEGIN/ {p=1} '

Exercise 19.2. The placement of the match with `BEGIN` and `END` may seem strange. Rear-range the awk program, test it out, and explain the results you get.

**19.10 Review questions**

Exercise 19.3. Write a shell script for making backups. When you call this script as `./backup somefile` it should test whether `somefile.bak` exists, and give a warning if it does. In either case, it should copy the original file to a backup.
# Chapter 20

## Compilers and libraries

### 20.1 An introduction to binary files

**Purpose.** In this section you will be introduced to the different types of binary files that you encounter while programming.

One of the first things you become aware of when you start programming is the distinction between the readable source code, and the unreadable, but executable, program code. In this tutorial you will learn about a couple more file types:

* A source file can be compiled to an *object file*, which is a bit like a piece of an executable: by itself it does nothing, but it can be combined with other object files to form an executable.

* A *library* is a bundle of object files that can be used to form an executable. Often, libraries are written by an expert and contain code for specialized purposes such as linear algebra manipulations. Libraries are important enough that they can be commercial, to be bought if you need expert code for a certain purpose.

You will now learn how these types of files are created and used.

### 20.2 Simple compilation

**Purpose.** In this section you will learn about executables and object files.

Let's start with a simple program that has the whole source in one file.

One file: `hello.c`

```c
#include <stdlib.h>
#include <stdio.h>

int main() {
    printf("hello world\n");
    return 0;
}
```
Compile this program with your favourite compiler; we will use `gcc` in this tutorial, but substitute your own as desired. As a result of the compilation, a file `a.out` is created, which is the executable.

```
%% gcc hello.c
%% ./a.out
hello world
```

You can get a more sensible program name with the `-o` option:

```
%% gcc -o helloprog hello.c
%% ./helloprog
hello world
```

Now we move on to a program that is in more than one source file.

### Main program: fooprog.c

```c
#include <stdlib.h>
#include <stdio.h>

extern bar(char*);

int main() {
    bar("hello world\n");
    return 0;
}
```

### Subprogram: fooprog.c

```c
#include <stdlib.h>
#include <stdio.h>

void bar(char *s) {
    printf("%s",s);
    return;
}
```

As before, you can make the program with one command.

```
%% gcc -o foo fooprog.c foosub.c
%% ./foo
hello world
```

However, you can also do it in steps, compiling each file separately and then linking them together.

```
%% gcc -c fooprog.c
%% gcc -c foosub.c
%% gcc -o foo fooprog.o foosub.o
%% ./foo
```
hello world

The `-c` option tells the compiler to compile the source file, giving an *object file*. The third command than acts as the *linker*, tieing together the object files into an executable. (With programs that are spread over several files there is always the danger of editing a subroutine definition and then forgetting to update all the places it is used. See the 'make' tutorial, section 21, for a way of dealing with this.)

## 20.3 Libraries

**Purpose.** In this section you will learn about libraries.

If you have written some subprograms, and you want to share them with other people (perhaps by selling them), then handing over individual object files is inconvenient. Instead, the solution is to combine them into a library. First we look at *static libraries*, for which the *archive utility* `ar` is used. A static library is linked into your executable, becoming part of it. This may lead to large executables; you will learn about shared libraries next, which do not suffer from this problem.

Create a directory to contain your library (depending on what your library is for this can be a system directory such as `/usr/bin`), and create the library file there.

```bash
%% mkdir ../lib
%% ar cr ../lib/libfoo.a foosub.o
```

The `nm` command tells you what's in the library:

```bash
%% nm ../lib/libfoo.a
../lib/libfoo.a(foosub.o):
00000000 T _bar
U _printf
```

Line with `T` indicate functions defined in the library file; a `U` indicates a function that is used.

The library can be linked into your executable by explicitly giving its name, or by specifying a library path:

```bash
%% gcc -o foo fooprog.o ../lib/libfoo.a
# or
%% gcc -o foo fooprog.o -L../lib -lfoo
%% ./foo
hello world
```

A third possibility is to use the `LD_LIBRARY_PATH` shell variable. Read the man page of your compiler for its use, and give the commandlines that create the `foo` executable, linking the library through this path.

Although they are somewhat more complicated to use, *shared libraries* have several advantages. For instance, since they are not linked into the executable but only loaded at runtime, they lead to (much) smaller executables. They are not created with `ar`, but through the compiler. For instance:
```
%% gcc -dynamiclib -o ../lib/libfoo.so foosub.o
%% nm ../lib/libfoo.so
```

```
../lib/libfoo.so(single module):
00000fc4 t __dyld_func_lookup
00000000 t __mh Dylib_header
00000fd2 T _bar
U _printf
00001000 d dyld__mach_header
00000fb0 t dyldstub_bindinG_helper
```

Shared libraries are not actually linked into the executable; instead, the executable will contain the information where the library is to be found at execution time:

```
%% gcc -o foo fooprog.o -L../lib -Wl,-rpath,`pwd`/../lib -lfoo
%% ./foo
hello world
```

The link line now contains the library path twice:

1. once with the `-L` directive so that the linker can resolve all references, and

2. once with the linker directive `-Wl, -rpath, `pwd`/../lib` which stores the path into the executable so that it can be found at runtime.

Build the executable again, but without the `-Wl` directive. Where do things go wrong and why? You can also fix this problem by using `LD_LIBRARY_PATH`. Explore this.

Use the command `ldd` to get information about what shared libraries your executable uses. (On Mac OS X, use `otool -L` instead.)
# Chapter 21

## Managing projects with Make

The *Make* utility helps you manage the building of projects: its main task is to facilitate rebuilding only those parts of a multi-file project that need to be recompiled or rebuilt. This can save lots of time, since it can replace a minutes-long full installation by a single file compilation. *Make* can also help maintaining multiple installations of a program on a single machine, for instance compiling a library with more than one compiler, or compiling a program in debug and optimized mode.

*Make* is a Unix utility with a long history, and traditionally there are variants with slightly different behaviour, for instance on the various flavours of Unix such as HP-UX, AUX, IRIX. These days, it is advisable, no matter the platform, to use the GNU version of *Make* which has some very powerful extensions; it is available on all Unix platforms (on Linux it is the only available variant), and it is a *de facto* standard. The manual is available at http://www.gnu.org/software/make/manual/make.html, or you can read the book [121].

There are other build systems, most notably *Scons* and *Bjam*. We will not discuss those here. The examples in this tutorial will be for the C and Fortran languages, but *Make* can work with any language, and in fact with things like TeX that are not really a language at all; see section 21.6.

### 21.1 A simple example

**Purpose.** In this section you will see a simple example, just to give the flavour of *Make*.

The files for this section can be downloaded from http://tinyurl.com/ISTC-make-tutorial.

#### 21.1.1 C

Make the following files:

foo.c
```c
#include "bar.h"
int c=3;
int d=4;
int main()
{
    int a=2;
    return (bar(a*c*d));
}
```

### bar.c

```c
#include "bar.h"
int bar(int a)
{
    int b=10;
    return (b*a);
}
```

### bar.h

```c
extern int bar(int);
```

and a makefile:

### Makefile

```
fooprog : foo.o bar.o
        cc -o fooprog foo.o bar.o

foo.o : foo.c
        cc -c foo.c

bar.o : bar.c
        cc -c bar.c

clean :
    rm -f *.o fooprog
```

The makefile has a number of rules like

```
foo.o : foo.c
<TAB>cc -c foo.c
```

which have the general form

```
target : prerequisite(s)
<TAB>rule(s)
```

where the rule lines are indented by a `TAB` character.

A rule, such as above, states that a 'target' file `foo.o` is made from a 'prerequisite' `foo.c`, namely by executing the command `cc -c foo.c`. The precise definition of the rule is:
* if the target `foo.o` does not exist or is older than the prerequisite `foo.c`,

* then the command part of the rule is executed: `cc -c foo.c`

* If the prerequisite is itself the target of another rule, then that rule is executed first.

Probably the best way to interpret a rule is:

* if any prerequisite has changed,

* then the target needs to be remade,

* and that is done by executing the commands of the rule;

* checking the prerequisite requires a recursive application of make:

- if the prerequisite does not exist, find a rule to create it;

- if the prerequisite already exists, check applicable rules to see if it needs to be remade.

If you call `make` without any arguments, the first rule in the makefile is evaluated. You can execute other rules by explicitly invoking them, for instance `make foo.o` to compile a single file.

**Exercise.** Call `make`.

*Expected outcome.* The above rules are applied: `make` without arguments tries to build the first target, `fooprog`. In order to build this, it needs the prerequisites `foo.o` and `bar.o`, which do not exist. However, there are rules for making them, which `make` recursively invokes. Hence you see two compilations, for `foo.o` and `bar.o`, and a link command for `fooprog`.

*Caveats.* Typos in the makefile or in file names can cause various errors. In particular, `make` sure you use tabs and not spaces for the rule lines. Unfortunately, debugging a makefile is not simple. *Make's* error message will usually give you the line number in the make file where the error was detected.

**Exercise.** Do `make clean`, followed by `mv foo.c boo.c` and `make` again. Explain the error message. Restore the original file name.

*Expected outcome.* *Make* will complain that there is no rule to make `foo.c`. This error was caused when `foo.c` was a prerequisite for making `foo.o`, and was found not to exist. *Make* then went looking for a rule to make it and no rule for creating `.c` files exists.

Now add a second argument to the function `bar`. This requires you to edit `bar.c` and `bar.h`: go ahead and make these edits. However, it also requires you to edit `foo.c`, but let us for now 'forget' to do that. We will see how *Make* can help you find the resulting error.

**Exercise.** Call `make` to recompile your program. Did it recompile `foo.c`?

*Expected outcome.* Even through conceptually `foo.c` would need to be recompiled since it uses the `bar` function, *Make* did not do so because the makefile had no rule that forced it.

In the makefile, change the line

`foo.o : foo.c`

to
foo.o : foo.c bar.h

which adds `bar.h` as a prerequisite for `foo.o`. This means that, in this case where `foo.o` already exists, `Make` will check that `foo.o` is not older than any of its prerequisites. Since `bar.h` has been edited, it is younger than `foo.o`, so `foo.o` needs to be reconstructed.

**Exercise.** Confirm that the new makefile indeed causes `foo.o` to be recompiled if `bar.h` is changed. This compilation will now give an error, since you 'forgot' to edit the use of the `bar` function.

### 21.1.2 Fortran

Make the following files:

```
foomain.F
```

```
program test
use testmod

call func(1,2)

end program
```

```
foomod.F
```

```
module testmod

contains

subroutine func(a,b)
integer a,b
print *,a,b,c
end subroutine func

end module
```

and a makefile:

```
Makefile
```

```
fooprog : fomain.o foomod.o
gfortran -o fooprog foo.o foomod.o
foomain.o : fomain.F
gfortran -c fomain.F
foomod.o : foomod.F
gfortran -c foomod.F
clean :
rm -f *.o fooprog
```
If you call `make`, the first rule in the makefile is executed. Do this, and explain what happens.

**Exercise.** Call `make`.

*Expected outcome.* The above rules are applied: `make` without arguments tries to build the first target, `foomain`. In order to build this, it needs the prerequisites `foomain.o` and `foomod.o`, which do not exist. However, there are rules for making them, which `make` recursively invokes. Hence you see two compilations, for `foomain.o` and `foomod.o`, and a link command for `fooprog`.

*Caveats.* Typos in the makefile or in file names can cause various errors. Unfortunately, debugging a makefile is not simple. You will just have to understand the errors, and make the corrections.

**Exercise.** Do `make clean`, followed by `mv foomod.c boomod.c` and `make` again. Explain the error message. Restore the original file name.

*Expected outcome.* `Make` will complain that there is no rule to make `foomod.c`. This error was caused when `foomod.c` was a prerequisite for `foomod.o`, and was found not to exist. `Make` then went looking for a rule to make it, and no rule for making `.F` files exists.

Now add an extra parameter to `func` in `foomod.F` and recompile.

**Exercise.** Call `make` to recompile your program. Did it recompile `foomain.F`?

*Expected outcome.* Even through conceptually `foomain.F` would need to be recompiled, `Make` did not do so because the makefile had no rule that forced it.

Change the line

`foomain.o : foomain.F`

to

`foomain.o : foomain.F foomod.o`

which adds `foomod.o` as a prerequisite for `foomain.o`. This means that, in this case where `foomain.o` already exists, `Make` will check that `foomain.o` is not older than any of its prerequisites. Recursively, `Make` will then check if `foomode.o` needs to be updated, which is indeed the case. After recompiling `foomode.F`, `foomode.o` is younger than `foomain.o`, so `foomain.o` will be reconstructed.

**Exercise.** Confirm that the corrected makefile indeed causes `foomain.F` to be recompiled.

### 21.1.3 About the make file

The make file needs to be called `makefile` or `Makefile`; it is not a good idea to have files with both names in the same directory. If you want `Make` to use a different file as make file, use the syntax `make -f My_Makefile`.
## 21.2 Variables and template rules

**Purpose.** In this section you will learn various work-saving mechanisms in *Make*, such as the use of variables and of template rules.

### 21.2.1 Makefile variables

It is convenient to introduce variables in your makefile. For instance, instead of spelling out the compiler explicitly every time, introduce a variable in the makefile:

```
CC = gcc
FC = gfortran
```

and use `${CC}` or `${FC}` on the compile lines:

```
foo.o : foo.c
       ${CC} -c foo.c

foomain.o : foomain.F
           ${FC} -c foomain.F
```

**Exercise.** Edit your makefile as indicated. First do `make clean`, then `make foo (C)` or `make fooprog (Fortran)`.

**Expected outcome.** You should see the exact same compile and link lines as before.

**Caveats.** Unlike in the shell, where braces are optional, variable names in a makefile have to be in braces or parentheses. Experiment with what happens if you forget the braces around a variable name.

One advantage of using variables is that you can now change the compiler from the commandline:

```
make CC="icc -O2"
make FC="gfortran -g"
```

**Exercise.** Invoke *Make* as suggested (after `make clean`). Do you see the difference in your screen output?

**Expected outcome.** The compile lines now show the added compiler option `-O2` or `-g`.

*Make* also has *automatic variables*:

* `$@` The target. Use this in the link line for the main program.

* `$^` The list of prerequisites. Use this also in the link line for the program.

* `$<` The first prerequisite. Use this in the compile commands for the individual object files.

* `$*` The stem of the target by stripping any extension.

Using these variables, the rule for `fooprog` becomes
```
fooprog : foo.o bar.o
${CC} -o $@ $^
```

and a typical compile line becomes

```
foo.o : foo.c bar.h
${CC} -c $<
```

You can also declare a variable

```
THEPROGRAM = fooprog
```

and use this variable instead of the program name in your makefile. This makes it easier to change your mind about the name of the executable later.

**Exercise.** Edit your makefile to add this variable definition, and use it instead of the literal program name. Construct a commandline so that your makefile will build the executable `fooprog_v2`.

*Expected outcome.* You need to specify the `THEPROGRAM` variable on the commandline using the syntax `make VAR=value`.

*Caveats.* Make sure that there are no spaces around the equals sign in your commandline.

### 21.2.2 Template rules

So far, you wrote a separate rule for each file that needed to be compiled. However, the rules for the various `.c` files are very similar:

* the rule header (`foo.o : foo.c`) states that a source file is a prerequisite for the object file with the same base name;

* and the instructions for compiling (`${CC} -c $<`) are even character-for-character the same, now that you are using `Make`'s built-in variables;

* the only rule with a difference is

```
foo.o : foo.c bar.h
${CC} -c $<
```

where the object file depends on the source file and another file.

We can take the commonalities and summarize them in one template rule<fnref n="1" />:

```
%.o : %.c
${CC} -c $<
```

```
%.o : %.F
${FC} -c $<
```

<footnote n="1">1. This mechanism is the first instance you'll see that only exists in GNU make, though in this particular case there is a similar mechanism in standard make. That will not be the case for the wildcard mechanism in the next section.</footnote>
This states that any object file depends on the C or Fortran file with the same base name. To regenerate the object file, invoke the C or Fortran compiler with the `-c` flag. These template rules can function as a replacement for the multiple specific targets in the makefiles above, except for the rule for `foo.o`.

The dependence of `foo.o` on `bar.h`, or `foomain.o` on `foomod.o`, can be handled by adding a rule

```
# C
foo.o : bar.h

# Fortran
foomain.o : foomod.o
```

with no further instructions. This rule states, 'if file `bar.h` or `foomod.o` changed, file `foo.o` or `foomain.o` needs updating' too. `Make` will then search the makefile for a different rule that states how this updating is done, and it will find the template rule.

**Exercise.** Change your makefile to incorporate these ideas, and test.

### 21.2.3 Wildcards

Your makefile now uses one general rule for compiling any source file. Often, your source files will be all the `.c` or `.F` files in your directory, so is there a way to state 'compile everything in this directory'? Indeed there is.

Add the following lines to your makefile, and use the variable `COBJECTS` or `FOBJECTS` wherever appropriate. The command `wildcard` gives the result of `ls`, and you can manipulate file names with `patsubst`.

```
# wildcard: find all files that match a pattern
CSOURCES := ${wildcard *.c}

# pattern substitution: replace one pattern string by another
COBJECTS := ${patsubst %.*c,%.o,${SRC}}

FSOURCES := ${wildcard *.F}
FOBJECTS := ${patsubst %.*F,%.o,${SRC}}
```

### 21.2.4 Conditionals

There are various ways of making the behaviour of a makefile dynamic. You can for instance put a shell conditional in an action line. However, this can make for a cluttered makefile; an easier way is to use makefile conditionals. There are two types of conditionals: tests on string equality, and tests on environment variables.

The first type looks like
```c
ifeq "${HOME}" "/home/thisisme"
# case where the executing user is me
else
# case where it's someone else
endif
```

and in the second case the test looks like

```c
ifdef SOME_VARIABLE
```

The text in the true and false part can be most any part of a makefile. For instance, it is possible to let one of the action lines in a rule be conditionally included. However, most of the times you will use conditionals to make the definition of variables dependent on some condition.

**Exercise.** Let's say you want to use your makefile at home and at work. At work, your employer has a paid license to the Intel compiler `icc`, but at home you use the open source Gnu compiler `gcc`. Write a makefile that will work in both places, setting the appropriate value for `CC`.

## 21.3 Miscellania

### 21.3.1 What does this makefile do?

Above you learned that issuing the `make` command will automatically execute the first rule in the makefile. This is convenient in one sense<fnref n="2" />, and inconvenient in another: the only way to find out what possible actions a makefile allows is to read the makefile itself, or the – usually insufficient – documentation.

A better idea is to start the makefile with a target

```make
info :
    @echo "The following are possible:"
    @echo " make"
    @echo " make clean"
```

Now `make` without explicit targets informs you of the capabilities of the makefile.

If your makefile gets longer, you might want to document each section like this. This runs into a problem: you can not have two rules with the same target, `info` in this case. However, if you use a double colon it is possible. Your makefile would have the following structure:

```make
info ::
    @echo "The following target are available:"
    @echo " make install"
```

<footnote n="2">2. There is a convention among software developers that a package can be installed by the sequence `./configure`; `make`; `make install`, meaning: Configure the build process for this computer, Do the actual build, Copy files to some system directory such as `/usr/bin`.</footnote>
```make
install :
    # ..... instructions for installing

info ::
    @echo " make clean"

clean :
    # ..... instructions for cleaning
```

### 21.3.2 Phony targets

The example makefile contained a target `clean`. This uses the `Make` mechanisms to accomplish some actions that are not related to file creation: calling `make clean` causes `Make` to reason 'there is no file called `clean`, so the following instructions need to be performed'. However, this does not actually cause a file `clean` to spring into being, so calling `make clean` again will make the same instructions being executed.

To indicate that this rule does not actually make the target, declare

`.PHONY : clean`

One benefit of declaring a target to be phony, is that the `Make` rule will still work, even if you have a file named `clean`.

### 21.3.3 Predefined variables and rules

Calling `make -p yourtarget` causes `make` to print out all its actions, as well as the values of all variables and rules, both in your makefile and ones that are predefined. If you do this in a directory where there is no makefile, you'll see that `make` actually already knows how to compile `.c` or `.F` files. Find this rule and find the definition of the variables in it.

You see that you can customize `make` by setting such variables as `CFLAGS` or `FFLAGS`. Confirm this with some experimentation. If you want to make a second makefile for the same sources, you can call `make -f othermakefile` to use this instead of the default `Makefile`.

Note, by the way, that both `makefile` and `Makefile` are legitimate names for the default makefile. It is not a good idea to have both `makefile` and `Makefile` in your directory.

### 21.3.4 Using the target as prerequisite

Suppose you have two different targets that are treated largely the same. You would want to write:

```make
PROGS = myfoo other
${PROGS} : $@.o # this is wrong!!
${CC} -o $@ $@.o ${list of libraries goes here}
```

and saying `make myfoo` would cause
```
cc -c myfoo.c
cc -o myfoo myfoo.o ${list of libraries}
```

and likewise for `make other`. What goes wrong here is the use of `$@.o` as prerequisite. In Gnu Make, you can repair this as follows<fnref n="3" />:

```
.SECONDEXPANSION:
${PROGS} : $$@.o
${CC} -o $@ $@.o ${list of libraries goes here}
```

**Exercise.** Write a second main program `foosecond.c` or `foosecond.F`, and change your makefile so that the calls `make foo` and `make foosecond` both use the same rule.

## 21.4 Shell scripting in a Makefile

**Purpose.** In this section you will see an example of a longer shell script appearing in a makefile rule.

In the makefiles you have seen so far, the command part was a single line. You can actually have as many lines there as you want. For example, let us make a rule for making backups of the program you are building.

Add a `backup` rule to your makefile. The first thing it needs to do is make a backup directory:

```
.PHONY : backup
backup :
    if [ ! -d backup ] ; then
        mkdir backup
    fi
```

Did you type this? Unfortunately it does not work: every line in the command part of a makefile rule gets executed as a single program. Therefore, you need to write the whole command on one line:

```
backup :
    if [ ! -d backup ] ; then mkdir backup ; fi
```

or if the line gets too long:

```
backup :
    if [ ! -d backup ] ; then \  
        mkdir backup ; \ 
    fi
```

<footnote n="3">3. Technical explanation: Make will now look at lines twice: the first time `$$` gets converted to a single `$`, and in the second pass `$@` becomes the name of the target.</footnote>
Next we do the actual copy:

```
backup :
    if [ ! -d backup ] ; then mkdir backup ; fi
cp myprog backup/myprog
```

But this backup scheme only saves one version. Let us make a version that has the date in the name of the saved program.

The Unix `date` command can customize its output by accepting a format string. Type the following: `date` This can be used in the makefile.

**Exercise.** Edit the `cp` command line so that the name of the backup file includes the current date.

**Expected outcome.** Hint: you need the backquote. Consult the Unix tutorial, section 19.3.2.4, if you do not remember what backquotes do.

If you are defining shell variables in the command section of a makefile rule, you need to be aware of the following. Extend your `backup` rule with a loop to copy the object files:

```
backup :
    if [ ! -d backup ] ; then mkdir backup ; fi
cp myprog backup/myprog
for f in ${OBJS} ; do \
cp $f backup ; \
done
```

(This is not the best way to copy, but we use it for the purpose of demonstration.) This leads to an error message, caused by the fact that *Make* interprets `$f` as an environment variable of the outer process. What works is:

```
backup :
    if [ ! -d backup ] ; then mkdir backup ; fi
cp myprog backup/myprog
for f in ${OBJS} ; do \
cp $$f backup ; \
done
```

(In this case *Make* replaces the double dollar by a single one when it scans the commandline. During the execution of the commandline, `$f` then expands to the proper filename.)

## 21.5 Practical tips for using Make

Here are a couple of practical tips.

* Debugging a makefile is often frustratingly hard. Just about the only tool is the `-p` option, which prints out all the rules that Make is using, based on the current makefile.
```
| wc -l` -gt 0 ] ; do \
pdflatex ${basename $@} ; \
done
```

The minus sign at the start of the line means that *Make* should not abort if these commands fail.

Finally, we would like to use *Make's* facility for taking dependencies into account. We could write a makefile that has the usual rules

```
mainfile.pdf : mainfile.tex includefile.tex
```

but we can also discover the include files explicitly. The following makefile is invoked with

```
make pdf TEXTFILE=mainfile
```

The `pdf` rule then uses some shell scripting to discover the include files (but not recursively), and it calls *Make* again, invoking another rule, and passing the dependencies explicitly.

```
pdf :
    export includes=`grep "^input" ${TEXFILE}.tex \
                      awk '{v=v FS $$2".tex"} END {print v}' ; \
                      ${MAKE} ${TEXFILE}.pdf INCLUDES="$$includes"

%.pdf : %.tex ${INCLUDES}
    pdflatex $< ; \
    while [ `cat ${basename $@}.log \
                      grep "Rerun to get" | wc -l` -gt 0 ] ; do \
        pdflatex $< ; \
    done
```

This shell scripting can also be done outside the makefile, generating the makefile dynamically.
# Chapter 22

## Source code control

Source code control systems, also called revision control or version control systems, are a way of storing software, where not only the current version is stored, but also all previous versions. This is done by maintaining a repository for all versions, while one or more users work on a 'checked out' copy of the latest version. Those of the users that are developers can then commit their changes to the repository. Other users then update their local copy. The repository typically resides on a remote machine that is reliably backup up.

There are various reasons for keeping your source in a repository.

* If you work in a team, it is the best way to synchronize your work with your colleagues. It is also a way to document what changes were made, by whom, and why.

* It will allow you to roll back a defective code to a version that worked.

* It allows you to have branches, for instance for customizations that need to be kept out of the main development line. If you are working in a team, a branch is a way to develop a major feature, stay up to date with changes your colleagues make, and only add your feature to the main development when it is sufficiently tested.

* If you work alone, it is a way to synchronize between more than one machine. (You could even imagine traveling without all your files, and installing them from the repository onto a borrowed machine as the need arises.)

* Having a source code repository is one way to backup your work.

There are various source code control systems; in this tutorial you can learn the basics of Subversion (also called `svn`), which is probably the most popular of the traditional source code control systems, and Mercurial (or `hg`), which is an example of the new generation of distributed source code control systems.

### 22.1 Workflow in source code control systems

Source code control systems are built around the notion of *repository*: a central store of the files of a project, together with their whole history. Thus, a repository allows you to share files with multiple people, but also to roll back changes, apply patches to old version, et cetera.

The basic actions on a repository are:
* Creating the repository; this requires you to have space and write permissions on some server. Maybe your sysadmin has to do it for you.

* Checking out the repository, that is, making a local copy of its contents in your own space.

* Adding your changes to the repository, and

* Updating your local copy with someone else's changes.

Adding your own changes is not always possible: there are many projects where the developer allows you to check out the repository, but not to incorporate changes. Such a repository is said to be read-only.

Figure 22.1 illustrates these actions for the Subversion system. Users who have checked out the repository

<figure><img src="image_138.png" /><figcaption>Figure 22.1: Workflow in traditional source code control systems such as Subversion</figcaption></figure>

can edit files, and check in the new versions with the `commit` command; to get the changes committed by other users you use `update`.

One of the uses of committing is that you can roll your code back to an earlier version if you realize you made a mistake or introduced a bug. It also allows you to easily see the difference between different code version. However, committing many small changes may be confusing to other developers, for instance if they come to rely on something you introduce which you later remove again. For this reason, *distributed source code control* systems use two levels of repositories.

There is still a top level that is authoritative, but now there is a lower level, typically of local copies, where you can commit your changes and accumulate them until you finally add them to the central repository. This also makes it easier to contribute to a read-only repository: you make your local changes, and when you are finished you tell the developer to inspect your changes and pull them into the top level repository. This structure is illustrated in figure 22.2.
<figure><img src="image_139.png" /><figcaption>Figure 22.2: Workflow in distributed source code control systems such as Mercurial</figcaption></figure>

## 22.2 Subversion or SVN

This lab should be done two people, to simulate a group of programmers working on a joint project. You can also do this on your own by using two copies of the repository.

### 22.2.1 Create and populate a repository

**Purpose.** In this section you will create a repository and make a local copy to work on.

First we need to have a repository. In practice, you will often use one that has been set up by a sysadmin, but there are several ways to set up a repository yourself.

* There are commercial and free hosting services such as Google code http://code.google.com/projecthosting (open source only) or BitBucket https://bitbucket.org/ (this has favourable academic licenses). Once you create a repository there, you can make a local copy on your computer:

```bash
%% svn co http://yourservice.com//yourproject/ project
Checked out revision 0.
```

where `co` is short for 'check out'<fnref n="1" />.

<footnote n="1">1. Alternatively, you can create a repository in the local file system with `svnadmin create ./repository --fs-type=fsfs` and check it out with `svn co file://path.../`; we will not go into this.</footnote>
You now have an empty directory `project`.

**Exercise.** Go into the project directory and see if it is really empty.

*Expected outcome.* There is a hidden directory `.svn`

The project directory knows where the master copy of the repository is stored.

**Exercise.** Do `svn info`

*Expected outcome.* You will see the URL of the repository and the revision number
*Caveats.* Make sure you are in the right directory. In the outer directories, you will get `svn: warning: '.' is not a working copy`

### 22.2.2 New files

**Purpose.** In this section you will make some simple changes: editing an existing file and creating a new file.

One student now makes a file to add to the repository<fnref n="2" />.

```bash
%% cat > firstfile
a
b
c
d
e
f
^D
```

This file is unknown to `svn`:

```bash
%% svn status
?          firstfile
```

We need to declare the file as belonging to the repository; a subsequent `svn commit` command then copies it into the master repository.

```bash
%% svn add firstfile
A          firstfile
%% svn status
A          firstfile
%% svn commit -m "made a first file"
Adding              firstfile
Transmitting file data .
Committed revision 1.
```

<footnote n="2">2. It is also possible to use `svn import` to create a whole repository from an existing directory tree.</footnote>
**Exercise.** The second student can now do `svn update` to update their copy of the repository

*Expected outcome.* `svn` should report `A firstfile` and `Updated to revision 1..` Check that the contents of the file are correct.

*Caveats.* In order to do the update command, you have to be in a checked-out copy of the repository. Do `svn info` to make sure that you are in the right place.

**Exercise.** Let both students create a new directory with a few files. Declare the directory and commit it. Do `svn update` to obtain the changes the other mde.

*Expected outcome.* You can do `svn add` on the directory, this will also add the files contained in it.

*Caveats.* Do not forget the `commit`.

In order for `svn` to keep track of your files, you should never do `cp` or `mv` on files that are in the repository. Instead, do `svn cp` or `svn mv`. Likewise, there are commands `svn rm` and `svn mkdir`.

### 22.2.3 Conflicts

**Purpose.** In this section you will learn about how do deal with conflicting edits by two users of the same repository.

Now let's see what happens when two people edit the same file. Let both students make an edit to `firstfile`, but one to the top, the other to the bottom. After one student commits the edit, the other will see

```
%% emacs firstfile # make some change
%% svn commit -m "another edit to the first file"
Sending firstfile
svn: Commit failed (details follow):
svn: Out of date: 'firstfile' in transaction '5-1'
```

The solution is to get the other edit, and commit again. After the update, `svn` reports that it has resolved a conflict successfully.

```
%% svn update
G firstfile
Updated to revision 5.
%% svn commit -m "another edit to the first file"
Sending firstfile
Transmitting file data.
Committed revision 6.
```

The `G` at the start of the line indicates that `svn` has resolved a conflicting edit.

If both students make edits on the same part of the file, `svn` can no longer resolve the conflicts. For instance, let one student insert a line between the first and the second, and let the second student edit the second line. Whoever tries to commit second, will get messages like this:
```
%% svn commit -m "another edit to the first file"

svn: Commit failed (details follow):

svn: Aborting commit: '/share/home/12345/yourname/myproject/firstfile'

remains in conflict

%% svn update

C firstfile

Updated to revision 7.
```

Subversion will give you several options. For instance, you can type `e` to open the file in an editor. You can also type `p` for 'postpone' and edit it later. Opening the file in an editor, it will look like:

```
aa
<<<<<<<<< .mine
bb
========
123
b
>>>>>>>> .r7
cc
```

indicating the difference between the local version ('mine') and the remote. You need to edit the file to resolve the conflict.

After this, you tell `svn` that the conflict was resolved, and you can commit:

```
%% svn resolved firstfile

Resolved conflicted state of 'firstfile'

%% svn commit -m "another edit to the first file"

Sending firstfile

Transmitting file data .

Committed revision 8.
```

The other student then needs to do another update to get the correction.

### 22.2.4 Inspecting the history

**Purpose.** In this section, you will learn how to get information about the repository.

You've already seen `svn info` as a way of getting information about the repository. To get the history, do `svn log` to get all log messages, or `svn log 2:5` to get a range.

To see differences in various revisions of individual files, use `svn diff`. First do `svn commit` and `svn update` to make sure you are up to date. Now do `svn diff firstfile`. No output, right? Now make an edit in `firstfile` and do `svn diff firstfile` again. This gives you the difference between the last committed version and the working copy.
You can also ask for differences between committed versions with `svn diff -r 4:6 firstfile`.

The output of this `diff` command is a bit cryptic, but you can understand it without too much trouble. There are also fancy GUI implementations of `svn` for every platform that show you differences in a much nicer way.

If you simply want to see what a file used to look like, do `svn cat -r 2 firstfile`. To get a copy of a certain revision of the repository, do `svn export -r 3 ../rev3`, which exports the repository at the current directory ('dot') to the directory `../rev3`.

If you save the output of `svn diff`, it is possible to apply it with the Unix `patch` command. This is a quick way to send patches to someone without them needing to check out the repository.

### 22.2.5 Shuffling files around

We now realize that we really wanted all these files in a subdirectory in the repository. First we create the directory, putting it under `svn` control:

<table><tr><td>%%</td><td>svn mkdir trunk</td></tr><tr><td>A</td><td>trunk</td></tr></table>

Then we move all files there, again prefixing all commands with `svn`:

<table><tr><td>%%</td><td>for</td><td>f</td><td>in</td><td>firstfile</td><td>otherfile</td><td>myfile</td><td>mysecondfile</td><td>do</td><td>\</td></tr><tr><td></td><td></td><td></td><td></td><td>svn mv</td><td>$f</td><td>trunk/</td><td>;</td><td>done</td><td></td></tr><tr><td>A</td><td></td><td></td><td></td><td>trunk/firstfile</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>D</td><td></td><td></td><td></td><td>firstfile</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>A</td><td></td><td></td><td></td><td>trunk/otherfile</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>D</td><td></td><td></td><td></td><td>otherfile</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>A</td><td></td><td></td><td></td><td>trunk/myfile</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>D</td><td></td><td></td><td></td><td>myfile</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>A</td><td></td><td></td><td></td><td>trunk/mysecondfile</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>D</td><td></td><td></td><td></td><td>mysecondfile</td><td></td><td></td><td></td><td></td><td></td></tr></table>

Finally, we commit these changes:

<table><tr><td>%%</td><td>svn commit</td><td>-m</td><td>"trunk created"</td></tr><tr><td>Deleting</td><td></td><td>firstfile</td><td></td></tr><tr><td>Adding</td><td></td><td>trunk/firstfile</td><td></td></tr><tr><td>Deleting</td><td></td><td>myfile</td><td></td></tr><tr><td>Deleting</td><td></td><td>mysecondfile</td><td></td></tr><tr><td>Deleting</td><td></td><td>otherfile</td><td></td></tr><tr><td>Adding</td><td></td><td>trunk</td><td></td></tr><tr><td>Adding</td><td></td><td>trunk/myfile</td><td></td></tr><tr><td>Adding</td><td></td><td>trunk/mysecondfile</td><td></td></tr><tr><td>Adding</td><td></td><td>trunk/otherfile</td><td></td></tr></table>
You probably now have picked up on the message that you always use `svn` to do file manipulations. Let's pretend this has slipped your mind.

**Exercise.** Create a file `somefile` and commit it to the repository. Then do `rm somefile`, thereby deleting a file without `svn` knowing about it. What is the output of `svn status`?

*Expected outcome.* `svn` indicates with an exclamation point that the file has disappeared.

You can fix this situation in a number of ways:

* `svn revert` restores the file to the state in which it was last restored. For a deleted file, this means that it is brought back into existence from the repository. This command is also useful to undo any local edits, if you change your mind about something.

* `svn rm firstfile` is the official way to delete a file. You can do this even if you have already deleted the file outside `svn`.

* Sometimes `svn` will get confused about your attempts to delete a file. You can then do `svn --force rm yourfile`.

### 22.2.6 Branching and merging

Suppose you want to tinker with the repository, while still staying up to date with changes that other people make.

```bash
%% svn copy \
    <the URL of your repository>/trunk \
    <the URL of your repository>/onebranch \
    -m "create a branch"
```

Committed revision 11.

You can check this out as before:

```bash
%% svn co <the URL of your repository>/onebranch \
    projectbranch
```

* `projectbranch/mysecondfile`

* `projectbranch/otherfile`

* `projectbranch/myfile`

* `projectbranch/firstfile`

Checked out revision 11.

Now, if you make edits in this branch, they will not be visible in the trunk:

```bash
%% emacs firstfile # do some edits here
```

```bash
%% svn commit -m "a change in the branch"
```

```
Sending firstfile
```

```
Transmitting file data.
```

Committed revision 13.
```bash
%% (cd ../myproject/trunk/ ; svn update )
At revision 13.
```

On the other hand, edits in the main trunk can be pulled into this branch:

```bash
%% svn merge ^/trunk
--- Merging r13 through r15 into '.':
U      secondfile
```

When you are done editing, the branch edits can be added back to the trunk. For this, it is best to have a clean checkout of the branch:

```bash
%% svn co file://`pwd`/repository/trunk mycleanproject
A      # all the current files
```

and then do a special merge:

```bash
%% cd mycleanproject
%% svn merge --reintegrate ^/branch
--- Merging differences between repository URLs into '.':
U      firstfile
U      .
%% svn info
Path: .
URL: <the URL of your repository>/trunk
Repository Root: <the URL of your repository>
Repository UUID: dc38b821-b9c6-4a1a-a194-a894fba1d7e7
Revision: 16
Node Kind: directory
Schedule: normal
Last Changed Author: build
Last Changed Rev: 14
Last Changed Date: 2009-05-18 13:34:55 -0500 (Mon, 18 May 2009)
%% svn commit -m "trunk updates from the branch"
Sending .
Sending firstfile
Transmitting file data .
Committed revision 17.
```

### 22.2.7 Repository browsers

The `svn` command can give you some amount of information about a repository, but there are graphical tools that are easier to use and that give a better overview of a repository. For the common platforms, several such tools exist, free or commercial. Here is a browser based tool: http://www.websvn.info.
## 22.3 Mercurial or hg

Mercurial is one of a new generation of *distributed source code control* systems. Many of its commands are the same as for subversion, but there are some new ones, corresponding to the new level of sophistication.

A good introduction to Mercurial can be found at http://www.hginit.com/.

This lab should be done two people, to simulate a group of programmers working on a joint project. You can also do this on your own by using two clones of the repository, preferably opening two windows on your computer.

### 22.3.1 Create and populate a repository

**Purpose.** In this section you will create a repository and make a local copy to work on.

First we need to have a repository. In practice, you will often use one that has been previously set up, but there are several ways to set up a repository yourself. There are commercial and free hosting services such as http://bitbucket.org. (Academic users can have more private repositories.)

Let's assume that one student has created a repository `your-project` on Bitbucket. Both students can then clone it:

```bash
% hg clone https://YourName@bitbucket.org/YourName/your-project
updating to branch default
0 files updated, 0 files merged, 0 files removed,
0 files unresolved
```

You now have an empty directory `your-project`.

**Exercise.** Go into the project directory and see if it is really empty.

*Expected outcome.* There is a hidden directory `.hg`

### 22.3.2 New files

**Purpose.** In this section you will make some simple changes: editing an existing file and creating a new file.

One student now makes a file to add to the repository:

```bash
% cat > firstfile
a
b
c
d
e
```
f

^D

(where ^D stands for control-D, which terminates the input.) This file is unknown to hg:

%% hg status

? firstfile

We need to declare the file as belonging to the repository; a subsequent hg commit command then copies
it into the repository.

%% hg add firstfile

%% hg status

A firstfile

%% hg commit -m "made a first file"

Unlike with Subversion, the file has now only been copied into the local repository, so that you can, for
instance, roll back your changes. If you want this file added to the master repository, you need the hg
push command:

%% hg push https://YourName@bitbucket.org/YourName/your-project

pushing to https://YourName@bitbucket.org/YourName/your-project

searching for changes

remote: adding changesets

remote: adding manifests

remote: adding file changes

remote: added 1 changesets with 1 changes to 1 files

remote: bb/acl: YourName is allowed. accepted payload.

In the push step you were probably asked for your password. You can prevent that by having some lines
in your $HOME/.hgrc file:

[paths]

projectrepo = https://YourName:yourpassword@bitbucket.org/YourName/my-project

[ui]

username=Your Name <you@somewhere.youruniversity.edu>

Now the command hg push projectrepo will push the local changes to the global repository without
asking for your password. Of course, now you have a file with a cleartext password, so you should set the
permissions of this file correctly.

The second student now needs to update their repository. Just like the upload took two commands, this pass
also takes two. First you do hg pull to update your local repository. This does not update the local files
you have: for that you need to do hg update.

**Exercise.** Do this and check that the contents of the file are correct.

*Expected outcome.* In order to do the update command, you have to be in a checked-out copy
of the repository.
Caveats.

**Exercise.** Let both students create a new directory with a few files. Declare the directory and commit it. Pull and update to obtain the changes the other mde.

*Expected outcome.* You can do `hg add` on the directory, this will also add the files contained in it.

In order for Mercurial to keep track of your files, you should never do the shell commands `cp` or `mv` on files that are in the repository. Instead, do `hg cp` or `hg mv`. Likewise, there is a command `hg rm`.

### 22.3.3 Conflicts

**Purpose.** In this section you will learn about how do deal with conflicting edits by two users of the same repository.

Now let's see what happens when two people edit the same file. Let both students make an edit to firstfile, but one to the top, the other to the bottom. After one student commits the edit, the other can commit changes, after all, these only affect the local repository. However, trying to push that change gives an error:

```
%% emacs firstfile # make some change
%% hg commit -m ``first again''
%% hg push test
abort: push creates new remote head b0d31ea209b3!
(you should pull and merge or use push -f to force)
```

The solution is to get the other edit, and commit again. This takes a couple of commands:

```
%% hg pull myproject
searching for changes
adding changesets
adding manifests
adding file changes
added 1 changesets with 1 changes to 1 files (+1 heads)
 run 'hg heads' to see heads, 'hg merge' to merge

%% hg merge
merging firstfile
0 files updated, 1 files merged, 0 files removed, 0 files unresolved
 branch merge, don't forget to commit

%% hg status
M firstfile
%% hg commit -m ``my edit again''
%% hg push test
pushing to https://VictorEijkhout:***@bitbucket.org/
VictorEijkhout/my-project
searching for changes
```
remote: adding changesets
remote: adding manifests
remote: adding file changes
remote: added 2 changesets with 2 changes to 1 files
remote: bb/acl: VictorEijkhout is allowed. accepted payload.

This may seem complicated, but you see that Mercurial prompts you for what commands to execute, and the workflow is clear, if you refer to figure 22.2.

**Exercise.** Do a `cat` on the file that both of you have been editing. You should find that both edits are incorporated. That is the 'merge' that Mercurial referred to.

If both students make edits on the same part of the file, Mercurial can no longer resolve the conflicts. For instance, let one student insert a line between the first and the second, and let the second student edit the second line. Whoever tries to push second, will get messages like this:

```
%% hg pull test
added 3 changesets with 3 changes to 1 files (+1 heads)
 run('hg heads' to see heads, 'hg merge' to merge)
%% hg merge
merging firstfile
warning: conflicts during merge.
merging firstfile incomplete!
 (edit conflicts, then use 'hg resolve --mark')
0 files updated, 0 files merged, 0 files removed, 1 files unresolved
use 'hg resolve' to retry unresolved file merges
or 'hg update -C .' to abandon
```

Mercurial will give you several options. It is easiest to resolve the conflict with a text editor. If you open the file that has the conflict you'll see something like:

```
<<<<<<<<< local
aa
bbbb
========
aaa
a2
b
>>>>>>>>> other
c
```

indicating the difference between the local version ('mine') and the other, that is the version that you pulled and tried to merge. You need to edit the file to resolve the conflict.

After this, you tell hg that the conflict was resolved:

```
hg resolve --mark
%% hg status
```
```
M firstfile
? firstfile.orig
```

After this you can commit and push again. The other student then needs to do another update to get the correction.

Not all files can be merged: for binary files Mercurial will ask you:

```
%% hg merge
merging proposal.tex
merging summary.tex
merking references.tex
no tool found to merge proposal.pdf
keep (l)ocal or take (o)ther? o
```

This means that the only choices are to keep your local version (type `l` and hit return) or take the other version (type `o` and hit return). In the case of a binary file that was obvious generated automatically, some people argue that they should not be in the repository to begin with.

### 22.3.4 Inspecting the history

**Purpose.** In this section, you will learn how to get information about the repository.

If you want to know where you cloned a repository from, look in the file `.hg/hgrc`.

The main sources of information about the repository are `hg log` and `hg id`. The latter gives you global information, depending on what option you use. For instance, `hg id -n` gives the local revision number.

`hg log` gives you a list of all changesets so far, with the comments you entered.

`hg log -v` tells you what files were affected in each changeset.

`hg log -r 5` or `hg log -r 6:8` gives information on one or more changesets.

To see differences in various revisions of individual files, use `hg diff`. First make sure you are up to date. Now do `hg diff firstfile`. No output, right? Now make an edit in `firstfile` and do `hg diff firstfile` again. This gives you the difference between the last committed version and the working copy.

Check for yourself what happens when you do a commit but no push, and you issue the above `diff` command.

You can also ask for differences between committed versions with `hg diff -r 4:6 firstfile`.

The output of this `diff` command is a bit cryptic, but you can understand it without too much trouble. There are also fancy GUI implementations of `hg` for every platform that show you differences in a much nicer way.

If you simply want to see what a file used to look like, do `hg cat -r 2 firstfile`. To get a copy of a certain revision of the repository, do `hg export -r 3 ../rev3`, which exports the repository at the current directory ('dot') to the directory `../rev3`.

If you save the output of `hg diff`, it is possible to apply it with the Unix `patch` command. This is a quick way to send patches to someone without them needing to check out the repository.
## Chapter 23

### Scientific Data Storage

There are many ways of storing data, in particular data that comes in arrays. A surprising number of people stores data in spreadsheets, then exports them to ascii files with comma or tab delimiters, and expects other people (or other programs written by themselves) to read that in again. Such a process is wasteful in several respects:

* The ascii representation of a number takes up much more space than the internal binary representation. Ideally, you would want a file to be as compact as the representation in memory.

* Conversion to and from ascii is slow; it may also lead to loss of precision.

For such reasons, it is desirable to have a file format that is based on binary storage. There are a few more requirements on a useful file format:

* Since binary storage can differ between platforms, a good file format is platform-independent. This will, for instance, prevent the confusion between *big-endian* and *little-endian* storage, as well as conventions of 32 versus 64 bit floating point numbers.

* Application data can be heterogeneous, comprising integer, character, and floating point data. Ideally, all this data should be stored together.

* Application data is also structured. This structure should be reflected in the stored form.

* It is desirable for a file format to be *self-documenting*. If you store a matrix and a right-hand side vector in a file, wouldn't it be nice if the file itself told you which of the stored numbers are the matrix, which the vector, and what the sizes of the objects are?

This tutorial will introduce the `HDF5` library, which fulfills these requirements. `HDF5` is a large and complicated library, so this tutorial will only touch on the basics. For further information, consult http://www.hdfgroup.org/HDF5/. While you do this tutorial, keep your browser open on http://www.hdfgroup.org/HDF5/doc/ or http://www.hdfgroup.org/HDF5/RM/RM_H5Front.html for the exact syntax of the routines.

#### 23.1 Introduction to HDF5

As described above, `HDF5` is a file format that is machine-independent and self-documenting. Each `HDF5` file is set up like a directory tree, with subdirectories, and leaf nodes which contain the actual data. This means that data can be found in a file by referring to its name, rather than its location in the file. In this
Expected outcome. A file `file.h5` should be created.

Caveats. Be sure to add `HDF5` include and library directories:

```
cc -c create.c -I. -I/opt/local/include
and
```

```
cc -o create create.o -L/opt/local/lib -lhdf5. The include and lib di-
rectories will be system dependent.
```

```
/*
 * File: create.c
 * Author: Victor Eijkhout
 */
#include "myh5defs.h"
#define FILE "file.h5"
```

```c
main() {
```

```c
hid_t      file_id;  /* file identifier */
herr_t     status;
```

```c
/* Create a new file using default properties. */
file_id = H5Fcreate(FILE, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
H5REPORT(file_id);
```

```c
/* Terminate access to the file. */
status = H5Fclose(file_id);
```

```c
}
```

You can display the created file on the commandline:

```
%% h5dump file.h5
HDF5 "file.h5" {
  GROUP "/" {
  }
}
```

Note that an empty file corresponds to just the root of the directory tree that will hold the data.

## 23.3 Datasets

Next we create a dataset, in this example a 2D grid. To describe this, we first need to construct a `dataspace`:

```
dims[0] = 4; dims[1] = 6;
dataspace_id = H5Screate_simple(2, dims, NULL);
dataset_id = H5Dcreate(file_id, "/dset", dataspace_id, .... );
....
status = H5Dclose dataset_id);
```
status = H5Close(dataspace_id);

Note that datasets and dataspaces need to be closed, just like files.

**Exercise.** Create a dataset by compiling and running the `dataset.c` code below

*Expected outcome.* This creates a file `dset.h` that can be displayed with `h5dump`.

```c
/*
 * File: dataset.c
 * Author: Victor Eijkhout
 */
#include "myh5defs.h"
#define FILE "dset.h5"

main() {
    hid_t      file_id, dataset_id, dataspace_id; /* identifiers */
    hsize_t    dims[2];
    herr_t     status;

    /* Create a new file using default properties. */
    file_id = H5Fcreate(FILE, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);

    /* Create the data space for the dataset. */
    dims[0] = 4;
    dims[1] = 6;
    dataspace_id = H5Screate_simple(2, dims, NULL);

    /* Create the dataset. */
    dataset_id = H5Dcreate(file_id, "/dset", H5T_native_INT,
                             dataspace_id, H5P_DEFAULT);

    /*H5T_STD_I32BE*/

    /* End access to the dataset and release resources used by it. */
    status = H5Dclose dataset_id);

    /* Terminate access to the data space. */
    status = H5Sclose dataspace_id);

    /* Close the file. */
    status = H5Fclose(file_id);
}
```

We again view the created file online:

```
%% h5dump dset.h5
HDF5 "dset.h5" {
  GROUP "/" {
```
```c
DATASET "dset" {
    DATATYPE H5T_STD_I32BE
    DATASPACE SIMPLE { ( 4, 6 ) / ( 4, 6 ) }
    DATA {
        (0,0): 0, 0, 0, 0, 0, 0,
        (1,0): 0, 0, 0, 0, 0, 0,
        (2,0): 0, 0, 0, 0, 0, 0,
        (3,0): 0, 0, 0, 0, 0, 0
    }
}
}
}
```

The datafile contains such information as the size of the arrays you store. Still, you may want to add related scalar information. For instance, if the array is output of a program, you could record with what input parameter was it generated.

```c
parmspace = H5Screate(H5S_SCALAR);
parm_id = H5Dcreate
    (file_id, "/parm", H5T_native_INT, parmspace, H5P_DEFAULT);
```

**Exercise.** Add a scalar dataspace to the HDF5 file, by compiling and running the `parmwrite.c` code below.

*Expected outcome.* A new file `wdset.h5` is created.

```c
/*
 * File: parmdataset.c
 * Author: Victor Eijkhout
 */
#include "myh5defs.h"
#define FILE "pdset.h5"

main() {
    hid_t        file_id, dataset_id, dataspace_id; /* identifiers */
    hid_t        parm_id, parmspace;
    hsize_t      dims[2];
    herr_t       status;

    /* Create a new file using default properties. */
    file_id = H5Fcreate(FILE, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);

    /* Create the data space for the dataset. */
    dims[0] = 4;
    dims[1] = 6;
    dataspace_id = H5Screate_simple(2, dims, NULL);
```
/* Create the dataset. */
dataset_id = H5Dcreate
(file_id, "/dset", H5T_STD_I32BE, dataspace_id, H5P_DEFAULT);

/* Add a descriptive parameter */
parmspace = H5Screate(H5S_SCALAR);
parm_id = H5Dcreate
(file_id, "/parm", H5T_native_INT, parmspace, H5P_DEFAULT);

/* End access to the dataset and release resources used by it. */
status = H5Dclose_id);
status = H5Dclose(parm_id);

/* Terminate access to the data space. */
status = H5Sclose dataspace_id);
status = H5Sclose(parmspace);

/* Close the file. */
status = H5Fclose(file_id);
}

%% h5dump wdset.h5
HDF5 "wdset.h5" {
GROUP "/" {
    DATASET "dset" {
        DATATYPE H5T_IEEE_F64LE
        DATASPACE SIMPLE { ( 4, 6 ) / ( 4, 6 ) }
        DATA {
            (0,0): 0.5, 1.5, 2.5, 3.5, 4.5, 5.5,
            (1,0): 6.5, 7.5, 8.5, 9.5, 10.5, 11.5,
            (2,0): 12.5, 13.5, 14.5, 15.5, 16.5, 17.5,
            (3,0): 18.5, 19.5, 20.5, 21.5, 22.5, 23.5
        }
    }
    DATASET "parm" {
        DATATYPE H5T_STD_I32LE
        DATASPACE SCALAR
        DATA {
            (0): 37
        }
    }
}
}
## 23.4 Writing the data

The datasets you created allocate the space in the `hdf5` file. Now you need to put actual data in it. This is done with the `H5Dwrite` call.

```c
/* Write floating point data */
for (i=0; i<24; i++) data[i] = i+.5;
status = H5Dwrite
    (dataset, H5T_native_double, H5S_ALL, H5S_ALL, H5P_DEFAULT,
     data);
/* write parameter value */
parm = 37;
status = H5Dwrite
    (parmset, H5T_native_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
     &parm);

/*
 * File: parmwrite.c
 * Author: Victor Eijkhout
 */
#include "myh5defs.h"
#define FILE "wdset.h5"

main() {
    hid_t      file_id, dataset, dataspace; /* identifiers */
    hid_t      parmset,parmspace;
    hsize_t    dims[2];
    herr_t     status;
    double data[24]; int i,parm;

    /* Create a new file using default properties. */
    file_id = H5Fcreate(FILE, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);

    /* Create the dataset. */
    dims[0] = 4; dims[1] = 6;
    dataspace = H5Screate_simple(2, dims, NULL);
    dataset = H5Dcreate
        (file_id, "/dset", H5T_native_double, dataspace, H5P_DEFAULT);

    /* Add a descriptive parameter */
    parmspace = H5Screate(H5S_SCALAR);
    parmset = H5Dcreate
        (file_id, "/parm", H5T_native_INT, parmspace, H5P_DEFAULT);

    /* Write data to file */
    for (i=0; i<24; i++) data[i] = i+.5;
    status = H5Dwrite
        (dataset, H5T_native_double, H5S_ALL, H5S_ALL, H5P_DEFAULT,
         data); H5REPORT(status);
```
```c
/* write parameter value */
parm = 37;
status = H5Dwrite
    (parmset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
     &parm); H5REPORT(status);

/* End access to the dataset and release resources used by it. */
status = H5Dclose dataset;
status = H5Dclose (parmset);

/* Terminate access to the data space. */
status = H5Sclose (dataspace);
status = H5Sclose (parmspace);

/* Close the file. */
status = H5Fclose (file_id);
}
```

```c
%% h5dump wdset.h5
HDF5 "wdset.h5" {
  GROUP "/" {
    DATASET "dset" {
      DATATYPE H5T_IEEE_F64LE
      DATASPACE SIMPLE { ( 4, 6 ) / ( 4, 6 ) }
      DATA {
        (0,0): 0.5, 1.5, 2.5, 3.5, 4.5, 5.5,
        (1,0): 6.5, 7.5, 8.5, 9.5, 10.5, 11.5,
        (2,0): 12.5, 13.5, 14.5, 15.5, 16.5, 17.5,
        (3,0): 18.5, 19.5, 20.5, 21.5, 22.5, 23.5
      }
    }
    DATASET "parm" {
      DATATYPE H5T_STD_I32LE
      DATASPACE SCALAR
      DATA {
        (0): 37
      }
    }
  }
}
```

If you look closely at the source and the dump, you see that the data types are declared as 'native', but rendered as LE. The 'native' declaration makes the datatypes behave like the built-in C or Fortran data types. Alternatively, you can explicitly indicate whether data is little-endian or big-endian. These terms describe how the bytes of a data item are ordered in memory. Most architectures use little endian, as you can see in the dump output, but, notably, IBM uses big endian.
## 23.5 Reading

Now that we have a file with some data, we can do the mirror part of the story: reading from that file. The essential commands are

```R
h5file = H5Fopen(....)
....
H5Dread dataset, .... data ....)
```

where the `H5Dread` command has the same arguments as the corresponding `H5Dwrite`.

**Exercise.** Read data from the `wdset.h5` file that you create in the previous exercise, by compiling and running the `allread.c` example below.

*Expected outcome.* Running the `allread` executable will print the value 37 of the parameter, and the value 8.5 of the (1, 2) data point of the array.

*Caveats.* Make sure that you run `parmwrite` to create the input file.

```R
/*
 * File: allread.c
 * Author: Victor Eijkhout
 */
#include "myh5defs.h"
#define FILE "wdset.h5"

main() {
    hid_t      file_id, dataset, parmset;
    herr_t     status;
    double data[24]; int parm;

    /* Open an existing file */
    file_id = H5Fopen(FILE, H5F_ACC_RDONLY, H5P_DEFAULT);
    H5REPORT(file_id);

    /* Locate the datasets. */
    dataset = H5Dopen(file_id, "/dset"); H5REPORT;
    parmset = H5Dopen(file_id, "/parm"); H5REPORT(parmset);

    /* Read data back */
    status = H5Dread
        (parmset, H5T_native_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
         &parm); H5REPORT(status);
    printf("parameter value: %d\n", parm);

    status = H5Dread
        (dataset, H5T_native_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,
         data); H5REPORT(status);
    printf("arbitrary data point [1,2]: %e\n", data[1*6+2]);
```
/* Terminate access to the datasets */
status = H5Dclose(dataset); H5REPORT(status);
status = H5Dclose(parmset); H5REPORT(status);

/* Close the file. */
status = H5Fclose(file_id);

%% ./allread
parameter value: 37
arbitrary data point [1,2]: 8.500000e+00
# Chapter 24

## Scientific Libraries

There are many libraries for scientific computing. Some are specialized to certain application areas, others are quite general. In this section we will take a brief look at the PETSc library for sparse matrix computations, and the BLAS/Lapack libraries for dense computations.

### 24.1 The Portable Extendable Toolkit for Scientific Computing

PETSc, the Portable Extendable Toolkit for Scientific Computation [], is a large powerful library, mostly concerned with linear and nonlinear system of equations that arise from discretized PDEs. Since it has many hundreds of routines (and a good manual already exists) we limit this tutorial to going through a few simple, yet representative, PETSc programs.

#### 24.1.1 What is in PETSc?

PETSc can be used as a library in the traditional sense, where you use some high level functionality, such as solving a nonlinear system of equations, in your program. However, it can also be used as a toolbox, to compose your own numerical applications using low-level tools.

* Linear system solvers (sparse/dense, iterative/direct)

* Nonlinear system solvers

* Tools for distributed matrices

* Support for profiling, debugging, graphical output

The basic functionality of PETSc can be extended through external packages:

* Dense linear algebra: `Scalapack`, `Plapack`

* Grid partitioning software: `ParMetis`, `Jostle`, `Chaco`, `Party`

* ODE solvers: `PVODE`

* Eigenvalue solvers (including SVD): `SLEPc`

* Optimization: `TAO`
### 24.1.2 Design of PETSc

PETSc has an object-oriented design, even though it is implemented in C, a not object oriented language. PETSc is also parallel: all objects in `Petsc` are defined on an `MPI` communicator (section 2.7.3.3). For an object such as a matrix this means that it is stored distributed over the processors in the communicator; for objects such as solvers it means that their operation is distributed over the communicator. PETSc objects can only interact if they are based on the same communicator. Most of the time objects will be defined on the `MPI_COMM WORLD` communicator, but subcommunicators can be used too. This way, you can define a matrix or a linear system on all your available processors, or on a subset.

Parallelism is handled through `MPI`. You may need to have occasional `MPI` calls in a PETSc code, but the vast bulk of communication is done behind the scenes inside PETSc calls. Shared memory parallel (such as through OpenMP; section 2.7.2) is not used explicitly, but the user can incorporate it without problems.

The object-oriented design means that a call such as matrix-vector multiplication

```c
MATMult(A,x,y); // y <- A x
```

looks the same, regardless whether `A` is sparse or dense, sequential or parallel.

One implication of this is that the actual data are usually hidden from the user. Data can be accessed through routines such as

```c
double *array;
VecGetArray(myvector,&array);
```

but most of the time the application does not explicitly maintain the data, only the PETSc objects containing the data. As the ultimate consequence of this, the user usually does not allocate data; rather, matrix and vector arrays get created through the PETSc create calls, and subsequent values are inserted through PETSc calls:

```c
MatSetValue(A,i,j,v,INSERT_VALUES); // A[i,j] <- v
```

This may feel like the user is giving up control, but it actually makes programming a lot easier, since the user does not have to worry about details of storage schemes, especially in parallel.

### 24.1.3 Small examples

In this section we will go through a number of successively more complicated examples of the use of PETSc. The files can be downloaded from the repository of this book. While doing these examples it is a good idea to keep the manual page open:

http://tinyurl.com/PETSc-man-page.

When you do these examples, make sure to use a version of PETSc that has debug mode enabled!
### 24.1.3.1 Program structure

The first example (we only list C sources in this book; the download includes their Fortran equivalents) illustrates the basic structure of a PETSc program: the include file at the top and the calls to `PetscInitialize`, `PetscFinalize`. Furthermore it uses two routines:

* `PetscOptionsGetInt` is used to check if the program was invoked with a `-n` option and a numerical value: if you ran your program as `./init -n 45`, the variable `n` should have the value 45.

* `Petscuatef` functions like the `printf` function, except that only one processor executes the print command. If you ran your program as `mpiexec -np 4 init` and you used `printf`, you would get four copies of the output.

```c
// init.c
int main(int argc, char **argv)
{
    MPI_Comms       comm;
    PetscInt        n = 20;
    PetscErrorCode  ierr;

    PetscFunctionBegin;
    ierr = PetscInitialize(&argc, &argv, 0, 0); CHKERRQ(ierr);
    comm = PETSC_COMM WORLD;
    ierr = PetscOptionsGetInt(PETSC_NULL, "-n", &n, PETSC_NULL); CHKERRQ(ierr);
    ierr = Petscuatef comm, "Input parameter: %d\n", n); CHKERRQ(ierr);

    ierr = PetscFinalize(); CHKERRQ(ierr);
    PetscFunctionReturn(0);
}
```

Just about the only lines of MPI that you need to know when programming PETSc are:

```
C:
    ierr = MPI_Comms_size(comm, &ntids);
    ierr = MPI_Comms_rank(comm, &mytid);

F:
    call MPI_Comms_size(comm, &ntids, ierr);
    call MPI_Comms_rank(comm, &mytid, ierr);
```

The first line gives the size of the communicator, meaning how many processes there are; the second one gives the rank of the current process as a number from zero to the number of processes minus one.

**Exercise 24.1.** Add these two lines to your program. Look up the routine `PetscSynchronized PRINTf` in the documentation and use it to let each process print out a line like `Process 3 out of 7`. You may also need to use `PetscSynchronizedFlush`.
### 24.1.3.2 Vectors

Next we start making PETSc objects, first a vector.

vec.c

```c
#include "petsc.h"

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc, char **argv)
{
    MPI_Comms        comm;
    int ntids, mytid, myfirst, mylast;
    Vec              x;
    PetscInt         n = 20;
    PetscReal        one = 1.0;
    PetscErrorCode    ierr;

    PetscFunctionBegin;
    ierr = PetscInitialize(&argc, &argv, 0, 0); CHKERRQ(ierr);
    comm = PETSC_COMM WORLD;
    ierr = PetscOptionsGetInt(PETSC_NULL, "-n", &n, PETSC_NULL);
    CHKERRQ(ierr);

    MPI_Comms_size( comm, &ntids); MPI_Comms_rank( comm, &mytid);

    ierr = VecCreate( comm, &x); CHKERRQ(ierr);
    ierr = VecSetSizes( x, PETSC_DECIDE, n); CHKERRQ(ierr);
    ierr = VecSetType( x, VECMPI); CHKERRQ(ierr);

    ierr = VecGetOwnershipRange( x, &myfirst, &mylast); CHKERRQ(ierr);
    ierr = PetscSynchronized( comm, "Proc %d, range %d--%d\n",
                                   mytid, myfirst, mylast); CHKERRQ(ierr);

    ierr = VecSet( x, one); CHKERRQ(ierr);
    ierr = VecAssemblyBegin( x); CHKERRQ(ierr);
    ierr = VecAssemblyEnd( x); CHKERRQ(ierr);
    ierr = VecView( x, 0); CHKERRQ(ierr);

    ierr = VecDestroy( &x); CHKERRQ(ierr);
    ierr = PetscFinalize(); CHKERRQ(ierr);
    PetscFunctionReturn( 0);
}
```

Note how it takes several calls to fully create the vector object:

* The type VECMPI means that the vector will be distributed.

* The routine setting the size has two size parameters; the second specifies that the global size is n, and the first one says that you leave the distribution for PETSc to decide.
At the end of the program there is a to `VecDestroy`, which deallocates the memory for the vector. While this is strictly speaking not necessary for the functioning of the program, it is a good idea to issue `Destroy` calls for each `Create` call, since it can prevent potential *memory leaks*.

**Exercise 24.2.** Comment out the `VecDestroy` call, and run the program with the option `-malloc_dump`. PETSc will now report on all memory that had not been freed at the time of the `PetscFinalize` call.

If you run the program in parallel the vector will be created distributed over the processors. The program contains a call to `VecGetOwnershipRange` to discover what part of the vector lives on what processor. You see that the `VecView` calls display for each processor precisely that part of the vector.

**Exercise 24.3.** The listing has a call to `VecSet` to set the vector elements to a constant value. Remove this line. Use the `myfirst`, `mylast` variables and the PETSc routine `VecSetValue` to let every process set its local vector elements to the processor number. (The last argument of `VecSetValue` should be `INSERT_VALUES`.) That is, if the vector is six elements long and there are three processors, the resulting vector should be `(0, 0, 1, 1, 2, 2)`.

**Exercise 24.4.** Change the code from the previous exercise. Now every vector element should be the sum of the processor number and the previous processor numbers; in the above example the result should be `(0, 0, 1, 1, 3, 3)`. Read the man page for `VecSetValue`!

Run the code from the previous two exercises again, adding a commandline argument `-log_summary`. Observe that the first code has no messages in the `VecAssemblyBegin/End` calls, but the second one does.

### 24.1.3.3 Matrices

Let's move up to matrices. Creating a matrix is similar to creating a vector. In this case the type is `MPIAIJ` which stands for a distributed sparse matrix.

mat.c

/*
* This program file is part of the textbook and course
* 'Introduction to High-Performance Scientific Computing'
* by Victor Eijkhout, copyright 2012/3/4/5
*
* This book is distributed under a Creative Commons Attribution 3.0
* Unported (CC BY 3.0) license and made possible by funding from
* The Saylor Foundation \url{http://www.saylor.org}.
*
* include "petsc.h"
*/
```c
#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc, char **argv)
{
    MPI_Comms       comm;
    int ntids, mytid, localsize, myfirst, mylast, i;
    Mat             A;
    PetscInt        n = 20;
    PetscErrorCode  Perr;

    PetscFunctionBegin;
    Perr = PetscInitialize(&argc, &argv, 0, 0); CHKERRQ(Perr);
    comm = PETSC_COMM_WORLD;
    Perr = PetscOptionsGetInt(PETSC_NULL, "-n", &n, PETSC_NULL);
    CHKERRQ(Perr);

    MPI_Comms_size(comm, &ntids); MPI_Comms_rank(comm, &mytid);
    localsize = PETSC_DECIDE;
    Perr = PetscSplitOwnership(comm, &localsize, &n); CHKERRQ(Perr);

    Perr = MatCreate(comm, &A); CHKERRQ(Perr);
    Perr = MatSetType(A, MATMPIAIJ); CHKERRQ(Perr);
    Perr = MatSetSizes(A, localsize, localsize,
                        PETSC_DECIDE, PETSC_DECIDE); CHKERRQ(Perr);
    Perr = MatGetOwnershipRange(A, &myfirst, &mylast); CHKERRQ(Perr);

    for (i=myfirst; i<mylast; i++) {
        PetscReal v=1.0*mytid;
        Perr = MatSetValues(A, 1, &i, 1, &i, &v, INSERT_VALUES); CHKERRQ(Perr);
        /*
         * Exercise: add some off-diagonal matrix elements
         */
    }

    Perr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); CHKERRQ(Perr);
    Perr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY); CHKERRQ(Perr);
    Perr = MatView(A, 0); CHKERRQ(Perr);

    Perr = MatDestroy(&A); CHKERRQ(Perr);
    Perr = PetscFinalize(); CHKERRQ(Perr);
    PetscFunctionReturn(0);
}
```

In this example a diagonal matrix is constructed, with each processor setting only its locally stored elements.

**Exercise 24.5.** In addition to setting the diagonal, also set the first subdiagonal and the first superdiagonal, that is, the elements $(i, i - 1)$ and $(i, i + 1)$. To set all three elements in a row with one call you can do this:

```
vm = 1.*mytid;
```
```c
for (i=myfirst; i<mylast; i++) {
    PetscInt j[3];      PettscReal v[3];
    j[0] = i-1; j[1] = i; j[2] = i+1;
    v[0] = --vm-1; v[1] = 2*vm; v[2] = -vm+1;
    ierr = MatSetValues
        (A,1,&i,3,j,v,INSERT_VALUES); CHKERRQ(ierr);
}
```

However, this code is not entirely correct. Edit the program using this fragment and run it. Diagnose the problem and fix it.

### 24.1.3.4 Matrix-vector operations

Next we will create a file `mul.c` based on `mat.c` and multiply the matrix and the vector. Make sure that the size declarations of the matrix and the vector are compatible. You also need a second vector to store the result of the multiplication. This is easiest done by

```c
ierr = VecDuplicate(x,&y); CHKERRQ(ierr);
```

**Exercise 24.6.** Look up the `MatMult` routine in the documentation and use it your program. Use `VecView` to inspect the result. Note that no size parameters or anything pertaining to parallelism appears in the calling sequence.

### 24.1.3.5 Solvers

Copy your `mat.c` file to `sys.c`: we are going to explore linear system solving. First you need to create a solver object and give it the matrix as operator:

```c
ierr = KSPCreate comm, &solver); CHKERRQ(ierr);
ierr = KSPSetOperators(solver, A, A, 0); CHKERRQ(ierr);
ierr = KSPSolve(solver, x, y); CHKERRQ(ierr);
ierr = VecView(y, 0); CHKERRQ(ierr);
```

**Exercise 24.7.** Add these lines to your code. Make sure you know what the correct solution is by using `MatMult` to obtain the right hand side.

You have just used the default linear system solver. Run the program again, but with the option `-ksp_view`. This will tell you all the details of what solver was used.

Solving the linear system is a one line call to `KSPSolve`. The story would end there if it weren't for some complications:

* Iterative methods can fail, and the `solve` call does not tell us whether that happened.

* If the system was solved successfully, we would like to know in how many iterations.

* There can be other reason for the iterative method to halt, such as reaching its maximum number of iterations without converging.
Exercise 24.8. Use the routine `KSPGetConvergedReason` to inspect the status of the solution vector. Use `KSPGetIterationNumber` to see how many iterations it took.

Exercise 24.9. Add code to your program to compute the residual and its norm. For the residual, look up the routines `VecDuplicate` and `VecAXPY`; for computing its norm look up `VecNorm`.

Exercise 24.10. Add a call to `KSPSetFromOptions` to your code. Use the option `-ksp_monitor` to observe the convergence behaviour.

The workings of a PETSc program can be customized to a great degree through the use of commandline options. This includes setting the type of the solver. In order for such options to be obeyed, you first need to put a command `KSPSetFromOptions` before the `KSPSolve` call.

Exercise 24.11. The `-ksp_view` option told you what solver and preconditioner were used. Look up the routines `KSPSetType` and `PCSetType` and use those to change the iterative method to CG and the preconditioner to Jacobi. Do this first by using commandline options, and then by editing the code.

### 24.1.4 A realistic program

This section will give you some further help towards solving a realistic PDE problem.

#### 24.1.4.1 Construction of the coefficient matrix

In the examples above you used a commandline argument to determine the matrix size directly. Here we construct the matrix of 5-point stencil for the Poisson operator (see section 4.2.3 and in particular figure 4.1). Determining its size takes two steps: you need to read the domain size $n = 1/h - 1$ and compute the matrix size from it.

C:

```c
int domain_size, matrix_size;
PetscOptionsGetInt
    (PETSC_NULL, "-n", &domain_size, &flag);
matrix_size = domain_size*domain_size;
```

Fortran:

```fortran
integer :: domain_size, matrix_size
call PetscOptionsGetInt(PETSC_NULL_CHARACTER,
    > "-n", domain_size, flag)
matrix_size = domain_size*domain_size;
```

Now you use the `matrix_size` parameter for constructing the matrix.
### 24.1.4.2 Filling in matrix elements

Just like in the examples above, you want each processor to set only its local rows. The easiest way to iterate over those is to iterate over all variables / matrix rows and select only the local ones.

We will now set matrix elements (refer to the full domain, but only inserting those elements that are in its matrix block row.

C:

```c
MatGetOwnershipRange(A, &myfirst, &mylast);
for ( i=0; i<domain_size; i++ ) {
    for ( j=0; j<domain_size; j++ ) {
        I = j + matrix_size*i;
        if (I>=myfirst && I<mylast) {
            J = I; // for the diagonal element
            MatSetValues
                (A, 1, &I, 1, &J, &v, INSERT_VALUES);
            J = ... // for the other points
            J = ....
        }
    }
}

MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
```

Fortran:

```fortran
call MatGetOwnershipRange(A, myfirst, mylast)
do i=0, matrix_size-1
    do j=0, domain_size-1
        ii = j + domain_size*i
        if (ii>=myfirst .and. ii<mylast) then
            jj = ii; for the diagonal element
            call MatSetValues
        >     (A, 1, ii, 1, jj, v, INSERT_VALUES)
            jj = ii...; for the other elements
            jj = ii...
        end if
    end do
end do

call MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)
call MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)
```

**Exercise 24.12.** Construct the matrix from equation (4.14) in section 4.2.3. Compute and output the product with the identity vector (meaning that all elements are 1), and check
that the result is correct. Make sure to test your program in parallel.

### 24.1.4.3 Finite Element Matrix assembly

PETSc's versatility in dealing with Finite Element matrices (see sections 4.2.5 and 6.6.2), where elements are constructed by adding together contributions, sometimes from different processors. This is no problem in PETSc: any processor can set (or add to) any matrix element. The assembly calls will move data to their eventual location on the correct processors.

```c
for (e=myfirstelement; e<mylastelement; e++) {
    for (i=0; i<nlocalnodes; i++) {
        I = localtoglobal(e, i);
        for (j=0; j<nlocalnodes; j++) {
            J = localtoglobal(e, j);
            v = integration(e, i, j);
            MatSetValues
                (mat, 1, &I, 1, &J, &v, ADD_VALUES);
            ...
        }
    }
}

MatAssemblyBegin (mat, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd (mat, MAT_FINAL_ASSEMBLY);
```

### 24.1.4.4 Linear system solver

We have covered the basics of setting up the solver and solving the system above.

As an illustration of the toolbox nature of PETSc, you can now use routines you have already seen to compute the residual and its norm.

Exercise 24.13. Create a new vector z (use `VecDuplicate`) and store the product of A and the computed solution y (use `MatMult`) in it. If you solved the system accurately, z should now be equal to x. To see how close it is, use

```c
PetscReal norm;
VecAXPY(z, -1, x);
VecNorm(z, NORM_2, &norm);
```

to subtract x from z and compute the norm of the result.

## 24.1.5 Quick experimentation

Reading a parameter from the commandline above is actually a special case of a general mechanism for influencing PETSc's workings through commandline options.

Here is an example of setting the iterative solver and preconditioner from the commandline:
```
yourprog -ksp_type gmres -ksp_gmres_restart 25
-pc_type ilu -pc_factorLevels 3
```

In order for this to work, your code needs to call

```
KSPSetFromOptions(solver);
```

before the system solution. This mechanism is very powerful, and it obviates the need for much code recompilation.

## 24.2 Libraries for dense linear algebra: BLAS, Lapack, Scalapack

In this section we will discuss libraries for dense linear algebra operations.

Dense linear algebra, that is linear algebra on matrices that are stored as two-dimensional arrays (as opposed to sparse linear algebra; see section 5.4, as well as the tutorial on PETSc 24) has been standardized for a considerable time. The basic operations are defined by the three levels of Basic Linear Algebra Subprograms (BLAS):

* Level 1 defines vector operations that are characterized by a single loop [111].

* Level 2 defines matrix vector operations, both explicit such as the matrix-vector product, and implicit such as the solution of triangular systems [45].

* Level 3 defines matrix-matrix operations, most notably the matrix-matrix product [44].

The name 'BLAS' suggests a certain amount of generality, but the original authors were clear [111] that these subprograms only covered dense linear algebra. Attempts to standardize sparse operations have never met with equal success.

Based on these building blocks libraries have been built that tackle the more sophisticated problems such as solving linear systems, or computing eigenvalues or singular values. Linpack<fnref n="1" /> and Eispack were the first to formalize these operations involved, using Blas Level 1 and Blas Level 2 respectively. A later development, Lapack uses the blocked operations of Blas Level 3. As you saw in section 1.5.1, this is needed to get high performance on cache-based CPUs. (Note: the reference implementation of the BLAS [13] will not give good performance with any compiler; most platforms have vendor-optimized implementations, such as the MKL library from Intel.)

With the advent of parallel computers, several projects arose that extended the Lapack functionality to distributed computing, most notably Scalapack [26, 12], PLapack [156, 155], and most recently Elemental [137]. These packages are harder to use than Lapack because of the need for a two-dimensional cyclic distribution; sections 6.2.2 and 6.3.2. We will not go into the details here.

<footnote n="1">1. The linear system solver from this package later became the Linpack benchmark; see section 2.11.5.</footnote>
### 24.2.1 Some general remarks

#### 24.2.1.1 The Fortran heritage

The original BLAS routines were written in Fortran, and the reference implementation is still in Fortran. For this reason you will see the routine definitions first in Fortran in this tutorial. It is possible to use the Fortran routines from a C/C++ program:

* You typically need to append an underscore to the Fortran name;

* You need to include a prototype file in your source, for instance `mkl.h`;

* Every argument needs to be a 'star'-argument, so you can not pass literal constants: you need to pass the address of a variable.

* You need to create a column-major matrix.

There are also C/C++ interfaces:

* The C routine names are formed by prefixing the original name with `cblas_`; for instance `dasum` becomes `cblas_dasum`.

* Fortran character arguments have been replaced by enumerated constants, for instance `CblasNoTrans` instead of the 'N' parameter.

* The `Cblas` interface can accommodate both row-major and column-major storage.

* Array indices are 1-based, rather than 0-based; this mostly becomes apparent in error messages and when specifying pivot indices.

#### 24.2.1.2 Routine naming

Routines conform to a general naming scheme: `XYYZZZ` where

* `X` precision: `S`, `D`, `C`, `Z` stand for single and double, single complex and double complex, respectively.

* `YY` storage scheme: general rectangular, triangular, banded.

* `ZZZ` operation. See the manual for a list.

#### 24.2.1.3 Data formats

Lapack and Blas use a number of data formats, including

* `GE` General matrix: stored two-dimensionally as `A` (`LDA`, `*`)

* `SY/HE` Symmetric/Hermitian: general storage; `UPLO` parameter to indicate upper or lower (e.g. `SPOTRF`)

* `GB/SB/HB` General/symmetric/Hermitian band; these formats use column-major storage; in `SGBTRF` over-allocation needed because of pivoting

* `PB` Symmetric of Hermitian positive definite band; no overallocation in `SPDTRF`

#### 24.2.1.4 Lapack operations

For Lapack, we can further divide the routines into an organization with three levels:

* Drivers. These are powerful top level routine for problems such as solving linear systems or computing an SVD. There are simple and expert drivers; the expert ones have more numerical sophistication.
* Computational routines. These are the routines that drivers are built up out of<fnref n="2" />. A user may have occasion to call them by themselves.

* Auxiliary routines.

Expert driver names end on 'X'.

* Linear system solving. Simple drivers: `-SV` (e.g., DGESV) Solve $AX = B$, overwrite A with LU (with pivoting), overwrite B with X.

Expert driver: `-SVX` Also transpose solve, condition estimation, refinement, equilibration

* Least squares problems. Drivers:

* `xGELS` using QR or LQ under full-rank assumption

* `xGELSY` "complete orthogonal factorisation"

* `xGELSS` using SVD

* `xGELSD` using divide-conquer SVD (faster, but more workspace than `xGELSS`)

Also: LSE & GLM linear equality constraint & general linear model

* Eigenvalue routines. Symmetric/Hermitian: `xSY` or `xHE` (also `SP`, `SB`, `ST`)

* simple driver `-EV`

* expert driver `-EVX`

* divide and conquer `-EVD`

* relative robust representation `-EVR`

* General (only `xGE`)

* Schur decomposition `-ES` and `-ESX`

* eigenvalues `-EV` and `-EVX`

* SVD (only `xGE`)

* simple driver `-SVD`

* divide and conquer `SDD`

* Generalized symmetric (`SY` and `HE`; `SP`, `SB`)

* simple driver `GV`

* expert `GVX`

* divide-conquer `GVD`

* Nonsymmetric:

* Schur: simple `GGES`, expert `GGESX`

* eigen: simple `GGEV`, expert `GGEVX`

* svd: `GGSVD`

### 24.2.1.5 BLAS matrix storage

There are a few points to bear in mind about the way matrices are stored in the BLAS and LAPACK<fnref n="3" />:

#### 24.2.1.5.1 Array indexing
Since these libraries originated in a Fortran environment, they use 1-based indexing. Users of languages such as C/C++ are only affected by this when routines use index arrays, such as the location of pivots in LU factorizations.

<footnote n="2">2. Ha! Take that, Winston.</footnote>

<footnote n="3">3. We are not going into band storage here.</footnote>
<figure><img src="image_140.png" /><figcaption>Figure 24.1: Column-major storage of an array in Fortran</figcaption></figure>

### 24.2.1.5.2 Fortran column-major ordering
Since computer memory is one-dimensional, some conversion is needed from two-dimensional matrix coordinates to memory locations. The Fortran language uses column-major storage, that is, elements in a column are stored consecutively; see figure 24.1. This is also described informally as ‘the leftmost index varies quickest’.

Arrays in C, on the other hand, are laid out in row-major order. How to create a C array that can be handled by Blas routines is discussed in section 28.2.

### 24.2.1.6 Submatrices and the LDA parameter

Using the storage scheme described above, it is clear how to store an $m \times n$ matrix in $mn$ memory locations. However, there are many cases where software needs access to a matrix that is a subblock of another, larger, matrix. As you see in figure 24.2 such a subblock is no longer contiguous in memory. The way to describe this is by introducing a third parameter in addition to M, N: we let LDA be the ‘leading dimension of A’, that is, the allocated first dimension of the surrounding array. This is illustrated in figure 24.3. To pass the subblock to a routine, you would specify it as

call routine( A(3,2), /* M= */ 2, /* N= */ 3, /* LDA= */ Mbig, ... )

### 24.2.2 Performance issues

The collection of BLAS and LAPACK routines are a de facto standard: the Application Programmer Interface (API) is fixed, but the implementation is not. You can find reference implementations on the netlib website (netlib.org), but these will be very low in performance.

On the other hand, many LAPACK routines can be based on the matrix-matrix product (BLAS routine gemm), which you saw in section 1.6.11 has the potential for a substantial fraction of peak performance. To achieve this, you should use an optimized version, such as
<figure><img src="image_141.png" /><figcaption>Figure 24.2: A subblock out of a larger matrix</figcaption></figure>

* MKL, the Intel math-kernel library;

* OpenBlas (http://www.openblas.net/), an open source version of the original Goto BLAS; or

* blis (https://code.google.com/p/blis/), a BLAS replacement and extension project.

### 24.2.3 Some simple examples

Let's look at some simple examples.

The routine `xscal` scales a vector in place.

```fortran
! Fortran
subroutine dscal (integer N, double precision DA,
                  double precision, dimension(*) DX, integer INCX )
  // C
void cblas_dscal (const MKL_INT n, const double a,
                  double *x, const MKL_INT incx);
```

A simple example:

```fortran
// example1.F90
  do i=1,n
    xarray(i) = 1.d0
  end do
  call dscal(n,scale,xarray,1)
  do i=1,n
    if (.notenk assert_equal(xarray(i),scale)) print *, "Error in index",i
  end do
```
<figure><img src="image_142.png" /><figcaption>Figure 24.3: A subblock out of a larger matrix, using LDA</figcaption></figure>

The same in C:

```c
// example1c.cxx
xarray = new double[n]; yarray = new double[n];

for (int i=0; i<n; i++)
    xarray[i] = 1.;
cblas_dscal(n,scale,xarray,1);
for (int i=0; i<n; i++)
    if (!assert_equal(xarray[i],scale))
        printf("Error in index %d",i);
```

Many routines have an increment parameter. For `xscale` that's the final parameter:

```c
// example2.F90
integer :: inc=2
call dscal(n/inc,scale,xarray,inc)
do i=1,n
    if (mod(i,inc)==1) then
        if (.not assert_equal(xarray(i),scale)) print *,"Error in index",
            else
                if (.not assert_equal(xarray(i),1.d0)) print *,"Error in index",
                    end if
    end do
```

The matrix-vector product `xgemv` computes $y \leftarrow \alpha Ax + \beta y$, rather than $y \leftarrow Ax$. The specification of
the matrix takes the M, N size parameters, and a character argument 'N' to indicate that the matrix is not
transposed. Both of the vectors have an increment argument.

subroutine dgemv (character TRANS,
integer M, integer N,
double precision ALPHA,
double precision, dimension (lda, *) A, integer LDA,
double precision, dimension (* ) X, integer INCX,
double precision BETA, double precision, dimension (* ) Y, integer INCY
)

An example of the use of this routine:

// example3.F90
do j=1,n
  xarray(j) = 1.d0
do i=1,m
  matrix(i,j) = 1.d0
end do
end do

alpha = 1.d0; beta = 0.d0
call dgemv( 'N',M,N, alpha,matrix,M, xarray,1, beta,yarray,1)
do i=1,m
  if (.not assert_equal( yarray(i), dble(n) )) &
    print *, "Error in index",i, ":", yarray(i)
end do

The same example in C has an extra parameter to indicate whether the matrix is stored in row or column
major storage:

// example3c.cxx
for (int j=0; j<n; j++) {
  xarray[j] = 1.;
  for (int i=0; i<m; i++)
    matrix[i+j*m] = 1.;
}

alpha = 1.; beta = 0.;
cblas_dgemv(CblasColMajor,
             CblasNoTrans,m,n, alpha,matrix,m, xarray,1, beta,yarray,1);

for (int i=0; i<m; i++)
  if (!assert_equal( yarray[i], (double)n ))
    printf("Error in index %d",i);
# Chapter 25

## Plotting with GNUplot

The *gnuplot* utility is a simple program for plotting sets of points or curves. This very short tutorial will show you some of the basics. For more commands and options, see the manual http://www.gnuplot.info/docs/gnuplot.html.

### 25.1 Usage modes

The two modes for running *gnuplot* are *interactive* and *from file*. In interactive mode, you call *gnuplot* from the command line, type commands, and watch output appear; you terminate an interactive session with *quit*. If you want to save the results of an interactive session, do *save "name.plt"*. This file can be edited, and loaded with *load "name.plt"*.

Plotting non-interactively, you call *gnuplot <your file>*.

The output of *gnuplot* can be a picture on your screen, or drawing instructions in a file. Where the output goes depends on the setting of the *terminal*. By default, *gnuplot* will try to draw a picture. This is equivalent to declaring

```
set terminal x11
```

or *aqua*, *windows*, or any choice of graphics hardware.

For output to file, declare

```
set terminal pdf
```

or *fig*, *latex*, *pbm*, et cetera. Note that this will only cause the pdf commands to be written to your screen: you need to direct them to file with

```
set output "myplot.pdf"
```

or capture them with

```
gnuplot my.plt > myplot.pdf
```
## 25.2 Plotting

The basic plot commands are `plot` for 2D, and `splot` ('surface plot') for 3D plotting.

### 25.2.1 Plotting curves

By specifying

`plot x**2`

you get a plot of $f(x) = x^2$; `gnuplot` will decide on the range for $x$. With

`set xrange [0:1]`

`plot 1-x title "down", x**2 title "up"`

you get two graphs in one plot, with the $x$ range limited to $[0, 1]$, and the appropriate legends for the graphs. The variable $x$ is the default for plotting functions.

Plotting one function against another – or equivalently, plotting a parametric curve – goes like this:

`set parametric`

`plot [t=0:1.57] cos(t),sin(t)`

which gives a quarter circle.

To get more than one graph in a plot, use the command `set multiplot`.

### 25.2.2 Plotting data points

It is also possible to plot curves based on data points. The basic syntax is `plot 'datafile'`, which takes two columns from the data file and interprets them as $(x, y)$ coordinates. Since data files can often have multiple columns of data, the common syntax is `plot 'datafile' using 3:6` for columns 3 and 6. Further qualifiers like `with lines` indicate how points are to be connected.

Similarly, `splot "datafile3d.dat" 2:5:7` will interpret three columns as specifying $(x, y, z)$ coordinates for a 3D plot.

If a data file is to be interpreted as level or height values on a rectangular grid, do `splot "matrix.dat" matrix` for data points; connect them with

`split "matrix.dat" matrix with lines`

### 25.2.3 Customization

Plots can be customized in many ways. Some of these customizations use the `set` command. For instance,
```matlab
set xlabel "time"
set ylabel "output"
set title "Power curve"
```

You can also change the default drawing style with

```matlab
set style function dots
```

(dots, lines, dots, points, et cetera), or change on a single plot with

```matlab
plot f(x) with points
```

## 25.3 Workflow

Imagine that your code produces a dataset that you want to plot, and you run your code for a number of inputs. It would be nice if the plotting can be automated. Gnuplot itself does not have the facilities for this, but with a little help from shell programming this is not hard to do.

Suppose you have data files

```
data1.dat data2.dat data3.dat
```

and you want to plot them with the same gnuplot commands. You could make a file `plot.template`:

```matlab
set term pdf
set output "FILENAME.pdf"
plot "FILENAME.dat"
```

The string `FILENAME` can be replaced by the actual file names using, for instance `sed`:

```bash
for d in data1 data2 data3 ; do
  cat plot.template | sed s/FILENAME/$d/ > plot.cmd
  gnuplot plot.cmd
done
```

Variations on this basic idea are many.
# Chapter 26

## Good coding practices

Sooner or later, and probably sooner than later, every programmer is confronted with code not behaving as intended. In this section you will learn some techniques of dealing with this problem. At first we will see a number of techniques for preventing errors; in the next chapter we will discuss debugging, the process of finding the inevitable errors in a program, once they have occurred.

### 26.1 Defensive programming

In this section we will discuss a number of techniques that are aimed at preventing the likelihood of programming errors, or increasing the likelihood of them being found at runtime. We call this *defensive programming*.

Scientific codes are often large and involved, so it is a good practice to code knowing that you are going to make mistakes and prepare for them. Another good coding practice is the use of tools: there is no point in reinventing the wheel if someone has already done it for you. Some of these tools are be described in other sections:

* Build systems, such as Make, Scons, Bjam; see section 21.

* Source code management with SVN, Git; see section 22.

* Regression testing and designing with testing in mind (unit testing)

First we will have a look at runtime sanity checks, where you test for things that can not or should not happen.

#### 26.1.1 Assertions

In the things that can go wrong with a program we can distinguish between errors and bugs. Errors are things that legitimately happen but that should not. File systems are common sources of errors: a program wants to open a file but the file doesn't exist because the user mistyped the name, or the program writes to a file but the disk is full. Other errors can come from arithmetic, such as *overflow errors*.

On the other hand, a *bug* in a program is an occurrence that cannot legitimately occur. Of course, 'legitimately' here means 'according to the programmer's intentions'. Bugs can often be described as 'the computer always does what you ask, not necessarily what you want'.