from the problem of figuring out which operations are independent, the main problem is that the problem of locating data in a parallel context is very hard. A compiler would need to consider the whole code, rather than a subroutine at a time. Even then, results have been disappointing.

More productive is the approach where the user writes mostly a sequential program, but gives some indications about what computations can be parallelized, and how data should be distributed. Indicating parallelism of operations explicitly is done in OpenMP (section 2.7.2); indicating the data distribution and leaving parallelism to the compiler and runtime is the basis for PGAS languages (section 2.7.5). Such approaches work best with shared memory.

By far the hardest way to program in parallel, but with the best results in practice, is to expose the parallelism to the programmer and let the programmer manage everything explicitly. This approach is necessary in the case of distributed memory programming. We will have a general discussion of distributed programming in section 2.7.3.1; section 2.7.3.3 will discuss the MPI library.

### 2.7.1 Thread parallelism

As a preliminary to OpenMP (section 2.7.2), we will briefly go into 'threads'.

To explain what a *thread* is, we first need to get technical about what a *process* is. A unix process corresponds to the execution of a single program. Thus, it has in memory:

* The program code, in the form of machine language instructions;

* A *heap*, containing for instance arrays that were created with `malloc`;

* A stack with quick-changing information, such as the *program counter* that indicates what instruction is currently being executed, and data items with local scope, as well as intermediate results from computations.

This process can have multiple threads; these are similar in that they see the same program code and heap, but they have their own stack. Thus, a thread is an independent 'strand' of execution through a process.

Processes can belong to different users, or be different programs that a single user is running concurrently, so they have their own data space. On the other hand, threads are part of one process and therefore share the process heap. Threads can have some private data, for instance by have their own data stack, but their main characteristic is that they can collaborate on the same data.

#### 2.7.1.1 The fork-join mechanism

Threads are dynamic, in the sense that they can be created during program execution. (This is different from the MPI model, where every processor run one process, and they are all created and destroyed at the same time.) When a program starts, there is one thread active: the *master thread*. Other threads are created by *thread spawning*, and the master thread can wait for their completion. This is known as the *fork-join model*; it is illustrated in figure 2.10. A group of threads that is forked from the same thread and active simultaneously is known as a *thread team*.
<figure><img src="image_41.png" /><figcaption>Figure 2.10: Thread creation and deletion during parallel execution</figcaption></figure>

### 2.7.1.2 Threads and parallel computing

Since threads share data, they are one possible mechanism for parallel programming. Of course, for this to be faster than sequential single-threaded execution, this assumes that the hardware can work on more than one thread at the same time. This can happen in two main circumstances: a multicore processor can accommodate one thread per core, and secondly, some cores have hardware support for executing more than one thread.

The thread mechanism has long existed, even on a single processor. By having more than one thread on a single processor, a higher processor utilization can result, since the instructions of one thread can be processed while another thread is waiting for data. On traditional CPUs, switching between threads is somewhat expensive (an exception is the *hyper-threading* mechanism) but on GPUs it is not, and in fact they *need* many threads to attain high performance.

In the context of parallel processing we are also interested in threads, since in a shared memory context multiple threads running on multiple processors or processor cores can be an easy way to parallelize a computation. The shared memory allows the threads to all see the same data. This can also lead to problems; see section 2.7.1.5.

### 2.7.1.3 Threads example

The following example<fnref n="6" /> is a clear illustration of the *fork-join* model. It uses the `pthread` library to spawn a number of tasks that all update a global counter. Since threads share the same memory space, they indeed see and update the same memory location.

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

int sum=0;

void adder() {
    sum = sum+1;
    return;
}
```

<footnote n="6">6. This is strictly Unix-centric and will not work on Windows.</footnote>
```c
#define NTHREADS 50
int main() {
    int i;
    pthread_t threads[NTHREADS];
    printf("forking\n");
    for (i=0; i<NTHREADS; i++)
        if (pthread_create(threads+i, NULL, &adder, NULL) !=0) return i+1;
    printf("joining\n");
    for (i=0; i<NTHREADS; i++)
        if (pthread_join(threads[i], NULL) !=0) return NTHREADS+i+1;
    printf("Sum computed: %d\n", sum);

    return 0;
}
```

The fact that this code gives the right result is a coincidence: it only happens because updating the variable is so much quicker than creating the thread. (On a multicore processor the chance of errors will greatly increase.) If we artificially increase the time for the update, we will no longer get the right result:

```c
void adder() {
    int t = sum; sleep(1); sum = t+1;
    return;
}
```

Now all threads read out the value of `sum`, wait a while (presumably calculating something) and then update.

This can be fixed by having a *lock* on the code region that should be 'mutually exclusive':

```c
pthread_mutex_t lock;

void adder() {
    int t,r;
    pthread_mutex_lock(&lock);
    t = sum; sleep(1); sum = t+1;
    pthread_mutex_unlock(&lock);
    return;
}

int main() {
    ...
    pthread_mutex_init(&lock, NULL);
```
The lock and unlock commands guarantee that no two threads can interfere with each other's update.

For more information on pthreads, see for instance https://computing.llnl.gov/tutorials/pthreadsguide.

### 2.7.1.4 Contexts

In the above example and its version with the `sleep` command we glanced over the fact that there were two types of data involved. First of all, the variable `s` was created outside the thread spawning part. Thus, this variable was *shared*.

On the other hand, the variable `t` was created once in each spawned thread. We call this *private data*.

The totality of all data that a thread can access is called its *context*. It contains private and shared data, as well as temporary results of computations that the thread is working on[^7].

It is quite possible to create more threads than a processor has cores, so a processor may need to switch between the execution of different threads. This is known as a *context switch*.

Context switches are not for free on regular CPUs, so they only pay off if the *granularity* of the threaded work is high enough. The exceptions to this story are:

* CPUs that have hardware support for multiple threads, for instance through *hyperthreading* (section 2.7.1.10), or as in the *Intel Xeon Phi* (section 2.10);

* GPUs, which in fact rely on fast context switching (section 2.10.3);

* certain other 'exotic' architectures such as the *Cray XMT* (section 2.9).

### 2.7.1.5 Data races, thread safety, sequential consistency, and atomic operations

Shared memory makes life easy for the programmer, since every processor has access to all of the data: no explicit data traffic between the processor is needed. On the other hand, multiple processes/processors can also write to the same variable, which is a source of potential problems.

Suppose that two processes both try to increment an integer variable `I`:

process 1: `I=I+2`

process 2: `I=I+3`

This is a legitimate activity if the variable is an accumulator for values computed by independent processes. The result of these two updates depends on the sequence in which the processors read and write the variable. Here are three scenarios:

<footnote>[^7]: It also contains the program counter and stack pointer. If you don't know what those are, don't worry.</footnote>
<table><thead><tr><th colspan="2">scenario 1.</th><th colspan="2">scenario 2.</th><th colspan="2">scenario 3.</th></tr></thead><tbody><tr><td colspan="6" style="text-align: center;">I = 0</td></tr><tr><td>read I = 0</td><td>read I = 0</td><td>read I = 0</td><td>read I = 0</td><td>read I = 0</td><td>read I = 0</td></tr><tr><td>compute I = 2</td><td>compute I = 3</td><td>compute I = 2</td><td>compute I = 3</td><td>compute I = 2</td><td>compute I = 2</td></tr><tr><td>write I = 2</td><td></td><td>write I = 3</td><td></td><td>write I = 2</td><td></td></tr><tr><td></td><td></td><td></td><td>write I = 2</td><td></td><td>read I = 2</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>compute I = 5</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>write I = 5</td></tr><tr><td>I = 3</td><td></td><td>I = 2</td><td></td><td>I = 5</td><td></td></tr></tbody></table>

Such a scenario, where the final result depends on which thread executes first, is known as a *race condition* or *data race*.

A very practical example of such conflicting updates is the inner product calculation:

```c
for (i=0; i<1000; i++)
    sum = sum + a[i] * b[i];
```

Here the products are truly independent, so we could choose to have the loop iterations do them in parallel, for instance by their own threads. However, all threads need to update the same variable `sum`.

Code that behaves the same whether it's executed sequentially or threaded is called *thread safe*. As you can see from the above examples, a lack of thread safety is typically due to the treatment of shared data. This implies that the more your program uses local data, the higher the chance that it is thread safe. Unfortunately, sometimes the threads need to write to shared/global data, for instance when the program does a *reduction*.

There are essentially two ways of solving this problem. One is that we declare such updates of a shared variable a *critical section* of code. This means that the instructions in the critical section (in the inner product example 'read sum from memory, update it, write back to memory') can be executed by only one thread at a time. In particular, they need to be executed entirely by one thread before any other thread can start them so the ambiguity problem above will not arise. Of course, the above code fragment is so common that systems like OpenMP (section 2.7.2) have a dedicated mechanism for it, by declaring it a *reduction* operation.

Critical sections can for instance be implemented through the *semaphore* mechanism [39]. Surrounding each critical section there will be two atomic operations controlling a semaphore, a sign post. The first process to encounter the semaphore will lower it, and start executing the critical section. Other processes see the lowered semaphore, and wait. When the first process finishes the critical section, it executes the second instruction which raises the semaphore, allowing one of the waiting processes to enter the critical section.

The other way to resolve common access to shared data is to set a temporary *lock* on certain memory areas. This solution may be preferable, if common execution of the critical section is likely, for instance if it implements writing to a database or hash table. In this case, one process entering a critical section would prevent any other process from writing to the data, even if they might be writing to different locations; locking the specific data item being accessed is then a better solution.
The problem with locks is that they typically exist on the operating system level. This means that they are relatively slow. Since we hope that iterations of the inner product loop above would be executed at the speed of the floating point unit, or at least that of the memory bus, this is unacceptable.

One implementation of this is *transactional memory*, where the hardware itself supports atomic operations; the term derives from database transactions, which have a similar integrity problem. In *transactional memory*, a process will perform a normal memory update, unless the processor detects a conflict with an update from another process. In that case, the updates ('transactions') are aborted and retried with one processor locking the memory and the other waiting for the lock. This is an elegant solution; however, aborting the transaction may carry a certain cost of *pipeline flushing* (section 1.2.5) and cache line invalidation (section 1.4.1).

### 2.7.1.6 Memory models and sequential consistency

The above signaled phenomenon of a *race condition* means that the result of some programs can be non-deterministic, depending on the sequence in which instructions are executed. There is a further factor that comes into play, and which is called the *memory model* that a processor uses. The memory model controls how the activity of one thread or core is seen by other threads or cores.

As an example, consider

initially: $A=B=0$; then

process 1: $A=1$; $x = B$;

process 2: $B=1$; $y = A$;

As above, we have three scenarios, which we describe by giving a global sequence of statements:

<table><thead><tr><th>scenario 1.</th><th>scenario 2.</th><th>scenario 3.</th></tr></thead><tbody><tr><td>A &leftarrow; 1</td><td>A &leftarrow; 1</td><td>B &leftarrow; 1</td></tr><tr><td>x &leftarrow; B</td><td>B &leftarrow; 1</td><td>y &leftarrow; A</td></tr><tr><td>B &leftarrow; 1</td><td>x &leftarrow; B</td><td>A &leftarrow; 1</td></tr><tr><td>y &leftarrow; A</td><td>y &leftarrow; A</td><td>x &leftarrow; B</td></tr><tr><td>x = 0, y = 1</td><td>x = 1, y = 1</td><td>x = 1, y = 0</td></tr></tbody></table>

(In the second scenario, statements 1,2 can be reversed, as can 3,4, without change in outcome.)

The three different outcomes can be characterized as being computed by a global ordering on the statements that respects the local orderings. This is known as *sequential consistency*: the parallel outcome is consistent with a sequential execution that interleaves the parallel computations, respecting their local statement orderings.

Maintaining sequential consistency is expensive: it means that any change to a variable immediately needs to be visible on all other threads, or that any access to a variable on a thread needs to consult all other threads. We discussed this in section 1.4.1.

In a *relaxed memory model* it is possible to get a result that is not sequentially consistent. Suppose, in the above example, the instructions $A \leftarrow 1$, $B \leftarrow 1$ are executed first, but these changes are not immediately visible to the other threads. Then subsequent instructions $x \leftarrow B$, $y \leftarrow A$ can receive the result $x = 0, y = 0$, which was not possible under the sequentially consistent model above.
Sequential consistency implies that

```
integer n
n = 0
!$omp parallel shared(n)
n = n + 1
!$omp end parallel
```

should have the same effect as

```
n = 0
n = n+1 ! for processor 0
n = n+1 ! for processor 1
! et cetera
```

With sequential consistency it is no longer necessary to declare atomic operations or critical sections; however, this puts strong demands on the implementation of the model, so it may lead to inefficient code.

**Exercise 2.15.** In section 2.7.1.5 you saw an example that needed a critical section to get the right final result. Argue that having a critical section is not enough for sequential consistency.

* Write a piece of sequential code that, when executed in parallel, corresponds to the example in section 2.7.1.5.

* Show that, using a critical section, there are two execution orderings that give the correct result.

* Show that one of these orderings is not sequentially consistent.

### 2.7.1.7 Affinity

Thread programming is very flexible, effectively creating parallelism as needed. However, a large part of this book is about the importance of data movement in scientific computations, and that aspect can not be ignored in thread programming.

In the context of a multicore processor, any thread can be scheduled to any core, and there is no immediate problem with this. The problem in cases where we have two subsequent regions that are handled with thread-level parallelism. If the first region computes data, and the second uses that data, then we have to make sure that the thread producing a particular data item and the thread consuming (or updating) it are scheduled to the same core.

We call *affinity* the maintaining of a fixed mapping between threads (*thread affinity*) or processes (*process affinity*) and cores.

Maintaining thread affinity is easy in some cases. If the loop structure that is being parallelized stays the same, a fixed mapping of threads to cores makes sense:

```
for (i=0; i<ndata; i++) // this loop will be done by threads
    x[i] = ....
for (i=0; i<ndata; i++) // as will this one
```
... = ... x[i] ...

In other cases a fixed mapping is not the right solution:

```c
for (i=0; i<ndata; i++) // produces loop
    x[i] = ....
for (i=0; i<ndata; i+=2) // use even indices
    ... = ... x[i] ...
for (i=1; i<ndata; i+=2) // use odd indices
    ... = ... x[i] ...
```

In this second example, either the program has to be transformed, or the programmer has to maintain in effect a *task queue*.

It is natural to think of affinity to think in terms of 'put the execution where the data is'. However, in practice the opposite view sometimes makes sense. For instance, figure 2.7 showed how the shared memory of a cluster node can actually be distributed. Thus, a thread can be attached to a socket, but data can be allocated by the Operating System (OS) on any of the sockets. The mechanism that is often used by the OS is called the *first-touch policy*:

* When the program allocates data, the OS does not actually create it;

* instead, the memory area for the data is created the first time a thread accesses it;

* thus, the first thread to touch the area in effect causes the data to be allocated on the memory of its socket.

Exercise 2.16. Explain the problem with the following code:

```c
// serial initialization
for (i=0; i<N; i++)
    a[i] = 0.;

#pragma omp parallel for
for (i=0; i<N; i++)
    a[i] = b[i] + c[i];
```

### 2.7.1.8 Concurrency

In this book, we discuss threads in the context of scientific computing. However, threads also facilitate parallelism in other contexts, such as Operating Systems (OSs). In this case, the term *concurrency* is often used, rather than 'parallelism'. Concurrent tasks are tasks that are not dependent on each other; therefore, they can be executed in any order, or indeed in parallel. Parallel tasks, on the other hand, are tasks that are indeed running simultaneously.

Concurrency has been studied for a long time. Since a lot of this research came from OS research, the question addressed is that of *resource contention*: what if two processes, for instance users on the same system, request access to the same resource, for instance a printer. This means that the system call to print is another example of a critical section of code; see section 2.7.1.5 above. Another problem with concurrency is that two processes can both request access to two devices with exclusive use, for instance
a printer and some hardware input device. If one process gets hold of the one resource first, and the other process the other resource, we have another case of deadlock.

### 2.7.1.9 Cilk Plus

Other programming models based on threads exist. For instance, Intel Cilk Plus (http://www.cilkplus.org/) is a set of extensions of C/C++ with which a programmer can create threads.

```
Sequential code:
    int fib(int n) {
        if (n<2) return 1;
        else {
            int rst=0;
            rst += fib(n-1);
            rst += fib(n-2);
            return rst;
        }
    }

Cilk code:
    cilk int fib(int n) {
        if (n<2) return 1;
        else {
            int rst=0;
            rst += cilk_spawn fib(n-1);
            rst += cilk_spawn fib(n-2);
            cilk_sync;
            return rst;
        }
    }
```

In this example, the variable `rst` is updated by two, potentially independent threads. The semantics of this update, that is, the precise definition of how conflicts such as simultaneous writes are resolved, is defined by sequential consistency; see section 2.7.1.5.

### 2.7.1.10 Hyperthreading versus multi-threading

In the above examples you saw that the threads that are spawned during one program run essentially execute the same code, and have access to the same data. Thus, at a hardware level, a thread is uniquely determined by a small number of local variables, such as its location in the code (the *program counter*) and intermediate results of the current computation it is engaged in.

Hyperthreading is an Intel technology to let multiple threads use the processor truly simulateneously, so that part of the processor would be optimally used.

If a processor switches between executing one thread and another, it saves this local information of the one thread, and loads the information of the other. The cost of doing this is modest compared to running a whole program, but can be expensive compared to the cost of a single instruction. Thus, hyperthreading may not always give a performance improvement.

Certain architectures have support for *multi-threading*. This means that the hardware actually has explicit storage for the local information of multiple threads, and switching between the threads can be very fast. This is the case on GPUs (section 2.10.3), and on the *Intel Xeon Phi* architecture, where each core can support up to four threads.

## 2.7.2 OpenMP

OpenMP is an extension to the programming languages C and Fortran. Its main approach to parallelism is the parallel execution of loops: based on *compiler directives*, a preprocessor can schedule the parallel
execution of the loop iterations.

Since OpenMP is based on *threads*, it features *dynamic parallelism*: the number of execution streams operating in parallel can vary from one part of the code to another. Parallelism is declared by creating parallel regions, for instance indicating that all iterations of a loop nest are independent, and the runtime system will then use whatever resources are available.

OpenMP is not a language, but an extension to the existing C and Fortran languages. It mostly operates by inserting directives into source code, which are interpreted by the compiler. It also has a modest number of library calls, but these are not the main point, unlike in MPI (section 2.7.3.3). Finally, there is a runtime system that manages the parallel execution.

OpenMP has an important advantage over MPI in its programmability: it is possible to start with a sequential code and transform it by *incremental parallelization*. By contrast, turning a sequential code into a distributed memory MPI program is an all-or-nothing affair.

Many compilers, such as GCC or the Intel compiler, support the OpenMP extensions. In Fortran, OpenMP directives are placed in comment statements; in C, they are placed in `#pragma` CPP directives, which indicate compiler specific extensions. As a result, OpenMP code still looks like legal C or Fortran to a compiler that does not support OpenMP. Programs need to be linked to an OpenMP runtime library, and their behaviour can be controlled through environment variables.

For more information about OpenMP, see [25] and http://openmp.org/wp/.

### 2.7.2.1 OpenMP examples

The simplest example of OpenMP use is the parallel loop.

```c
#pragma omp for
for (i=0; i<ProblemSize; i++) {
    a[i] = b[i];
}
```

Clearly, all iterations can be executed independently and in any order. The pragma CPP directive then conveys this fact to the compiler.

Some loops are fully parallel conceptually, but not in implementation:

```c
for (i=0; i<ProblemSize; i++) {
    t = b[i]*b[i];
    a[i] = sin(t) + cos(t);
}
```

Here it looks as if each iteration writes to, and reads from, a shared variable *t*. However, *t* is really a temporary variable, local to each iteration. Code that should be parallelizable, but is not due to such constructs, is called not *thread safe*.

OpenMP indicates that the temporary is private to each iteration as follows:
<figure><img src="image_42.png" /><figcaption>Figure 2.11: Static or round-robin (left) vs dynamic (right) thread scheduling; the task numbers are indicated.</figcaption></figure>

```cpp
#pragma parallel for shared(a,b), private(t)
for (i=0; i<ProblemSize; i++) {
    t = b[i]*b[i];
    a[i] = sin(t) + cos(t);
}
```

If a scalar *is* indeed shared, OpenMP has various mechanisms for dealing with that. For instance, shared variables commonly occur in *reduction operations*:

```cpp
s = 0;

#pragma parallel for reduction(+:sum)
for (i=0; i<ProblemSize; i++) {
    s = s + a[i]*b[i];
}
```

As you see, a sequential code can be easily parallelized this way.

The assignment of iterations to threads is done by the runtime system, but the user can guide this assignment. We are mostly concerned with the case where there are more iterations than threads: if there are *P* threads and *N* iterations and *N* > *P*, how is iteration *i* going to be assigned to a thread?

The simplest assignment uses *round-robin task scheduling*, a *static scheduling* strategy where thread *p* get iterations *p* × (N/P), ..., (*p* + 1) × (N/P) − 1. This has the advantage that if some data is reused between iterations, it will stay in the data cache of the processor executing that thread. On the other hand, if the iterations differ in the amount of work involved, the process may suffer from *load unbalance* with static scheduling. In that case, a *dynamic scheduling* strategy would work better, where each thread starts work on the next unprocessed iteration as soon as it finishes its current iteration. The example in figure 2.11 shows static versus dynamic scheduling of a number of tasks that gradually decrease in individual running time. In static scheduling, the first thread gets tasks 1 and 4, the second 2 and 5, etcetera. In dynamic scheduling, any thread that finishes its task gets the next task. This clearly gives a better running time in this particular example. On the other hand, dynamic scheduling is likely to have a higher overhead.

You can control OpenMP scheduling of loop iterations with the `schedule` keyword; its values include *static* and *dynamic*. It is also possible to indicate a *chunksize*, which controls the size of the block of iterations that gets assigned together to a thread. If you omit the chunksize, OpenMP will divide the iterations into as many blocks as there are threads.
Exercise 2.17. Let's say there are t threads, and your code looks like

```c
for (i=0; i<N; i++) {
    a[i] = // some calculation
}
```

If you specify a chunksize of 1, iterations 0, t, 2t, ... go to the first thread, 1, 1+t, 1+2t, ... to the second, etcetera. Discuss why this is a bad strategy from a performance point of view. Hint: look up the definition of false sharing. What would be a good chunksize?

### 2.7.3 Distributed memory programming through message passing

While OpenMP programs, and programs written using other shared memory paradigms, still look very much like sequential programs, this does not hold true for message passing code. Before we discuss the Message Passing Interface (MPI) library in some detail, we will take a look at this shift the way parallel code is written.

#### 2.7.3.1 The global versus the local view in distributed programming

There can be a marked difference between how a parallel algorithm looks to an observer, and how it is actually programmed. Consider the case where we have an array of processors {$P_i$}$_{i=0..p-1}$, each containing one element of the arrays x and y, and $P_i$ computes

$$ \begin{cases} y_i \leftarrow y_i + x_{i-1} & i > 0 \\ y_i \text{ unchanged} & i = 0 \end{cases} \qquad (2.4) $$

The global description of this could be

* Every processor $P_i$ except the last sends its x element to $P_{i+1}$;

* every processor $P_i$ except the first receive an x element from their neighbour $P_{i-1}$, and

* they add it to their y element.

However, in general we can not code in these global terms. In the SPMD model (section 2.4.2) each processor executes the same code, and the overall algorithm is the result of these individual behaviours. The local program has access only to local data – everything else needs to be communicated with send and receive operations – and the processor knows its own number.

One possible way of writing this would be

* If I am processor 0 do nothing, otherwise receive a y element from the left, add it to my x element.

* If I am the last processor do nothing, otherwise send my y element to the right.

At first we look at the case where sends and receives are so-called blocking communication instructions: a send instruction does not finish until the sent item is actually received, and a receive instruction waits for the corresponding send. This means that sends and receives between processors have to be carefully paired. We will now see that this can lead to various problems on the way to an efficient code.

The above solution is illustrated in figure 2.12, where we show the local timelines depicting the local
<figure><img src="image_43.png" /><figcaption>Figure 2.12: Local and resulting global view of an algorithm for sending data to the right</figcaption></figure>

<figure><img src="image_44.png" /><figcaption>Figure 2.13: Local and resulting global view of an algorithm for sending data to the right</figcaption></figure>

processor code, and the resulting global behaviour. You see that the processors are not working at the same time: we get serialized execution.

What if we reverse the send and receive operations?

* If I am not the last processor, send my *x* element to the right;

* If I am not the first processor, receive an *x* element from the left and add it to my *y* element.

This is illustrated in figure 2.13 and you see that again we get a serialized execution, except that now the processor are activated right to left.

If the algorithm in equation 2.4 had been cyclic:

$$ \begin{cases} y_i \leftarrow y_i + x_{i-1} & i = 1 \dots n-1 \\ y_0 \leftarrow y_0 + x_{n-1} & i = 0 \end{cases} \qquad (2.5) $$
<figure><img src="image_45.png" /><figcaption>Figure 2.14: Local and resulting global view of an algorithm for sending data to the right</figcaption></figure>

the problem would be even worse. Now the last processor can not start its receive since it is blocked
sending $x_{n-1}$ to processor 0. This situation, where the program can not progress because every processor
is waiting for another, is called deadlock.

The solution to getting an efficient code is to make as much of the communication happen simultaneously
as possible. After all, there are no serial dependencies in the algorithm. Thus we program the algorithm as
follows:

* If I am an odd numbered processor, I send first, then recieve;

* If I am an even numbered processor, I receive first, then send.

This is illustrated in figure 2.14, and we see that the execution is now parallel.

Exercise 2.18. Take another look at figure 2.3 of a parallel reduction. The basic actions are:

* receive data from a neighbour

* add it to your own data

* send the result on.

As you see in the diagram, there is at least one processor who does not send data on,
and others may do a variable number of receives before they send their result on.

Write node code so that an SPMD program realizes the distributed reduction. Hint: write each processor number in binary. The algorithm uses a number of steps that is equal to the length of this bitstring.

* Assuming that a processor receives a message, express the distance to the origin of that message in the step number.

* Every processor sends at most one message. Express the step where this happens in terms of the binary processor number.

### 2.7.3.2 Blocking and non-blocking communication

The reason for blocking instructions is to prevent accumulation of data in the network. If a send instruction were to complete before the corresponding receive started, the network would have to store the data somewhere in the mean time. Consider a simple example:

```
buffer = ... ; // generate some data
send(buffer, 0); // send to processor 0
```
```c
buffer = ...; // generate more data
send(buffer, 1); // send to processor 1
```

After the first send, we start overwriting the buffer. If the data in it hasn't been received, the first set of values would have to be buffered somewhere in the network, which is not realistic. By having the send operation block, the data stays in the sender's buffer until it is guaranteed to have been copied to the recipient's buffer.

One way out of the problem of sequentialization or deadlock that arises from blocking instruction is the use of non-blocking communication instructions, which include explicit buffers for the data. With non-blocking send instruction, the user needs to allocate a buffer for each send, and check when it is safe to overwrite the buffer.

```c
buffer0 = ...; // data for processor 0
send(buffer0, 0); // send to processor 0
buffer1 = ...; // data for processor 1
send(buffer1, 1); // send to processor 1
...
// wait for completion of all send operations.
```

### 2.7.3.3 The MPI library

If OpenMP is the way to program shared memory, Message Passing Interface (MPI) [148] is the standard solution for programming distributed memory. MPI ('Message Passing Interface') is a specification for a library interface for moving data between processes that do not otherwise share data. The MPI routines can be divided roughly in the following categories:

* Process management. This includes querying the parallel environment and constructing subsets of processors.

* Point-to-point communication. This is a set of calls where two processes interact. These are mostly variants of the send and receive calls.

* Collective calls. In these routines, all processors (or the whole of a specified subset) are involved. Examples are the *broadcast* call, where one processor shares its data with every other processor, or the *gather* call, where one processor collects data from all participating processors.

Let us consider how the OpenMP examples can be coded in MPI<fnref n="8" />. First of all, we no longer allocate

```c
double a[ProblemSize];
```

but

```c
double a[LocalProblemSize];
```

where the local size is roughly a 1/P fraction of the global size. (Practical considerations dictate whether you want this distribution to be as evenly as possible, or rather biased in some way.)

<footnote n="8">8. This is not a course in MPI programming, and consequently the examples will leave out many details of the MPI calls. If you want to learn MPI programming, consult for instance [75].</footnote>
The parallel loop is trivially parallel, with the only difference that it now operates on a fraction of the arrays:

```c
for (i=0; i<LocalProblemSize; i++) {
    a[i] = b[i];
}
```

However, if the loop involves a calculation based on the iteration number, we need to map that to the global value:

```c
for (i=0; i<LocalProblemSize; i++) {
    a[i] = b[i] + f(i + MyFirstVariable);
}
```

(We will assume that each process has somehow calculated the values of `LocalProblemSize` and `MyFirstVariable`.) Local variables are now automatically local, because each process has its own instance:

```c
for (i=0; i<LocalProblemSize; i++) {
    t = b[i] * b[i];
    a[i] = sin(t) + cos(t);
}
```

However, shared variables are harder to implement. Since each process has its own data, the local accumulation has to be explicitly assembled:

```c
for (i=0; i<LocalProblemSize; i++) {
    s = s + a[i] * b[i];
}
MPI_Allreduce(s, globals, 1, MPI_DOUBLE, MPI_SUM);
```

The 'reduce' operation sums together all local values `s` into a variable `globals` that receives an identical value on each processor. This is known as a *collective operation*.

Let us make the example slightly more complicated:

```c
for (i=0; i<ProblemSize; i++) {
    if (i==0)
        a[i] = (b[i]+b[i+1])/2
    else if (i==ProblemSize-1)
        a[i] = (b[i]+b[i-1])/2
    else
        a[i] = (b[i]+b[i-1]+b[i+1])/3
}
```

If we had shared memory, we could write the following parallel code:

```c
for (i=0; i<LocalProblemSize; i++) {
    bleft = b[i-1];
    bright = b[i+1];
```
a[i] = (b[i]+bleft+bright)/3

To turn this into valid distributed memory code, first we account for the fact that `bleft` and `bright` need to be obtained from a different processor for `i==0` (`bleft`), and for `i==LocalProblemSize-1` (`bright`). We do this with a exchange operation with our left and right neighbour processor:

```c
// get bfromleft and bfromright from neighbour processors, then
for (i=0; i<LocalProblemSize; i++) {
    if (i==0) bleft=bfromleft;
    else bleft = b[i-1]
    if (i==LocalProblemSize-1) bright=bfromright;
    else bright = b[i+1];
    a[i] = (b[i]+bleft+bright)/3
```

Obtaining the neighbour values is done as follows. First we need to ask our processor number, so that we can start a communication with the processor with a number one higher and lower.

```c
MPI_Comm_rank(MPI_COMM WORLD, &myTaskID);
MPI Sendrecv
    (/* to be sent: */ &b[LocalProblemSize-1],
     /* destination */ myTaskID+1,
     /* to be recvd: */ &bfromleft,
     /* source: */ myTaskID-1,
     /* some parameters omitted */
    );
MPI Sendrecv(&b[0], myTaskID-1,
    &bfromright, /* ... */ );
```

There are still two problems with this code. First, the `sendrecv` operations need exceptions for the first and last processors. This can be done elegantly as follows:

```c
MPI.Comm_rank(MPI COMM WORLD, &myTaskID);
MPI.Comm_size(MPI COMM WORLD, &nTasks);
if (myTaskID==0) leftproc = MPI_PROC_NULL;
else leftproc = myTaskID-1;
if (myTaskID==nTasks-1) rightproc = MPI_PROC_NULL;
else rightproc = myTaskID+1;
MPI Sendrecv( &b[LocalProblemSize-1], &bfromleft, rightproc );
MPI Sendrecv( &b[0], &bfromright, leftproc);
```

**Exercise 2.19.** There is still a problem left with this code: the boundary conditions from the original, global, version have not been taken into account. Give code that solves that problem.

MPI gets complicated if different processes need to take different actions, for example, if one needs to send data to another. The problem here is that each process executes the same executable, so it needs to contain
both the send and the receive instruction, to be executed depending on what the rank of the process is.

```c
if (myTaskID==0) {
    MPI Send(myInfo, 1, MPI_INT, /* to: */ 1, /* labeled: */ 0,
              MPI_COMM WORLD);
} else {
    MPI Recv(myInfo, 1, MPI_INT, /* from: */ 0, /* labeled: */ 0,
              /* not explained here: */ &status, MPI_COMM WORLD);
}
```

### 2.7.3.4 Blocking

Although MPI is sometimes called the ‘assembly language of parallel programming’, for its perceived difficulty and level of explicitness, it is not all that hard to learn, as evinced by the large number of scientific codes that use it. The main issues that make MPI somewhat intricate to use, are buffer management and blocking semantics.

These issues are related, and stem from the fact that, ideally, data should not be in two places at the same time. Let us briefly consider what happens if processor 1 sends data to processor 2. The safest strategy is for processor 1 to execute the send instruction, and then wait until processor 2 acknowledges that the data was successfully received. This means that processor 1 is temporarily blocked until processor 2 actually executes its receive instruction, and the data has made its way through the network. This is the standard behaviour of the `MPI Send` and `MPI Recv` calls, which are said to use *blocking communication*.

Alternatively, processor 1 could put its data in a buffer, tell the system to make sure that it gets sent at some point, and later checks to see that the buffer is safe to reuse. This second strategy is called *non-blocking communication*, and it requires the use of a temporary buffer.

### 2.7.3.5 Collective operations

In the above examples, you saw the `MPI_Allreduce` call, which computed a global sum and left the result on each processor. There is also a local version `MPI_Reduce` which computes the result only on one processor. These calls are examples of *collective operations* or *collectives*. The collectives are:

**reduction**: each processor has a data item, and these items need to be combined arithmetically with an addition, multiplication, max, or min operation. The result can be left on one processor, or on all, in which case we call this an **allreduce** operation.

**broadcast**: one processor has a data item that all processors need to receive.

**gather**: each processor has a data item, and these items need to be collected in an array, without combining them in an operations such as an addition. The result can be left on one processor, or on all, in which case we call this an **allgather**.

**scatter**: one processor has an array of data items, and each processor receives one element of that array.

**all-to-all**: each processor has an array of items, to be scattered to all other processors.

Collective operations are blocking (see section 2.7.3.4), although MPI 3.0 (which is currently only a draft) will have non-blocking collectives. We will analyze the cost of collective operations in detail in section 6.1.
### 2.7.3.6 Non-blocking communication

In a simple computer program, each instruction takes some time to execute, in a way that depends on what goes on in the processor. In parallel programs the situation is more complicated. A send operation, in its simplest form, declares that a certain buffer of data needs to be sent, and program execution will then stop until that buffer has been safely sent and received by another processor. This sort of operation is called a *non-local operation* since it depends on the actions of other processes, and a *blocking communication* operation since execution will halt until a certain event takes place.

Blocking operations have the disadvantage that they can lead to deadlock. In the context of message passing this describes the situation that a process is waiting for an event that never happens; for instance, it can be waiting to receive a message and the sender of that message is waiting for something else. You can easily recognize that deadlock occurs if two processes are waiting for each other, or more generally, if you have a cycle of processes where each is waiting for the next process in the cycle. Example:

```c
if ( /* this is process 0 */ )
    // wait for message from 1
else if ( /* this is process 1 */ )
    // wait for message from 0
```

A block receive here leads to deadlock. Even without deadlock, they can lead to considerable *idle time* in the processors, as they wait without performing any useful work. On the other hand, they have the advantage that it is clear when the buffer can be reused: after the operation completes, there is a guarantee that the data has been safely received at the other end.

The blocking behaviour can be avoided, at the cost of complicating the buffer semantics, by using *non-blocking communication* operations. A non-blocking send (`MPI_Isend`) declares that a data buffer needs to be sent, but then does not wait for the completion of the corresponding receive. There is a second operation `MPI_Wait` that will actually block until the receive has been completed. The advantage of this decoupling of sending and blocking is that it now becomes possible to write:

```c
MPI_ISend(somebuffer, &handle); // start sending, and
                                   // get a handle to this particular communication
{ ... } // do useful work on local data
MPI_Wait(handle); // block until the communication is completed;
{ ... } // do useful work on incoming data
```

With a little luck, the local operations take more time than the communication, and you have completely eliminated the communication time.

In addition to non-blocking sends, there are non-blocking receives. A typical piece of code then looks like

```c
MPI_ISend(sendbuffer, &sendhandle);
MPI_IReceive(recvbuffer, &recvhandle);
{ ... } // do useful work on local data
MPI_Wait(sendhandle); Wait(recvhandle);
{ ... } // do useful work on incoming data
```
Exercise 2.20. Take another look at equation (2.5) and give pseudocode that solves the problem using non-blocking sends and receives. What is the disadvantage of this code over a blocking solution?

### 2.7.3.7 MPI version 1 and 2

The first MPI standard [131] had a number of notable omissions, which are included in the MPI 2 standard [76]. One of these concerned parallel input/output: there was no facility for multiple processes to access the same file, even if the underlying hardware would allow that. A separate project MPI-I/O has now been rolled into the MPI-2 standard. We will discuss parallel I/O in this book.

A second facility missing in MPI, though it was present in PVM [42, 62] which predates MPI, is process management: there is no way to create new processes and have them be part of the parallel run.

Finally, MPI-2 has support for one-sided communication: one process put data into the memory of another, without the receiving process doing an actual receive instruction. We will have a short discussion in section 2.7.3.8 below.

### 2.7.3.8 One-sided communication

The MPI way of writing matching send and receive instructions is not ideal for a number of reasons. First of all, it requires the programmer to give the same data description twice, once in the send and once in the receive call. Secondly, it requires a rather precise orchestration of communication if deadlock is to be avoided; the alternative of using asynchronous calls is tedious to program, requiring the program to manage a lot of buffers. Lastly, it requires a receiving processor to know how many incoming messages to expect, which can be tricky in irregular applications. Life would be so much easier if it was possible to pull data from another processor, or conversely to put it on another processor, without that other processor being explicitly involved.

This style of programming is further encouraged by the existence of Remote Direct Memory Access (RDMA) support on some hardware. An early example was the Cray T3E. These days, one-sided communication is widely available through its incorporation in the MPI-2 library; section 2.7.3.7.

Let us take a brief look at one-sided communication in MPI-2, using averaging of array values as an example:

$$ \forall i: a_i \leftarrow (a_i + a_{i-1} + a_{i+1})/3. $$

The MPI parallel code will look like

```c
// do data transfer
a_local = (a_local+left+right)/3
```

It is clear what the transfer has to accomplish: the `a_local` variables needs to become the `left` variable on the processor with the next higher rank, and the `right` variables on the one with the next lower rank.

First of all, processors need to declare explicitly what memory area is available for one-sided transfer, the so-called 'window'. In this example, that consists of the `a_local`, `left`, and `right` variables on the processors:
MPI_Win_create(&a_local, ..., &data_windows);
MPI_Win_create(&left, ..., &left_windows);
MPI_Win_create(&right, ..., &right_windows);

The code now has two options: it is possible to push data out

target = my_tid-1;
MPI_Put(&a_local, ..., target, right_windows);
target = my_tid+1;
MPI_Put(&a_local, ..., target, left_windows);

or to pull it in

data_windows = a_local;
source = my_tid-1;
MPI_Get(&right, ..., data_windows);
source = my_tid+1;
MPI_Get(&left, ..., data_windows);

The above code will have the right semantics if the Put and Get calls are blocking; see section 2.7.3.4. However, part of the attraction of one-sided communication is that it makes it easier to express communication, and for this, a non-blocking semantics is assumed.

The problem with non-blocking one-sided calls is that it becomes necessary to ensure explicitly that communication is successfully completed. For instance, if one processor does a one-sided *put* operation on another, the other processor has no way of checking that the data has arrived, or indeed that transfer has begun at all. Therefore it is necessary to insert a global barrier in the program, for which every package has its own implementation. In MPI-2 the relevant call is the `MPI_Win_fence` routine. These barriers in effect divide the program execution in *supersteps*; see section 2.7.8.

Another form of one-sided communication is used in the Charm++ package; see section 2.7.7.

### 2.7.4 Hybrid shared/distributed memory computing

Modern architectures are often a mix of shared and distributed memory. For instance, a cluster will be distributed on the level of the nodes, but sockets and cores on a node will have shared memory. One level up, each socket can have a shared L3 cache but separate L2 and L1 caches. Intuitively it seems clear that a mix of shared and distributed programming techniques would give code that is optimally matched to the architecture. In this section we will discuss such hybrid programming models, and discuss their efficacy.

A common setup of clusters uses distributed memory *nodes*, where each node contains several *sockets*, that share memory. This suggests using MPI to communicate between the nodes (*inter-node communication*) and OpenMP for parallelism on the node (*intra-node communication*).

In practice this is realized as follows:

* On each node a single MPI process is started (rather than one per socket or core);
* This one MPI process than uses OpenMP (or another threading protocol) to spawn as many threads are there are independent sockets or cores on the node.

* The OpenMP threads can then easily access the same shared memory.

The alternative would be to have an MPI process on each core or socket, and do all communication through message passing, even between processes that can see the same shared memory.

This hybrid strategy may sound like a good idea but the truth is complicated.

* Message passing between MPI processes sounds like it's more expensive than communicating through shared memory. However, optimized versions of MPI can typically detect when processes are on the same node, and they will replace the message passing by a simple data copy. The only argument against using MPI is then that each process has its own data space, so there is memory overhead because each process has to allocate space for buffers and duplicates of the data that is copied.

* Threading is more flexible: if a certain part of the code needs more memory per process, an OpenMP approach could limit the number of threads on that part. On the other hand, flexible handling of threads incurs a certain amount of OS overhead that MPI does not have with its fixed processes.

* Shared memory programming is conceptually simple, but there can be unexpected performance pitfalls. For instance, the performance of two processes now can impeded by the need for maintaining *cache coherence* and by *false sharing*.

On the other hand, the hybrid approach offers some advantage since it bundles messages. For instance, if two MPI processes on one node send messages to each of two processes on another node there would be four messages; in the hybrid model these would be bundled into one message.

**Exercise 2.21.** Analyze the discussion in the last item above. Assume that the bandwidth between the two nodes is only enough to sustain one message at a time. What is the cost savings of the hybrid model over the purely distributed model? Hint: consider bandwidth and latency separately.

This bundling of MPI processes may have an advantage for a deeper technical reason. In order to support a *handshake protocol*, each MPI process needs a small amount of buffer space for each other process. With a larger number of processes this can be a limitation, so bundling is attractive on high core count processors such as the *Intel Xeon Phi*.

The MPI library is explicit about what sort of threading it supports: you can query whether multi-threading is supported at all, whether all MPI calls have to originate from one thread or one thread at-a-time, or whether there is complete freedom in making MPI calls from threads.

### 2.7.5 Parallel languages

One approach to mitigating the difficulty of parallel programming is the design of languages that offer explicit support for parallelism. There are several approaches, and we will see some examples.

* Some languages reflect the fact that many operations in scientific computing are data parallel (section 2.6.1). Languages such as *High Performance Fortran (HPF)* (section 2.7.5.3) have an *array syntax*, where operations such addition of arrays can be expressed $A = B+C$. This syntax
simplifies programming, but more importantly, it specifies operations at an abstract level, so that a lower level can make specific decision about how to handle parallelism. However, the data parallelism expressed in HPF is only of the simplest sort, where the data are contained in regular arrays. Irregular data parallelism is harder; the Chapel language (section 2.7.5.5) makes an attempt at addressing this.

* Another concept in parallel languages, not necessarily orthogonal to the previous, is that of Partitioned Global Address Space (PGAS) model: there is only one address space (unlike in the MPI model), but this address space is partitioned, and each partition has affinity with a thread or process. Thus, this model encompasses both SMP and distributed shared memory. A typical PGAS language, *Unified Parallel C (UPC)*, allows you to write programs that for the most part looks like regular C code. However, by indicating how the major arrays are distributed over processors, the program can be executed in parallel.

### 2.7.5.1 Discussion

Parallel languages hold the promise of making parallel programming easier, since they make communication operations appear as simple copies or arithmetic operations. However, by doing so they invite the user to write code that may not be efficient, for instance by inducing many small messages.

<figure><img src="image_46.png" /><figcaption>Figure 2.15: Data shift that requires communication</figcaption></figure>

As an example, consider arrays a, b that have been horizontally partitioned over the processors, and that are shifted (see figure 2.15):

```
for (i=0; i<N; i++)
    for (j=0; j<N/np; j++)
        a[i][j+joffset] = b[i][j+1+joffset]
```
If this code is executed on a shared memory machine, it will be efficient, but a naive translation in the distributed case will have a single number being communicated in each iteration of the *i* loop. Clearly, these can be combined in a single buffer send/receive operation, but compilers are usually unable to make this transformation. As a result, the user is forced to, in effect, re-implement the blocking that needs to be done in an MPI implementation:

```c
for (i=0; i<N; i++)
    t[i] = b[i][N/np+joffset]
for (i=0; i<N; i++)
    for (j=0; j<N/np-1; j++) {
        a[i][j] = b[i][j+1]
        a[i][N/np] = t[i]
    }
```

On the other hand, certain machines support direct memory copies through global memory hardware. In that case, PGAS languages can be more efficient than explicit message passing, even with physically distributed memory.

### 2.7.5.2 Unified Parallel C

Unified Parallel C (UPC) [153] is an extension to the C language. Its main source of parallelism is *data parallelism*, where the compiler discovers independence of operations on arrays, and assigns them to separate processors. The language has an extended array declaration, which allows the user to specify whether the array is partitioned by blocks, or in a *round-robin* fashion.

The following program in UPC performs a vector-vector addition.

```c
//vect_add.c
#include <upc_relaxed.h>
#define N 100*THREADS
shared int v1[N], v2[N], v1plusv2[N];
void main() {
    int i;
    for (i=MYTHREAD; i<N; i+=THREADS)
        v1plusv2[i]=v1[i]+v2[i];
}
```

The same program with an explicitly parallel loop construct:

```c
//vect_add.c
#include <upc_relaxed.h>
#define N 100*THREADS
shared int v1[N], v2[N], v1plusv2[N];
void main()
{
    int i;
```
```c
upc_forall (i=0; i<N; i++; i)
    v1plusv2[i]=v1[i]+v2[i];
}
```

is comparable to UPC in spirit, but based on Java rather than on C.

### 2.7.5.3 High Performance Fortran

High Performance Fortran<fnref n="9" /> (HPF) is an extension of Fortran 90 with constructs that support parallel computing, published by the High Performance Fortran Forum (HPFF). The HPFF was convened and chaired by Ken Kennedy of Rice University. The first version of the HPF Report was published in 1993.

Building on the array syntax introduced in Fortran 90, HPF uses a data parallel model of computation to support spreading the work of a single array computation over multiple processors. This allows efficient implementation on both SIMD and MIMD style architectures. HPF features included:

* New Fortran statements, such as FORALL, and the ability to create PURE (side effect free) procedures;

* The use of *compiler directives* for recommended distributions of array data;

* Extrinsic procedure interface for interfacing to non-HPF parallel procedures such as those using message passing;

* Additional library routines, including environmental inquiry, parallel prefix/suffix (e.g., 'scan'), data scattering, and sorting operations.

Fortran 95 incorporated several HPF capabilities. While some vendors did incorporate HPF into their compilers in the 1990s, some aspects proved difficult to implement and of questionable use. Since then, most vendors and users have moved to OpenMP-based parallel processing. However, HPF continues to have influence. For example the proposed BIT data type for the upcoming Fortran-2008 standard contains a number of new intrinsic functions taken directly from HPF.

### 2.7.5.4 Co-array Fortran

Co-array Fortran (CAF) is an extension to the Fortran 95/2003 language. The main mechanism to support parallelism is an extension to the array declaration syntax, where an extra dimension indicates the parallel distribution. For instance, in

```fortran
Real, dimension(100), codimension[*] :: X
Real :: Y(100)[*]
Real :: Z(100,200)[10,0:9,*]
```

arrays X, Y have 100 elements on each processor. Array Z behaves as if the available processors are on a three-dimensional grid, with two sides specified and the third adjustable to accommodate the available processors.

Communication between processors is now done through copies along the (co-)dimensions that describe the processor grid. The Fortran 2008 standard includes co-arrays.

<footnote n="9">9. This section quoted from Wikipedia</footnote>
### 2.7.5.5 Chapel

Chapel [24] is a new parallel programming language<fnref n="10" /> being developed by Cray Inc. as part of the DARPA-led High Productivity Computing Systems program (HPCS). Chapel is designed to improve the productivity of high-end computer users while also serving as a portable parallel programming model that can be used on commodity clusters or desktop multicore systems. Chapel strives to vastly improve the programmability of large-scale parallel computers while matching or beating the performance and portability of current programming models like MPI.

Chapel supports a multithreaded execution model via high-level abstractions for data parallelism, task parallelism, concurrency, and nested parallelism. Chapel's locale type enables users to specify and reason about the placement of data and tasks on a target architecture in order to tune for locality. Chapel supports global-view data aggregates with user-defined implementations, permitting operations on distributed data structures to be expressed in a natural manner. In contrast to many previous higher-level parallel languages, Chapel is designed around a multiresolution philosophy, permitting users to initially write very abstract code and then incrementally add more detail until they are as close to the machine as their needs require. Chapel supports code reuse and rapid prototyping via object-oriented design, type inference, and features for generic programming.

Chapel was designed from first principles rather than by extending an existing language. It is an imperative block-structured language, designed to be easy to learn for users of C, C++, Fortran, Java, Perl, Matlab, and other popular languages. While Chapel builds on concepts and syntax from many previous languages, its parallel features are most directly influenced by ZPL, High-Performance Fortran (HPF), and the Cray MTA's extensions to C and Fortran.

Here is vector-vector addition in Chapel:

```pascal
const BlockDist = newBlock1D(bbox=[1..m], tasksPerLocale=...);
const ProblemSpace: domain(1, 64)) distributed BlockDist = [1..m];
var A, B, C: [ProblemSpace] real;
forall (a, b, c) in (A, B, C) do
  a = b + alpha * c;
```

### 2.7.5.6 Fortress

Fortress [58] is a programming language developed by Sun Microsystems. Fortress<fnref n="11" /> aims to make parallelism more tractable in several ways. First, parallelism is the default. This is intended to push tool design, library design, and programmer skills in the direction of parallelism. Second, the language is designed to be more friendly to parallelism. Side-effects are discouraged because side-effects require synchronization to avoid bugs. Fortress provides transactions, so that programmers are not faced with the task of determining lock orders, or tuning their locking code so that there is enough for correctness, but not so much that performance is impeded. The Fortress looping constructions, together with the library, turns "iteration" inside out; instead of the loop specifying how the data is accessed, the data structures specify how the loop is run,

<footnote n="10">10. This section quoted from the Chapel homepage.</footnote>

<footnote n="11">11. This section quoted from the Fortress homepage.</footnote>
and aggregate data structures are designed to break into large parts that can be effectively scheduled for parallel execution. Fortress also includes features from other languages intended to generally help productivity – test code and methods, tied to the code under test; contracts that can optionally be checked when the code is run; and properties, that might be too expensive to run, but can be fed to a theorem prover or model checker. In addition, Fortress includes safe-language features like checked array bounds, type checking, and garbage collection that have been proven-useful in Java. Fortress syntax is designed to resemble mathematical syntax as much as possible, so that anyone solving a problem with math in its specification can write a program that can be more easily related to its original specification.

### 2.7.5.7 X10

X10 is an experimental new language currently under development at IBM in collaboration with academic partners. The X10 effort is part of the IBM PERCS project (Productive Easy-to-use Reliable Computer Systems) in the DARPA program on High Productivity Computer Systems. The PERCS project is focused on a hardware-software co-design methodology to integrate advances in chip technology, architecture, operating systems, compilers, programming language and programming tools to deliver new adaptable, scalable systems that will provide an order-of-magnitude improvement in development productivity for parallel applications by 2010.

X10 aims to contribute to this productivity improvement by developing a new programming model, combined with a new set of tools integrated into Eclipse and new implementation techniques for delivering optimized scalable parallelism in a managed runtime environment. X10 is a type-safe, modern, parallel, distributed object-oriented language intended to be very easily accessible to Java(TM) programmers. It is targeted to future low-end and high-end systems with nodes that are built out of multi-core SMP chips with non-uniform memory hierarchies, and interconnected in scalable cluster configurations. A member of the Partitioned Global Address Space (PGAS) family of languages, X10 highlights the explicit reification of locality in the form of places; lightweight activities embodied in async, future, foreach, and ateach constructs; constructs for termination detection (finish) and phased computation (clocks); the use of lock-free synchronization (atomic blocks); and the manipulation of global arrays and data structures.

### 2.7.5.8 Linda

As should be clear by now, the treatment of data is by far the most important aspect of parallel programming, far more important than algorithmic considerations. The programming system *Linda* [63, 64], also called a *coordination language*, is designed to address the data handling explicitly. Linda is not a language as such, but can, and has been, incorporated into other languages.

The basic concept of Linda is *tuple space*: data is added to a pool of globally accessible information by adding a label to it. Processes then retrieve data by their label, and without needing to know which processes added them to the tuple space.

Linda is aimed primarily at a different computation model than is relevant for *High-Performance Computing (HPC)*: it addresses the needs of asynchronous communicating processes. However, is has been used for scientific computation [38]. For instance, in parallel simulations of the heat equation (section 4.3), processors can write their data into tuple space, and neighbouring processes can retrieve their *ghost region* without having to know its provenance. Thus, Linda becomes one way of implementing *one-sided communication*.
### 2.7.5.9 The Global Arrays library

The Global Arrays library (http://www.emsl.pnl.gov/docs/global/) is another example of one-sided communication, and in fact it predates MPI. This library has as its prime data structure cartesian product arrays<fnref n="12" />, distributed over a processor grid of the same or lower dimension. Through library calls, any processor can access any sub-brick out of the array in either a put or get operation. These operations are non-collective. As with any one-sided protocol, a barrier sync is necessary to ensure completion of the sends/receives.

### 2.7.6 OS-based approaches

It is possible to design an architecture with a shared address space, and let the data movement be handled by the operating system. The Kendall Square computer [98] had an architecture name ‘all-cache’, where no data was natively associated with any processor. Instead, all data was considered to be cached on a processor, and moved through the network on demand, much like data is moved from main memory to cache in a regular CPU. This idea is analogous to the NUMA support in current SGI architectures.

### 2.7.7 Active messages

The MPI paradigm (section 2.7.3.3) is traditionally based on two-sided operations: each data transfer requires an explicit send and receive operation. This approach works well with relatively simple codes, but for complicated problems it becomes hard to orchestrate all the data movement. One of the ways to simplify consists of using active messages. This is used in the package Charm++[94].

With active messages, one processor can send data to another, without that second processor doing an explicit receive operation. Instead, the recipient declares code that handles the incoming data, a ‘method’ in objective orientation parlance, and the sending processor calls this method with the data that it wants to send. Since the sending processor in effect activates code on the other processor, this is also known as remote method invocation. A big advantage of this method is that overlap of communication and computation becomes easier to realize.

As an example, consider the matrix-vector multiplication with a tridiagonal matrix

$$ \forall i: y_i \leftarrow 2x_i - x_{i+1} - x_{i-1}. $$

See section 4.2.2 for an explanation of the origin of this problem in PDEs. Assuming that each processor has exactly one index $i$, the MPI code could look like:

```c
if ( /* I am the first or last processor */ )
    n_neighbours = 1;
else
    n_neighbours = 2;
/* do the MPI_Isend operations on my local data */
```

<footnote n="12">12. This means that if the array is three-dimensional, it can be described by three integers $n_1, n_2, n_3$, and each point has a coordinate $(i_1, i_2, i_3)$ with $1 \le i_1 \le n_1$ et cetera.</footnote>
```c
sum = 2*local_x_data;
received = 0;
for (neighbour=0; neighbour<n_neighbours; neighbour++) {
    MPI_WaitAny( /* wait for any incoming data */ )
    sum = sum - /* the element just received */
    received++
    if (received==n_neighbours)
        local_y_data = sum
}
```

With active messages this looks like:

```c
void incorporate_neighbour_data(x) {
    sum = sum-x;
    if (received==n_neighbours)
        local_y_data = sum
}
sum = 2*local_xdata;
received = 0;
alls[myid+1].incorporate_neighbour_data(local_x_data);
allung[myid-1].incorporate_neighbour_data(local_x_data);
```

### 2.7.8 Bulk synchronous parallelism

The MPI library (section 2.7.3.3) can lead to very efficient code. The price for this is that the programmer needs to spell out the communication in great detail. On the other end of the spectrum, PGAS languages (section 2.7.5) ask very little of the programmer, but give not much performance in return. One attempt to find a middle ground is the *Bulk Synchronous Parallel (BSP)* model [154, 147]. Here the programmer needs to spell out the communications, but not their ordering.

The BSP model orders the program into a sequence of *supersteps*, each of which ends with a *barrier* synchronization. The communications that are started in one superstep are all asynchronous and rely on the barrier for their completion. This makes programming easier and removes the possibility of deadlock. Moreover, all communication are of the *one-sided communication* type.

**Exercise 2.22.** Consider the parallel summing example in section 2.1. Argue that a BSP implementation needs $\log_2 n$ supersteps.

Because of its synchronization of the processors through the barriers concluding the supersteps the BSP model can do a simple cost analysis of parallel algorithms.

Another aspect of the BSP model is its use of *overdecomposition* of the problem, where multiple processes are assigned to each processor, as well as *random placement* of data and tasks. This is motivated with a statistical argument that shows it can remedy *load imbalance*. If there are *p* processors and if in a superstep *p* remote accesses are made, with high likelihood some processor receives $\log p / \log \log p$ accesses, while
others receive none. Thus, we have a load imbalance that worsens with increasing processor count. On the other hand, if $p \log p$ accesses are made, for instance because there are $\log p$ processes on each processor, the maximum number of accesses is $3 \log p$ with high probability. This means the load balance is within a constant factor of perfect.

The BSP model is implemented in BSPlib [87]. Other system can be said to be BSP-like in that they use the concept of supersteps; for instance Google's Pregel [119].

### 2.7.9 Program design for parallelism

A long time ago it was thought that some magic combination of compiler and runtime system could transform an existing sequential program into a parallel one. That hope has long evaporated, so these days a parallel program will have been written from the ground up as parallel. Of course there are different types of parallelism, and they each have their own implications for precisely how you design your parallel program. In this section we will briefly look into some of the issues.

#### 2.7.9.1 Parallel data structures

One of the issues in parallel program design is the use of Array-Of-Structures (AOS) vs Structure-Of-Arrays (SOA). In normal program design you often define a structure

```c
struct { int number; double xcoord,ycoord; } _Node;
struct { double xtrans,ytrans} _Vector;
typedef struct _Node* Node;
typedef struct _Vector* Vector;
```

and if you need a number of them you create an array of such structures.

```c
Node *nodes = (node) malloc( n_nodes*sizeof(struct _Node) );
```

This is the AOS design.

Now suppose that you want to parallelize an operation

```c
void shift(node the_point, vector by) {
    the_point->xcoord += by->xtrans;
    the_point->ycoord += by->ytrans;
}
```

which is done in a loop

```c
for (i=0; i<n_nodes; i++) {
    shift nodes[i], shift_vector);
}
```
This code has the right structure for MPI programming (section 2.7.3.3), where every processor has its own local array of nodes. This loop is also easily parallelizable with OpenMP (section 2.7.2).

However, in the 1980s codes had to be substantially rewritten as it was realized that the AOS design was not good for vector computers. In that case you operands need to be contiguous, and so codes had to go to a SOA design:

```c
node_numbers = (int*) malloc( n_nodes*sizeof(int) );
node_xcoordinates = // et cetera
node_ycoordinates = // et cetera
```

and you would iterate

```c
for (i=0; i<n_nodes; i++) {
    node_xcoords[i] += shift_vector->xtrans;
    node_ycoords[i] += shift_vector->ytrans;
}
```

Oh, did I just say that the original SOA design was best for distributed memory programming? That meant that 10 years after the vector computer era everyone had to rewrite their codes again for clusters. And of course nowadays, with increasing SIMD width, we need to go part way back to the AOS design. (There is some experimental software support for this transformation in the Intel ispc project, http://ispc.github.io/, which translates SPMD code to SIMD.)

### 2.7.9.2 Latency hiding

Communication between processors is typically slow, slower than data transfer from memory on a single processor, and much slower than operating on data. For this reason, it is good to think about the relative volumes of network traffic versus ‘useful’ operations when designing a parallel program. There has to be enough work per processor to offset the communication.

Another way of coping with the relative slowness of communication is to arrange the program so that the communication actually happens while some computation is going on. This is referred to as *overlapping computation with communication or latency hiding*.

For example, consider the parallel execution of a matrix-vector product $y = Ax$ (there will be further discussion of this operation in section 6.2.1). Assume that the vectors are distributed, so each processor $p$ executes

$$ \forall_{i \in I_p} : y_i = \sum_j a_{ij} x_j. $$

Since $x$ is also distributed, we can write this as

$$ \forall_{i \in I_p} : y_i = \left( \sum_{j \text{ local}} + \sum_{j \text{ not local}} \right) a_{ij} x_j. $$

This scheme is illustrated in figure 2.16. We can now proceed as follows:
<figure><img src="image_47.png" /><figcaption>Figure 2.16: The parallel matrix-vector product with a blockrow distribution.</figcaption></figure>

* Start the transfer of non-local elements of $x$;

* Operate on the local elements of $x$ while data transfer is going on;

* Make sure that the transfers are finished;

* Operate on the non-local elements of $x$.

**Exercise 2.23.** How much can you gain from overlapping computation and communication?
Hint: consider the border cases where computation takes zero time and there is only communication, and the reverse. Now consider the general case.

Of course, this scenario presupposes that there is software and hardware support for this overlap. MPI allows for this (see section 2.7.3.6), through so-called *asynchronous communication* or *non-blocking communication* routines. This does not immediately imply that overlap will actually happen, since hardware support is an entirely separate question.

## 2.8 Topologies

If a number of processors are working together on a single task, most likely they need to communicate data. For this reason there needs to be a way for data to make it from any processor to any other. In this section we will discuss some of the possible schemes to connect the processors in a parallel machine. Such a scheme is called a (processor) *topology*.

In order to get an appreciation for the fact that there is a genuine problem here, consider two simple schemes that do not 'scale up':

* Ethernet is a connection scheme where all machines on a network are on a single cable<sup>13</sup>. If one machine puts a signal on the wire to send a message, and another also wants to send a

<footnote n="13">13. We are here describing the original design of Ethernet. With the use of switches, especially in an HPC context, this description does not really apply anymore.</footnote>
message, the latter will detect that the sole available communication channel is occupied, and it will wait some time before retrying its send operation. Receiving data on ethernet is simple: messages contain the address of the intended recipient, so a processor only has to check whether the signal on the wire is intended for it.

The problems with this scheme should be clear. The capacity of the communication channel is finite, so as more processors are connected to it, the capacity available to each will go down. Because of the scheme for resolving conflicts, the average delay before a message can be started will also increase<fnref n="14" />.

* In a fully connected configuration, each processor has one wire for the communications with each other processor. This scheme is perfect in the sense that messages can be sent in the minimum amount of time, and two messages will never interfere with each other. The amount of data that can be sent from one processor is no longer a decreasing function of the number of processors; it is in fact an increasing function, and if the network controller can handle it, a processor can even engage in multiple simultaneous communications.

The problem with this scheme is of course that the design of the network interface of a processor is no longer fixed: as more processors are added to the parallel machine, the network interface gets more connecting wires. The network controller similarly becomes more complicated, and the cost of the machine increases faster than linearly in the number of processors.

In this section we will see a number of schemes that can be increased to large numbers of processors.

### 2.8.1 Some graph theory

The network that connects the processors in a parallel computer can conveniently be described with some elementary graph theory concepts. We describe the parallel machine with a graph where each processor is a node, and two nodes are connected<fnref n="15" /> if there is a direct connection between them.

We can then analyze two important concepts of this graph.

First of all, the degree of a node in a graph is the number of other nodes it is connected to. With the nodes representing processors, and the edges the wires, it is clear that a high degree is not just desirable for efficiency of computing, but also costly from an engineering point of view. We assume that all processors have the same degree.

Secondly, a message traveling from one processor to another, through one or more intermediate nodes, will most likely incur some delay at each stage of the path between the nodes. For this reason, the diameter of the graph is important. The diameter is defined as the maximum shortest distance, counting numbers of links, between any two nodes:

$$d(G) = \max_{i,j} \text{shortest path between } i \text{ and } j.$$

If $d$ is the diameter, and if sending a message over one wire takes unit time, this means a message will always arrive in at most time $d$.

<footnote n="14">14. It was initially thought that ethernet would be inferior to other solutions such as IBM's 'token ring'. It takes fairly sophisticated statistical analysis to prove that it works a lot better than was naively expected.</footnote>

<footnote n="15">15. We assume that connections are symmetric, so that the network is an *undirected graph*.</footnote>
Exercise 2.24. Find a relation between the number of processors, their degree, and the diameter of the connectivity graph.

In addition to the question ‘how long will a message from processor A to processor B take’, we often worry about conflicts between two simultaneous messages: is there a possibility that two messages, under way at the same time, will need to use the same network link? In figure 2.17 we illustrate what happens if every processor $p_i$ with $i < n/2$ send a message to $p_{i+n/2}$: there will be $n/2$ messages trying to get through the wire between $p_{n/2-1}$ and $p_{n/2}$. This sort of conflict is called *congestion* or *contention*. Clearly, the more

<figure><img src="image_48.png" /><figcaption>Figure 2.17: Contention for a network link due to simultaneous messages</figcaption></figure>

links a parallel computer has, the smaller the chance of congestion.

A precise way to describe the likelihood of congestion, is to look at the *bisection width*. This is defined as the minimum number of links that have to be removed to partition the processor graph into two unconnected graphs. For instance, consider processors connected as a linear array, that is, processor $P_i$ is connected to $P_{i-1}$ and $P_{i+1}$. In this case the bisection width is 1.

The bisection width $w$ describes how many messages can, guaranteed, be under way simultaneously in a parallel computer. Proof: take $w$ sending and $w$ receiving processors. The $w$ paths thus defined are disjoint: if they were not, we could separate the processors into two groups by removing only $w-1$ links.

In practice, of course, more than $w$ messages can be under way simultaneously. For instance, in a linear array, which has $w = 1$, $P/2$ messages can be sent and received simultaneously if all communication is between neighbours, and if a processor can only send or receive, but not both, at any one time. If processors can both send and receive simultaneously, $P$ messages can be under way in the network.

Bisection width also describes redundancy in a network: if one or more connections are malfunctioning, can a message still find its way from sender to receiver?

While bisection width is a measure express as a number of wires, in practice we care about the capacity through those wires. The relevant concept here is *bisection bandwidth*: the bandwidth across the bisection width, which is the product of the bisection width, and the capacity (in bits per second) of the wires. Bisection bandwidth can be considered as a measure for the bandwidth that can be attained if an arbitrary half of the processors communicates with the other half. Bisection bandwidth is a more realistic measure than the *aggregate bandwidth* which is sometimes quoted and which is defined as the total data rate if every processor is sending: the number of processors times the bandwidth of a connection times the number of simultaneous sends a processor can perform. This can be quite a high number, and it is typically not representative of the communication rate that is achieved in actual applications.
### 2.8.2 Busses

The first interconnect design we consider is to have all processors on the same memory bus. This design connects all processors directly to the same memory pool, so it offers a UMA or SMP model.

The main disadvantage of using a bus is the limited scalability, since only one processor at a time can do a memory access. To over this, we need to assume that processors are slower than memory, or that the processors have cache or other local memory to operate out of. In the latter case, maintaining cache coherence is easy with a bus by letting processors listen to all the memory traffic on the bus – a process known as snooping.

### 2.8.3 Linear arrays and rings

A simple way to hook up multiple processors is to connect them in a linear array: every processor has a number $i$, and processor $P_i$ is connected to $P_{i-1}$ and $P_{i+1}$. The first and last processor are possible exceptions: if they are connected to each other, we call the architecture a *ring network*.

This solution requires each processor to have two network connections, so the design is fairly simple.

**Exercise 2.25.** What is the bisection width of a linear array? Of a ring?

**Exercise 2.26.** With the limited connections of a linear array, you may have to be clever about how to program parallel algorithms. For instance, consider a 'broadcast' operation: processor 0 has a data item that needs to be sent to every other processor.
We make the following simplifying assumptions:

* a processor can send any number of messages simultaneously,

* but a wire can carry only one message at a time; however,

* communication between any two processors takes unit time, regardless the number of processors in between them.

In a fully connected network or a star network you can simply write
for $i = 1...N-1$:

send the message to processor $i$

With the assumption that a processor can send multiple messages, this means that the operation is done in one step.

Now consider a linear array. Show that, even with this unlimited capacity for sending, the above algorithm runs into trouble because of congestion.

Find a better way to organize the send operations. Hint: pretend that your processors are connected as a binary tree. Assume that there are $N = 2^n - 1$ processors. Show that the broadcast can be done in $\log N$ stages, and that processors only need to be able to send a single message simultaneously.

This exercise is an example of embedding a 'logical' communication pattern in a physical one.

### 2.8.4 2D and 3D arrays

A popular design for parallel computers is to organize the processors in a two-dimensional or three-dimensional cartesian mesh. This means that every processor has a coordinate $(i, j)$ or $(i, j, k)$, and it is
connected to its neighbours in all coordinate directions. The processor design is still fairly simple: the num-
ber of network connections (the degree of the connectivity graph) is twice the number of space dimensions
(2 or 3) of the network.

It is a fairly natural idea to have 2D or 3D networks, since the world around us is three-dimensional, and
computers are often used to model real-life phenomena. If we accept for now that the physical model
requires *nearest neighbour* type communications (which we will see is the case in section 4.2.3), then a
mesh computer is a natural candidate for running physics simulations.

Exercise 2.27. What is the diameter of a 3D cube of $n \times n \times n$ processors? What is the bisection width? How does that change if you add wraparound torus connections?

Exercise 2.28. Your parallel computer has its processors organized in a 2D grid. The chip manufacturer comes out with a new chip with same clock speed that is dual core instead of single core, and that will fit in the existing sockets. Critique the following argument: "the amount work per second that can be done (that does not involve communication) doubles; since the network stays the same, the bisection bandwidth also stays the same, so I can reasonably expect my new machine to become twice as fast<fnref n="16" />."

Grid-based designs often have so-called *wrap-around* or *torus* connections, which connect the left and right sides of a 2D grid, as well as the top and bottom. This is illustrated in figure 2.18.

<figure><img src="image_49.png" /><figcaption>Figure 2.18: A 2D grid with torus connections</figcaption></figure>

Some computer designs claim to be a grid of high dimensionality, for instance 5D, but not all dimensional are equal here. For instance a 3D grid where each *node* is a quad-socket quad-core can be considered as a 5D grid. However, the last two dimensions are fully connected.

**2.8.5 Hypercubes**

Above we gave a hand-waving argument for the suitability of mesh-organized processors, based on the prevalence of nearest neighbour communications. However, sometimes sends and receives between arbitrary processors occur. One example of this is the above-mentioned broadcast. For this reason, it is desirable

<footnote n="16">16. With the numbers one and two replaced by higher numbers, this is actually not a bad description of current trends in processor design.</footnote>
<figure><img src="image_50.png" /><figcaption>Figure 2.20: Hypercubes</figcaption></figure>

to have a network with a smaller diameter than a mesh. On the other hand we want to avoid the complicated design of a fully connected network.

A good intermediate solution is the *hypercube* design. An *n*-dimensional hypercube computer has $2^n$ processors, with each processor connected to one other in each dimension; see figure 2.20.

An easy way to describe this is to give each processor an address consisting of *d* bits: we give each node of a hypercube a number that is the bit pattern describing its location in the cube; see figure 2.19.

With this numbering scheme, a processor is then connected to all others that have an address that differs by exactly one bit. This means that, unlike in a grid, a processor's neighbours do not have numbers that differ by 1 or $\sqrt{p}$, but by 1, 2, 4, 8, ...

The big advantages of a hypercube design are the small diameter and large capacity for traffic through the network.

**Exercise 2.29.** What is the diameter of a hypercube? What is the bisection width?

One disadvantage is the fact that the processor design is dependent on the total machine size. In practice, processors will be designed with a maximum number of possible connections, and someone buying a smaller machine then will be paying for unused capacity. Another disadvantage is the fact that extending a given machine can only be done by doubling it: other sizes than $2^p$ are not possible.

**Exercise 2.30.** Consider the parallel summing example of section 2.1, and give the execution

<figure><img src="image_51.png" /><figcaption>Figure 2.19: Numbering of the nodes of a hypercube</figcaption></figure>
time of a parallel implementation on a hypercube. Show that the theoretical speedup from the example is attained (up to a factor) for the implementation on a hypercube.

### 2.8.5.1 Embedding grids in a hypercube

Above we made the argument that mesh-connected processors are a logical choice for many applications that model physical phenomena. Hypercubes do not look like a mesh, but they have enough connections that they can simply pretend to be a mesh by ignoring certain connections.

Let's say that we want the structure of a 1D array: we want processors with a numbering so that processor $i$ can directly send data to $i-1$ and $i+1$. We can not use the obvious numbering of nodes as in figure 2.19. For instance, node 1 is directly connected to node 0, but has a distance of 2 to node 2. The right neighbour of node 3 in a ring, node 4, even has the maximum distance of 3 in this hypercube. Clearly we need to renumber the nodes in some way.

What we will show is that it's possible to walk through a hypercube, touching every corner exactly once, which is equivalent to embedding a 1D mesh in the hypercube.

<figure><img src="image_52.png" /><figcaption>Figure 2.21: Gray codes</figcaption></figure>

The basic concept here is a (binary reflected) *Gray code* [73]. This is a way of ordering the binary numbers $0 \dots 2^d - 1$ as $g_0, \dots, g_{2^d-1}$ so that $g_i$ and $g_{i+1}$ differ in only one bit. Clearly, the ordinary binary numbers do not satisfy this: the binary representations for 1 and 2 already differ in two bits. Why do Gray codes help us? Well, since $g_i$ and $g_{i+1}$ differ only in one bit, it means they are the numbers of nodes in the hypercube that are directly connected.

Figure 2.21 illustrates how to construct a Gray code. The procedure is recursive, and can be described informally as 'divide the cube into two subcubes, number the one subcube, cross over to the other subcube, and number its nodes in the reverse order of the first one'. The result for a 2D cube is in figure 2.22.

<figure><img src="image_53.png" /><figcaption>Figure 2.22: Gray code numbering of the nodes of a hypercube</figcaption></figure>
Since a Gray code offers us a way to embed a one-dimensional ‘mesh’ into a hypercube, we can now work our way up.

**Exercise 2.31.** Show how a square mesh of $2^{2d}$ nodes can be embedded in a hypercube by appending the bit patterns of the embeddings of two $2^d$ node cubes. How would you accommodate a mesh of $2^{d_1+d_2}$ nodes? A three-dimensional mesh of $2^{d_1+d_2+d_3}$ nodes?

### 2.8.6 Switched networks

Above, we briefly discussed fully connected processors. They are impractical if the connection is made by making a large number of wires between all the processors. There is another possibility, however, by connecting all the processors to a *switch* or switching network. Some popular network designs are the *crossbar*, the *butterfly exchange*, and the *fat tree*.

Switching networks are made out of switching elements, each of which have a small number (up to about a dozen) of inbound and outbound links. By hooking all processors up to some switching element, and having multiple stages of switching, it then becomes possible to connect any two processors by a path through the network.

#### 2.8.6.1 Cross bar

The simplest switching network is a cross bar, an arrangement of $n$ horizontal and vertical lines, with a switch element on each intersection that determines whether the lines are connected; see figure 2.23. If we designate the horizontal lines as inputs the vertical as outputs, this is clearly a way of having $n$ inputs be mapped to $n$ outputs. Every combination of inputs and outputs (sometimes called a ‘permutation’) is allowed.

One advantage of this type of network is that no connection blocks another. The main disadvantage is that the number of switching elements is $n^2$, a fast growing function of the number of processors $n$.

<figure><img src="image_54.png" /><figcaption>Figure 2.23: A simple cross bar connecting 6 inputs to 6 outputs</figcaption></figure>

#### 2.8.6.2 Butterfly exchange

Butterfly exchanges are typically built out of small switching elements, and they have multiple stages; as the number of processors grows, the number of stages grows with it. Figure 2.24 shows a one-stage butterfly used to implement both shared memory and distributed (or distributed shared) memory, and the extension of the latter to two stages.

As you can see in figure 2.25, butterfly exchanges allow several processors to access memory simultaneously. Also, their access times are identical, so exchange networks are a way of implementing a *UMA* architecture; see section 2.5.1. One computer that was based on a Butterfly exchange network was the *BBN Butterfly*<fnref n="17" />.

<footnote n="17">17. http://en.wikipedia.org/wiki/BBN_Butterfly</footnote>
<figure><img src="image_55.png" /><figcaption>Figure 2.24: A butterfly exchange network for two and four processors/memories</figcaption></figure>

<figure><img src="image_56.png" /><figcaption>Figure 2.25: Two independent routes through a butterfly exchange network</figcaption></figure>

**Exercise 2.32.** For both the simple cross bar and the butterfly exchange, the network needs to be expanded as the number of processors grows. Give the number of wires (of some unit length) and the number of switching elements that is needed in both cases to connect *n* processors and memories. What is the time that a data packet needs to go from memory to processor, expressed in the unit time that it takes to traverse a unit length of wire and the time to traverse a switching element?

Routing through a butterfly network is done based on considering the bits in the destination address. On the *i*-th level the *i*-th digit is considered; if this is 1, the left exit of the switch is taken, if 0, the right exit. This is illustrated in figure 2.26. If we attach the memories to the processors, as in figure 2.25, we need only two bits (to the last switch) but a further three bits to describe the reverse route.

### 2.8.6.3 Fat-trees

If we were to connect switching nodes like a tree, there would be a big problem with congestion close to the root since there are only two wires attached to the root note. Say we have a *k*-level tree, so there are $2^k$ leaf nodes. If all leaf nodes in the left subtree try to communicate with nodes in the right subtree, we have $2^{k-1}$ messages going through just one wire into the root, and similarly out through one wire. A fat-tree
<figure><img src="image_57.png" /><figcaption>Figure 2.26: Routing through a three-stage butterfly exchange network</figcaption></figure>

is a tree network where each level has the same total bandwidth, so that this congestion problem does not occur: the root will actually have $2^{k-1}$ incoming and outgoing wires attached [74]. Figure 2.27 shows this structure on the left; on the right is shown a cabinet of the Stampede cluster, with a *leaf switch* for top and bottom half of the cabinet.

The first successful computer architecture based on a fat-tree was the Connection Machines CM5.

In fat-trees, as in other switching networks, each message carries its own routing information. Since in a fat-tree the choices are limited to going up a level, or switching to the other subtree at the current level, a message needs to carry only as many bits routing information as there are levels, which is $\log_2 n$ for n processors.

**Exercise 2.33.** Show that the bisection width of a fat tree is $P/2$ where P is the number of processor leaf nodes. Hint: show that there is only one way of splitting a fat tree-connected set of processors into two connected subsets.

The theoretical exposition of fat-trees in [113] shows that fat-trees are optimal in some sense: it can deliver messages as fast (up to logarithmic factors) as any other network that takes the same amount of space to build. The underlying assumption of this statement is that switches closer to the root have to connect more wires, therefore take more components, and correspondingly are larger. This argument, while theoretically interesting, is of no practical significance, as the physical size of the network hardly plays a role in the biggest currently available computers that use fat-tree interconnect. For instance, in the Ranger supercomputer of The University of Texas at Austin, the fat-tree switch connects 60,000 cores, yet takes less than 10 percent of the floor space.
<figure><img src="image_58.png" /><figcaption>Figure 2.27: A fat tree with a three-level interconnect (left); the leaf switches in a cabinet of the Stampede cluster (right)</figcaption></figure>

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

<figure><img src="image_60.png" /><figcaption>Figure 2.28: A fat-tree built from simple switching elements</figcaption></figure>

A fat tree, as sketched above, would be costly to build, since for every next level a new, bigger, switch would have to be designed. In practice, therefore, a network with the characteristics of a fat-tree is constructed from simple switching elements; see figure 2.28. This network is equivalent in its bandwidth and routing possibilities to a fat-tree. Routing algorithms will be slightly more complicated: in a fat-tree, a data packet can go up in only one way, but here a packet has to know to which of the two higher switches to route.

This type of switching network is one case of a *Clos network* [29].

### 2.8.7 Cluster networks

The above discussion was somewhat abstract, but in real-life clusters you can actually see the network designs reflected. For instance, *fat tree cluster networks* will have a central cabinet corresponding to the top level in the tree. Figure 2.29 shows the switches of the TACC Ranger (no longer in service) and *Stampede*
<figure><img src="image_61.png" /></figure>

<figure><img src="image_62.png" /><figcaption>Figure 2.29: Networks switches for the TACC Ranger and Stampede clusters</figcaption></figure>

clusters. In the second picture it can be seen that there are actually multiple redundant fat-tree networks.

On the other hand, clusters such as the *IBM BlueGene*, which is based on a *torus-based cluster*, will look like a collection of identical cabinets, since each contains an identical part of the network; see figure 2.30.

### 2.8.8 Bandwidth and latency

The statement above that sending a message can be considered a unit time operation, is of course unrealistic. A large message will take longer to transmit than a short one. There are two concepts to arrive at a more realistic description of the transmission process; we have already seen this in section 1.3.2 in the context of transferring data between cache levels of a processor.

**latency** Setting up a communication between two processors takes an amount of time that is independent of the message size. The time that this takes is known as the *latency* of a message. There are various causes for this delay.

* The two processors engage in 'hand-shaking', to make sure that the recipient is ready, and that appropriate buffer space is available for receiving the message.

* The message needs to be encoded for transmission by the sender, and decoded by the receiver.

* The actual transmission may take time: parallel computers are often big enough that, even at lightspeed, the first byte of a message can take hundreds of cycles to traverse the distance between two processors.
<figure><img src="image_63.png" /><figcaption>Figure 2.30: A BlueGene computer</figcaption></figure>

**bandwidth** After a transmission between two processors has been initiated, the main number of interest is the number of bytes per second that can go through the channel. This is known as the *bandwidth*. The bandwidth can usually be determined by the *channel rate*, the rate at which a physical link can deliver bits, and the *channel width*, the number of physical wires in a link. The channel width is typically a multiple of 16, usually 64 or 128. This is also expressed by saying that a channel can send one or two 8-byte words simultaneously.

Bandwidth and latency are formalized in the expression

$$T(n) = \alpha + \beta n$$

for the transmission time of an *n*-byte message. Here, $\alpha$ is the latency and $\beta$ is the time per byte, that is, the inverse of bandwidth. Sometimes we consider data transfers that involve communication, for instance in the case of a *collective operation*; see section 6.1. We then extend the transmission time formula to

$$T(n) = \alpha + \beta n + \gamma n$$

where $\gamma$ is the time per operation, that is, the inverse of the *computation rate*.

It would also be possible to refine this formulas as

$$T(n, p) = \alpha + \beta n + \delta p$$

where *p* is the number of network 'hops' that is traversed. However, on most networks the value of $\delta$ is far lower than of $\alpha$, so we will ignore it here. Also, in fat-tree networks (section 2.8.6.3) the number of hops is of the order of $\log P$, where *P* is the total number of processors, so it can never be very large anyway.

### 2.8.9 Locality in parallel computing

In section 1.5.2 you found a discussion of locality concepts in single-processor computing. The concept of *locality in parallel computing* comprises all this and a few levels more.
**Between cores; private cache** Cores on modern processors have private coherent caches. This means that it looks like you don't have to worry about locality, since data is accessible no matter what cache it is in. However, maintaining coherence costs bandwidth, so it is best to keep access localized.

**Between cores; shared cache** The cache that is shared between cores is one place where you don't have to worry about locality: this is memory that is truly symmetric between the processing cores.

**Between sockets** The sockets on a node (or motherboard) appear to the programmer to have shared memory, but this is really NUMA access (section 2.5.2) since the memory is associated with specific sockets.

**Through the network structure** Some networks have clear locality effects. You saw a simple example in section 2.8.1, and in general it is clear that any grid-type network will favour communication between 'nearby' processors. Networks based on fat-trees seem free of such contention issues, but the levels induce a different form of locality. One level higher than the locality on a node, small groups of nodes are often connected by a *leaf switch*, which prevents data from going to the central switch.

## 2.9 Multi-threaded architectures

The architecture of modern CPUs is largely dictated by the fact that getting data from memory is much slower than processing it. Hence, a hierarchy of ever faster and smaller tries to keep data as close to the processing unit as possible, mitigating the long latency and small bandwidth of main memory. The ILP in the processing unit also helps to hide the latency and more fully utilize the available bandwidth.

However, finding ILP is a job for the compiler and there is a limit to what it can practically find. On the other hand, scientific codes are often very *data parallel* in a sense that is obvious to the programmer, though not necessarily to the compiler. Would it be possible for the programmer to specify this parallelism explicitly and for the processor to use it?

In section 2.4.1 you saw that SIMD architectures can be programmed in an explicitly data parallel way. What if we have a great deal of data parallelism but not that many processing units? In that case, we could turn the parallel instruction streams into threads (see section 2.7.1) and have multiple threads be executed on each processing unit. Whenever a thread would stall because of an outstanding memory request, the processor could switch to another thread for which all the necessary inputs are available. This is called *multi-threading*. While it sounds like a way of preventing the processor from waiting for memory, it can also be viewed as a way of keeping memory maximally occupied.

**Exercise 2.34.** If you consider the long latency and limited bandwidth of memory as two separate problems, does multi-threading address them both?

The problem here is that most CPUs are not good at switching quickly between threads. A *context switch* (switching between one thread and another) takes a large number of cycles, comparable to a wait for data from main memory. In a so-called *Multi-Threaded Architecture (MTA)* a context-switch is very efficient,
sometimes as little as a single cycle, which makes it possible for one processor to work on many threads simultaneously.

The multi-threaded concept was explored in the Tera Computer MTA machine, which evolved into the current Cray XMT<sup>18</sup>.

The other example of an MTA is the GPU, where the processors work as SIMD units, while being themselves multi-threaded; see section 2.10.3.

## 2.10 Co-processors

Current CPUs are built to be moderately efficient at just about any conceivable computation. This implies that by restricting the functionality of a processor it may be possible to raise its efficiency, or lower its power consumption at similar efficiency. Thus, the idea of incorporating a *co-processor* attached to the *host* process has been explored many times. For instance, Intel's 8086 chip, which powered the first generation of IBM PCs, could have a numerical co-processor, the 80287, added to it. This processor was very efficient at transcendental functions and it also incorporated SIMD technology. Using separate functionality for graphics has also been popular, leading to the SSE instructions for the x86 processor, and separate GPU units to be attached to the PCI-X bus.

Further examples are the use of co-processors for Digital Signal Processing (DSP) instructions, as well as FPGA boards which can be reconfigured to accomodate specific needs. Early *array processors* such as the ICL DAP were also co-processors.

In this section we look briefly at some modern incarnations of this idea, in particular GPUs.

### 2.10.1 A little history

Co-processors can be programmed in two different ways: sometimes it is seamlessly integrated, and certain instructions are automatically executed there, rather than on the 'host' processor. On the other hand, it is also possible that co-processor functions need to be explicitly invoked, and it may even be possible to overlap co-processor functions with host functions. The latter case may sound attractive from an efficiency point of view, but it raises a serious problem of programmability. The programmer now needs to identify explicitly two streams of work: one for the host processor and one for the co-processor.

Some notable parallel machines with co-processors are:

* The *Intel Paragon* (1993) had two processors per node, one for communication and the other for computation. These were in fact identical, the *Intel i860* Intel i860 processor. In a later revision, it became possible to pass data and function pointers to the communication processors.

* The *IBM Roadrunner* at Los Alamos was the first machine to reach a PetaFlop<sup>19</sup>. It achieved this speed through the use of Cell co-processors. Incidentally, the Cell processor is in essence the engine of the Sony Playstation3, showing again the commoditization of supercomputers (section 2.4.3).

<footnote n="18">18. Tera Computer renamed itself Cray Inc. after acquiring Cray Research from SGI.</footnote>

<footnote n="19">19. The *Grape computer* had reached this point earlier, but that was a special purpose machine for molecular dynamics calculations.</footnote>
* The Chinese *Tianhe-1A* topped the Top 500 list in 2010, reaching about 2.5 PetaFlops through the use of NVidia GPUs.

* The *Tianhe-2* and the *TACC Stampede cluster* use *Intel Xeon Phi* co-processors.

The Roadrunner and Tianhe-1A are examples of co-processors that are very powerful, and that need to be explicitly programmed independently of the host CPU. For instance, code running on the GPUs of the Tianhe-1A is programmed in *CUDA* and compiled separately.

In both cases the programmability problem is further exacerbated by the fact that the co-processor can not directly talk to the network. To send data from one co-processor to another it has to be passed to a host processor, from there through the network to the other host processor, and only then moved to the target co-processor.

### 2.10.2 Bottlenecks

Co-processors often have their own memory, and the *Intel Xeon Phi* can run programs independently, but more often there is the question of how to access the memory of the host processor. A popular solution is to connect the co-processor through a *PCI bus*. Accessing host memory this way is slower than the direct connection that the host processor has. For instance, the *Intel Xeon Phi* has a *bandwidth* of 512-bit wide at 5.5GT per second (we will get to that 'GT' in a second), while its connection to host memory is 5.0GT/s, but only 16-bit wide.

**GT measure** We are used to seeing bandwidth measured in gigabits per second. For a *PCI bus* one often see the *GT* measure. This stands for giga-transfer, and it measures how fast the bus can change state between zero and one. Normally, every state transition would correspond to a bit, except that the bus has to provide its own clock information, and if you would send a stream of identical bits the clock would get confused. Therefore, every 8 bits are encoded in 10 bits, to prevent such streams. However, this means that the effective bandwidth is lower than the theoretical number, by a factor of 4/5 in this case.

And since manufacturers like to give a positive spin on things, they report the higher number.

### 2.10.3 GPU computing

A *Graphics Processing Unit (GPU)* (or sometimes *General Purpose Graphics Processing Unit (GPGPU)*) is a special purpose processor, designed for fast graphics processing. However, since the operations done for graphics are a form of arithmetic, GPUs have gradually evolved a design that is also useful for non-graphics computing. The general design of a GPU is motivated by the 'graphics pipeline': identical operations are performed on many data elements in a form of *data parallelism* (section 2.6.1), and a number of such blocks of data parallelism can be active at the same time.

The basic limitations that hold for a CPU hold for a GPU: accesses to memory incur a long latency. The solution to this problem in a CPU is to introduce levels of cache; in the case of a GPU a different approach is taken (see also section 2.9). GPUs are concerned with *throughput computing*, delivering large amounts of data with high average rates, rather than any single result as quickly as possible. This is made possible by supporting many threads (section 2.7.1) and having very fast switching between them. While one thread is waiting for data from memory, another thread that already has its data can proceed with its computations.
### 2.10.3.1 SIMD-type programming with kernels

Present day GPUs<fnref n="20" /> have an architecture that combines SIMD and SPMD parallelism. Threads are not completely independent, but are ordered in *thread blocks*, where all threads in the block execute the same instruction, making the execution SIMD. It is also possible to schedule the same instruction stream (a 'kernel' in Cuda terminology) on more than one thread block. In this case, thread blocks can be out of sync, much like processes in an SPMD context. However, since we are dealing with threads here, rather than processes, the term *Single Instruction Multiple Thread (SIMT)* is used.

This software design is apparent in the hardware; for instance, an NVidia GPU has 16–30 Streaming Multiprocessors (SMs), and a SMs consists of 8 Streaming Processors (SPs), which correspond to processor cores; see figure 2.31. The SPs act in true SIMD fashion. The number of cores in a GPU is typically larger than in traditional multi-core processors, but the cores are more limited. Hence, the term *manycore* is used here.

<figure><img src="image_64.png" /><figcaption>Figure 2.31: Diagram of a GPU</figcaption></figure>

The SIMD, or *data parallel*, nature of GPUs becomes apparent in the way *CUDA* starts processes. A *kernel*, that is, a function that will be executed on the GPU, is started on *mn* cores by:

`KernelProc<<< m, n >>> (args)`

The collection of *mn* cores executing the kernel is known as a *grid*, and it is structured as *m* *thread blocks* of *n* threads each. A thread block can have up to 512 threads.

Recall that threads share an address space (see section 2.7.1), so they need a way to identify what part of the data each thread will operate on. For this, the blocks in a thread are numbered with *x*, *y* coordinates, and the threads in a block are numbered with *x*, *y*, *z* coordinates. Each thread knows its coordinates in the block, and its block's coordinates in the grid.

<footnote n="20">20. The most popular GPUs today are made by NVidia, and are programmed in *CUDA*, an extension of the C language.</footnote>
We illustrate this with a vector addition example:

```c
// Each thread performs one addition
__global__ void vecAdd(float* A, float* B, float* C)
{
    int i = threadIdx.x + blockDim.x * blockIdx.x;
    C[i] = A[i] + B[i];
}

int main()
{
    // Run grid of N/256 blocks of 256 threads each
    vecAdd<<< N/256, 256 >>>(d_A, d_B, d_C);
}
```

This shows the SIMD nature of GPUs: every thread executes the same scalar program, just on different data.

Threads in a thread block are truly data parallel: if there is a conditional that makes some threads take the *true* branch and others the *false* branch, then one branch will be executed first, with all threads in the other branch stopped. Subsequently, *and not simultaneously*, the threads on the other branch will then execute their code. This may induce a severe performance penalty.

GPUs rely on large amount of data parallelism and the ability to do a fast *context switch*. This means that they will thrive in graphics and scientific applications, where there is lots of data parallelism. However they are unlikely to do well on 'business applications' and operating systems, where the parallelism is of the Instruction Level Parallelism (ILP) type, which is usually limited.

### 2.10.3.2 GPUs versus CPUs

These are some of the differences between GPUs and regular CPUs:

* First of all, as of this writing (late 2010), GPUs are attached processors, for instance over a PCI-X bus, so any data they operate on has to be transferred from the CPU. Since the memory bandwidth of this transfer is low, at least 10 times lower than the memory bandwidth in the GPU, sufficient work has to be done on the GPU to overcome this overhead.

* Since GPUs are graphics processors, they put an emphasis on *single precision* floating point arithmetic. To accomodate the scientific computing community, *double precision* support is increasing, but double precision speed is typically half the single precision flop rate. This discrepancy is likely to be addressed in future generations.

* A CPU is optimized to handle a single stream of instructions that can be very heterogeneous in character; a GPU is made explicitly for data parallelism, and will perform badly on traditional codes.

* A CPU is made to handle one *thread*, or at best a small number of threads. A GPU *needs* a large number of threads, far larger than the number of computational cores, to perform efficiently.
### 2.10.3.3 Expected benefit from GPUs

GPUs have rapidly gained a reputation for achieving high performance, highly cost effectively. Stories abound of codes that were ported with minimal effort to CUDA, with a resulting speedup of sometimes 400 times. Is the GPU really such a miracle machine? Were the original codes badly programmed? Why don't we use GPUs for everything if they are so great?

The truth has several aspects.

First of all, a GPU is not as general-purpose as a regular CPU: GPUs are very good at doing *data parallel* computing, and CUDA is good at expressing this fine-grained parallelism elegantly. In other words, GPUs are suitable for a certain type of computation, and will be a poor fit for many others.

Conversely, a regular CPU is not necessarily good at *data parallelism*. Unless the code is very carefully written, performance can degrade from optimal by approximately the following factors:

* Unless directives or explicit parallel constructs are used, compiled code will only use 1 out of the available cores, say 4.

* If instructions are not pipelined, the latency because of the floating point pipeline adds another factor of 4.

* If the core has independent add and multiply pipelines, another factor of 2 will be lost if they are not both used simultaneously.

* Failure to use SIMD registers can add more to the slowdown with respect to peak performance.

Writing the optimal CPU implementation of a computational kernel often requires writing in assembler, whereas straightforward CUDA code will achieve high performance with comparatively little effort, provided of course the computation has enough data parallelism.

### 2.10.4 Intel Xeon Phi

Recently, Intel has released a co-processor, the *Intel Xeon Phi* (also known by its architecture design as Many Integrated Cores (MIC)) which is expressly designed for numerical computing. It has both differences and similarities with GPUs.

* Both are connected through a *PCI-X* bus, which means that operations on the device have a considerable latency in starting up.

* The Xeon Phi has general purpose cores, so that it can run whole programs; GPUs has that only to a limited extent (see section 2.10.3.1).

* The Xeon Phi accepts ordinary C code.

* Both architectures require a large amount of SIMD-style parallelism, in the case of the Xeon Phi because of the 8-wide AVX instructions.

* Both devices work, or can work, through *offloading* from a host program. In the case of the Xeon Phi this can happen with OpenMP constructs and *MKL* calls.
### 2.11 Remaining topics

#### 2.11.1 Load balancing

In much of this chapter, we assumed that a problem could be perfectly divided over processors, that is, a processor would always be performing useful work, and only be *idle* because of latency in communication. In practice, however, a processor may be idle because it is waiting for a message, and the sending processor has not even reached the send instruction in its code. Such a situation, where one processor is working and another is idle, is described as *load unbalance*: there is no intrinsic reason for the one processor to be idle, and it could have been working if we had distributed the work load differently.

There is an asymmetry between processors having too much work and having not enough work: it is better to have one processor that finishes a task early, than having one that is overloaded so that all others wait for it.

**Exercise 2.35.** Make this notion precise. Suppose a parallel task takes time 1 on all processors but one.

* Let $0 < \alpha < 1$ and let one processor take time $1 + \alpha$. What is the speedup and efficiency as function of the number of processors? Consider this both in the Amdahl and Gustafsson sense (section 2.3.3).

* Answer the same questions if one processor takes time $1 - \alpha$.

Load balancing is often expensive since it requires moving relatively large amounts of data. For instance, section 6.5 has an analysis showing that the data exchanges during a sparse matrix-vector product is of a lower order than what is stored on the processor. However, we will not go into the actual cost of moving: our main concerns here are to balance the workload, and to preserve any locality in the original load distribution.

##### 2.11.1.1 Load balancing versus data distribution

There is a duality between work and data: in many applications the distribution of data implies a distribution of work and the other way around. If an application updates a large array, each element of the array typically 'lives' on a uniquely determined processor, and that processor takes care of all the updates to that element. This strategy is known as *owner computes*.

Thus, there is a direct relation between data and work, and, correspondingly, data distribution and load balancing go hand in hand. For instance, in section 6.2 we will talk about how data distribution influences the efficiency, but this immediately translates to concerns about load distribution:

* Load needs to be evenly distributed. This can often be done by evenly distributing the data, but sometimes this relation is not linear.

* Tasks need to be placed to minimize the amount of traffic between them. In the matrix-vector multiplication case this means that a two-dimensional distribution is to be preferred over a one-dimensional one; the discussion about *space-filling curves* is similarly motivated.

As a simple example of how the data distribution influences the load balance, consider a linear array where each point undergoes the same computation, and each computation takes the same amount of time. If the length of the array, $N$, is perfectly divisible by the number of processors, $p$, the work is perfectly evenly distributed. If the data is not evenly divisible, we start by assigning $\lfloor N/p \rfloor$ points to each processor, and the remaining $N - p\lfloor N/p \rfloor$ points to the last processors.
**Exercise 2.36.** In the worst case, how unbalanced does that make the processors’ work? Compare this scheme to the option of assigning $\lceil N/p \rceil$ points to all processors except one, which gets fewer; see the exercise above.

It is better to spread the surplus $r = N - p\lceil N/p \rceil$ over $r$ processors than one. This could be done by giving one extra data point to the first or last $r$ processors. While this scheme is decently balanced, computing for instance to what processor a given point belongs is tricky. The following scheme makes such computations easier: let $f(i) = \lceil iN/p \rceil$, then processor $i$ gets points $f(i)$ up to $f(i+1)$.

**Exercise 2.37.** Show that $\lceil N/p \rceil \le f(i+1) - f(i) \le \lceil N/p \rceil$.

Under this scheme, the processor that owns index $i$ is $\lfloor (p(i+1) - 1)/N \rfloor$.

### 2.11.1.2 Load balancing of independent tasks

In other cases, work load is not directly determined by data. This can happen if there is a pool of work to be done, and the processing time for each work item is not or easily computed from its description. In such cases we may want some flexibility in assigning work to processes.

Let us first consider the case of a job that can be partitioned into independent tasks that do not communicate. An example would be computing the pixels of a *Mandelbrot set*, where each pixel is set according to a mathematical function that does not depend on surrounding pixels. If we could predict the time it would take to draw an arbitrary part of the picture, we could make a perfect division of the work and assign it to the processors. This is known as *static load balancing*.

More realistically, we can not predict the running time of a part of the job perfectly, and we use an *overde-composition* of the work: we divide the work in more tasks than there are processors. These tasks are then assigned to a *work pool*, and processors take the next job from the pool whenever they finish a job. This is known as *dynamic load balancing*. Many graph and combinatorial problems can be approached this way; see section 2.6.3. For task assignment in a multicore context, see section 6.12.

There are results that show that randomized assignment of tasks to processors is statistically close to optimal [96], but this ignores the aspect that in scientific computing tasks typically communicate frequently.

**Exercise 2.38.** Suppose you have tasks $\{T_i\}_{i=1,...,N}$ with running times $t_i$, and an unlimited number of processors. Look up *Brent's theorem* in section 2.3.1.2, and derive from it that the fastest execution scheme for the tasks can be characterized as follows: there is one processor that only executes the task with maximal $t_i$ value.

### 2.11.1.3 Load balancing as graph problem

Next let us consider a parallel job where the parts do communicate. In this case we need to balance both the scalar workload and the communication.

A parallel computation can be formulated as a graph (see Appendix 16 for an introduction to graph theory) where the processors are the vertices, and there is an edge between two vertices if their processors need to communicate at some point. Such a graph is often derived from an underlying graph of the problem being solved. As an example consider the matrix-vector product $y = Ax$ where $A$ is a sparse matrix, and look
in detail at the processor that is computing $y_i$ for some $i$. The statement $y_i \leftarrow y_i + A_{ij}x_j$ implies that this processor will need the value $x_j$, so, if this variable is on a different processor, it needs to be sent over.

We can formalize this: Let the vectors $x$ and $y$ be distributed disjointly over the processors, and define uniquely $P(i)$ as the processor that owns index $i$. Then there is an edge $(P, Q)$ if there is a nonzero element $a_{ij}$ with $P = P(i)$ and $Q = P(j)$. This graph is undirected in the case of a *structurally symmetric matrix*, that is $a_{ij} \neq 0 \Leftrightarrow a_{ji} \neq 0$.

The distribution of indices over the processors now gives us vertex and edge weights: a processor has a vertex weight that is the number of indices owned by it; an edge $(P, Q)$ has a weight that is the number of vector components that need to be sent from $Q$ to $P$, as described above.

The load balancing problem can now be formulated as follows:

Find a partitioning $\mathbb{P} = \cup_i \mathbb{P}_i$, such the variation in vertex weights is minimal, and simultaneously the edge weights are as low as possible.

Minimizing the variety in vertex weights implies that all processor have approximately the same amount of work. Keeping the edge weights low means that the amount of communication is low. These two objectives need not be satisfiable at the same time: some trade-off is likely.

**Exercise 2.39.** Consider the limit case where processors are infinitely fast and bandwidth between processors is also unlimited. What is the sole remaining factor determining the runtime? What graph problem do you need to solve now to find the optimal load balance? What property of a sparse matrix gives the worst case behaviour?

An interesting approach to load balancing comes from spectral graph theory (section 16.5): if $A_G$ is the adjacency matrix of an undirected graph and $D_G - A_G$ the *graph Laplacian*, then the eigenvector $u_1$ to the smallest eigenvalue zero is positive, and the eigenvector $u_2$ to the next eigenvalue is orthogonal to it. Therefore $u_2$ has to have elements of alternating sign; further analysis shows that the elements with positive sign are connected, as are the negative ones. This leads to a natural bisection of the graph.

### 2.11.1.4 Load redistributing

In certain applications an initial load distribution is clear, but later adjustments are needed. A typical example is in Finite Element Method (FEM) codes, where load can be distributed by a partitioning of the physical domain; see section 6.5.3. If later the discretization of the domain changes, the load has to be *rebalanced* or *redistributed*. In the next subsection we will see one technique for load balancing and rebalancing that is aimed at preserving locality.

### 2.11.1.5 Load balancing with space-filling curves

In the previous sections we considered two aspects of load balancing: making sure all processors have an approximately equal amount of work, and letting the distribution reflect the structure of the problem so that communication is kept within reason. We can phrase the second point trying to preserve the locality of the problem when distributed over a parallel machine: points in space that are close together are likely to interact, so they should be on the same processor, or at least one not too far away.

Striving to preserve locality is not obviously the right strategy. In BSP (see section 2.7.8) a statistical argument is made that *random placement* will give a good load balance as well as balance of communication.
Exercise 2.40. Consider the assignment of processes to processors, where the structure of the problem is such that each processes only communicates with its nearest neighbours, and let processors be ordered in a two-dimensional grid. If we do the obvious assignment of the process grid to the processor grid, there will be no contention. Now write a program that assigns processes to random processors, and evaluate how much contention there will be.

In the previous section you saw how graph partitioning techniques can help with the second point of preserving problem locality. In this section you will see a different technique that is attractive both for the initial load assignment and for subsequent *load rebalancing*. In the latter case, a processor's work may increase or decrease, necessitating moving some of the load to a different processor.

For instance, some problems are adaptively refined<sup>21</sup>. This is illustrated in figure 2.32. If we keep track

<figure><img src="image_65.png" /><figcaption>Figure 2.32: Adaptive refinement of a domain in subsequent levels</figcaption></figure>

of these refinement levels, the problem gets a tree structure, where the leaves contain all the work. Load balancing becomes a matter of partitioning the leaves of the tree over the processors; figure 2.33. Now we observe that the problem has a certain locality: the subtrees of any non-leaf node are physically close, so there will probably be communication between them.

* Likely there will be more subdomains than processors; to minimize communication between processors, we want each processor to contain a simply connected group of subdomains. Moreover, we want each processor to cover a part of the domain that is 'compact' in the sense that it has low aspect ratio, and low surface-to-volume ratio.

* When a subdomain gets further subdivided, part of the load of its processor may need to be shifted to another processor. This process of *load redistributing* should preserve locality.

To fulfill these requirements we use Space-Filling Curves (SFCs). A Space-Filling Curve (SFC) for the load balanced tree is shown in figure 2.34. We will not give a formal discussion of SFCs; instead we will let

<footnote n="21">21. For a detailed discussion, see [21].</footnote>
<figure><img src="image_66.png" /><figcaption>Figure 2.33: Load distribution of an adaptively refined domain</figcaption></figure>

<figure><img src="image_67.png" /><figcaption>Figure 2.34: A space filling curve for the load balanced tree</figcaption></figure>

figure 2.35 stand for a definition: a SFC is a recursively defined curve that touches each subdomain once<fnref n="22" />. The SFC has the property that domain elements that are close together physically will be close together on the curve, so if we map the SFC to a linear ordering of processors we will preserve the locality of the problem.

More importantly, if the domain is refined by another level, we can refine the curve accordingly. Load can then be redistributed to neighbouring processors on the curve, and we will still have locality preserved.

<footnote n="22">22. Space-Filling Curves (SFCs) were introduced by Peano as a mathematical device for constructing a continuous surjective function from the line segment [0, 1] to a higher dimensional cube [0, 1]¹⁰. This upset the intuitive notion of dimension that 'you can not stretch and fold a line segment to fill up the square'. A proper treatment of the concept of dimension was later given by Brouwer.</footnote>
<figure><img src="image_68.png" /><figcaption>Figure 2.35: Space filling curves, regularly and irregularly refined</figcaption></figure>

## 2.11.2 Distributed computing, grid computing, cloud computing

In this section we will take a short look at terms such as *cloud computing*, and an earlier term *distributed computing*. These are concepts that have a relation to parallel computing in the scientific sense, but that differ in certain fundamental ways.

Distributed computing can be traced back as coming from large database servers, such as airline reservations systems, which had to be accessed by many travel agents simultaneously. For a large enough volume of database accesses a single server will not suffice, so the mechanism of *remote procedure call* was invented, where the central server would call code (the procedure in question) on a different (remote) machine. The remote call could involve transfer of data, the data could be already on the remote machine, or there would be some mechanism that data on the two machines would stay synchronized. This gave rise to the *Storage Area Network (SAN)*. A generation later than distributed database systems, web servers had to deal with the same problem of many simultaneous accesses to what had to act like a single server.

We already see one big difference between distributed computing and high performance parallel computing. Scientific computing needs parallelism because a single simulation becomes too big or slow for one machine; the business applications sketched above deal with many users executing small programs (that is, database or web queries) against a large data set. For scientific needs, the processors of a parallel machine (the nodes in a cluster) have to have a very fast connection to each other; for business needs no such network is needed, as long as the central dataset stays coherent.

Both in *HPC* and in business computing, the server has to stay available and operative, but in distributed computing there is considerably more liberty in how to realize this. For a user connecting to a service such as a database, it does not matter what actual server executes their request. Therefore, distributed computing can make use of *virtualization*: a virtual server can be spawned off on any piece of hardware.

An analogy can be made between remote servers, which supply computing power wherever it is needed, and the electric grid, which supplies electric power wherever it is needed. This has led to *grid computing* or *utility computing*, with the Teragrid, owned by the US National Science Foundation, as an example. Grid computing was originally intended as a way of hooking up computers connected by a *Local Area Network (LAN)* or *Wide Area Network (WAN)*, often the Internet. The machines could be parallel themselves, and were often owned by different institutions. More recently, it has been viewed as a way of sharing resources, both datasets, software resources, and scientific instruments, over the network.

The notion of utility computing as a way of making services available, which you recognize from the above description of distributed computing, went mainstream with Google's search engine, which indexes the
whole of the Internet. Another example is the GPS capability of Android mobile phones, which combines GIS, GPS, and mashup data. The computing model by which Google's gathers and processes data has been formalized in MapReduce [33]. It combines a data parallel aspect (the 'map' part) and a central accumulation part ('reduce'). Neither involves the tightly coupled neighbour-to-neighbour communication that is common in scientific computing. An open source framework for MapReduce computing exists in Hadoop [78]. Amazon offers a commercial Hadoop service.

The concept of having a remote computer serve user needs is attractive even if no large datasets are involved, since it absolves the user from the need of maintaining software on their local machine. Thus, Google Docs offers various 'office' applications without the user actually installing any software. This idea is sometimes called *Software As a Service (SAS)*, where the user connects to an 'application server', and accesses it through a client such as a web browser. In the case of Google Docs, there is no longer a large central dataset, but each user interacts with their own data, maintained on Google's servers. This of course has the large advantage that the data is available from anywhere the user has access to a web browser.

The SAS concept has several connections to earlier technologies. For instance, after the mainframe and workstation eras, the so-called *thin client* idea was briefly popular. Here, the user would have a workstation rather than a terminal, yet work on data stored on a central server. One product along these lines was Sun's *Sun Ray* (circa 1999) where users relied on a smartcard to establish their local environment on an arbitrary, otherwise stateless, workstation.

### 2.11.2.1 Usage scenarios

The model where services are available on demand is attractive for businesses, which increasingly are using cloud services. The advantages are that it requires no initial monetary and time investment, and that no decisions about type and size of equipment have to be made. At the moment, cloud services are mostly focused on databases and office applications, but scientific clouds with a high performance interconnect are under development.

The following is a broad classification of usage scenarios of cloud resources<fnref n="23" />.

* Scaling. Here the cloud resources are used as a platform that can be expanded based on user demand. This can be considered Platform-as-a-Service (PAS): the cloud provides software and development platforms, eliminating the administration and maintenance for the user. We can distinguish between two cases: if the user is running single jobs and is actively waiting for the output, resources can be added to minimize the wait time for these jobs (capability computing). On the other hand, if the user is submitting jobs to a queue and the time-to-completion of any given job is not crucial (capacity computing), resources can be added as the queue grows. In HPC applications, users can consider the cloud resources as a cluster; this falls under Infrastructure-as-a-Service (IAS): the cloud service is a computing platform allowing customization at the operating system level.

* Multi-tenancy. Here the same software is offered to multiple users, giving each the opportunity for individual customizations. This falls under Software-as-a-Service (SAS): software is provided on demand; the customer does not purchase software, but only pays for its use.

<footnote n="23">23. Based on a blog post by Ricky Ho: http://blogs.globallogic.com/five-cloud-computing-patterns.</footnote>
* Batch processing. This is a limited version of one of the Scaling scenarios above: the user has a large amount of data to process in batch mode. The cloud then becomes a batch processor. This model is a good candidate for MapReduce computations; section 2.11.4.

* Storage. Most cloud providers offer database services, so this model absolves the user from maintaining their own database, just like the Scaling and Batch processing models take away the user's concern with maintaining cluster hardware.

* Synchronization. This model is popular for commercial user applications. Netflix and Amazon's Kindle allow users to consume online content (streaming movies and ebooks respectively); after pausing the content they can resume from any other platform. Apple's recent iCloud provides synchronization for data in office applications, but unlike Google Docs the applications are not 'in the cloud' but on the user machine.

The first Cloud to be publicly accessible was Amazon's Elastic Compute cloud (EC2), launched in 2006. EC2 offers a variety of different computing platforms and storage facilities. Nowadays more than a hundred companies provide cloud based services, well beyond the initial concept of computers-for-rent.

The infrastructure for cloud computing can be interesting from a computer science point of view, involving distributed file systems, scheduling, virtualization, and mechanisms for ensuring high reliability.

An interesting project, combining aspects of grid and cloud computing is the Canadian Advanced Network For Astronomical Research[144]. Here large central datasets are being made available to astronomers as in a grid, together with compute resources to perform analysis on them, in a cloud-like manner. Interestingly, the cloud resources even take the form of user-configurable virtual clusters.

### 2.11.2.2 Characterization

Summarizing<fnref n="24" /> we have three cloud computing service models:

**Software as a Service** The consumer runs the provider's application, typically through a thin client such as a browser; the consumer does not install or administer software. A good example is Google Docs

**Platform as a Service** The service offered to the consumer is the capability to run applications developed by the consumer, who does not otherwise manage the processing platform or data storage involved.

**Infrastructure as a Service** The provider offers to the consumer both the capability to run software, and to manage storage and networks. The consumer can be in charge of operating system choice and network components such as firewalls.

These can be deployed as follows:

**Private cloud** The cloud infrastructure is managed by one organization for its own exclusive use.

**Public cloud** The cloud infrastructure is managed for use by a broad customer base.

One could also define hybrid models such as community clouds.

The characteristics of cloud computing are then:

<footnote n="24">24. The remainder of this section is based on the NIST definition of cloud computing http://csrc.nist.gov/publications/nistpubs/800-145/SP800-145.pdf.</footnote>
**On-demand and self service** The consumer can quickly request services and change service levels, without requiring human interaction with the provider.

**Rapid elasticity** The amount of storage or computing power appears to the consumer to be unlimited, subject only to budgeting constraints. Requesting extra facilities is fast, in some cases automatic.

**Resource pooling** Virtualization mechanisms make a cloud appear like a single entity, regardless its underlying infrastructure. In some cases the cloud remembers the 'state' of user access; for instance, Amazon's Kindle books allow one to read the same book on a PC, and a smartphone; the cloud-stored book 'remembers' where the reader left off, regardless the platform.

**Network access** Clouds are available through a variety of network mechanisms, from web browsers to dedicated portals.

**Measured service** Cloud services are typically 'metered', with the consumer paying for computing time, storage, and bandwidth.

### 2.11.3 Capability versus capacity computing

Large parallel computers can be used in two different ways. In later chapters you will see how scientific problems can be scaled up almost arbitrarily. This means that with an increasing need for accuracy or scale, increasingly large computers are needed. The use of a whole machine for a single problem, with only time-to-solution as the measure of success, is known as *capability computing*.

On the other hand, many problems need less than a whole supercomputer to solve, so typically a computing center will set up a machine so that it serves a continuous stream of user problems, each smaller than the full machine. In this mode, the measure of success is the sustained performance per unit cost. This is known as *capacity computing*, and it requires a finely tuned *job scheduling* strategy.

A popular scheme is *fair-share scheduling*, which tries to allocate resources equally between users, rather than between processes. This means that it will lower a user's priority if the user had recent jobs, and it will give higher priority to short or small jobs. Examples of schedulers on this principle are *SGE* and *Slurm*.

Jobs can have dependencies, which makes scheduling harder. In fact, under many realistic conditions scheduling problems are *NP-complete*, so in practice heuristics will be used. This topic, while interesting, is not further discussed in this book.

### 2.11.4 MapReduce

MapReduce [33] is a programming model for certain parallel operations. One of its distinguishing characteristics is that it is implemented using *functional programming*. The MapReduce model handles computations of the following form:

* For all available data, select items that satisfy a certain criterium;

* and emit a key-value pair for them. This is the mapping stage.

* Optionally there can be a combine/sort stage where all pairs with the same key value are grouped together.

* Then do a global reduction on the keys, yielding one or more of the corresponding values. This is the reduction stage.

We will now give a few examples of using MapReduce, and present the functional programming model that underlies the MapReduce abstraction.
### 2.11.4.1 Expressive power of the MapReduce model

The reduce part of the MapReduce model makes it a prime candidate for computing global statistics on a dataset. One example would be to count how many times each of a set of words appears in some set of documents. The function being mapped knows the set of words, and outputs for each document a pair of document name and a list with the occurrence counts of the words. The reduction then does a component-wise sum of the occurrence counts.

The combine stage of MapReduce makes it possible to transform data. An example is a 'Reverse Web-Link Graph': the map function outputs target-source pairs for each link to a target URL found in a page named "source". The reduce function concatenates the list of all source URLs associated with a given target URL and emits the pair target-list(source).

A less obvious example is computing PageRank (section 9.3) with MapReduce. Here we use the fact that the PageRank computation relies on a distributed sparse matrix-vector product. Each web page corresponds to a column of the web matrix $W$; given a probability $p_j$ of being on page $j$, that page can then compute tuples $(i, w_{ij}p_j)$. The combine stage of MapReduce then sums together $(Wp)_i = \sum_j w_{ij}p_j$.

Database operations can be implemented with MapReduce but since it has a relatively large latency, it is unlikely to be competitive with standalone databases, which are optimized for fast processing of a single query, rather than bulk statistics. For other applications see http://horicky.blogspot.com/2010/08/designing-algorithmis-for-map-reduce.html.

### 2.11.4.2 Mapreduce software

The implementation of MapReduce by Google was released under the name *Hadoop*. While it suited the Google model of single-stage reading and processing of data, it had considerable disadvantages for many other users:

* Hadoop would flush all its data back to disc after each MapReduce cycle, so for operations that take more than a single cycle the file system and bandwidth demands are too great.

* In computing center environments, where a user's data is not continuously online, the time required for loading data into *Hadoop File System (HDFS)* would like overwhelm the actual analysis.

For these reasons, further projects such as *Apache Spark* (https://spark.apache.org/) offer caching of data.

### 2.11.4.3 Functional programming

The mapping and reduction operations are easily implemented on any type of parallel architecture, using a combination of threading and message passing. However, at Google where this model was developed traditional parallelism was not attractive for two reasons. First of all, processors could fail during the computation, so a traditional model of parallelism would have to be enhanced with *fault tolerance* mechanisms. Secondly, the computing hardware could already have a load, so parts of the computation may need to be migrated, and in general any type of synchronization between tasks would be very hard.
MapReduce is one way to abstract from such details of parallel computing, namely through adopting a functional programming model. In such a model the only operation is the evaluation of a function, applied to some arguments, where the arguments are themselves the result of a function application, and the result of the computation is again used as argument for another function application. In particular, in a strict functional model there are no variables, so there is no static data.

A function application, written in Lisp style as `(f a b)` (meaning that the function `f` is applied to arguments `a` and `b`) would then be executed by collecting the inputs from whereven they are to the processor that evaluates the function `f`. The mapping stage of a MapReduce process is denoted

`(map f (some list of arguments))`

and the result is a list of the function results of applying `f` to the input list. All details of parallelism and of guaranteeing that the computation successfully finishes are handled by the `map` function.

Now we are only missing the reduction stage, which is just as simple:

`(reduce g (map f (the list of inputs)))`

The `reduce` function takes a list of inputs and performs a reduction on it.

The attractiveness of this functional model lies in the fact that functions can not have *side effects*: because they can only yield a single output result, they can not change their environment, and hence there is no coordination problem of multiple tasks accessing the same data.

Thus, MapReduce is a useful abstraction for programmers dealing with large amounts of data. Of course, on an implementation level the MapReduce software uses familiar concepts such as decomposing the data space, keeping a work list, assigning tasks to processors, retrying failed operations, et cetera.

### 2.11.5 The top500 list

There are several informal ways of measuring just 'how big' a computer is. The most popular is the TOP500 list, maintained at http://www.top500.org/, which records a computer's performance on the Linpack benchmark. Linpack is a package for linear algebra operations, and no longer in use, since it has been superseded by Lapack for shared memory and Scalapack for distritubed memory computers. The benchmark operation is the solution of a (square, nonsingular, dense) linear system through LU factorization with partial pivoting, with subsequent forward and backward solution.

The LU factorization operation is one that has great opportunity for cache reuse, since it is based on the matrix-matrix multiplication kernel discussed in section 1.5.1. It also has the property that the amount of work outweighs the amount of communication: $O(n^3)$ versus $O(n^2)$. As a result, the Linpack benchmark is likely to run at a substantial fraction of the peak speed of the machine. Another way of phrasing this is to say that the Linpack benchmark is a *CPU-bound* or *compute-bound* algorithm.

Typical efficiency figures are between 60 and 90 percent. However, it should be noted that many scientific codes do not feature the dense linear solution kernel, so the performance on this benchmark is not indicative of the performance on a typical code. Linear system solution through iterative methods (section 5.5), for
instance, is much less efficient in a flops-per-second sense, being dominated by the bandwidth between CPU and memory (a *bandwidth bound algorithm*).

One implementation of the Linpack benchmark that is often used is ‘High-Performance LINPACK’ (http://www.netlib.org/benchmark/hpl/), which has several parameters such as blocksize that can be chosen to tune the performance.

### 2.11.5.1 The top500 list as a recent history of supercomputing

The top500 list offers a history of almost 20 years of supercomputing. In this section we will take a brief look at historical developments<fnref n="25" />. First of all, figure 2.36 shows the evolution of architecture types by

<figure><img src="image_69.png" /><figcaption>Figure 2.36: Evolution of the architecture types on the top500 list</figcaption></figure>

charting what portion of the aggregate peak performance of the whole list is due to each type.

* Vector machines feature a relatively small number of very powerful vector pipeline processors (section 2.4.1.1). This type of architecture has largely disappeared; the last major machine of this type was the Japanese *Earth Simulator* which is seen as the spike in the graph around 2002, and which was at the top of the list for two years.

* Micro-processor based architectures get their power from the large number of processors in one machine. The graph distinguishes between x86 (*Intel* and *AMD* processors with the exception of the *Intel Itanium*) processors and others; see also the next graph.

<footnote n="25">25. The graphs contain John McCalpin’s analysis of the top500 data.</footnote>
* A number of systems were designed as highly scalable architectures: these are denoted MPP for 'massively parallel processor'. In the early part of the timeline this includes architectures such as the *Connection Machine*, later it is almost exclusively the *IBM BlueGene*.

* In recent years 'accelerated systems' are the upcoming trend. Here, a processing unit such as a GPU is attached to the networked main processor.

Next, figure 2.37 shows the dominance of the x86 processor type relative to other micro-processors. (Since

<figure><img src="image_70.png" /><figcaption>Figure 2.37: Evolution of the architecture types on the top500 list</figcaption></figure>

we classified the *IBM BlueGene* as an MPP, its processors are not in the 'Power' category here.)

Finally, figure 2.38 shows the gradual increase in core count. Here we can make the following observations:

* In the 1990s many processors consisted of more than one chip. In the rest of the graph, we count the number of cores per 'package', that is, per *socket*. In some cases a socket can actually contain two separate dies.

* With the advent of multi-core processors, it is remarkable how close to vertical the section in the graph are. This means that new processor types are very quickly adopted, and the lower core counts equally quickly completely disappear.

* For accelerated systems (mostly systems with GPUs) the concept of 'core count' is harder to define; the graph merely shows the increasing importance of this type of architecture.
<figure><img src="image_71.png" /><figcaption>Figure 2.38: Evolution of the architecture types on the top500 list</figcaption></figure>

### 2.11.6 Heterogeneous computing

You have now seen several computing models: single core, shared memory multicore, distributed memory clusters, GPUs. These models all have in common that, if there is more than one instruction stream active, all streams are interchangable. With regard to GPUs we need to refine this statement: all instruction stream on the GPU are interchangable. However, a GPU is not a standalone device, but can be considered a co-processor to a host processor.

If we want to let the host perform useful work while the co-processor is active, we now have two different instruction streams or types of streams. This situation is known as heterogeneous computing. In the GPU case, these instruction streams are even programmed by a slightly different mechanisms – using CUDA for the GPU – but this need not be the case: the Intel Many Integrated Cores (MIC) architecture is programmed in ordinary C.
# Chapter 3

## Computer Arithmetic

Of the various types of data that one normally encounters, the ones we are concerned with in the context of scientific computing are the numerical types: integers (or whole numbers) ..., -2, -1, 0, 1, 2, ..., real numbers 0, 1, -1.5, 2/3, √2, log 10, ..., and complex numbers 1 + 2i, √3 - √5i, .... Computer memory is organized to give only a certain amount of space to represent each number, in multiples of *bytes*, each containing 8 *bits*. Typical values are 4 bytes for an integer, 4 or 8 bytes for a real number, and 8 or 16 bytes for a complex number.

Since only a certain amount of memory is available to store a number, it is clear that not all numbers of a certain type can be stored. For instance, for integers only a range is stored. In the case of real numbers, even storing a range is not possible since any interval [*a*, *b*] contains infinitely many numbers. Therefore, any *representation of real numbers* will cause gaps between the numbers that are stored. Calculations in a computer are sometimes described as *finite precision arithmetic*. Since many results are not representable, any computation that results in such a number will have to be dealt with by issuing an error or by approximating the result. In this chapter we will look at the ramifications of such approximations of the ‘true’ outcome of numerical calculations.

For detailed discussions, see the book by Overton [132]; it is easy to find online copies of the essay by Goldberg [68]. For extensive discussions of round-off error analysis in algorithms, see the books by Higham [86] and Wilkinson [164].

### 3.1 Integers

In scientific computing, most operations are on real numbers. Computations on integers rarely add up to any serious computation load<fnref n="1" />. It is mostly for completeness that we start with a short discussion of integers.

Integers are commonly stored in 16, 32, or 64 bits, with 16 becoming less common and 64 becoming more and more so. The main reason for this increase is not the changing nature of computations, but the fact that integers are used to index arrays. As the size of data sets grows (in particular in parallel computations), larger indices are needed. For instance, in 32 bits one can store the numbers zero through $2^{32} - 1 \approx 4 \cdot 10^9$.

<footnote n="1">1. Some computations are done on bit strings. We will not mention them at all.</footnote>
In other words, a 32 bit index can address 4 gigabytes of memory. Until recently this was enough for most purposes; these days the need for larger data sets has made 64 bit indexing necessary.

When we are indexing an array, only positive integers are needed. In general integer computations, of course, we need to accommodate the negative integers too. We will now discuss several strategies for implementing negative integers. Our motivation here will be that arithmetic on positive and negative integers should be as simple as on positive integers only: the circuitry that we have for comparing and operating on bitstrings should be usable for (signed) integers.

There are several ways of implementing negative integers. The simplest solution is to reserve one bit as a sign bit, and use the remaining 31 (or 15 or 63; from now on we will consider 32 bits the standard) bits to store the absolute magnitude. By comparison, we will call the straightforward interpretation of bitstring unsigned integers.

<table><tr><td>bitstring</td><td>00 ... 0</td><td>...</td><td>01 ... 1</td><td>10 ... 0</td><td>...</td><td>11 ... 1</td></tr><tr><td>interpretation as unsigned int</td><td>0</td><td>...</td><td>2<sup>31</sup> - 1</td><td>2<sup>31</sup></td><td>...</td><td>2<sup>32</sup> - 1</td></tr><tr><td>interpretation as signed integer</td><td>0</td><td>...</td><td>2<sup>31</sup> - 1</td><td>-0</td><td>...</td><td>-(2<sup>31</sup> - 1)</td></tr></table>

This scheme has some disadvantages, one being that there is both a positive and negative number zero. This means that a test for equality becomes more complicated than simply testing for equality as a bitstring. More importantly, in the second half of the bitstrings, the interpretation as signed integer decreases, going to the right. This means that a test for greater-than becomes complex; also adding a positive number to a negative number now has to be treated differently from adding it to a positive number.

Another solution would be to let an unsigned number $n$ be interpreted as $n - B$ where $B$ is some plausible base, for instance $2^{31}$.

<table><tr><td>bitstring</td><td>00 ... 0</td><td>...</td><td>01 ... 1</td><td>10 ... 0</td><td>...</td><td>11 ... 1</td></tr><tr><td>interpretation as unsigned int</td><td>0</td><td>...</td><td>2<sup>31</sup> - 1</td><td>2<sup>31</sup></td><td>...</td><td>2<sup>32</sup> - 1</td></tr><tr><td>interpretation as shifted int</td><td>-2<sup>31</sup></td><td>...</td><td>-1</td><td>0</td><td>...</td><td>2<sup>31</sup> - 1</td></tr></table>

This shifted scheme does not suffer from the $\pm 0$ problem, and numbers are consistently ordered. However, if we compute $n - n$ by operating on the bitstring that represents $n$, we do not get the bitstring for zero. To get we rotate the number line to put the pattern for zero back at zero.

The resulting scheme, which is the one that is used most commonly, is called 2's complement. Using this scheme, the representation of integers is formally defined as follows.

* If $0 \le m \le 2^{31} - 1$, the normal bit pattern for $m$ is used.

* For $-2^{31} \le n \le -1$, $n$ is represented by the bit pattern for $2^{32} - |n|$.

The following diagram shows the correspondence between bitstrings and their interpretation as 2's complement integer:

<table><tr><td>bitstring</td><td>00 ... 0</td><td>...</td><td>01 ... 1</td><td>10 ... 0</td><td>...</td><td>11 ... 1</td></tr><tr><td>interpretation as unsigned int</td><td>0</td><td>...</td><td>2<sup>31</sup> - 1</td><td>2<sup>31</sup></td><td>...</td><td>2<sup>32</sup> - 1</td></tr><tr><td>interpretation as 2's comp. integer</td><td>0</td><td>...</td><td>2<sup>31</sup> - 1</td><td>-2<sup>31</sup></td><td>...</td><td>-1</td></tr></table>
Some observations:

* There is no overlap between the bit patterns for positive and negative integers, in particular, there is only one pattern for zero.

* The positive numbers have a leading bit zero, the negative numbers have the leading bit set.

Exercise 3.1. For the ‘naive’ scheme and the 2’s complement scheme for negative numbers, give pseudocode for the comparison test $m < n$, where $m$ and $n$ are integers. Be careful to distinguish between all cases of $m$, $n$ positive, zero, or negative.

Adding two numbers with the same sign, or multiplying two numbers of any sign, may lead to a result that is too large or too small to represent. This is called *overflow*; see section 3.2.3 for details.

Exercise 3.2. Investigate what happens when you perform such a calculation. What does your compiler say if you try to write down a nonrepresentible number explicitly, for instance in an assignment statement?

In exercise 3.1 above you explored comparing two integers. Let us know how subtracting numbers in two’s complement is implemented. Consider $0 \le m \le 2^{31} - 1$ and $1 \le n \le 2^{31}$ and let us see what happens in the computation of $m - n$.

Suppose we have an algorithm for adding and subtracting unsigned 32-bit numbers. Can we use that to subtract two’s complement integers? We start by observing that the integer subtraction $m - n$ becomes the unsigned addition $m + (2^{32} - n)$.

* Case: $m < n$. In this case, $m - n$ is negative and $1 \le |m - n| \le 2^{31}$, so the bit pattern for $m - n$ is that of $2^{32} - (n - m)$. Now, $2^{32} - (n - m) = m + (2^{32} - n)$, so we can compute $m - n$ in 2’s complement by adding the bit patterns of $m$ and $-n$ as unsigned integers.

* Case: $m > n$. Here we observe that $m + (2^{32} - n) = 2^{32} + m - n$. Since $m - n > 0$, this is a number $> 2^{32}$ and therefore not a legitimate representation of a negative number. However, if we store this number in 33 bits, we see that it is the correct result $m - n$, plus a single bit in the 33-rd position. Thus, by performing the unsigned addition, and ignoring the *overflow bit*, we again get the correct result.

In both cases we conclude that we can perform the subtraction $m - n$ by adding the unsigned number that represent $m$ and $-n$ and ignoring overflow if it occurs.

## 3.2 Real numbers

In this section we will look at how real numbers are represented in a computer, and the limitations of various schemes. The next section will then explore the ramifications of this for arithmetic involving computer numbers.

### 3.2.1 They're not really real numbers

In the mathematical sciences, we usually work with real numbers, so it's convenient to pretend that computers can do this too. However, since numbers in a computer have only a finite number of bits, most real numbers can not be represented exactly. In fact, even many fractions can not be represented exactly, since they repeat; for instance, $1/3 = 0.333...$. An illustration of this is given in appendix 31.6.
**Exercise 3.3.** Some programming languages allow you to write loops with not just an integer, but also with a real number as 'counter'. Explain why this is a bad idea. Hint: when is the upper bound reached?

Whether a fraction repeats depends on the number system. (How would you write 1/3 in ternary, base 3, arithmetic?) In binary computers this means that fractions such as 1/10, which in decimal arithmetic are terminating, are repeating. Since decimal arithmetic is important in financial calculations, some people care about this kind of arithmetic right; see section 3.2.2.2 for what was done about it.

### 3.2.2 Representation of real numbers

Real numbers are stored using a scheme that is analogous to what is known as 'scientific notation', where a number is represented as a *significant* and an *exponent*, for instance $6.022 \cdot 10^{23}$, which has a significant 6022 with a *radix point* after the first digit, and an exponent 23. This number stands for

$$6.022 \cdot 10^{23} = [6 \times 10^{0} + 0 \times 10^{-1} + 2 \times 10^{-2} + 2 \times 10^{-3}] \cdot 10^{23}.$$

We introduce a *base*, a small integer number, 10 in the preceding example, and 2 in computer numbers, and write numbers in terms of it as a sum of *t* terms:

$$x = \pm 1 \times [d_1\beta^0 + d_2\beta^{-1} + d_3\beta^{-2} + \dots + d_t\beta^{-t+1}b] \times \beta^e = \pm \sum_{i=1}^{t} d_i\beta^{1-i} \times \beta^e \quad (3.1)$$

where the components are

* the *sign bit*: a single bit storing whether the number is positive or negative;

* $\beta$ is the base of the number system;

* $0 \le d_i \le \beta - 1$ the digits of the *mantissa* or *significant* – the location of the radix point (decimal point in decimal numbers) is implicitly assumed to the immediately following the first digit;

* *t* is the length of the mantissa;

* $e \in [L, U]$ exponent; typically $L < 0 < U$ and $L \approx -U$.

Note that there is an explicit sign bit for the whole number; the sign of the exponent is handled differently. For reasons of efficiency, *e* is not a signed number; instead it is considered as an unsigned number in excess of a certain minimum value. For instance, the bit pattern for the number zero is interpreted as $e = L$.

#### 3.2.2.1 Some examples

Let us look at some specific examples of floating point representations. Base 10 is the most logical choice for human consumption, but computers are binary, so base 2 predominates there. Old IBM mainframes grouped bits to make for a base 16 representation.

<table><thead><tr><th></th><th>β</th><th>t</th><th>L</th><th>U</th></tr></thead><tbody><tr><td>IEEE single precision (32 bit)</td><td>2</td><td>24</td><td>-126</td><td>127</td></tr><tr><td>IEEE double precision (64 bit)</td><td>2</td><td>53</td><td>-1022</td><td>1023</td></tr><tr><td>Old Cray 64 bit</td><td>2</td><td>48</td><td>-16383</td><td>16384</td></tr><tr><td>IBM mainframe 32 bit</td><td>16</td><td>6</td><td>-64</td><td>63</td></tr><tr><td>packed decimal</td><td>10</td><td>50</td><td>-999</td><td>999</td></tr><tr><td>Setun</td><td>3</td><td></td><td></td><td></td></tr></tbody></table>

Of these, the single and double precision formats are by far the most common. We will discuss these in section 3.2.7 and further.
### 3.2.2.2 Binary coded decimal

Decimal numbers are not relevant in scientific computing, but they are useful in financial calculations, where computations involving money absolutely have to be exact. Binary arithmetic is at a disadvantage here, since numbers such as 1/10 are repeating fractions in binary. With a finite number of bits in the mantissa, this means that the number 1/10 can not be represented exactly in binary. For this reason, binary-coded-decimal schemes were used in old IBM mainframes, and are in fact being standardized in revisions of IEEE754 [90]; see also section 3.2.7.

In BCD schemes, one or more decimal digits are encoded in a number of bits. The simplest scheme would encode the digits 0...9 in four bits. This has the advantage that in a BCD number each digit is readily identified; it has the disadvantage that about 1/3 of all bits are wasted, since 4 bits can encode the numbers 0...15. More efficient encodings would encode 0...999 in ten bits, which could in principle store the numbers 0...1023. While this is efficient in the sense that few bits are wasted, identifying individual digits in such a number takes some decoding. For this reason, BCD arithmetic needs hardware support from the processor, which is rarely found these days; one example is the IBM Power architecture, starting with the IBM Power6.

### 3.2.2.3 Other number bases for computer arithmetic

There have been some experiments with ternary arithmetic (see http://en.wikipedia.org/wiki/Ternary_computer and http://www.computer-museum.ru/english/setun.htm), however, no practical hardware exists.

## 3.2.3 Limitations

Since we use only a finite number of bits to store floating point numbers, not all numbers can be represented. The ones that can not be represented fall into two categories: those that are too large or too small (in some sense), and those that fall in the gaps. Numbers can be too large or too small in the following ways.

**Overflow** The largest number we can store is

$$ (\beta - 1) \cdot 1 + (\beta - 1) \cdot \beta^{-1} + \cdots + (\beta - 1) \cdot \beta^{-(t-1)} = \beta - 1 \cdot \beta^{-(t-1)}, $$

and the smallest number (in an absolute sense) is $-(\beta - \beta^{-(t-1)}$; anything larger than the former or smaller than the latter causes a condition called overflow.

**Underflow** The number closest to zero is $\beta^{-(t-1)} \cdot \beta^L$. A computation that has a result less than that (in absolute value) causes a condition called underflow. (See section 3.2.4 for gradual underflow.)

The fact that only a small number of real numbers can be represented exactly is the basis of the field of round-off error analysis. We will study this in some detail in the following sections.

The occurrence of overflow or underflow means that your computation will be 'wrong' from that point on. Underflow will make the computation proceed with a zero where there should have been a nonzero; overflow is represented as $Inf$, short of 'infinite'. Computing with $Inf$ is possible to an extent: adding two of those quantities will again give $Inf$. However, subtracting them gives $NaN$: 'not a number'.
In none of these cases will the computation abort: the processor will continue, unless you tell it otherwise. The ‘otherwise’ consists of you telling the compiler to generate an *interrupt*, which halts the computation with an error message. In C, you specify this with a library call:

```c
#include <fenv.h>
int main() {
    ...
    feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
```

Interrupts can also be enabled by compiler options.

### 3.2.4 Normalized and unnormalized numbers

The general definition of floating point numbers, equation (3.1), leaves us with the problem that numbers have more than one representation. For instance, $.5 \times 10^2 = .05 \times 10^3$. Since this would make computer arithmetic needlessly complicated, for instance in testing equality of numbers, we use *normalized floating point numbers*. A number is normalized if its first digit is nonzero. The implies that the mantissa part is $\beta > x_m \ge 1$.

A practical implication in the case of binary numbers is that the first digit is always 1, so we do not need to store it explicitly. In the IEEE 754 standard, this means that every floating point number is of the form

$$1.d_1 d_2 \dots d_t \times 2^{exp}$$

and only the digits $d_1 d_2 \dots d_t$ are stored.

Another implication of this scheme is that we have to amend the definition of underflow (see section 3.2.3 above): instead of having a smallest number $-(\beta - \beta^{-(t-1)})$, any number less than $1 \cdot \beta^L$ now causes underflow. Trying to compute a number less than that in absolute value is sometimes handled by using *unnormalized floating point numbers*, a process known as *gradual underflow*. In this case, a special value of the exponent indicates that the number is no longer normalized. In the case IEEE standard arithmetic (section 3.2.7) this is done through a zero exponent field.

However, this is typically tens or hundreds of times slower than computing with regular floating point numbers<fnref n="2" />. At the time of this writing, only the IBM Power6 has hardware support for gradual underflow.

### 3.2.5 Representation error

Let us consider a real number that is not representable in a computer’s number system.

An unrepresentable number is approximated either by ordinary *rounding*, rounding up or down, or *truncation*. This means that a machine number $\tilde{x}$ is the representation for all $x$ in an interval around it. With $t$

<footnote n="2">2. In real-time applications such as audio processing this phenomenon is especially noticeable; see http://phonophunk.com/articles/pentium4-denormalization.php?pg=3.</footnote>
digits in the mantissa, this is the interval of numbers that differ from $\tilde{x}$ in the $t+1$st digit. For the mantissa part we get:

$$ \begin{cases} x \in [\tilde{x}, \tilde{x} + \beta^{-t+1}) & \text{truncation} \\ x \in [\tilde{x} - \frac{1}{2}\beta^{-t+1}, \tilde{x} + \frac{1}{2}\beta^{-t+1}) & \text{rounding} \end{cases} $$

If $x$ is a number and $\tilde{x}$ its representation in the computer, we call $x - \tilde{x}$ the *representation error* or *absolute representation error*, and $\frac{x-\tilde{x}}{x}$ the *relative representation error*. Often we are not interested in the sign of the error, so we may apply the terms error and relative error to $|x - \tilde{x}|$ and $|\frac{x-\tilde{x}}{x}|$ respectively.

Often we are only interested in bounds on the error. If $\epsilon$ is a bound on the error, we will write

$$ \tilde{x} = x \pm \epsilon \stackrel{\text{D}}{=} |x - \tilde{x}| \le \epsilon \Leftrightarrow \tilde{x} \in [x - \epsilon, x + \epsilon] $$

For the relative error we note that

$$ \tilde{x} = x(1 + \epsilon) \Leftrightarrow \left| \frac{\tilde{x} - x}{x} \right| \le \epsilon $$

Let us consider an example in decimal arithmetic, that is, $\beta = 10$, and with a 3-digit mantissa: $t = 3$. The number $x = 1.256$ has a representation that depends on whether we round or truncate: $\tilde{x}_{\text{round}} = 1.26$, $\tilde{x}_{\text{truncate}} = 1.25$. The error is in the 4th digit: if $\epsilon = x - \tilde{x}$ then $|\epsilon| < \beta^{t-1}$.

**Exercise 3.4.** The number in this example had no exponent part. What are the error and relative error if there had been one?

### 3.2.6 Machine precision

Often we are only interested in the order of magnitude of the representation error, and we will write $\tilde{x} = x(1 + \epsilon)$, where $|\epsilon| \le \beta^{-t}$. This maximum relative error is called the *machine precision*, or sometimes *machine epsilon*. Typical values are:

$$ \begin{cases} \epsilon \approx 10^{-7} & 32-bit \text{ single precision} \\ \epsilon \approx 10^{-16} & 64-bit \text{ double precision} \end{cases} $$

Machine precision can be defined another way: $\epsilon$ is the smallest number that can be added to 1 so that $1 + \epsilon$ has a different representation than 1. A small example shows how aligning exponents can shift a too small operand so that it is effectively ignored in the addition operation:

$$ \begin{array}{c} \begin{array}{cc} + & 1.0000 \times 10^0 \\ & 1.0000 \times 10^{-5} \end{array} \quad \Rightarrow \quad \begin{array}{cc} + & 1.0000 \times 10^0 \\ & 0.00001 \times 10^0 \\ - & 1.0000 \times 10^0 \end{array} \end{array} $$

Yet another way of looking at this is to observe that, in the addition $x + y$, if the ratio of $x$ and $y$ is too large, the result will be identical to $x$.

The machine precision is the maximum attainable accuracy of computations: it does not make sense to ask for more than 6-or-so digits accuracy in single precision, or 15 in double.
<figure><img src="image_72.png" /><figcaption>Figure 3.1: Single precision arithmetic</figcaption></figure>

Exercise 3.5. Write a small program that computes the machine epsilon. Does it make any difference if you set the *compiler optimization levels* low or high?

Exercise 3.6. The number $e \approx 2.72$, the base for the natural logarithm, has various definitions. One of them is

$$e = \lim_{n \to \infty} (1 + 1/n)^n.$$

Write a single precision program that tries to compute $e$ in this manner. Evaluate the expression for $n = 10^k$ with $k = 1, \dots, 10$. Explain the output for large $n$. Comment on the behaviour of the error.

### 3.2.7 The IEEE 754 standard for floating point numbers

Some decades ago, issues like the length of the mantissa and the rounding behaviour of operations could differ between computer manufacturers, and even between models from one manufacturer. This was obviously a bad situation from a point of portability of codes and reproducibility of results. The IEEE standard 754³<fnref n="4" /> codified all this, for instance stipulating 24 and 53 bits for the mantissa in single and double precision arithmetic, using a storage sequence of sign bit, exponent, mantissa.

The standard also declared the rounding behaviour to be *correct rounding*: the result of an operation should be the rounded version of the exact result. There will be much more on the influence of rounding (and truncation) on numerical computations, below.

<footnote n="3">3. IEEE 754 is a standard for binary arithmetic; there is a further standard, IEEE 854, that allows decimal arithmetic.</footnote>

<footnote n="4">4. “It was remarkable that so many hardware people there, knowing how difficult p754 would be, agreed that it should benefit the community at large. If it encouraged the production of floating-point software and eased the development of reliable software, it would help create a larger market for everyone’s hardware. This degree of altruism was so astonishing that MATLAB’s creator Dr. Cleve Moler used to advise foreign visitors not to miss the country’s two most awesome spectacles: the Grand Canyon, and meetings of IEEE p754.” W. Kahan, http://www.cs.berkeley.edu/~wkahan/ieee754status/754story.html.</footnote>
Above (section 3.2.3), we have seen the phenomena of overflow and underflow, that is, operations leading to unrepresentable numbers. There is a further exceptional situation that needs to be dealt with: what result should be returned if the program asks for illegal operations such as $\sqrt{-4}$? The IEEE 754 standard has two special quantities for this: $Inf$ and $NaN$ for ‘infinity’ and ‘not a number’. Infinity is the result of overflow or dividing by zero, not-a-number is the result of, for instance, subtracting infinity from infinity. If $NaN$ appears in an expression, the whole expression will evaluate to that value. The rule for computing with $Inf$ is a bit more complicated [68].

An inventory of the meaning of all bit patterns in IEEE 754 double precision is given in figure 3.1. Recall from section 3.2.4 above that for normalized numbers the first nonzero digit is a 1, which is not stored, so the bit pattern $d_1d_2...d_t$ is interpreted as $1.d_1d_2...d_t$.

**Exercise 3.7.** Every programmer, at some point in their life, makes the mistake of storing a real number in an integer or the other way around. This can happen for instance if you call a function differently from how it was defined.

```c
void a(double x) {...}
int main() {
    int i;
    .... a(i) ....
}
```

What happens when you print $x$ in the function? Consider the bit pattern for a small integer, and use the table in figure 3.1 to interpret it as a floating point number. Explain that it will be an unnormalized number<fnref n="5" />.

These days, almost all processors adhere to the IEEE 754 standard. Early generations of the NVidia Tesla GPUs were not standard-conforming in single precision. The justification for this was that single precision is more likely used for graphics, where exact compliance matters less. For many scientific computations, double precision is necessary, since the precision of calculations gets worse with increasing problem size or runtime. This is true for the sort of calculations in chapter 4, but not for others such as Lattice Boltzmann Method (LBM)

### 3.2.7.1 Setting the rounding behaviour

The 754 standard also declares that a processor should be able to switch its rounding behaviour between ordinary rounding, rounding up or down (sometimes phrased as ‘towards plus infinity’ and ‘towards minus infinity’ respectively) or truncation. In C99, the API for this is contained in `fenv.h`:

```c
#include <fenv.h>

int roundings[] ={FE_TONEAREST, FE_UPWARD, FE_DOWNWARD, FE_TOWARDZERO};
rchoice = ....
int status = fesetround(roundings[rchoice]);
```

<footnote n="5">5. This is one of those errors you won't forget after you make it. In the future, whenever you see a number on the order of $10^{-305}$ you'll recognize that you probably made this error.</footnote>
In Fortran2003 the function IEEE_SET_ROUNTING_MODE is available in the IEEE_ARITHMETIC module.

Setting the rounding behaviour can serve as a quick test for the stability of an algorithm: if the result changes appreciably between two different rounding strategies, the algorithm is likely not stable.

## 3.3 Round-off error analysis

Numbers that are too large or too small to be represented, leading to overflow and underflow, are uncommon: usually computations can be arranged so that this situation will not occur. By contrast, the case that the result of a computation between computer numbers (even something as simple as a single addition) is not representable is very common. Thus, looking at the implementation of an algorithm, we need to analyze the effect of such small errors propagating through the computation. This is commonly called *round-off error analysis*.

### 3.3.1 Correct rounding

The IEEE 754 standard, mentioned in section 3.2.7, does not only declare the way a floating point number is stored, it also gives a standard for the accuracy of operations such as addition, subtraction, multiplication, division. The model for arithmetic in the standard is that of *correct rounding*: the result of an operation should be as if the following procedure is followed:

* The exact result of the operation is computed, whether this is representable or not;

* This result is then rounded to the nearest computer number.

In short: the representation of the result of an operation is the rounded exact result of that operation. (Of course, after two operations it no longer needs to hold that the computed result is the exact rounded version of the exact result.)

If this statement sounds trivial or self-evident, consider subtraction as an example. In a decimal number system with two digits in the mantissa, the computation $1.0 - 9.4 \cdot 10^{-1} = 1.0 - 0.94 = 0.06 = 0.6 \cdot 10^{-2}$. Note that in an intermediate step the mantissa .094 appears, which has one more digit than the two we declared for our number system. The extra digit is called a *guard digit*.

Without a guard digit, this operation would have proceeded as $1.0 - 9.4 \cdot 10^{-1}$, where $9.4 \cdot 10^{-1}$ would be rounded to 0.9, giving a final result of 0.1, which is almost double the correct result.

**Exercise 3.8.** Consider the computation $1.0 - 9.5 \cdot 10^{-1}$, and assume again that numbers are rounded to fit the 2-digit mantissa. Why is this computation in a way a lot worse than the example?

One guard digit is not enough to guarantee correct rounding. An analysis that we will not reproduce here shows that three extra bits are needed [67].
### 3.3.1.1 Mul-Add operations

In 2008, the IEEE 754 standard was revised to include the behaviour of Fused Multiply-Add (FMA) operations, that is, operations of the form $a*b+c$. The standard here defines correct rounding to be that the result of this combined computation should be the rounded correct result. A naive implementation of this operations would involve two roundings: one after the multiplication and one after the addition.

**Exercise 3.9.** Can you come up with an example where correct rounding of an FMA is considerably more accurate than rounding the multiplication and addition seperately?

Hint: let the $c$ term be of opposite sign as $a*b$, and try to force cancellation in the subtraction.

Mul-add operations sometimes use the 80-bit extended precision format for great precision of the intermediate result.

## 3.3.2 Addition

Addition of two floating point numbers is done in a couple of steps. First the exponents are aligned: the smaller of the two numbers is written to have the same exponent as the larger number. Then the mantissas are added. Finally, the result is adjusted so that it again is a normalized number.

As an example, consider $1.00 + 2.00 \times 10^{-2}$. Aligning the exponents, this becomes $1.00 + 0.02 = 1.02$, and this result requires no final adjustment. We note that this computation was exact, but the sum $1.00 + 2.55 \times 10^{-2}$ has the same result, and here the computation is clearly not exact: the exact result is 1.0255, which is not representable with three digits to the mantissa.

In the example $6.15 \times 10^1 + 3.98 \times 10^1 = 10.13 \times 10^1 = 1.013 \times 10^2 \rightarrow 1.01 \times 10^2$ we see that after addition of the mantissas an adjustment of the exponent is needed. The error again comes from truncating or rounding the first digit of the result that does not fit in the mantissa: if $x$ is the true sum and $\tilde{x}$ the computed sum, then $\tilde{x} = x(1 + \epsilon)$ where, with a 3-digit mantissa $|\epsilon| < 10^{-3}$.

Formally, let us consider the computation of $s = x_1+x_2$, and we assume that the numbers $x_i$ are represented as $\tilde{x}_i = x_i(1 + \epsilon_i)$. Then the sum $s$ is represented as

$$ 
\begin{aligned}
\tilde{s} &= (\tilde{x}_1 + \tilde{x}_2)(1 + \epsilon_3) \\
&= x_1(1 + \epsilon_1)(1 + \epsilon_3) + x_2(1 + \epsilon_2)(1 + \epsilon_3) \\
&\approx x_1(1 + \epsilon_1 + \epsilon_3) + x_2(1 + \epsilon_2 + \epsilon_3) \\
&\approx s(1 + 2\epsilon)
\end{aligned}
$$

under the assumptions that all $\epsilon_i$ are small and of roughly equal size, and that both $x_i > 0$. We see that the relative errors are added under addition.

## 3.3.3 Multiplication

Floating point multiplication, like addition, involves several steps. In order to multiply two numbers $m_1 \times \beta^{e_1}$ and $m_2 \times \beta^{e_2}$, the following steps are needed.

* The exponents are added: $e \leftarrow e_1 + e_2$.
• The mantissas are multiplied: $m \leftarrow m_1 \times m_2$.

• The mantissa is normalized, and the exponent adjusted accordingly.

For example: $1.23 \cdot 10^0 \times 5.67 \cdot 10^1 = 0.69741 \cdot 10^1 \rightarrow 6.9741 \cdot 10^0 \rightarrow 6.97 \cdot 10^0$.

**Exercise 3.10.** Analyze the relative error of multiplication.

### 3.3.4 Subtraction

Subtraction behaves very differently from addition. Whereas in addition errors are added, giving only a gradual increase of overall roundoff error, subtraction has the potential for greatly increased error in a single operation.

For example, consider subtraction with 3 digits to the mantissa: $1.24-1.23 = 0.01 \rightarrow 1.00 \cdot 10^{-2}$. While the result is exact, it has only one significant digit<fnref n="6" />. To see this, consider the case where the first operand 1.24 is actually the rounded result of a computation that should have resulted in 1.235. In that case, the result of the subtraction should have been $5.00 \cdot 10^{-3}$, that is, there is a 100% error, even though the relative error of the inputs was as small as could be expected. Clearly, subsequent operations involving the result of this subtraction will also be inaccurate. We conclude that subtracting almost equal numbers is a likely cause of numerical roundoff.

There are some subtleties about this example. Subtraction of almost equal numbers is exact, and we have the correct rounding behaviour of IEEE arithmetic. Still, the correctness of a single operation does not imply that a sequence of operations containing it will be accurate. While the addition example showed only modest decrease of numerical accuracy, the cancellation in this example can have disastrous effects. You'll see an example in section 3.3.6.1.

### 3.3.5 Associativity

Another effect of the way floating point numbers are treated is on the *associativity* of mathematical operations. Operations such as addition and multiplication, that ought to be associative, are no longer so in computer arithmetic.

Let's consider a simple example, showing how this can be caused by the rounding behaviour of floating point numbers. Let floating point numbers be stored as a single digit for the mantissa, one digit for the exponent, and one *guard digit*; now consider the computation of $4+6+7$. Evaluation left-to-right gives:

$$ 
\begin{aligned}
(4 \cdot 10^0 + 6 \cdot 10^0) + 7 \cdot 10^0 &\Rightarrow 10 \cdot 10^0 + 7 \cdot 10^0 && \text{addition} \\
&\Rightarrow 1 \cdot 10^1 + 7 \cdot 10^0 && \text{rounding} \\
&\Rightarrow 1.0 \cdot 10^1 + 0.7 \cdot 10^1 && \text{using guard digit} \\
&\Rightarrow 1.7 \cdot 10^1 \\
&\Rightarrow 2 \cdot 10^1 && \text{rounding}
\end{aligned}
$$

<footnote n="6">6. Normally, a number with 3 digits to the mantissa suggests an error corresponding to rounding or truncating the fourth digit. We say that such a number has 3 significant digits. In this case, the last two digits have no meaning, resulting from the normalization process.</footnote>
On the other hand, evaluation right-to-left gives:

$$
\begin{align*}
4 \cdot 10^0 + (6 \cdot 10^0 + 7 \cdot 10^0) &\Rightarrow 4 \cdot 10^0 + 13 \cdot 10^0 \quad \text{addition} \\
&\Rightarrow 4 \cdot 10^0 + 1 \cdot 10^1 \quad \text{rounding} \\
&\Rightarrow 0.4 \cdot 10^1 + 1.0 \cdot 10^1 \quad \text{using guard digit} \\
&\Rightarrow 1.4 \cdot 10^1 \\
&\Rightarrow 1 \cdot 10^1 \quad \text{rounding}
\end{align*}
$$

The conclusion is that the sequence in which rounding and truncation is applied to intermediate results
makes a difference. You can also observe that starting with the smaller numbers gives a more accurate
result. In section 3.3.6.2 you'll see a more elaborate example of this principle.

Usually, the evaluation order of expressions is determined by the definition of the programming language, or at least by the compiler. In section 3.3.7 we will see how in parallel computations the associativity is not so uniquely determined.

**3.3.6 Examples**

From the above, the reader may got the impression that roundoff errors only lead to serious problems in
exceptional circumstances. In this section we will discuss some very practical examples where the inex-
actness of computer arithmetic becomes visible in the result of a computation. These will be fairly simple
examples; more complicated examples exist that are outside the scope of this book, such as the instability
of matrix inversion. The interested reader is referred to [164, 86].

3.3.6.1 The ‘abc-formula’

As a practical example, consider the quadratic equation $ax^2 + bx + c = 0$ which has solutions $x = \frac{-b \pm \sqrt{b^2 - 4ac}}{2a}$. Suppose $b > 0$ and $b^2 \gg 4ac$ then $\sqrt{b^2 - 4ac} \approx b$ and the ‘+’ solution will be inaccurate. In this case it is better to compute $x_- = \frac{-b - \sqrt{b^2 - 4ac}}{2a}$ and use $x_+ \cdot x_- = -c/a$.

Exercise 3.11. Write a program that computes the roots of the quadratic equation, both the
'textbook' way, and as described above.

* Let $b = 1$ and $a = -c$, and $4ac \downarrow 0$ by taking progressively smaller values for $a$ and $c$.

* Print out the computed root, the root using the stable computation, and the value of $f(x) = ax^2 + bx + c$ in the computed root.

Now suppose that you don't care much about the actual value of the root: you want
to make sure the residual $f(x)$ is small in the computed root. Let $x^*$ be the exact
root, then

$$
f(x^* + h) \approx f(x^*) + hf'(x^*) = hf'(x^*).
$$

Now investigate separately the cases $a \downarrow 0, c = -1$ and $a = -1, c \downarrow 0$. Can you
explain the difference?
### 3.3.6.2 Summing series

The previous example was about preventing a large roundoff error in a single operation. This example shows that even gradual buildup of roundoff error can be handled in different ways.

Consider the sum $\sum_{n=1}^{10000} \frac{1}{n^2} = 1.644834$ and assume we are working with single precision, which on most computers means a machine precision of $10^{-7}$. The problem with this example is that both the ratio between terms, and the ratio of terms to partial sums, is ever increasing. In section 3.2.6 we observed that a too large ratio can lead to one operand of an addition in effect being ignored.

If we sum the series in the sequence it is given, we observe that the first term is 1, so all partial sums ($\sum_{n=1}^{N}$ where $N < 10000$) are at least 1. This means that any term where $1/n^2 < 10^{-7}$ gets ignored since it is less than the machine precision. Specifically, the last 7000 terms are ignored, and the computed sum is 1.644725. The first 4 digits are correct.

However, if we evaluate the sum in reverse order we obtain the exact result in single precision. We are still adding small quantities to larger ones, but now the ratio will never be as bad as one-to-$\epsilon$, so the smaller number is never ignored. To see this, consider the ratio of two terms subsequent terms:

$$ \frac{n^2}{(n-1)^2} = \frac{n^2}{n^2 - 2n + 1} = \frac{1}{1 - 2/n + 1/n^2} \approx 1 + \frac{2}{n} $$

Since we only sum $10^5$ terms and the machine precision is $10^{-7}$, in the addition $1/n^2 + 1/(n-1)^2$ the second term will not be wholly ignored as it is when we sum from large to small.

**Exercise 3.12.** There is still a step missing in our reasoning. We have shown that in adding two subsequent terms, the smaller one is not ignored. However, during the calculation we add partial sums to the next term in the sequence. Show that this does not worsen the situation.

The lesson here is that series that are monotone (or close to monotone) should be summed from small to large, since the error is minimized if the quantities to be added are closer in magnitude. Note that this is the opposite strategy from the case of subtraction, where operations involving similar quantities lead to larger errors. This implies that if an application asks for adding and subtracting series of numbers, and we know a priori which terms are positive and negative, it may pay off to rearrange the algorithm accordingly.

### 3.3.6.3 Unstable algorithms

We will now consider an example where we can give a direct argument that the algorithm can not cope with problems due to inexactely represented real numbers.

Consider the recurrence $y_n = \int_0^1 \frac{x^n}{x-5} dx = \frac{1}{n} - 5y_{n-1}$. This is easily seen to be monotonically decreasing; the first term can be computed as $y_0 = \ln 6 - \ln 5$.

Performing the computation in 3 decimal digits we get: