Assertions serve to detect bugs in your program: an *assertion* is a predicate that should be true at a certain point in your program. Thus, an assertion failing means that you didn't code what you intended to code. An assertion is typically a statement in your programming language, or a preprocessor macro; upon failure of the assertion, your program will take some abortive action.

Some examples of assertions:

* If a subprogram has an array argument, it is a good idea to test whether the actual argument is a null pointer before indexing into the array.

* Similarly, you could test a dynamically allocated data structure for not having a null pointer.

* If you calculate a numerical result for which certain mathematical properties hold, for instance you are writing a sine function, for which the result has to be in [-1, 1], you should test whether this property indeed holds for the result.

Assertions are often disabled in a program once it's sufficiently tested. The reason for this is that assertions can be expensive to execute. For instance, if you have a complicated data structure, you could write a complicated integrity test, and perform that test in an assertion, which you put after every access to the data structure.

Because assertions are often disabled in the 'production' version of a code, they should not affect any stored data. If they do, your code may behave differently when you're testing it with assertions, versus how you use it in practice without them. This is also formulated as 'assertions should not have side-effects'.

### 26.1.1.1 The C assert macro

The C standard library has a file assert.h which provides an assert() macro. Inserting assert(foo) has the following effect: if foo is zero (false), a diagnostic message is printed on standard error:

```
Assertion failed: foo, file filename, line line-number
```

which includes the literal text of the expression, the file name, and line number; and the program is subsequently aborted. Here is an example:

```c
#include assert.h>

void open_record(char *record_name)
{
    assert(record_name!=NULL);
    /* Rest of code */
}

int main(void)
{
    open_record(NULL);
}
```

The assert macro can be disabled by defining the NDEBUG macro.
26.1.1.2 An assert macro for Fortran

(Thanks to Robert Mclay for this code.)

#if (defined( GFORTRAN ) || defined( G95 ) || defined ( PGI ) )
# define MKSTR(x) "x"
#else
# define MKSTR(x) #x
#endif
#ifndef NDEBUG
# define EXPRESS(x, msg) if (.not. (x) ) \
    call assert( FILE , LINE ,MKSTR(x),msg)
#else
# define EXPRESS(x, msg)
#endif
subroutine assert(file, ln, testStr, msgIn)
implicit none
character(*) :: file, testStr, msgIn
integer :: ln
print *, "Assert: ",trim(testStr)," Failed at ",trim(file),":",ln
print *, "Msg:", trim(msgIn)
stop
end subroutine assert

which is used as

asser (nItemsSet .gt. arraySize, "Too many elements set")

26.1.2 Use of error codes

In some software libraries (for instance MPI or PETSc) every subprogram returns a result, either the func-
tion value or a parameter, to indicate success or failure of the routine. It is good programming practice to
check these error parameters, even if you think that nothing can possibly go wrong.

It is also a good idea to write your own subprograms in such a way that they always have an error parameter.
Let us consider the case of a function that performs some numerical computation.

float compute(float val)
{
    float result;
    result = ... /* some computation */
    return result;
}

float value result;
result = compute(value);
### 26.2.2 Memory leaks

We say that a program has a *memory leak*, if it allocates memory, and subsequently loses track of that memory. The operating system then thinks the memory is in use, while it is not, and as a result the computer memory can get filled up with allocated memory that serves no useful purpose.

In this example data is allocated inside a lexical scope:

```c
for (i=... ) {
    real *block = malloc( /* large number of bytes */ )
        /* do something with that block of memory */
        /* and forget to call "free" on that block */
}
```

The block of memory is allocated in each iteration, but the allocation of one iteration is no longer available in the next. A similar example can be made with allocating inside a conditional.

It should be noted that this problem is far less serious in Fortran, where memory is deallocated automatically as a variable goes out of scope.

There are various tools for detecting memory errors: Valgrind, DALLOC, Electric Fence. For valgrind, see section 27.3.

### 26.2.3 Roll-your-own malloc

Many programming errors arise from improper use of dynamically allocated memory: the program writes beyond the bounds, or writes to memory that has not been allocated yet, or has already been freed. While some compilers can do bound checking at runtime, this slows down your program. A better strategy is to write your own memory management. Some libraries such as PETSc already supply an enhanced malloc; if this is available you should certainly make use of it.

If you write in C, you will probably know the `malloc` and `free` calls:

```c
int *ip;

ip = (int*) malloc(500*sizeof(int));
if (ip==0) { /* could not allocate memory */ }
..... do stuff with ip ....
free(ip);
```

You can save yourself some typing by

```c
#define MYMALLOC(a,b,c) \
    a = (c*)malloc(b*sizeof(c)); \
    if (a==0) { /* error message and appropriate action */ }

int *ip;

MYMALLOC(ip,500,int);
```
Runtime checks on memory usage (either by compiler-generated bounds checking, or through tools like valgrind or Rational Purify) are expensive, but you can catch many problems by adding some functionality to your malloc. What we will do here is to detect memory corruption after the fact.

We allocate a few integers to the left and right of the allocated object (line 1 in the code below), and put
a recognizable value in them (line 2 and 3), as well as the size of the object (line 2). We then return the
pointer to the actually requested memory area (line 4).

```c
#define MEMCOOKIE 137
#define MYMALLOC(a,b,c) { \
    char *aa; int *ii; \
    aa = malloc(b*sizeof(c)+3*sizeof(int)); /* 1 */ \
    ii = (int*)aa; ii[0] = b*sizeof(c); \
        ii[1] = MEMCOOKIE;                 /* 2 */ \
    aa = (char*) (ii+2); a = (c*)aa;     /* 4 */ \
    aa = aa+b*sizesof(c); ii = (int*)aa; \
        ii[0] = MEMCOOKIE;                /* 3 */ \
}
```

Now you can write your own `free`, which tests whether the bounds of the object have not been written over.

```c
#define MYFREE(a) { \
    char *aa; int *ii,; ii = (int*)a; \
    if (* (--ii) != MEMCOOKIE) printf("object corrupted\n"); \
    n = * (--ii); aa = a+n; ii = (int*)aa; \
    if (*ii != MEMCOOKIE) printf("object corrupted\n"); \
}
```

You can extend this idea: in every allocated object, also store two pointers, so that the allocated memory areas become a doubly linked list. You can then write a macro `CHECKMEMORY` which tests all your allocated objects for corruption.

Such solutions to the memory corruption problem are fairly easy to write, and they carry little overhead. There is a memory overhead of at most 5 integers per object, and there is practically no performance penalty.

(Instead of writing a wrapper for `malloc`, on some systems you can influence the behaviour of the system routine. On linux, `malloc` calls hooks that can be replaced with your own routines; see http://www.gnu.org/s/libc/manual/html_node/Hooks-for-Malloc.html.)

### 26.2.4 Specific techniques: Fortran

Use `Implicit none`.

Put all subprograms in modules so that the compiler can check for missing arguments and type mismatches. It also allows for automatic dependency building with `fdepend`.

Use the C preprocessor for conditional compilation and such.
**26.3 Testing**

There are various philosophies for testing the correctness of a code.

* Correctness proving: the programmer draws up predicates that describe the intended behaviour of code fragments and proves by mathematical techniques that these predicates hold [88, 40].

* Unit testing: each routine is tested separately for correctness. This approach is often hard to do for numerical codes, since with floating point numbers there is essentially an infinity of possible inputs, and it is not easy to decide what would constitute a sufficient set of inputs.

* Integration testing: test subsystems

* System testing: test the whole code. This is often appropriate for numerical codes, since we often have model problems with known solutions, or there are properties such as bounds that need to hold on the global solution.

* Test-driven design: the program development process is driven by the requirement that testing is possible at all times.

With parallel codes we run into a new category of difficulties with testing. Many algorithms, when executed in parallel, will execute operations in a slightly different order, leading to different roundoff behaviour. For instance, the parallel computation of a vector sum will use partial sums. Some algorithms have an inherent damping of numerical errors, for instance stationary iterative methods (section 5.5.1), but others have no such built-in error correction (nonstationary methods; section 5.5.7). As a result, the same iterative process can take different numbers of iterations depending on how many processors are used.

**26.3.1 Test-driven design and development**

In test-driven design there is a strong emphasis on the code always being testable. The basic ideas are as follows.

* Both the whole code and its parts should always be testable.

* When extending the code, make only the smallest change that allows for testing.

* With every change, test before and after.

* Assure correctness before adding new features.
# Chapter 27

## Debugging

When a program misbehaves, *debugging* is the process of finding out *why*. There are various strategies of finding errors in a program. The crudest one is debugging by print statements. If you have a notion of where in your code the error arises, you can edit your code to insert print statements, recompile, rerun, and see if the output gives you any suggestions. There are several problems with this:

* The edit/compile/run cycle is time consuming, especially since

* often the error will be caused by an earlier section of code, requiring you to edit, compile, and rerun repeatedly. Furthermore,

* the amount of data produced by your program can be too large to display and inspect effectively, and

* if your program is parallel, you probably need to print out data from all proccessors, making the inspection process very tedious.

For these reasons, the best way to debug is by the use of an interactive *debugger*, a program that allows you to monitor and control the behaviour of a running program. In this section you will familiarize yourself with `gdb`, which is the open source debugger of the GNU project. Other debuggers are proprietary, and typically come with a compiler suite. Another distinction is that `gdb` is a commandline debugger; there are graphical debuggers such as `ddd` (a frontend to `gdb`) or `DDT` and `TotalView` (debuggers for parallel codes). We limit ourselves to `gdb`, since it incorporates the basic concepts common to all debuggers.

In this tutorial you will debug a number of simple programs with `gdb` and `valgrind`. The files can be downloaded from http://tinyurl.com/ISTC-debug-tutorial.

### 27.1 Invoking gdb

There are three ways of using `gdb`: using it to start a program, attaching it to an already running program, or using it to inspect a *core dump*. We will only consider the first possibility.

Here is an example of how to start `gdb` with a program that has no arguments (Fortran users, use `hello.F`):

tutorials/gdb/c/hello.c
```c
#include <stdlib.h>
#include <stdio.h>
int main() {
    printf("hello world\n");
    return 0;
}

%% cc -g -o hello hello.c
# regular invocation:
%% ./hello
hello world
# invocation from gdb:
%% gdb hello
GNU gdb 6.3.50-20050815 # ..... version info
Copyright 2004 Free Software Foundation, Inc. .... copyright info ....
(gdb) run
Starting program: /home/eijkhout/tutorials/gdb/hello
Reading symbols for shared libraries +. done
hello world

Program exited normally.
(gdb) quit
%%
```

Important note: the program was compiled with the *debug flag -g*. This causes the *symbol table* (that is, the translation from machine address to program variables) and other debug information to be included in the binary. This will make your binary larger than strictly necessary, but it will also make it slower, for instance because the compiler will not perform certain optimizations<fnref n="1" />.

To illustrate the presence of the symbol table do

```bash
%% cc -g -o hello hello.c
%% gdb hello
GNU gdb 6.3.50-20050815 # ..... version info
(gdb) list
```

and compare it with leaving out the -g flag:

```bash
%% cc -o hello hello.c
%% gdb hello
GNU gdb 6.3.50-20050815 # ..... version info
(gdb) list
```

<footnote n="1">1. Compiler optimizations are not supposed to change the semantics of a program, but sometimes do. This can lead to the nightmare scenario where a program crashes or gives incorrect results, but magically works correctly with compiled with debug and run in a debugger.</footnote>
For a program with commandline input we give the arguments to the run command (Fortran users use `say.F`):

tutorials/gdb/c/say.c

```c
#include <stdlib.h>
#include <stdio.h>
int main(int argc, char **argv) {
    int i;
    for (i=0; i<atoi(argv[1]); i++)
        printf("hello world\n");
    return 0;
}

%% cc -o say -g say.c
%% ./say 2
hello world
hello world
%% gdb say
.... the usual messages ...
(gdb) run 2
Starting program: /home/eijkhout/tutorials/gdb/c/say 2
Reading symbols for shared libraries +. done
hello world
hello world

Program exited normally.
```

## 27.2 Finding errors

Let us now consider some programs with errors.

### 27.2.1 C programs

```c
// square.c
    int nmax, i;
    float *squares, sum;

   [f]SCAN独独独独独独独独独独独独
    for (i=1; i<=nmax; i++) {
        squares[i] = 1./((i*i)/2.); sum += squares[i];
    }
    printf("Sum: %e\n", sum);
```
```
% cc -g -o square square.c
% ./square
5000
Segmentation fault
```

The segmentation fault (other messages are possible too) indicates that we are accessing memory that we are not allowed to, making the program abort. A debugger will quickly tell us where this happens:

```
% gdb square
(gdb) run
50000

Program received signal EXC_BAD_ACCESS, Could not access memory.
Reason: KERN_INVALID_ADDRESS at address: 0x000000000000eb4a
0x00007fff824295ca in __svfscanf_1 ()
```

Apparently the error occurred in a function `__svfscanf_1`, which is not one of ours, but a system function. Using the `backtrace` (or `bt`, also `where` or `w`) command we quickly find out how this came to be called:

```
(gdb) backtrace
#0  0x00007fff824295ca in __svfscanf_1 ()
#1  0x00007fff8244011b in fscanf ()
#2  0x0000000100000e89 in main (argc=1, argv=0x7fff5fbfc7c0) at square.c:7
```

We take a close look at line 7, and see that we need to change `nmax` to `&nmax`.

There is still an error in our program:

```
(gdb) run
50000

Program received signal EXC_BAD_ACCESS, Could not access memory.
Reason: KERN_PROTECTION_FAILURE at address: 0x000000010000f000
0x0000000100000ebe in main (argc=2, argv=0x7fff5fbfc7a8) at square1.c:9
9          squares[i] = 1./((i*i)); sum += squares[i];
```

We investigate further:

```
(gdb) print i
$1 = 11237
(gdb) print squares[i]
Cannot access memory at address 0x10000f000
(gdb) print squares
$2 = (float *) 0x0
```

and we quickly see that we forgot to allocate `squares`.

Memory errors can also occur if we have a legitimate array, but we access it outside its bounds.
```c
// up.c
int nlocal = 100,i;
double s, *array = (double*) malloc(nlocal*sizeof(double));
for (i=0; i<nlocal; i++) {
    double di = (double)i;
    array[i] = 1/(di*di);
}
s = 0.;
for (i=nlocal-1; i>=0; i++) {
    double di = (double)i;
    s += array[i];
}
```

```
Program received signal EXC_BAD_ACCESS, Could not access memory.
Reason: KERN_INVALID_ADDRESS at address: 0x0000000100200000
0x0000000100000f43 in main (argc=1, argv=0x7fff5fbfe2c0) at up.c:15
15          s += array[i];
(gdb) print array
$1 = (double *) 0x100104d00
(gdb) print i
$2 = 128608
```

### 27.2.2 Fortran programs

Compile and run the following program:

tutorials/gdb/f/square.F

```fortran
Program square
real squares(1)
integer i

do i=1,100
    squares(i) = sqrt(1.*i)
    sum = sum + squares(i)
end do
print *, "Sum:", sum

End
```

It should abort with a message such as 'Illegal instruction'. Running the program in gdb quickly tells you where the problem lies:

```
(gdb) run
Starting program: tutorials/gdb//fsquare
```
Reading symbols for shared libraries +++++. done

Program received signal EXC_BAD_INSTRUCTION,
Illegal instruction/operand.
0x0000000100000da3 in square () at square.F:7
7          sum = sum + squares(i)

We take a close look at the code and see that we did not allocate squares properly.

## 27.3 Memory debugging with Valgrind

Insert the following allocation of squares in your program:

squares = (float *) malloc(nmax*sizeof(float));

Compile and run your program. The output will likely be correct, although the program is not. Can you see the problem?

To find such subtle memory errors you need a different tool: a memory debugging tool. A popular (because open source) one is `valgrind`; a common commercial tool is `purify`.

tutorials/gdb/c/square1.c

```c
#include <stdlib.h>
#include <stdio.h>
int main(int argc, char **argv) {
    int nmax, i;
    float *squares, sum;

    fprintf(stderr, "%d", &nmax);
    squares = (float*) malloc(nmax*sizeof(float));
    for (i=1; i<=nmax; i++) {
        squares[i] = 1./((i*i));
        sum += squares[i];
    }
    printf("Sum: %e\n", sum);

    return 0;
}
```

Compile this program with `cc -o square1 square1.c` and run it with `valgrind square1` (you need to type the input value). You will lots of output, starting with:

```
% valgrind square1
==53695== Memcheck, a memory error detector
==53695== Copyright (C) 2002-2010, and GNU GPL'd, by Julian Seward et al.
==53695== Using Valgrind-3.6.1 and LibVEX; rerun with -h for copyright info
==53695== Command: a.out
```
```c
==53695==
10
==53695== Invalid write of size 4
==53695==  at 0x100000EB0: main (square1.c:10)
==53695== Address 0x10027e148 is 0 bytes after a block of size 40 alloc'd
==53695==  at 0x1000101EF: malloc (vg_replace_malloc.c:236)
==53695==  by 0x100000E77: main (square1.c:8)
==53695==
==53695== Invalid read of size 4
==53695==  at 0x100000EC1: main (square1.c:11)
==53695== Address 0x10027e148 is 0 bytes after a block of size 40 alloc'd
==53695==  at 0x1000101EF: malloc (vg_replace_malloc.c:236)
==53695==  by 0x100000E77: main (square1.c:8)
```

Valgrind is informative but cryptic, since it works on the bare memory, not on variables. Thus, these error messages take some exegesis. They state that a line 10 writes a 4-byte object immediately after a block of 40 bytes that was allocated. In other words: the code is writing outside the bounds of an allocated array. Do you see what the problem in the code is?

Note that valgrind also reports at the end of the program run how much memory is still in use, meaning not properly freed.

If you fix the array bounds and recompile and rerun the program, valgrind still complains:

```c
==53785== Conditional jump or move depends on uninitialized value(s)
==53785==  at 0x10006FC68: __dtoa (in /usr/lib/libSystem.B.dylib)
==53785==  by 0x10003199F: __vfprintf (in /usr/lib/libSystem.B.dylib)
==53785==  by 0x1000738AA: vfprintf_1 (in /usr/lib/libSystem.B.dylib)
==53785==  by 0x1000A1006: printf (in /usr/lib/libSystem.B.dylib)
==53785==  by 0x100000EF3: main (in ./square2)
```

Although no line number is given, the mention of `printf` gives an indication where the problem lies. The reference to an 'uninitialized value' is again cryptic: the only value being output is `sum`, and that is not uninitialized: it has been added to several times. Do you see why valgrind calls it uninitialized all the same?

## 27.4 Stepping through a program

Often the error in a program is sufficiently obscure that you need to investigate the program run in detail. Compile the following program

```
tutorials/gdb/c/roots.c
```

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

float root(int n)
{
```
```c
float r;
r = sqrt(n);
return r;
}

int main() {
    int i;
    float x=0;
    for (i=100; i>-100; i--)
        x += root(i+5);
    printf("sum: %e\n",x);
    return 0;
}
```

and run it:

```
%% ./roots
sum: nan
```

Start it in gdb as before:

```
%% gdb roots
GNU gdb 6.3.50-20050815
Copyright 2004 Free Software Foundation, Inc.
...
```

but before you run the program, you set a breakpoint at main. This tells the execution to stop, or 'break', in the main program.

```
(gdb) break main
Breakpoint 1 at 0x100000ea6: file root.c, line 14.
```

Now the program will stop at the first executable statement in main:

```
(gdb) run
Starting program: tutorials/gdb/c/roots
Reading symbols for shared libraries +. done
Breakpoint 1, main () at roots.c:14
14      float x=0;
```

If execution is stopped at a breakpoint, you can do various things, such as issuing the step command:

```
Breakpoint 1, main () at roots.c:14
14      float x=0;
(gdb) step
15      for (i=100; i>-100; i--)
(gdb)
```
16
(gdb)

x += root(i);

(if you just hit return, the previously issued command is repeated). Do a number of steps in a row by hitting return. What do you notice about the function and the loop?

Switch from doing step to doing next. Now what do you notice about the loop and the function?

Set another breakpoint: break 17 and do cont. What happens?

Rerun the program after you set a breakpoint on the line with the sqrt call. When the execution stops there do where and list.

## 27.5 Inspecting values

Run the previous program again in gdb: set a breakpoint at the line that does the sqrt call before you actually call run. When the program gets to line 8 you can do print n. Do cont. Where does the program stop?

If you want to repair a variable, you can do `set var=value`. Change the variable n and confirm that the square root of the new value is computed. Which commands do you do?

## 27.6 Breakpoints

If a problem occurs in a loop, it can be tedious keep typing cont and inspecting the variable with print. Instead you can add a condition to an existing breakpoint. First of all, you can make the breakpoint subject to a condition: with

```
condition 1 if (n<0)
```

breakpoint 1 will only obeyed if `n<0` is true.

You can also have a breakpoint that is only activated by some condition. The statement

```
break 8 if (n<0)
```

means that breakpoint 8 becomes (unconditionally) active after the condition `n<0` is encountered.

Another possibility is to use `ignore 1 50`, which will not stop at breakpoint 1 the next 50 times.

Remove the existing breakpoint, redefine it with the condition `n<0` and rerun your program. When the program breaks, find for what value of the loop variable it happened. What is the sequence of commands you use?

You can set a breakpoint in various ways:

* `break foo.c` to stop when code in a certain file is reached;
* `break 123` to stop at a certain line in the current file;

* `break foo` to stop at subprogram `foo`

* or various combinations, such as `break foo.c:123`.

* Finally,

* If you set many breakpoints, you can find out what they are with `info breakpoints`.

* You can remove breakpoints with `delete n` where `n` is the number of the breakpoint.

* If you restart your program with `run` without leaving gdb, the breakpoints stay in effect.

* If you leave gdb, the breakpoints are cleared but you can save them: `save breakpoints <file>`. Use `source <file>` to read them in on the next gdb run.

Finally, you can execute commands at a breakpoint:

```
break 45
command
print x
cont
end
```

This states that at line 45 variable `x` is to be printed, and execution should immediately continue.

If you want to run repeated gdb sessions on the same program, you may want to save an reload breakpoints. This can be done with

```
save-breakpoint filename
source filename
```

## 27.7 Parallel debugging

Debugging in parallel is harder than sequentially, because you will run errors that are only due to interaction of processes such as `deadlock`; see section 2.7.3.6.

As an example, consider this segment of MPI code:

```c
MPI_Init(0,0);
// set comm, ntids, mytid
for (int it=0; ; it++) {
    double randomnumber = ntids * ( rand() / (double)RAND_MAX );
    printf("[%d] iteration %d, random %e\n",mytid,it,randomnumber);
    if (randomnumber>mytid && randomnumber<mytid+1./(ntids+1))
        MPI_Finalize();
}
MPI_Finalize();
```

Each process computes random numbers until a certain condition is satisfied, then exits. However, consider introducing a barrier (or something that acts like it, such as a reduction):
<figure><img src="image_143.png" /><figcaption>Figure 27.1: Display of 16 processes in the DDT debugger</figcaption></figure>

```c
for (int it=0; ; it++) {
    double randomnumber = ntids * ( rand() / (double)RAND_MAX );
    printf("[%d] iteration %d, random %e\n", mytid, it, randomnumber);
    if (randomnumber>mytid && randomnumber<mytid+1./(ntids+1))
        MPI_Finalize();
    MPI_BARRIER comm;
}
MPI_Finalize();
```

Now the execution will hang, and this is not due to any particular process: each process has a code path from init to finalize that does not develop any memory errors or other runtime errors. However as soon as one process reaches the finalize call in the conditional it will stop, and all other processes will be waiting at the barrier.

Figure 27.1 shows the main display of the Allinea DDT debugger (http://www.allinea.com/products/ddt) at the point where this code stops. Above the source panel you see that there are 16 processes, and that the status is given for process 1. In the bottom display you see that out of 16 processes
15 are calling `MPI Barriers` on line 19, while one is at line 18. In the right display you see a listing of the local variables: the value specific to process 1. A rudimentary graph displays the values over the processors: the value of `ntids` is constant, that of `mytid` is linearly increasing, and `it` is constant except for one process.

Exercise 27.1. Make and run `ring_la`. The program does not terminate and does not crash. In the debugger you can interrupt the execution, and see that all processes are executing a receive statement. This is probably a case of deadlock. Diagnose and fix the error.

Exercise 27.2. The author of `ring_lc` was very confused about how MPI works. Run the program. While it terminates without a problem, the output is wrong. Set a breakpoint at the send and receive statements to figure out what is happening.

## 27.8 Further reading

A good tutorial: http://www.dirac.org/linux/gdb/.

Reference manual: http://www.ofb.net/gnu/gdb/gdb_toc.html.
# Chapter 28

## C/Fortran interoperability

Most of the time, a program is written in a single language, but in some circumstances it is necessary or desirable to mix sources in more than one language for a single executable. One such case is when a library is written in one language, but used by a program in another. In such a case, the library writer will probably have made it easy for you to use the library; this section is for the case that you find yourself in the place of the library writer. We will focus on the common case of interoperability between C/C++ and Fortran.

This issue is complicated by the fact that both languages have been around for a long time, and various recent language standards have introduced mechanisms to facilitate interoperability. However, there is still a lot of old code around, and not all compilers support the latest standards. Therefore, we discuss both the old and the new solutions.

### 28.1 Linker conventions

As explained above, a compiler turns a source file into a binary, which no longer has any trace of the source language: it contains in effect functions in machine language. The linker will then match up calls and definitions, which can be in different files. The problem with using multiple languages is then that compilers have different notions of how to translate function names from the source file to the binary file.

Let's look at codes (you can find example files in tutorials/linking):

```c
// C:
    Subroutine foo()
    Return
    End Subroutine

! Fortran
void foo() {
    return;
}
```

After compilation you can use *nm* to investigate the binary *object file*:
```fortran
%% nm fprog.o
0000000000000000 T _foo_

%% nm cprog.o
0000000000000000 T _foo

```

You see that internally the `foo` routine has different names: the Fortran name has an underscore appended. This makes it hard to call a Fortran routine from C, or vice versa. The possible name mismatches are:

* The Fortran compiler appends an underscore. This is the most common case.

* Sometimes it can append two underscores.

* Typically the routine name is lowercase in the object file, but uppercase is a possibility too.

Since C is a popular language to write libraries in, this means that the problem is often solved in the C library by:

* Appending an underscore to all C function names; or

* Including a simple wrapper call:

```c
int SomeCFunction(int i, float f)
{
    ...
}

int SomeCFunction_(int i, float f)
{
    return SomeCFunction(i, f);
}
```

### 28.1.1 C bindings in Fortran 2003

With the latest Fortran standard there are explicit C bindings, making it possible to declare the external name of variables and routines:

```fortran
module operator
    real, bind(C) :: x
contains
    subroutine s() bind(C,name='s')
    return
end subroutine
end module
```

```fortran
%% ifort -c fbind.F90
%% nm fbind.o
.... T _s
.... C _x
```
It is also possible to declare data types to be C-compatible:

```fortran
Program fdata

    use iso_c_bind

type, bind(C) :: c_comp
    real (c_float)  :: data
    integer (c_int) :: i
    type (c_ptr)    :: ptr
end type

end Program fdata
```

### 28.1.2 C++ linking

Libraries written in C++ offer further problems. The C++ compiler makes external symbols by combining the names a class and its methods, in a process known as *name mangling*. You can force the compiler to generate names that are intelligible to other languages by

```c
#ifdef __cplusplus
extern "C" {
#endif

.

place declarations here

.

#ifdef __cplusplus
}
#endif
```

Example: compiling

```c
#include <stdlib.h>

int foo(int x) {
    return x;
}
```

and inspecting the output with `nm` gives:

```
000000000000010 s EH_framel
0000000000000000 T _foo
```
On the other hand, the identical program compiled as C++ gives

```
0000000000000010 s EH_frame1
0000000000000000 T __Z3fooi
```

You see that the name for `foo` is something mangled, so you can not call this routine from a program in a different language. On the other hand, if you add the `extern` declaration:

```c
#include <stdlib.h>

#ifdef __cplusplus
extern "C" {
#endif
int foo(int x) {
    return x;
}

#ifdef __cplusplus
}
#endif
```

you again get the same linker symbols as for C, so that the routine can be called from both C and Fortran.

If your main program is in C, you can use the C++ compiler as linker. If the main program is in Fortran, you need to use the Fortran compiler as linker. It is then necessary to link in extra libraries for the C++ system routines. For instance, with the Intel compiler `-lstdc++ -lc` needs to be added to the link line.

## 28.2 Arrays

C and Fortran have different conventions for storing multi-dimensional arrays. You need to be aware of this when you pass an array between routines written in different languages.

Fortran stores multi-dimensional arrays in column-major order; see figure 28.1. For two dimensional arrays A(i,j) this means that the elements in each column are stored contiguously: a 2 × 2 array is stored as A(1,1), A(2,1), A(1,2), A(2,2). Three and higher dimensional arrays are an obvious extension: it is sometimes said that ‘the left index varies quickest’.

C arrays are stored in row-major order: elements in each row are stored contiguous, and columns are then placed sequentially in memory. A 2 × 2 array A[2][2] is stored as A[1][1], A[1][2], A[2][1], A[2][2].

A number of remarks about arrays in C.

* C (before the C99 standard) has multi-dimensional arrays only in a limited sense. You can declare them, but if you pass them to another C function, they no longer look multi-dimensional: they have become plain `float*` (or whatever type) arrays. That brings us to the next point.
<figure><img src="image_144.png" /><figcaption>Figure 28.1: Fortran and C array storage by columns and rows respectively</figcaption></figure>

* Multi-dimensional arrays in C look as if they have type `float**`, that is, an array of pointers that point to (separately allocated) arrays for the rows. While you could certainly implement this:

```c
float **A;
A = (float**)malloc(m*sizeof(float*));
for (i=0; i<n; i++)
    A[i] = (float*)malloc(n*sizeof(float));
```

careful reading of the standard reveals that a multi-dimensional array is in fact a single block of memory, no further pointers involved.

Given the above limitation on passing multi-dimensional arrays, and the fact that a C routine can not tell whether it's called from Fortran or C, it is best not to bother with multi-dimensional arrays in C, and to emulate them:

```c
float *A;
A = (float*)malloc(m*n*sizeof(float));
#define SUB(i,j,m,n) i+j*m
for (i=0; i<m; i++)
    for (j=0; j<n; j++)
        A SUB(i,j,m,n)]
```

where for interoperability we store the elements in column-major fashion.

## 28.3 Strings

Programming languages differ widely in how they handle strings.
* In C, a string is an array of characters; the end of the string is indicated by a null character, that is the ascii character zero, which has an all zero bit pattern. This is called null termination.

* In Fortran, a string is an array of characters. The length is maintained in a internal variable, which is passed as a hidden parameter to subroutines.

* In Pascal, a string is an array with an integer denoting the length in the first position. Since only one byte is used for this, strings can not be longer than 255 characters in Pascal.

As you can see, passing strings between different languages is fraught with peril. This situation is made even worse by the fact that passing strings as subroutine arguments is not standard.

Example: the main program in Fortran passes a string

```pascal
Program Fstring
    character(len=5) :: word = "Word"
    call cstring(word)
end Program Fstring
```

and the C routine accepts a character string and its length:

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

void cstring_(char *txt, int txtlen) {
    printf("length = %d\n", txtlen);
    printf("<<");
    for (int i=0; i<txtlen; i++)
        printf("%c", txt[i]);
    printf(">>\n");
}
```

which produces:

```
length = 5
<<Word >>
```

To pass a Fortran string to a C program you need to append a null character:

```c
call cfunction ('A string'//CHAR(0))
```

Some compilers support extensions to facilitate this, for instance writing

```c
DATA forstring /'This is a null-terminated string.'C/
```

Recently, the 'C/Fortran interoperability standard' has provided a systematic solution to this.
## 28.4 Subprogram arguments

In C, you pass a `float` argument to a function if the function needs its value, and `float*` if the function has to modify the value of the variable in the calling environment. Fortran has no such distinction: every variable is passed by reference. This has some strange consequences: if you pass a literal value 37 to a subroutine, the compiler will allocate a nameless variable with that value, and pass the address of it, rather than the value<fnref n="1" />.

For interfacing Fortran and C routines, this means that a Fortran routine looks to a C program like all its argument are ‘star’ arguments. Conversely, if you want a C subprogram to be callable from Fortran, all its arguments have to be star-this or that. This means on the one hand that you will sometimes pass a variable by reference that you would like to pass by value.

Worse, it means that C subprograms like

```c
void mysub(int **iarray) {
    *iarray = (int*)malloc(8*sizeof(int));
    return;
}
```

can not be called from Fortran. There is a hack to get around this (check out the Fortran77 interface to the `Petsc` routine `VecGetValues`) and with more cleverness you can use `POINTER` variables for this.

## 28.5 Input/output

Both languages have their own system for handling input/output, and it is not really possible to meet in the middle. Basically, if Fortran routines do I/O, the main program has to be in Fortran. Consequently, it is best to isolate I/O as much as possible, and use C for I/O in mixed language programming.

## 28.6 Fortran/C interoperability in Fortran2003

The latest version of Fortran, unsupported by many compilers at this time, has mechanisms for interfacing to C.

* There is a module that contains named kinds, so that one can declare `INTEGER, KIND(C_SHORT) :: i`

* Fortran pointers are more complicated objects, so passing them to C is hard; Fortran2003 has a mechanism to deal with C pointers, which are just addresses.

* Fortran derived types can be made compatible with C structures.

<footnote n="1">1. With a bit of cleverness and the right compiler, you can have a program that says `print *, 7` and prints `8` because of this.</footnote>
# Chapter 29

## LaTeX for scientific documentation

### 29.1 The idea behind LATEX, some history of TEX

TEX is a typesetting system that dates back to the late 1970s. In those days, graphics terminals where you could design a document layout and immediately view it, the way you can with for instance Microsoft Word, were rare. Instead, TEX uses a two-step workflow, where you first type in your document with formatting instructions in an ascii document, using your favourite text editor. Next, you would invoke the `latex` program, as a sort of compiler, to translate this document to a form that can be printed or viewed.

```
%% edit mydocument.tex
%% latex mydocument
%% # print or view the resulting output
```

The process is comparable to making web pages by typing HTML commands.

This way of working may seem clumsy, but it has some advantages. For instance, the `TEX` input files are plain ascii, so they can easily be generated automatically, for instance from a database. Also, you can edit them with whatever your favourite editor happens to be.

Another point in favour of `TEX` is the fact that the layout is specified by commands that are written in a sort of programming language. This has some important consequences:

* Separation of concerns: when you are writing your document, you do not have to think about layout. You give the 'chapter' command, and the implementation of that command will be decided independently, for instance by you choosing a document style.

* Changing the layout of a finished document is easily done by choosing a different realization of the layout commands in the input file: the same 'chapter' command is used, but by choosing a different style the resulting layout is different. This sort of change can be as simple as a one-line change to the document style declaration.

* If you have unusual typesetting needs, it is possible to write new `TEX` commands for this. For many needs such extensions have in fact already been written; see section 29.4.

The commands in `TEX` are fairly low level. For this reason, a number of people have written systems on top of `TEX` that offer powerful features, such as automatic cross-referencing, or generation of a table of contents. The most popular of these systems is `LATEX`. Since `TEX` is an interpreted system, all of its mechanisms are still available to the user, even though `LATEX` is loaded on top of it.
### 29.1.1 Installing LaTeX

The easiest way to install LaTeX on your system is by downloading the TeXlive distribution from http://tug.org/texlive. Apple users can also use fink or macports. Various front-ends to TeX exist, such as TeXshop on the Mac.

### 29.1.2 Running LaTeX

**Purpose.** In this section you will run the LaTeX compiler

Originally, the `latex` compiler would output a device independent file format, named `dvi`, which could then be translated to PostScript or PDF, or directly printed. These days, many people use the `pdflatex` program which directly translates `.tex` files to `.pdf` files. This has the big advantage that the generated PDF files have automatic cross linking and a side panel with table of contents. An illustration is found below.

Let us do a simple example.

``` Latex
\documentclass{article}
\begin{document}
Hello world!
\end{document}
```

<caption>Figure 29.1: A minimal LaTeX document</caption>

**Exercise.** Create a text file `minimal.tex` with the content as in figure 29.1. Try the command `pdflatex minimal` or `latex minimal`. Did you get a file `minimal.pdf` in the first case or `minimal.dvi` in the second case? Use a pdf viewer, such as Adobe Reader, or `dvips` respectively to view the output.

*Caveats.* If you make a typo, TeX can be somewhat unfriendly. If you get an error message and TeX is asking for input, typing `x` usually gets you out, or `Ctrl-C`. Some systems allow you to type `e` to go directly into the editor to correct the typo.

## 29.2 A gentle introduction to LaTeX

Here you will get a very brief run-through of LaTeX features. There are various more in-depth tutorials available, such as the one by Oetiker [128].

### 29.2.1 Document structure

Each LaTeX document needs the following lines:
\documentclass{....} % the dots will be replaced
\begin{document}
\end{document}

The 'documentclass' line needs a class name in between the braces; typical values are 'article' or 'book'.
Some organizations have their own styles, for instance 'ieeeproc' is for proceedings of the IEEE.

All document text goes between the \begin{document} and \end{document} lines. (Matched 'begin' and 'end' lines are said to denote an 'environment', in this case the document environment.)

The part before \begin{document} is called the 'preamble'. It contains customizations for this particular document. For instance, a command to make the whole document double spaced would go in the preamble. If you are using pdflatex to format your document, you want a line

\usepackage{hyperref}

here.

Have you noticed the following?

* The backslash character is special: it starts a LaTeX command.

* The braces are also special: they have various functions, such as indicating the argument of a command.

* The percent character indicates that everything to the end of the line is a comment.

### 29.2.2 Some simple text

**Purpose.** In this section you will learn some basics of text formatting.

**Exercise.** Create a file first.tex with the content of figure 29.1 in it. Type some text in the preamble, that is, before the \begin{document} line and run pdflatex on your file.

*Expected outcome.* You should get an error message because you are not allowed to have text in the preamble. Only commands are allowed there; all text has to go after \begin{document}.

**Exercise.** Edit your document: put some text in between the \begin{document} and \end{document} lines. Let your text have both some long lines that go on for a while, and some short ones. Put superfluous spaces between words, and at the beginning or end of lines. Run pdflatex on your document and view the output.

*Expected outcome.* You notice that the white space in your input has been collapsed in the output. TeX has its own notions about what space should look like, and you do not have to concern yourself with this matter.

**Exercise.** Edit your document again, cutting and pasting the paragraph, but leaving a blank line between the two copies. Paste it a third time, leaving several blank lines. Format, and view the output.
*Expected outcome*. TeX interprets one or more blank lines as the separation between paragraphs.

**Exercise.** Add \usepackage{pslatex} to the preamble and rerun pdflatex on your document. What changed in the output?

*Expected outcome*. This should have the effect of changing the typeface from the default to Times Roman.

*Caveats*. Typefaces are notoriously unstandardized. Attempts to use different typefaces may or may not work. Little can be said about this in general.

Add the following line before the first paragraph:

\section{This is a section}

and a similar line before the second. Format. You see that LaTeX automatically numbers the sections, and that it handles indentation different for the first paragraph after a heading.

**Exercise.** Replace *article* by *artikel3* in the documentclass declaration line and reformat your document. What changed?

*Expected outcome*. There are many documentclasses that implement the same commands as *article* (or another standard style), but that have their own layout. Your document should format without any problem, but get a better looking layout.

*Caveats*. The *artikel3* class is part of most distributions these days, but you can get an error message about an unknown documentclass if it is missing or if your environment is not set up correctly. This depends on your installation. If the file seems missing, download the files from http://tug.org/texmf-dist/tex/latex/ntgclass/ and put them in your current directory; see also section 29.2.9.

### 29.2.3 Math

**Purpose.** In this section you will learn the basics of math typesetting

One of the goals of the original TeX system was to facilitate the setting of mathematics. There are two ways to have math in your document:

* Inline math is part of a paragraph, and is delimited by dollar signs.

* Display math is, as the name implies, displayed by itself.

**Exercise.** Put $x+y$ somewhere in a paragraph and format your document. Put \[x+y\] somewhere in a paragraph and format.

*Expected outcome*. Formulas between single dollars are included in the paragraph where you declare them. Formulas between \[...\] are typeset in a display.

For display equations with a number, use an *equation* environment. Try this.

Here are some common things to do in math. Make sure to try them out.
* Subscripts and superscripts: $x_i^{2}$. If the sub or superscript is more than a single symbol, it needs to be grouped: $x_{i+1}^{2}$. If you need a brace in a formula, use $\{ \}$.

* Greek letters and other symbols: $\alpha\otimes\beta$.

* Combinations of all these: $\int_{t=0}^{\infty}$.

**Exercise.** Take the last example and typeset it as display math. Do you see a difference with inline math?

*Expected outcome.* TeX tries not to include the distance between text lines, even if there is math in a paragraph. For this reason it typesets the bounds on an integral sign differently from display math.

### 29.2.4 Referencing

**Purpose.** In this section you will see TeX's cross referencing mechanism in action.

So far you have not seen LATEX do much that would save you any work. The cross referencing mechanism of LATEX will definitely save you work: any counter that LATEX inserts (such as section numbers) can be referenced by a label. As a result, the reference will always be correct.

Start with an example document that has at least two section headings. After your first section heading, put the command \label{sec:first}, and put \label{sec:other} after the second section heading. These label commands can go on the same line as the section command, or on the next. Now put

As we will see in section~\ref{sec:other}.

in the paragraph before the second section. (The tilde character denotes a non-breaking space.)

**Exercise.** Make these edits and format the document. Do you see the warning about an undefined reference? Take a look at the output file. Format the document again, and check the output again. Do you have any new files in your directory?

*Expected outcome.* On a first pass through a document, the TeX compiler will gather all labels with their values in a .aux file. The document will display a double question mark for any references that are unknown. In the second pass the correct values will be filled in.

*Caveats.* If after the second pass there are still undefined references, you probably made a typo. If you use the bibtex utility for literature references, you will regularly need three passes to get all references resolved correctly.

Above you saw that the *equation* environment gives displayed math with an equation number. You can add a label to this environment to refer to the equation number.

**Exercise.** Write a formula in an *equation* environment, and add a label. Refer to this label anywhere in the text. Format (twice) and check the output.

*Expected outcome.* The \label and \ref command are used in the same way for formulas as for section numbers. Note that you must use \begin/end{equation} rather than \[...\] for the formula.
### 29.2.5 Lists

**Purpose.** In this section you will see the basics of lists.

Bulleted and numbered lists are provided through an environment.

```
\begin{itemize}
\item This is an item;
\item this is one too.
\end{itemize}
\begin{enumerate}
\item This item is numbered;
\item this one is two.
\end{enumerate}
```

**Exercise.** Add some lists to your document, including nested lists. Inspect the output.

*Expected outcome.* Nested lists will be indented further and the labeling and numbering style changes with the list depth.

**Exercise.** Add a label to an item in an *enumerate* list and refer to it.

*Expected outcome.* Again, the `\label` and `\ref` commands work as before.

### 29.2.6 Source code and algorithms

As a computer scientist, you will often want to include algorithms in your writings; sometimes even source code.

In this tutorial so far you have seen that some characters have special meaning to LaTeX, and just can not just type them and expect them to show up in the output. Since funny characters appear quite regularly in programming languages, we need a tool for this: the *verbatim mode*.

To display bits of code inside a paragraph, you use the `\verb` command. This command delimits its argument with two identical characters that can not appear in the verbatim text. For instance, the output `if (x%5>0) { ... }` is produced by `\verb+if (x%5>0) { ... } +`. (Exercise: how did the author of this book get that verbatim command in the text?)

For longer stretches of verbatim text, that need to be displayed by themselves you use

```
\begin{verbatim}
stuff
\end{verbatim}
```

Finally, in order to include a whole file as verbatim listing, use `.`.

Verbatim text is one way of displaying algorithms, but there are more elegant solutions. For instance, in this book the following is used:

```
\usepackage[algo2e, noline, noend]{algorithm2e}
```

The result can be seen, for instance, on page 271.
### 29.2.7 Graphics

Since you can not immediately see the output of what you are typing, sometimes the output may come as a surprise. That is especially so with graphics. LaTeX has no standard way of dealing with graphics, but the following is a common set of commands:

```
\usepackage{graphicx} % this line in the preamble
```

```
\includegraphics{myfigure} % in the body of the document
```

The figure can be in any of a number of formats, except that PostScript figures (with extension `.ps` or `.eps`) can not be used if you use `pdflatex`.

Since your figure is often not the right size, the `include` line will usually have something like:

```
\includegraphics[scale=.5]{myfigure}
```

A bigger problem is that figures can be too big to fit on the page if they are placed where you declare them. For this reason, they are usually treated as 'floating material'. Here is a typical declaration of a figure:

```
\begin{figure}[ht]
\includegraphics{myfigure}
\caption{This is a figure}
\label{fig:first}
\end{figure}
```

It contains the following elements:

* The `figure` environment is for 'floating' figures; they can be placed right at the location where they are declared, at the top or bottom of the next page, at the end of the chapter, et cetera.

* The `[ht]` argument of the `\begin{figure}` line states that your figure should be attempted to be placed here; it that does not work, it should go top of the next page. The remaining possible specifications are `b` for placement at the bottom of a page, or `p` for placement on a page by itself. For example

```
\begin{figure}[hbp]
```

declares that the figure has to be placed here if possible, at the bottom of the page if that's not possible, and on a page of its own if it is too big to fit on a page with text.

* A caption to be put under the figure, including a figure number;

* A label so that you can refer to the figure number by its label: `figure~\ref{fig:first}`.

* And of course the figure material. There are various ways to fine-tune the figure placement. For instance

```
\begin{center}
\includegraphics{myfigure}
\end{center}
```

gives a centered figure.
### 29.2.8 Bibliography references

The mechanism for citing papers and books in your document is a bit like that for cross referencing. There are labels involved, and there is a `\cite{thatbook}` command that inserts a reference, usually numeric. However, since you are likely to refer to a paper or book in more than one document you write, LaTeX allows you to have a database of literature references in a file by itself, rather than somewhere in your document.

Make a file `mybibliography.bib` with the following content:

```bibtex
@article{JoeDoe1985,
author = {Joe Doe},
title = {A framework for bibliography references},
journal = {American Library Assoc. Mag.},
year = {1985}
}
```

In your document `mydocument.tex`, put

```bibtex
For details, refer to Doe\cite{JoeDoe1985} % somewhere in the text
```

```bibtex
\bibliography{mybibliography} % at the end of the document
```

```bibtex
\bibliographystyle{plain}
```

Format your document, then type on the commandline

```bash
bibtex mydocument
```

and format your document two more times. There should now be a bibliography in it, and a correct citation. You will also see that files `mydocument.bbl` and `mydocument.blg` have been created.

### 29.2.9 Environment variables

On Unix systems, `TEX` investigates the `TEXINPUTS` environment variable when it tries to find an include file. Consequently, you can create a directory for your styles and other downloaded include files, and set this variable to the location of that directory. Similarly, the `BIBINPUTS` variable indicates the location of bibliography files for `bibtex` (section 29.2.8).

## 29.3 A worked out example

The following example `demo.tex` contains many of the elements discussed above.

```tex
\documentclass{article}
```

```tex
\usepackage{pslatex, graphicx, amsmath, amssymb}
```

```tex
\usepackage{pdflatex}
```

```tex
\newtheorem{theorem}{Theorem}
```
\newcounter{excounter}
\newenvironment{exercise}
{
  \refstepcounter{excounter}
  \begin{quotation}
    \textbf{Exercise \arabic{excounter}.}
  \end{quotation}
}{
  \end{quotation}
}

\begin{document}
\title{SSC 335: demo}
\author{Victor Eijkhout}
\date{today}
\maketitle

\section{This is a section}
\label{sec:intro}

This is a test document, used in \cite{latexdemo}. It contains a discussion in section \ref{sec:discussion}.

\begin{exercise}\label{easy-ex}
Left to the reader.
\end{exercise}
\begin{exercise}
Also left to the reader, just like in exercise \ref{easy-ex}
\end{exercise}

\begin{theorem}
This is cool.
\end{theorem}
\begin{equation*}
\begin{split}
& \text{This is a formula: } \$a\text{Leftarrow } b\$. \\
& \text{\label{eq:one}}
\end{split}
\end{equation*}
\begin{equation*}
x_i\text{leftarrow} y_{\{ij\}}\cdot\text{箭} x^{\{(k)\}}_j
\end{equation*}
\text{Text: } \$\int_0^1 \sqrt{x}\,dx\$
\begin{center}
\text{\int_0^1 \sqrt{x}\,,dx}
\end{center}
\begin{center}
\text{This is another section}
\label{sec:discussion}
\end{center}

\begin{table}[ht]
\centering
\begin{tabular}{|l|}
\hline
one&value \\
\hline
another&values \\
\hline
\end{tabular}
\caption{This is the only table in my demo}
\label{tab:thetable}
\end{table}
\end{document}
\begin{figure}[ht]
\centering
\includegraphics{graphics/caches}
\caption{this is the only figure}
\label{fig:thefigure}
\end{figure}
As I showed in the introductory section~\ref{sec:intro}, in the
paper~\cite{AdJo:colorblind}, it was shown that
equation~\eqref{eq:one}
\begin{center}
\item There is an item.
\item There is another item
\begin{itemize}
\item sub one
\item sub two
\end{itemize}
\end{center}
\end{itemize}
\begin{enumerate}
\item item one
\item item two
\begin{enumerate}
\item sub one
\item sub two
\end{enumerate}
\end{enumerate}
\end{listofcontents}
\bibliography{math}
\bibliographystyle{plain}
\end{document}

You also need the file math.bib:

@article{AdJo:colorblind,
author = {Loyce M. Adams and Harry F. Jordan},
title = {Is {SOR} color-blind?},
journal = {SIAM J. Sci. Stat. Comput.},
year = {1986},
volume = {7},
pages = {490--506},
abstract = {For what stencils do ordinary and multi-colour SOR have
the same eigenvalues.},
keywords = {SOR, colouring}
}

@misc{latexdemo,
author = {Victor Eijkhout},
title = {Short {\LaTeX}\ demo},
note = {SSC 335, oct 1, 2008}
}

The following sequence of commands

pdflatex demo
bibtex demo
pdflatex demo
pdflatex demo

gives
# SSC 335: demo

Victor Eijkhout

today

## 1 This is a section

This is a test document, used in [2]. It contains a discussion in section 2.

**Exercise 1.** Left to the reader.

**Exercise 2.** Also left to the reader, just like in exercise 1

**Theorem 1** *This is cool.*

This is a formula: $a \leftarrow b$.

Text: $\int_{0}^{1} \sqrt{x} dx$

$$x_i \leftarrow y_{ij} \cdot x_j^{(k)} \quad (1)$$

$$\int_{0}^{1} \sqrt{x} dx$$

## 2 This is another section

<table><caption>Table 1: This is the only table in my demo</caption><tr><td>one</td><td>value</td></tr><tr><td>another</td><td>values</td></tr></table>

<figure><img src="image_145.png" /><figcaption>Figure 1: this is the only figure</figcaption></figure>

As I showed in the introductory section 1, in the paper [1], it was shown that equation (1)

* There is an item.
* There is another item

- sub one

- sub two

1. item one

2. item two

(a) sub one

(b) sub two

## Contents

1. This is a section 1

2. This is another section 1

## List of Figures

1. this is the only figure 1

## References

[1] Loyce M. Adams and Harry F. Jordan. Is SOR color-blind? *SIAM J. Sci. Stat. Comput.*, 7:490–506, 1986.

[2] Victor Eijkhout. Short LATEX demo. SSC 335, oct 1, 2008.
29.4 Where to take it from here

This tutorial touched only briefly on some essentials of T E X and L A T E X. You can find longer intros on-
line [128], or read a book [106, 100, 126]. Macro packages and other software can be found on the
Comprehensive T E X Archive http://www.ctan.org. For questions you can go to the newsgroup
comp.text.tex, but the most common ones can often already be found on web sites [152].
# PART V

## PROJECTS, CODES, INDEX
# Chapter 30

## Class projects

### 30.1 Cache simulation and analysis

In this project you will build a cache simulator and analyze the cache hit/miss behaviour of code, either real or simulated.

#### 30.1.1 Cache simulation

A simulated cache is a simple data structure that records for each cache address what memory address it contains, and how long the data has been present. Design this data structure and write the access routines. (Use of an object oriented language is recommended.)

Write your code so that your cache can have various levels of associativity, and different replacement policies.

For simplicity, do not distinguish between read and write access. Therefore, the execution of a program becomes a stream of

```
cache.access_address( 123456 );
cache.access_address( 70543 );
cache.access_address( 12338383 );
.....
```

calls where the argument is the memory address. Your code will record whether the request can be satisfied from cache, or whether the data needs to be loaded from memory.

#### 30.1.2 Code simulation

Find some example codes, for instance from a scientific project you are involved in, and translate the code to a sequence of memory address.

You can also simulate codes by generating a stream of access instructions as above:

* Some access will be to random locations, corresponding to use of scalar variables;

* At other times access will be to regularly space addresses, corresponding to the use of an array;
* Array operations with indirect addressing will cause prolonged episodes of random address access.

Try out various mixes of the above instruction types. For the array operations, try smaller and larger arrays with various degrees of reuse.

Does your simulated code behave like a real code?

**30.1.3 Investigation**

First implement a single cache level and investigate the behaviour of cache hits and misses. Explore different associativity amounts and different replacement policies.

**30.1.4 Analysis**

Do a statistical analysis of the cache hit/miss behaviour. You can start with [140]<fnref n="1" />. Hartstein [80] found a power law behaviour. Are you finding the same?

<footnote n="1">1. Strictly speaking that paper is about page swapping out of virtual memory (section 1.3.8), but everything translates to cacheline swapping out of main memory.</footnote>
**30.2 Evaluation of Bulk Synchronous Programming**

In this project you're asked to analyze the Bulk Synchronous Parallel (BSP) model, both abstractly and by simulation.

**30.2.1 Discussion**

For your assignment you need to investigate Bulk Synchronous Programming. Read the wikipedia article http://en.wikipedia.org/wiki/Bulk_synchronous_parallel and section 2.7.8 in the book.

Consider the following questions.

1. Discuss the relation between the complexity model of BSP and the $\alpha, \beta, \gamma$ model you learned in section 6.1.

2. BSP uses barrier synchronisations. Are these necessary in practice? Consider separately two different classes of algorithms, namely PDE solving, and large scale graph computations.

**30.2.2 Simulation**

In the defining paper on BSP, Valiant [154] advocates making more tasks than processors, and distributing them randomly. Program a simulation and test if this strategy solves the load balance problem. Discuss the random placement strategy in the context of various algorithms.
## 30.3 Heat equation

In this project you will combine some of the theory and practical skills you learned in class to solve a real world problem. In addition to writing code that solves the program, you will use the following software practices:

* use source code control,

* give the program checkpoint/restart capability, and

* use basic profiling and visualization of the results.

The Heat Equation (see section 4.3) is given by

$$ \frac{\partial T(x, t)}{\partial t} = \begin{cases} \alpha \frac{\partial^2 T(x,y)}{\partial x^2} + q(x, t) & \text{1D} \\ \alpha \frac{\partial^2 T(x,y)}{\partial x^2} + \alpha \frac{\partial^2 T(x,y)}{\partial y^2} + q(x, t) & \text{2D} \\ \dots & \text{3D} \end{cases} $$

where $t \ge 0$ and $x \in [0, 1]$, subject to boundary conditions

$$ T(x, 0) = T_0(x), \quad \text{for } x \in [0, 1], $$

and similar in higher dimensions, and

$$ T(0, t) = T_a(t), \quad T(1, t) = T_b(t), \quad \text{for } t > 0. $$

You will solve this problem using the explicit and implicit Euler methods.

### 30.3.1 Software

Write your software using the PETSc library (see tutorial 24. In particular, use the MatMult routine for matrix-vector multiplication and KSPSolve for linear system solution. Exception: code the Euler methods yourself.

Be sure to use a Makefile for building your project (tutorial 21).

Add your source files, Makefile, and job scripts to an svn repository (tutorial 22); do not add binaries or output files. Make sure that there is a README file with instructions on how to build and run your code.

Implement a checkpoint/restart facility by writing vector data, size of the time step, and other necessary items, to an hdf5 file (tutorial 23). Your program should be able to read this file and resume execution.

### 30.3.2 Tests

Do the following tests on a single core.
Method stability

Run your program for the 1D case with

$$
\begin{cases}
q = \sin \ell \pi x \\
T_0(x) = e^x \\
T_a(t) = T_b(t) = 0 \\
\alpha = 1
\end{cases}
$$

Take a space discretization at least $h = 10^{-2}$ but especially in parallel do not be afraid to try large problem sizes. Try various time steps and show that the explicit method can diverge. What is the maximum time step for which it is stable?

For the implicit method, at first use a direct method to solve the system. This corresponds to PETSc options KSPPREONLY and PCLU (see section 5.5.10 and the PETSc tutorial, 24).

Now use an iterative method (for instance `KSPCG` and `PCJACOBI`); is the method still stable? Explore using a low convergence tolerance and large time steps.

Since the forcing function $q$ and the boundary conditions have no time dependence, the solution $u(\cdot, t)$ will converge to a steady state solution $u_\infty(x)$ as $t \to \infty$. What is the influence of the time step on the speed with which implicit method converges to this steady state?

Hint: the steady state is described by $u_t \equiv 0$. Substitute this in the PDE. Can you find explicitly what the steady state is?

Run these tests with various values for $\ell$.

Timing

If you run your code with the commandline option `-log_summary`, you will get a table of timings of the various PETSc routines. Use that to do the following timing experiments. Make sure you use a version of PETSc that was *not* compiled with debug mode.

Construct your coefficient matrix as a dense matrix, rather than sparse. Report on the difference in total memory, and the runtime and flop rate of doing one time step. Do this for both the explicit and implicit method and explain the results.

With a sparse coefficient matrix, report on the timing of a single time step. Discuss the respective flop counts and the resulting performance.

Restart

Implement a restart facility: every 10 iterations write out the values of the current iterate, together with values of $\Delta x$, $\Delta t$, and $\ell$. Add a flag `-restart` to your program that causes it to read the restart file and resume execution, reading all parameters from the restart file.

Run your program for 25 iterations, and restart, causing it to run again from iteration 20. Check that the values in iterations 20...25 match.
### 30.3.3 Parallelism

Do the following tests to determine the parallel scaling of your code.

At first test the explicit method, which should be perfectly parallel. Report on actual speedup attained. Try larger and smaller problem sizes and report on the influence of the problem size.

The above settings for the implicit method (`KSPPREONLY` and `PCLU`) lead to a runtime error. One way out is to let the system be solved by an iterative method. Read the PETSc manual and web pages to find out some choices of iterative method and preconditioner and try them. Report on their efficacy.

### 30.3.4 Comparison of solvers

In the implicit timestepping method you need to solve linear systems. In the 2D (or 3D) case it can make a big difference whether you use a direct solver or an iterative one.

Set up your code to run with $q \equiv 0$ and zero boundary conditions everywhere; start with a nonzero initial solution. You should now get convergence to a zero steady state, thus the norm of the current solution is the norm of the iterate.

Now do the following comparison; take several values for the time step.

#### Direct solver

If your PETSc installation includes direct solvers such as `MUMPS`, you can invoke them with

```
myprog -pc_type lu -ksp_type preonly \
       -pc_factor MatSolver_package mumps
```

Run your code with a direct solver, both sequentially and in parallel, and record how long it takes for the error to get down to $10^{-6}$.

#### Iterative Solver

Use an iterative solver, for instance `KSPCG` and `KSPBCGS`. Experiment with the convergence tolerance: how many timesteps does it take to get a $10^{-6}$ error if you set the iterative method tolerance to $10^{-12}$, how much if you take a lesser tolerance?

Compare timings between direct and iterative method.

### 30.3.5 Reporting

Write your report using `LATEX` (tutorial 29). Use both `tables` and `graphs` to report numerical results. Use `gnuplot` (tutorial 25) or a related utility for graphs.
## 30.4 The memory wall

In this project you explore the ramifications of the *memory wall*; see section 1.3. This project involves literature search, but little programming.

### 30.4.1 Background

Wulf and McKee [166] observed trends in the average latency in processors. Let $t_m$ be the latency from memory, $t_c$ the latency from cache, and $p$ the probability of a cache hit, then the average latency is

$$t_{\text{avg}} = pt_c + (1-p)t_m.$$

As the gap between processor speed and memory speed increases, this latency will grow, unless one can drive $p$ down, that is, reduce the number of cache misses or at least lessen their importance.

### 30.4.2 Assignment

Do a literature search and discuss the following topics.

* What have the trends been in processor speed and memory speed? What bearing does the introduction of multicore chips have on this balance?

* Section 1.3.4.3 discussed various types of cache misses. Some misses are more related to the algorithm and some more to the hardware. What strategies have hardware designers used to less then impact of cache misses?

* Compulsory cache misses seem unavoidable, since they are a property of the algorithm. However, if they can be hidden (see section 1.3.2 for ‘latency hiding’) their performance impact disappears. Research the topic of prefetch streams and their relation to latency hiding.

* How does the size of a cacheline interact with this behaviour?

* Discuss compiler techniques that lessen the incidence of cache misses.

* Can you find examples of algorithm options, that is, algorithms that compute the same result (not necessarily in the arithmetic sense; compare direct square root versus iterative approximation) but with different computational behaviour?

For all of these topics you are encouraged to write simulations.
# Chapter 31

## Codes

This section contains several simple codes that illustrate various issues relating to the performance of a single CPU. The explanations can be found in section 1.6.

### 31.1 Preliminaries

#### 31.1.1 Hardware event counting

The codes in this chapter make calls to a library named PAPI for 'Performance Application Programming Interface' [19, 135]. This is a portable set of calls to query the hardware counters that are built into most processors. Since these counters are part of the processor hardware, they can measure detailed events such as cache misses without this measurement process disturbing the phenomenon it is supposed to observe.

While using hardware counters is fairly straightforward, the question of whether what they are reporting is what you actually meant to measure is another matter altogether. For instance, the presence of hardware prefetch streams (section 1.3.5) implies that data can be loaded into cache without this load being triggered by a cache miss. Thus, the counters may report numbers that seem off, or even impossible, under a naive interpretation.

#### 31.1.2 Test setup

In the next sections you will see several codes that try to do precise timings on some phenomenon. This requires some care if you want to be sure you are actually timing what you're interested in, and not something else.

First of all, of course you perform your experiment multiple times to get a good average timing. If the individual timings are too far apart, you can decide to throw away outliers, or maybe rethink your experiment: there may be a variable that you haven't accounted for.

Let's assume that your timings are all within an acceptable spread; next you have to compensate for cache effects stemming from the fact that you are doing an experiment multiple times. If your problem sizes is small, the first timing will bring all data in cache, and the second run will find it there, leading to much faster runtimes. You can solve this problem by, in between any pair of timings, touching an array that is larger than the cache size. This has the effect of flushing all experiment data from the cache.
On the other hand, if you want to time something that happens on data in cache, you want to be sure
that your data is in cache, so you could for instance write to the array, bringing it in cache, prior to your
experiment. This is sometimes referred to as *cache warming*, and a cache that contains the problem data is
called a *hot cache*. Failure to warm the cache leads to an results where the first run takes appreciably longer
than the rest.

**31.1.3 Memory issues**

For performance reasons it is often desirable to enforce *cacheline boundary alignment* of data. The easiest solution is to use *static allocation*, which puts the data on cacheline boundaries.

To put dynamically aligned data on cacheline boundaries, the following code does the tric:

```c
double *a;
a = malloc( /* some number of bytes */ +8 );
if ( (int)a % 8 != 0 ) { /* it is not 8-byte aligned */
    a += 1; /* advance address by 8 bytes, then */
    /* either: */
    a = ( (a>>3) <<3 );
    /* or: */
    a = 8 * ( ( (int)a )/8 );
}
```

This code allocates a block of memory, and, if necessary, shifts it right to have a starting address that is a multiple of 8.

However, a better solution is to use *posix_memalign*:

```c
int posix_memalign(
    void **memptr, size_t alignment, size_t size);
```

which allocates *size* bytes, aligned to a multiple of *alignment* bytes. For example:

```c
double x;
posix_memalign( (void**) &x, 64, N*sizeof(double) );
```

will allocate 64 doubles for *x*, and align the array to cacheline boundaries if there are 8 words per cacheline.

**31.2 Cache size**

This code demonstrates the fact that operations are more efficient if data is found in L1 cache, than in L2, L3, or main memory. To make sure we do not measure any unintended data movement, we perform one iteration to bring data in the cache before we start the timers.
```c
printf("Run: data set size=%d\n",size);

/* clear the cache by dragging the whole array through it */
for (n=0; n<arraysize; n++) array[n] = 0.;
/* now load the data in the highest cache level that fits */
for (n=0; n<size; n++) array[n] = 0.;

retval = PAPI_start_counter events,NEVENTS); PCHECK (retval);
/* run the experiment */
for (i=0; i<NRUNS; i++) {
    for (j=0; j<size; j++) array[j] = 2.3*array[j]+1.2;
}
retval = PAPI_stop_counter (values,NEVENTS); PCHECK (retval);
printf("size=%d\nTot cycles: %d\n",size,values[0]);
printf("cycles per array loc: %9.5f\n",size,values[0]/(1.*NRUNS*size));
printf("L1 misses:\t%d\nfraction of L1 lines missed:\t%9.5f\n",
       values[1],values[1]/(size/8.));
printf("L2 misses:\t%d\nfraction of L2 lines missed:\t%9.5f\n",
       values[2],values[2]/(size/8.));
printf("\n");
}
free (array);

return 0;
}
```

## 31.3 Cachelines

This code illustrates the need for small strides in vector code. The main loop operates on a vector, progressing by a constant stride. As the stride increases, runtime will increase, since the number of cachelines transferred increases, and the bandwidth is the dominant cost of the computation.

There are some subtleties to this code: in order to prevent accidental reuse of data in cache, the computation is preceded by a loop that accesses at least twice as much data as will fit in cache. As a result, the array is guaranteed not to be in cache.

```c
/*
 * File:      line.c
 * Author:    Victor Eijkhout <eijkhout@tacc.utexas.edu>
 *
 * Usage:     line
 */

#include "papi_test.h"
extern int TESTS_QUET;      /* Declared in test_utils.c */
#define PCHECK(e) \n
```
```c
values[3], (1.*L1WORDS)/values[3]);
printf("L1 accesses:\t%d\naccesses per operation:\t%9.5f\n",
values[2], (1.*L1WORDS)/values[2]);
printf("\n");
}
free(array);

return 0;
}
```

Note that figure 1.17 in section 1.6.4 only plots up to stride 8, while the code computes to 16. In fact, at stride 12 the prefetch behaviour of the Opteron changes, leading to peculiarities in the timing, as shown in figure 31.1.

<figure><img src="image_146.png" /><figcaption>Figure 31.1: Run time in kcycles and L1 reuse as a function of stride</figcaption></figure>

## 31.4 Cache associativity

This code illustrates the effects of cache associativity; see sections 1.3.4.9 and 1.6.6 for a detailed explanation. A number of vectors (dependent on the inner loop variable *i*) is traversed simultaneously. Their lengths are chosen to induce cache conflicts. If the number of vectors is low enough, cache associativity will resolve these conflicts; for higher values of *m* the runtime will quickly increase. By allocating the vectors with a larger size, the cache conflicts go away.

/*
* File: assoc.c
* Author: Victor Eijkhout <eijkhout@tacc.utexas.edu>
```c
array = (double*) malloc(13*(MAXN+8)*sizeof(double));
#else
    array = (double*) malloc(13*MAXN*sizeof(double));
#endif

/* clear the array and bring in cache if possible */
for (m=1; m<12; m++) {
    for (n=2048; n<MAXN; n=2*n) {
        printf("Run: %d,%d\n",m,n);
    #if defined(SHIFT)
        printf("shifted\n");
    #endif

    for (i=0; i<=m; i++)
        for (j=0; j<n; j++)
            array[INDEX(i,j,m+1,n)] = 0.;

    /* access the rows in a way to cause cache conflicts */
    retval = PAPI_start_counter events,NEVENTS); PCHECK retval);
    for (j=0; j<n; j++)
        for (i=1; i<=m; i++)
            array[INDEX(0,j,m+1,n)] += array[INDEX(i,j,m+1,n)];
    retval = PAPI_stop_counter values,NEVENTS); PCHECK retval);
    printf("m,n=%d,%d\n#elements:\t%d\nTot cycles: %d\n",
           m,n,m*n,values[1]);
    printf("L1 misses:\t%d\nmisses per accumulation:\t%9.5f\n\n",
           values[0],values[0]/(1.*n));
}
}
free(array);

return 0;
}
```

## 31.5 TLB

This code illustrates the behaviour of a TLB; see sections 1.3.8 and 1.6.5 for a thorough explanation. A two-dimensional array is declared in column-major ordering (Fortran style). This means that striding through the data by varying the *i* coordinate will have a high likelihood of TLB hits, since all elements on a page are accessed consecutively. The number of TLB entries accessed equals the number of elements divided by the page size. Striding through the array by the *j* coordinate will have each next element hitting a new page, so TLB misses will ensue when the number of columns is larger than the number of TLB entries.

/*
 * File:     tlb.c
 * Author:   Victor Eijkhout <eijkhout@tacc.utexas.edu>
```c
{
    int events[NEVENTS] = {PAPI_TLB_DM, PAPI_TOT_CYC};
    long_long values[NEVENTS];
    int retval, order=COL;
    PAPI_event_info_t info, info1;
    const PAPI hw_info_t *hwinfo = NULL;
    int event_code;
    const PAPI_substrate_info_t *s = NULL;

    tests_quiet(argc, argv); /* Set TESTS_QUIET variable */
    if (argc==2 && !strcmp(argv[1], "row")) {
        printf("wrong way\n"); order=ROW;
    } else printf("right way\n");

    retval = PAPI_library_init(PAPI_VER_CURRENT);
    if (retval != PAPI_VER_CURRENT)
        test_fail(__FILE__, __LINE__, "PAPI_library_init", retval);

    {
        int i;
        for (i=0; i<NEVENTS; i++) {
            retval = PAPI_query_event events[i]); PCHECK retval);
        }
    }

#define M 1000
#define N 2000

{
    int m,n;
    m = M;
    array = (double*) malloc(M*N*sizeof(double));
    for (n=10; n<N; n+=10) {
        if (order==COL)
            clear_right(m,n);
        else
            clear_wrong(m,n);
        retval = PAPI_start_counter events, NEVENTS); PCHECK retval);
        if (order==COL)
            do_operations(m,n);
        else
            do_operations_wrong(m,n);
        retval = PAPI_stop_counter values, NEVENTS); PCHECK retval);
        printf("m,n=%d,%d\n#elements:\t%d\n",m,n,m*n);
        printf("Tot cycles: %d\nTLB misses:\t%d\nmisses per column:\t%9.5f\n\n",
               values[1], values[0], values[0] / (1.*n));
    }
    free(array);
}
```
```c
        return 0;
    }
```

## 31.6 Unrepresentible numbers

Section 3.2.1 mentioned that many fractions are not representible in binary. This code illustrates that by dividing a number first by 7, then multiplying it by 7 and dividing by 49. Compile and execute this code, and see that an input of 3.5 gives equality, but 3.6 doesn't.

```c
#include <stdlib.h>
#include <stdio.h>
int main() {
    double x, div1, div2;
    scanf("%lg", &x);
    div1 = x/7; div2 = (7*x)/49;
    printf("%e %2.17e %2.17e\n", x, div1, div2);
    if (div1==div2) printf("Lucky guess\n");
    else printf("Bad luck\n");
}
```
# Chapter 32

## Bibliography, index, and list of acronyms

### 32.1 Bibliography

[1] Alfred V. Aho, Brian W. Kernighan, and Peter J. Weinberger. *The Awk Programming Language*. Addison-Wesley Series in Computer Science. Addison-Wesley Publ., 1988. ISBN 020107981X, 9780201079814.

[2] G. Amdahl. The validity of the single processor approach to achieving large scale computing capabilities. In *Proceedings of the AFIPS Computing Conference*, volume 30, pages 483-485, 1967.

[3] O. Axelsson and A.V. Barker. *Finite element solution of boundary value problems. Theory and computation*. Academic Press, Orlando, Fl., 1984.

[4] Owe Axelsson and Ben Polman. Block preconditioning and domain decomposition methods II. *J. Comp. Appl. Math.*, 24:55-72, 1988.

[5] Josh Barnes and Piet Hut. A hierarchical $o(nlogn)$ force-calculation algorithm. *Nature*, 324:446-449, 1986.

[6] Richard Barrett, Michael Berry, Tony F. Chan, James Demmel, June Donato, Jack Dongarra, Victor Eijkhout, Roldan Pozo, Charles Romine, and Henk van der Vorst. *Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods*. SIAM, Philadelphia PA, 1994. http://www.netlib.org/templates/.

[7] K.E. Batcher. MPP: A high speed image processor. In *Algorithmically Specialized Parallel Computers*. Academic Press, New York, 1985.

[8] Gordon Bell. The outlook for scalable parallel processing. Decision Resources, Inc, 1994.

[9] Abraham Berman and Robert J. Plemmons. *Nonnegative Matrices in the Mathematical Sciences*. SIAM, 1994. originally published by Academic Press, 1979, New York.

[10] Petter E. Bjorstad, William Gropp, and Barry Smith. *Domain decomposition : parallel multilevel methods for elliptic partial differential equations*. Cambridge University Press, 1996.

[11] Fischer Black and Myron S Scholes. The pricing of options and corporate liabilities. *Journal of Political Economy*, 81(3):637-54, May-June 1973.

[12] L.S. Blackford, J. Choi, A. Cleary, E. D'Azevedo, J. Demmel, I. Dhillon, J. Dongarra, S. Hammerling, G. Henry, A. Petitet, K. Stanley, D. Walker, and R.C. Whaley. *ScaLAPACK Users' Guide*. SIAM, 1997.

[13] Netlib.org BLAS reference implementation. http://www.netlib.org/blas/

[14] Guy E. Blelloch, Michael A. Heroux, and Marco Zagha. Segmented operations for sparse matrix computation on vector multiprocessors. Technical Report CMU-CS-93-173, CMU, 1993.
[15] Mark Bohr. A 30 year retrospective on Dennard's MOSFET scaling paper. *Solid-State Circuits Newsletter*, IEEE, 12(1):11–13, winter 2007.

[16] Mark Bohr. The new era of scaling in an soc world. In *ISSCC*, pages 23–28, 2009.

[17] Jeff Bolz, Ian Farmer, Eitan Grinspun, and Peter Schröoder. Sparse matrix solvers on the gpu: conjugate gradients and multigrid. *ACM Trans. Graph.*, 22(3):917–924, July 2003.

[18] BOOST interval arithmetic library. http://www.boost.org/libs/numeric/interval/doc/interval.htm.

[19] S. Browne, J Dongarra, N. Garner, G. Ho, and P. Mucci. A portable programming interface for performance evaluation on modern processors. *International Journal of High Performance Computing Applications*, 14:189–204, Fall 2000.

[20] Alfredo Buttari, Victor Eijkhout, Julien Langou, and Salvatore Filippone. Performance optimization and modeling of blocked sparse kernels. *Int. J. High Perf. Comput. Appl.*, 21:467–484, 2007.

[21] P. M. Campbell, K. D. Devine, J. E. Flaherty, L. G. Gervasio, and J. D. Teresco. Dynamic octree load balancing using space-filling curves, 2003.

[22] Ernie Chan, Marcel Heimlich, Avi Purkayastha, and Robert van de Geijn. Collective communication: theory, practice, and experience. *Concurrency and Computation: Practice and Experience*, 19:1749–1783, 2007.

[23] A. P. Chandrakasan, R. Mehra, M. Potkonjak, J. Rabaey, and R. W. Brodersen. Optimizing power using transformations. *IEEE Transaction on Computer Aided Design of Integrated Circuits and Systems*, pages 13–32, January 1995.

[24] Chapel programming language homepage. http://chapel.cray.com/.

[25] Barbara Chapman, Gabriele Jost, and Ruud van der Pas. *Using OpenMP: Portable Shared Memory Parallel Programming*, volume 10 of *Scientific Computation Series*. MIT Press, ISBN 0262533022, 9780262533027, 2008.

[26] Yaeyoung Choi, Jack J. Dongarra, Roldan Pozo, and David W. Walker. Scalapack: a scalable linear algebra library for distributed memory concurrent computers. In *Proceedings of the fourth symposium on the frontiers of massively parallel computation (Frontiers '92)*, McLean, Virginia, Oct 19–21, 1992, pages 120–127, 1992.

[27] A. Chronopoulos and C.W. Gear. *s*-step iterative methods for symmetric linear systems. *Journal of Computational and Applied Mathematics*, 25:153–168, 1989.

[28] Barry A. Cipra. An introduction to the ising model. *The American Mathematical Monthly*, pages 937–959, 1987.

[29] Charles Clos. A study of non-blocking switching networks. *Bell System Technical Journal*, 32:406–242, 1953.

[30] E. Cuthill and J. McKee. Reducing the bandwidth of sparse symmetric matrices. In *ACM proceedings of the 24th National Conference*, 1969.

[31] Eduardo F. D'Azevedo, Mark R. Fahey, and Richard T. Mills. Vectorized sparse matrix multiply for compressed row storage format. *Lecture Notes in Computer Science, Computational Science ICCS 2005*, pages 99–106, 2005.

[32] E.F. D'Azevedo, V.L. Eijkhout, and C.H. Romine. A matrix framework for conjugate gradient methods and some variants of cg with less synchronization overhead. In *Proceedings of the Sixth SIAM Conference on Parallel Processing for Scientific Computing*, pages 644–646, Philadelphia, 1993. SIAM.
[33] Jeffrey Dean and Sanjay Ghemawat. MapReduce: Simplified data processing on large clusters. In *OSDI'04: Sixth Symposium on Operating System Design and Implementation*, 2004.

[34] J. Demmel, M. Heath, and H. Van der Vorst. Parallel numerical linear algebra. In *Acta Numerica 1993*. Cambridge University Press, Cambridge, 1993.

[35] James Demmel, Mark Hoemmen, Marghoob Mohiyuddin, and Katherine Yelick. Avoiding communication in sparse matrix computations. In *IEEE International Parallel and Distributed Processing Symposium*, 2008.

[36] Jim Demmel, Jack Dongarra, Victor Eijkhout, Erika Fuentes, Antoine Petitet, Rich Vuduc, R. Clint Whaley, and Katherine Yelick. Self adapting linear algebra algorithms and software. *Proceedings of the IEEE*, 93:293–312, February 2005.

[37] R.H. Dennard, F.H. Gaensslen, V.L. Rideout, E. Bassous, and A.R. LeBlanc. Design of ion-implanted mosfet's with very small physical dimensions. *Solid-State Circuits, IEEE Journal of*, 9(5):256–268, oct 1974.

[38] Ashish Deshpande and Martin Schultz. Efficient parallel programming with linda. In *In Supercomputing '92 Proceedings*, pages 238–244, 1992.

[39] E. W. Dijkstra. Cooperating sequential processes. http://www.cs.utexas.edu/users/EWD/transcriptions/EWD01xx/EWD123.html. Technological University, Eindhoven, The Netherlands, September 1965.

[40] Edsger W. Dijkstra. Programming as a discipline of mathematical nature. *Am. Math. Monthly*, 81:608–612, 1974.

[41] Edsger W. Dijkstra. Programming considered as a human activity. published as EWD:EWD117pub, n.d.

[42] J. Dongarra, A. Geist, R. Manchek, and V. Sunderam. Integrated PVM Framework Supports Heterogeneous Network Computing. *Computers in Physics*, 7(2):166–75, April 1993.

[43] J. J. Dongarra. *The LINPACK benchmark: An explanation*, volume 297, chapter Supercomputing 1987, pages 456–474. Springer-Verlag, Berlin, 1988.

[44] Jack J. Dongarra, Jeremy Du Croz, Sven Hammarling, and Iain Duff. A set of level 3 basic linear algebra subprograms. *ACM Transactions on Mathematical Software*, 16(1):1–17, March 1990.

[45] Jack J. Dongarra, Jeremy Du Croz, Sven Hammarling, and Richard J. Hanson. An extended set of FORTRAN basic linear algebra subprograms. *ACM Transactions on Mathematical Software*, 14(1):1–17, March 1988.

[46] Dale Dougherty and Arnold Robbins. *sed & awk*. O'Reilly Media, 2nd edition edition. Print ISBN: 978-1-56592-225-9, ISBN 10:1-56592-225-5; Ebook ISBN: 978-1-4493-8700-6, ISBN 10:1-4493-8700-4.

[47] Dr. Dobbs. *Complex arithmetic: in the intersection of C and C++*. http://www.ddj.com/cpp/184401628.

[48] I. S. Duff, R. G. Grimes, and J. G. Lewis. Users' guide for the Harwell-Boeing sparse matrix collection (release I). Technical Report RAL 92-086, Rutherford Appleton Laboratory, 1992.

[49] C. Edwards, P. Geng, A. Patra, and R. van de Geijn. Parallel matrix distributions: have we been doing it all wrong? Technical Report TR-95-40, Department of Computer Sciences, The University of Texas at Austin, 1995.

[50] Victor Eijkhout. A general formulation for incomplete blockwise factorizations. *Comm. Appl. Numer. Meth.*, 4:161–164, 1988.
[51] Victor Eijkhout, Paolo Bientinesi, and Robert van de Geijn. Towards mechanical derivation of Krylov solver libraries. *Procedia Computer Science*, 1(1):1805–1813, 2010. Proceedings of ICCS 2010, http://www.sciencedirect.com/science/publication?issn=18770509&volume=1&issue=1.

[52] Paul Erdös and A. Rnyi. On the evolution of random graphs. *Publications of the Mathematical Institute of the Hungarian Academy of Sciences*, 5:17–61, 1960.

[53] V. Faber and T. Manteuffel. Necessary and sufficient conditions for the existence of a conjugate gradient method. *SIAM J. Numer. Anal.*, 21:352–362, 1984.

[54] R.D. Falgout, J.E. Jones, and U.M. Yang. Pursuing scalability for hypre's conceptual interfaces. Technical Report UCRL-JRNL-205407, Lawrence Livermore National Lab, 2004. submitted to ACM Transactions on Mathematical Software.

[55] M. Fiedler. A property of eigenvectors of nonnegative symmetric matrices and its applications to graph theory. *Czechoslovak Mathematical Journal*, 25:618–633, 1975.

[56] D. C. Fisher. Your favorite parallel algorithms might not be as fast as you think. *IEEE Trans. Computers*, 37:211–213, 1988.

[57] M. Flynn. Some computer organizations and their effectiveness. *IEEE Trans. Comput.*, C-21:948, 1972.

[58] Project fortress homepage. http://projectfortress.sun.com/Projects/Community.

[59] D. Frenkel and B. Smit. Understanding molecular simulations: From algorithms to applications, 2nd edition. 2002.

[60] Roland W. Freund and Noël M. Nachtigal. QMR: a quasi-minimal residual method for non-Hermitian linear systems. *Numer. Math.*, 60:315–339, 1991.

[61] M. Frigo, Charles E. Leiserson, H. Prokop, and S. Ramachandran. Cache oblivious algorithms. In *Proc. 40th Annual Symposium on Foundations of Computer Science*, 1999.

[62] A. Geist, A. Beguelin, J. Dongarra, W. Jiang, R. Manchek, and V. Sunderam. *PVM: A Users' Guide and Tutorial for Networked Parallel Computing*. MIT Press, 1994. The book is available electronically, the url is ftp://www.netlib.org/pvm3/book/pvm-book.ps.

[63] David Gelernter. Generative communication in Linda. *ACM Transactions on Programming Languages and Systems*, 7:80–112, 1985.

[64] David Gelernter and Nicholas Carriero. Coordination languages and their significance. *Commun. ACM*, 35(2):97–107, 1992.

[65] GNU multiple precision library. http://gmplib.org/.

[66] Stefan Goedecker and Adolfy Hoisie. *Performance Optimization of Numerically Intensive Codes*. SIAM, 2001.

[67] D. Goldberg. Compuer arithmetic. Appendix in [84].

[68] David Goldberg. What every computer scientist should know about floating-point arithmetic. *Computing Surveys*, March 1991.

[69] G. H. Golub and D. P. O'Leary. Some history of the conjugate gradient and Lanczos algorithms: 1948-1976. 31:50–102, 1989.

[70] G. H. Golub and C. F. Van Loan. *Matrix Computations*. North Oxford Academic, Oxford, 1983.

[71] Gene H. Golub and Charles F. Van Loan. *Matrix Computations*. The Johns Hopkins University Press, Baltimore, second edition, 1989.
[72] Kazushige Goto and Robert A. van de Geijn. Anatomy of high-performance matrix multiplication. *ACM Trans. Math. Softw.*, 34(3):1–25, 2008.

[73] F. Gray. Pulse code communication. U.S. Patent 2,632,058, March 17, 1953 (filed Nov. 1947).

[74] Ronald I. Greenberg and Charles E. Leiserson. Randomized routing on fat-trees. In *Advances in Computing Research*, pages 345–374. JAI Press, 1989.

[75] W. Gropp, E. Lusk, and A. Skjellum. *Using MPI*. The MIT Press, 1994.

[76] William Gropp, Steven Huss-Lederman, Andrew Lumsdaine, Ewing Lusk, Bill Nitzberg, William Saphir, and Marc Snir. *MPI: The Complete Reference, Volume 2 - The MPI-2 Extensions*. MIT Press, 1998.

[77] William Gropp, Thomas Sterling, and Ewing Lusk. *Beowulf Cluster Computing with Linux, 2nd Edition*. MIT Press, 2003.

[78] Hadoop wiki. http://wiki.apache.org/hadoop/FrontPage.

[79] Louis A. Hageman and David M. Young. *Applied Iterative Methods*. Academic Press, New York, 1981.

[80] A. Hartstein, V. Srinivasan, T. R. Puzak, and P. G. Emma. Cache miss behavior: is it &#8730;2? In *Proceedings of the 3rd conference on Computing frontiers*, CF '06, pages 313–320, New York, NY, USA, 2006. ACM.

[81] Michael T. Heath. *Scientific Computing: an introductory survey; second edition*. McGraw Hill, 2002.

[82] Don Heller. A survey of parallel algorithms in numerical linear algebra. *SIAM Review*, 20:740–777, 1978.

[83] B. A. Hendrickson and D. E. Womble. The torus-wrap mapping for dense matrix calculations on massively parallel computers. *SIAM J. Sci. Comput.*, 15(5):1201–1226, 1994.

[84] John L. Hennessy and David A. Patterson. *Computer Architecture, A Quantitative Approach*. Morgan Kaufman Publishers, 3rd edition edition, 1990, 3rd edition 2003.

[85] M.R. Hestenes and E. Stiefel. Methods of conjugate gradients for solving linear systems. *Nat. Bur. Stand. J. Res.*, 49:409–436, 1952.

[86] Nicholas J. Higham. *Accuracy and Stability of Numerical Algorithms*. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, second edition, 2002.

[87] Jonathan Hill, Bill McColl, Dan C. Stefanescu, Mark W. Goudreau, Kevin Lang, Satish B. Rao, Torsten Suel, Thanasis Tsantilas, and Rob H. Bisseling. BSPlib: The BSP programming library. *Parallel Computing*, 24(14):1947–1980, 1998.

[88] C. A. R. Hoare. An axiomatic basis for computer programming. *Communications of the ACM*, pages 576–580, October 1969.

[89] Torsten Hoefler, Christian Siebert, and Andrew Lumsdaine. Scalable communication protocols for dynamic sparse data exchange. *SIGPLAN Not.*, 45(5):159–168, January 2010.

[90] IEEE 754: Standard for binary floating-point arithmetic. http://grouper.ieee.org/groups/754.

[91] Interval arithmetic. http://en.wikipedia.org/wiki/Interval_(mathematics).

[92] C.R. Jesshope and R.W. Hockney editors. The DAP approach, volume 2. pages 311–329. Infotech Intl. Ltd., Maidenhead, 1979.

[93] M. T. Jones and P. E. Plassmann. The efficient parallel iterative solution of large sparse linear
systems. In A. George, J.R. Gilbert, and J.W.H. Liu, editors, *Graph Theory and Sparse Matrix Computations*, IMA Vol 56. Springer Verlag, Berlin, 1994.

[94] L. V. Kale and S. Krishnan. Charm++: Parallel programming with message-driven objects. In *Parallel Programming using C++, G. V. Wilson and P. Lu, editors*, pages 175–213. MIT Press, 1996.

[95] Michael Karbo. PC architecture. http://www.karbosguide.com/books/pcarchitecture/chapter00.htm.

[96] R.M. Karp and Y. Zhang. A randomized parallel branch-and-bound procedure. In *Proceedings of the Twentieth Annual ACM Symposium on Theory of Computing*, Chicago, IL, USA, 2-4 May 1988, pages 290–300. ACM Press, 1988.

[97] J. Katzenelson. Computational structure of the n-body problem. *SIAM Journal of Scientific and Statistical Computing*, 10:787–815, July 1989.

[98] Kendall Square Research. http://en.wikipedia.org/wiki/Kendall_Square_Research.

[99] Donald Knuth. *The Art of Computer Programming, Volume 2: Seminumerical algorithms*. Addison-Wesley, Reading MA, 3rd edition edition, 1998.

[100] Helmut Kopka and Patrick W. Daly. *A Guide to LATEX*. Addison-Wesley, first published 1992.

[101] Ulrich Kulisch. Very fast and exact accumulation of products. *Computing*, 91(4):397–405, April 2011.

[102] Ulrich Kulisch and Van Snyder. The exact dot product as basic tool for long interval arithmetic. *Computing*, 91(3):307–313, 2011.

[103] Milind Kulkarni, Martin Burtscher, Rajasekhar Inkulu, Keshav Pingali, and Calin Cascaval. How much parallelism is there in irregular applications? In *Principles and Practices of Parallel Programming (PPoPP)*, 2009.

[104] Vipin Kumar, Ananth Grama, Anshul Gupta, and George Karypis. *Introduction to Parallel Computing*. Benjamin Cummings, 1994.

[105] U. Kung, Charalampos E. Tsourakakis, and Christos Faloutsos. Pegasus: A peta-scale graph mining system - implementation and observations. In *Proc. Intl. Conf. Data Mining*, pages 229–238, 2009.

[106] L. Lamport. *LATEX, a Document Preparation System*. Addison-Wesley, 1986.

[107] C. Lanczos. Solution of systems of linear equations by minimized iterations. *Journal of Research, Nat. Bu. Stand.*, 49:33–53, 1952.

[108] Rubin H Landau, Manual José Páez, and Cristian C. Bordeianu. *A Survey of Computational Physics*. Princeton University Press, 2008.

[109] J. Langou. *Iterative methods for solving linear systems with multiple right-hand sides*. Ph.D. dissertation, INSA Toulouse, June 2003. CERFACS TH/PA/03/24.

[110] Amy N. Langville and Carl D. Meyer. A survey of eigenvector methods for web information retrieval. *SIAM Review*, 47(1):135–161, 2005.

[111] C. L. Lawson, R. J. Hanson, D. R. Kincaid, and F. T. Krogh. Basic linear algebra subprograms for fortran usage. *ACM Trans. Math. Softw.*, 5(3):308–323, September 1979.

[112] P. L'Ecuyer. Combined multiple recursive generators. *Operations Research*, 44, 1996.

[113] Charles E. Leiserson. Fat-Trees: Universal networks for hardware-efficient supercomputing. *IEEE Trans. Comput*, C-34:892–901, 1985.

[114] Richard J. Lipton, Donald J. Rose, and Robert Endre Tarjan. Generalized nested dissection. *SIAM J. Numer. Anal.*, 16:346–358, 1979.
[115] Richard J. Lipton and Robert Endre Tarjan. A separator theorem for planar graphs. *SIAM J. Appl. Math.*, 36:177–189, 1979.

[116] J.D.C. Little. A proof of the queueing formula $L = \lambda W$. *Ope. Res.*, pages 383–387, 1961.

[117] Yongchao Liu, Douglas L. Maskell, and Bertil Schmidt. CUDASW++: optimizing Smith-Waterman sequence database searches for CUDA-enabled graphics processing units. *BMC Res Notes*, 2:73, 2009. PMID-19416548.

[118] M. Luby. A simple parallel algorithm for the maximal independent set problem. *SIAM Journal on Computing*, 4, 1986.

[119] Grzegorz Malewicz, Matthew H. Austern, Aart J.C. Bik, James C. Dehnert, Ilan Horn, Naty Leiser, and Grzegorz Czajkowski. Pregel: a system for large-scale graph processing. In *Proceedings of the 28th ACM symposium on Principles of distributed computing*, PODC '09, pages 6–6, New York, NY, USA, 2009. ACM.

[120] M. Mascagni and A. Srinivasan. Algorithm 806: Sprng: A scalable library for pseudorandom number generation. *ACM Transactions on Mathematical Software*, 26:436–461, 2000.

[121] Robert Mecklenburg. *Managing Projects with GNU Make*. O'Reilly Media, 3rd edition edition, 2004. Print ISBN:978-0-596-00610-5
Ebook ISBN:978-0-596-10445-0
ISBN 10:0-596-10445-6.

[122] J.A. Meijerink and H.A. van der Vorst. An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix. *Math Comp*, 31:148–162, 1977.

[123] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller. Equation of state calculations by fast computing machines. *Journal of Chemical Physics*, 21:1087–1092, June 1953.

[124] Gerard Meurant. Multitasking the conjugate gradient method on the CRAY X-MP/48. *Parallel Computing*, 5:267–280, 1987.

[125] Gérard Meurant. Domain decomposition methods for partial differential equations on parallel computers. *Int. J. Supercomputing Appl.*, 2:5–12, 1988.

[126] Frank Mittelbach, Michel Goossens, Johannes Braams, David Carlisle, and Chris Rowley. *The LATEX Companion*, 2nd edition. Addison-Wesley, 2004.

[127] Kenneth Moreland and Ron Oldfield. Formal metrics for large-scale parallel performance. In Julian M. Kunkel and Thomas Ludwig, editors, *High Performance Computing*, volume 9137 of *Lecture Notes in Computer Science*, pages 488–496. Springer International Publishing, 2015.

[128] Tobi Oetiker. The not so short introductino to LATEX. http://tobi.oetiker.ch/lshort/.

[129] National Institute of Standards and Technology. Matrix market. http://math.nist.gov/MatrixMarket.

[130] Andrew T. Ogielski and William Aiello. Sparse matrix computations on parallel processor arrays. *SIAM J. Sci. Stat. Comput.* in press.

[131] S. Otto, J. Dongarra, S. Hess-Lederman, M. Snir, and D. Walker. *Message Passing Interface: The Complete Reference*. The MIT Press, 1995.

[132] Michael L. Overton. *Numerical Computing with IEEE Floating Point Arithmetic*. SIAM, Philadelphia PA, 2001.

[133] Larry Page, Sergey Brin, R. Motwani, and T. Winograd. The pagerank citation ranking: Bringing order to the web, 1998.

[134] V. Ya. Pan. New combination of methods for the acceleration of matrix multiplication. *Comp. & Maths. with Appl.*, 7:73–125, 1981.
[135] Performance application programming interface. http://icl.cs.utk.edu/papi/.

[136] S. Plimpton. Fast parallel algorithms for short-range molecular dynamics. *J. Comput. Phys.*, 117:1–19, 1995.

[137] Jack Poulson, Bryan Marker, Jeff R. Hammond, and Robert van de Geijn. Elemental: a new framework for distributed memory dense matrix computations. *ACM Transactions on Mathematical Software*. submitted.

[138] M. Püschel, B. Singer, J. Xiong, J. M. F. Moura, J. Johnson, D. Padua, M. Veloso, and R. W. Johnson. SPIRAL: A generator for platform-adapted libraries of signal processing algorithms. *Int’l Journal of High Performance Computing Applications*, 18(1):21–45, 2004.

[139] Texas Advanced Computing Center: Sun Constellation Cluster: Ranger. http://www.tacc.utexas.edu/resources/hpc/constellation.

[140] Gururaj S. Rao. Performance analysis of cache memories. *J. ACM*, 25:378–395, July 1978.

[141] J.K. Reid. On the method of conjugate gradients for the solution of large sparse systems of linear equations. In J.K. Reid, editor, *Large sparse sets of linear equations*, pages 231–254. Academic Press, London, 1971.

[142] Y. Saad. Iterative methods for sparse linear systems. PWS Publishing Company, Boston, 1996.

[143] Tetsuya Sato. The earth simulator: Roles and impacts. *Nuclear Physics B - Proceedings Supplements*, 129-130:102 – 108, 2004. Lattice 2003.

[144] David Schade. Canfar: Integrating cyberinfrastructure for astronomy. https://wiki.bc.net/atl-conf/display/BCNETPUBLIC/CANFAR+-+Integrating+Cyberinfrastructure+for+Astronomy.

[145] R. Schreiber. Scalability of sparse direct solvers. In A. George, J.R. Gilbert, and J.W.-H. Liu, editors, *Sparse Matrix Computations: Graph Theory Issues and Algorithms (An IMA Workshop Volume)*. Springer-Verlag, New York, 1993, 1993. also: Technical Report RIACS TR 92.13, NASA Ames Research Center, Moffet Field, Calif., May 1992.

[146] D. E. Shaw. A fast, scalable method for the parallel evaluation of distance-limited pairwise particle interactions. *J. Comput. Chem.*, 26:1318–1328, 2005.

[147] D. B. Skillicorn, Jonathan M. D. Hill, and W. F. McColl. Questions and answers about BSP, 1996.

[148] Marc Snir, Steve Otto, Steven Huss-Lederman, David Walker, and Jack Dongarra. *MPI: The Complete Reference, Volume 1, The MPI-1 Core*. MIT Press, second edition edition, 1998.

[149] Dan Spielman. Spectral graph theory, fall 2009. http://www.cs.yale.edu/homes/spielman/561/.

[150] G. W. Stewart. Communication and matrix computations on large message passing systems. *Parallel Computing*, 16:27–40, 1990.

[151] V. Strassen. Gaussian elimination is not optimal. *Numer. Math.*, 13:354–356, 1969.

[152] TeX frequently asked questions.

[153] Universal Parallel C at George Washington University. http://upc.gwu.edu/.

[154] Leslie G. Valiant. A bridging model for parallel computation. *Commun. ACM*, 33:103–111, August 1990.

[155] R. van de Geijn, Philip Alpatov, Greg Baker, Almadena Chtchelkanova, Joe Eaton, Carter Edwards, Murthy Guddati, John Gunnels, Sam Guyer, Ken Klimkowski, Calvin Lin, Greg Morrow, Peter Nagel, James Overfelt, and Michelle Pal. Parallel linear algebra package (PLAPACK): Release r0.1 (beta) users' guide. 1996.
[156] Robert A. van de Geijn. *Using PLAPACK: Parallel Linear Algebra Package*. The MIT Press, 1997.

[157] Robert A. van de Geijn and Enrique S. Quintana-Ortí. *The Science of Programming Matrix Computations*. www.lulu.com, 2008.

[158] Henk van der Vorst. Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems. *SIAM J. Sci. Stat. Comput.*, 13:631–644, 1992.

[159] Richard S. Varga. *Matrix Iterative Analysis*. Prentice-Hall, Englewood Cliffs, NJ, 1962.

[160] R. Vuduc, J. Demmel, and K. Yelikk. Oski: A library of automatically tuned sparse matrix kernels. In *Proceedings of SciDAC 2005, Journal of Physics: Conference Series, to appear.*, 2005.

[161] Richard W. Vuduc. *Automatic Performance Tuning of Sparse Matrix Kernels*. PhD thesis, University of California Berkeley, 2003.

[162] R. Clint Whaley, Antoine Petitet, and Jack J. Dongarra. Automated empirical optimization of software and the ATLAS project. *Parallel Computing*, 27(1-2):3-35, 2001. Also available as University of Tennessee LAPACK Working Note #147, UT-CS-00-448, 2000 (www.netlib.org/lapack/lawns/lawn147.ps).

[163] O. Widlund. On the use of fast methods for separable finite difference equations for the solution of general elliptic problems. In D.J. Rose and R.A. Willoughby, editors, *Sparse matrices and their applications*, pages 121–134. Plenum Press, New York, 1972.

[164] J.H. Wilkinson. *Rounding Errors in Algebraic Processes*. Prentice-Hall, Englewood Cliffs, N.J., 1963.

[165] Samuel Williams, Andrew Waterman, and David Patterson. Roofline: an insightful visual performance model for multicore architectures. *Commun. ACM*, 52:65-76, April 2009.

[166] Wm. A. Wulf and Sally A. McKee. Hitting the memory wall: implications of the obvious. *SIGARCH Comput. Archit. News*, 23(1):20–24, March 1995.

[167] L. T. Yand and R. Brent. The improved bicgstab method for large and sparse unsymmetric linear systems on parallel distributed memory architectures. In *Proceedings of the Fifth International Conference on Algorithms and Architectures for Parallel Processing*. IEEE, 2002.

[168] Randy Yates. Fixed point: an introduction. http://www.digitalsignallabs.com/fp.pdf, 2007.

[169] Andy Yoo, Edmond Chow, Keith Henderson, William McLendon, Bruce Hendrickson, and Umit Catalyurek. A scalable distributed parallel breadth-first search algorithm on BlueGene/L. In *Proceedings of the 2005 ACM/IEEE conference on Supercomputing, SC '05*, pages 25–, Washington, DC, USA, 2005. IEEE Computer Society.

[170] David M. Young. *Iterative method for solving partial differential equations of elliptic type*. PhD thesis, Harvard Univ., Cambridge, MA, 1950.
32.2 List of acronyms

AMR Adaptive Mesh Refinement

AOS Array-Of-Structures

API Application Programmer Interface

AVX Advanced Vector Extensions

BEM Boundary Element Method

BFS Breadth-First Search

BLAS Basic Linear Algebra Subprograms

BSP Bulk Synchronous Parallel

BVP Boundary Value Problem

CAF Co-array Fortran

CCS Compressed Column Storage

CG Conjugate Gradients

CGS Classical Gram-Schmidt

COO Coordinate Storage

CRS Compressed Row Storage

DAG Directed Acyclic Graph

DRAM Dynamic Random-Access Memory

DSP Digital Signal Processing

FD Finite Difference

FMA Fused Multiply-Add

FDM Finite Difference Method

FEM Finite Element Method

FMM Fast Multipole Method

FOM Full Orthogonalization Method

FPU Floating Point Unit

FFT Fast Fourier Transform

FSA Finite State Automaton

FSB Front-Side Bus

FPGA Field-Programmable Gate Array

GMRES Generalized Minimum Residual

GPU Graphics Processing Unit

GPGPU General Purpose Graphics Processing Unit

GS Gram-Schmidt

HDFS Hadoop File System

HPC High-Performance Computing

HPF High Performance Fortran

IBVP Initial Boundary Value Problem

ILP Instruction Level Parallelism

ILU Incomplete LU

IVP Initial Value Problem

LAN Local Area Network

LBM Lattice Boltzmann Method

LRU Least Recently Used

MIC Many Integrated Cores

MIMD Multiple Instruction Multiple Data

MGS Modified Gram-Schmidt

MPI Message Passing Interface

MSI Modified-Shared-Invalid

MTA Multi-Threaded Architecture

NUMA Non-Uniform Memory Access

ODE Ordinary Differential Equation

OS Operating System

PGAS Partitioned Global Address Space

PDE Partial Differential Equation

PRAM Parallel Random Access Machine

RDMA Remote Direct Memory Access

SAN Storage Area Network

SAS Software As a Service

SFC Space-Filling Curve

SIMD Single Instruction Multiple Data

SIMT Single Instruction Multiple Thread

SM Streaming Multiprocessor

SMP Symmetric Multi Processing

SOA Structure-Of-Arrays

SOR Successive Over-Relaxation

SP Streaming Processor

SPMD Single Program Multiple Data

SPD symmetric positive definite

SRAM Static Random-Access Memory

SSE SIMD Streaming Extensions

TLB Translation Look-aside Buffer

UMA Uniform Memory Access

UPC Unified Parallel C

WAN Wide Area Network
## 32.3 Index

* .bashrc, see shell, startup files

* .profile, see shell, startup files

2's complement, 153

* α, see latency

* β, see bandwidth

* γ, see computation rate

* active messages, 115

* acyclic graph, 351

* Adaptive Mesh Refinement (AMR), 308

* address space, 81

shared, 81

* adjacency

matrix, 322

* adjacency graph, 209, 353

* adjacency matrix, 352

* Advanced Vector Extensions (AVX), 79, 137, 328

* affinity, 94–95

* allgather, 237–238, 242, 246

* Alliant FX/8, 37

* allocation

static, 500

* allreduce, 237

* AMD, 79, 149

Barcelona, 32, 34

Opteron, 35, 46, 51, 52

* Amdahl's law, 72–74

* AMR, see Adaptive Mesh Refinement

* analytics

graph, 315

* AOS, see Array-Of-Structures

* API, see Application Programmer Interface

* Apple

iCloud, 145

* archive utility, 390

* arithmetic

computer, see floating point arithmetic

finite precision, 152

intensity, 39

* array processors, 78, 133

* array syntax, 109

* Array-Of-Structures (AOS), 117

* assembly

inline, 33

language, 24

* assertion, 451

* assertions, 450–452

* associative, see also cache, associative
* associativity

of mathematical operations, 163

* asynchronous communication, 119

* atomic operation, 37, 91–93

* atomicity, see atomic operation

* automaton, 359–360

linear bounded, 360

* autotuning, 56

* AVX, see Advanced Vector Extensions

* axpy, 40

* background process, 378

* backwards stability, 167

* banded matrix, 182

storage, 205–206

* bandwidth, 22, 131, 136, 235

aggregate, 121

measure in GT/s, 134

of a matrix, see halfbandwidth

of a matrix, 205

* bandwidth bound algorithm, 149

* bandwidth-bound, 41, 261

* Barnes-Hut algorithm, 326–327

* barrier, 116

* base, 155

* bash, 367

* Basic Linear Algebra Subprograms (BLAS), 440

* BBN

Butterfly, 126

* BEM, see Boundary Element Method

benchmarking, 28, 40

* BFS, see Breadth-First Search

bidirectional exchange, 238

* big-endian, 420, 427

* binary-coded-decimal, 156
bisection bandwidth, 121

bisection width, 121

BitBucket, 408

bitonic sequence, 312

Bitonic sort, 309

bitonic sort, 75, 312–313
sequential complexity of, 313

bits, 152

Bjam, 392

Black-Scholes model, 331

BLAS, see Basic Linear Algebra Subprograms
data format, 442–443

blis, 444

block Jacobi, 263–264, 279

block matrix, 184, 252, 272

block tridiagonal, 184, 215, 275

blocking communication, 99, 106

blocking for cache reuse, 49

Boundary Element Method (BEM), 256

Boundary Value Problem (BVP), 179–186

branch misprediction, 21

branch penalty, 21

breakpoint, 465, 466–467

Brent's theorem, 70, 139

broadcast, 102, 235–236

BSP, see Bulk Synchronous Parallel

bubble sort, 309

bucket brigade algorithm, 236

buffering, 240

bug, 450

Bulk Synchronous Parallel (BSP), 116

bus, 21

memory, 122

speed, 23

width, 23

butterfly exchange, 126

BVP, see Boundary Value Problem

by reference, 476

bytes, 152

C

array layout, 473–474

language standard, 168

C++

language standard, 168

linking to, 472–473

name mangling, 472

cache, 20, 21, 25–32

associative, 31

block, 28

blocking, 49, 54, 203

coherence, 27, 36–39, 46, 83, 122, 286

hierarchy, 26, 57

hit, 26

hot, 500

line, 28–29, 43

mapping, 29

miss, 26, 33

miss, capacity, 27

miss, compulsory, 27

miss, conflict, 27

miss, invalidation, 27

oblivious programming, 57–58

replacement policies, 28

shared, 82

warming, 500

cacheline

boundary alignment, 29, 500

invalidation, 37

CAF, see Co-array Fortran

Cannon's algorithm, see matrix-matrix product,

Cannon's algorithm

capability computing, 146

capacity computing, 146

cartesian mesh, 122

cartesian product, 115

cat, 368

Cayley-Hamilton theorem, 225

ccNUMA, 83

CCS, see Compressed Column Storage

cd, 370

CDC

Cyber205, 80

Cell processor, 133

CG, see Conjugate Gradients

CGS, see Classical Gram-Schmidt

channel rate, 131
channel width, 131
Chapel, 110, 113
characteristic polynomial, 225
Charm++, 115
checkerboard ordering, 184
Cheeger's
  constant, 357
  inequality, 357
chess, 86
chgrp, 373
chmod, 372
Cholesky factorization, 198, 284
Chomsky hierarchy, 360
Cilk Plus, 96
cleanup code, 48
clique, 353
clock speed, 19
Clos network, 129
cloud computing, 143–146
  service models, 145
cluster
  node, see node
clusters, 79, 80
  Beowulf, 80
Co-array Fortran (CAF), 112
co-processor, 133–134, 151
coherence
  cache, 109
collective communication, 102, 105, 235–238
collective operation, 103, 131, 235–238
  long vector, 236
  short vector, 236
colour number, 274, 352
colouring, 273
column-major, 61, 443, 473
communication
  blocking, 105
  overhead, 72
  overlapping computation with, 118, 283
compare-and-swap, 309, 313
compiler, 56, 85, 111
directives, 96, 112
flags, 169

optimization, 24
  vs round-off, 168–169
optimization levels, 159
parallelizing, 87
complexity, 344
  computational, 202–203, 344
  of iterative methods, 233
  space, 212–214, 344
Compressed Column Storage (CCS), 208
Compressed Row Storage (CRS), 207–209, 256
Compressed Row Storage (CRS)
  performance of the matrix-vector product, 261
computation rate, 131, 235
compute-bound, 41, 148
concurrency, 95
condition number, 167, 338
conditionally stable, 176
congestion, 121
Conjugate Gradients (CG), 232, 283
Connection Machine, 78, 150
contention, 121
context, 91
  switch, 91, 132, 136
control flow, 12, 20, 77, 285
conveniently parallel, 86
COO, see Coordinate Storage
coordinate storage, 208
coordination language, 114
core, 36, 14–36
  vs processor, 36
core dump, 458
correct rounding, 161
Courant-Friedrichs-Lewy condition, 179
CPU-bound, see compute-bound
cpu-bound, 14
Cramer's rule, 192
Crank-Nicolson method, 190
Cray, 252
  Cray-1, 80
  Cray-2, 80
  T3E, 107
  X/MP, 80
  XE6, 83
XMT, 91, 133
Y/MP, 80
Cray Inc., 133
Cray Research, 133
critical path, 70, 286
critical section, 92
crossbar, 82, 126
CRS, see Compressed Row Storage
csh, 367
CUDA, 78, 134, 135, 151
cut, 375
Cuthill-McKee ordering, 215–216, 277
cycle (in graph), 351
cyclic distribution, 250
Cygwin, 367
DAG, see Directed Acyclic Graph
data decomposition, 239
data flow, 12, 20, 77, 285
data parallel, 132, 135, 137
data parallelism, 67, 134
data race, see race condition
data reuse, 21, 39
ddd, 458
DDT, 458, 468
deadlock, 96, 101, 106, 467
debug flag, 459
dibugger, 458
debugging, 458–469
  in parallel, 467–469
DEC
  Alpha, 17, 80
defensive programming, 450
degree, 120, 350
Delauney mesh refinement, 85
Dennard scaling, 61
denormal, see floating point numbers, unnormalized
Dense linear algebra, 238–251
dependency, 19
dgemm, 59
diagonal dominance, 201
diagonal storage, 205–207, 253
diameter, 120, 351
die, 14, 26

difference stencil, **185**
differential operator, **180**
Dijkstra's shortest path algorithm, 319
direct mapping, 30
direct methods
  for linear systems, 261
Directed Acyclic Graph (DAG), 285, 351
directives, see compiler, directives
directories, 367
Dirichlet boundary condition, 179, 181
discretization, 175
distributed computing, 143–146
divide-and-conquer, 57
domain decomposition, 266
double precision, 136
DRAM, see Dynamic Random-Access Memory
DSP, see Digital Signal Processing
dynamic programming, 317
Dynamic Random-Access Memory (DRAM), 25
Earth Simulator, 149
ebook, 145
edges, 350
efficiency, 68
eigenvector
  dominant, 341
Eispack, 440
Elliptic PDEs, 180
embarrassingly parallel, 68, 86
embedding, 122
environment variable, 376, 379–381
Eratosthenes
  sieve of, 313
escape, 374, 385
ETA-10, 80
Euler
  explicit, 174–176
  implicit, 176–178
evicted, 38
excess, 155
executable, 367
exponent, 155
export, 380
extended precision, 162, 169, 170
Facebook, 315
factorization, see LU factorization
false sharing, 39, 99, 109
Fast Fourier Transform (FFT), 51
Fast Multipole Method (FMM), 327
fast solvers, 222
fat tree, 126–129
  bisection width of a, 128
  clusters based on, 129
fault tolerance, 147
FD, see Finite Difference
FDM, see Finite Difference Method
FEM, see Finite Element Method
fenv.h, 160
FFT, see Fast Fourier Transform
Fiedler vector, 357
Fiedler’s theorem, 357
field scaling, 61
files, 367
fill locations, 212
fill-in, 211–216, 275
  estimates, 212–214
  reduction, 214–216
finger, 385
finite difference, 175
Finite Difference Method (FDM), 181
Finite Element Method (FEM), 186, 260
Finite State Automaton (FSA), 38, 359–360
finite volume method, 186
first-touch policy, 95
floating point
  unit, 15–16
floating point arithmetic
  associativity of, 167, 361
floating point arithmetic, 152–171
floating point numbers
  normalized, 157
  representation, 155–156
  unnormalized, 157
floating point pipeline, 47
flops, 19
Floyd-Warshall algorithm, 317–319
  parallelization of the, 323

flushed, 28
flushing
  pipeline, 93
Flynn’s taxonomy, 77
FMA, see Fused Multiply-Add
FMM, see Fast Multipole Method
FOM, see Full Orthogonalization Method
foreground process, 378
fork-join, 88, 89
Fortran, 443
  array layout, 473–474
  declarations in, 169
  iso C bindings, 471
  language standard, 168
Fortress, 113–114
FPGA, see Field-Programmable Gate Array
FPU, see Floating Point Unit
Front-Side Bus (FSB), 21
FSA, see Finite State Automaton
FSB, see Front-Side Bus
Full Orthogonalization Method (FOM), 228
fully associative, 31
fully connected, 120
functional parallelism, 67
functional programming, 146–148
Fused Multiply-Add (FMA), 15, 162, 170
gather, 102, 167, 237
Gauss-Seidel, 220
GCC, 97
gdb, 458–467
General Purpose Graphics Processing Unit
  (GPGPU), 134
Generalized Minimum Residual (GMRES), 233
ghost region, 114, 254, 281, 282, 334
Global Arrays, 115
GMRES, see Generalized Minimum Residual
GNU, 458
  gdb, see gdb
  gnuplot, see gnuplot
  Make, see Make
gnuplot, 447
Goodyear
  MPP, 78
Google, 117, 256, 320
code, 408
Google Docs, 144, 145

Goto
Kazushige, 59

GPGPU, see General Purpose Graphics Processing Unit

GPU, see Graphics Processing Unit

Gram-Schmidt, 228, 338–340
modified, 228, 339

granularity, 87, 91

Grape computer, 133, 327

graph
adjacency, 215, see adjacency graph
colouring, 264, 274, 352
directed, 350
Laplacian, 140, 356
random, 320
social, 315
theory, 350–358
of parallel computers, 120
undirected, 120, 210, 350

Graphics Processing Unit (GPU), 80, 134–137, 150, 160, 257, 277

Gray code, 125

grep, 373

grid (CUDA), 135

grid computing, 143

groups, 373

GS, see Gram-Schmidt

guard digit, 161, 163

Gustafson's law, 73

Hadoop, 147

Hadoop File System (HDFS), 147

halfbandwidth, 213
left, 206
right, 206

halo, see ghost region

handshake protocol, 109

hardware prefetch, 33

Harwell-Boeing matrix format, 208

HDFS, see Hadoop File System

heap, 88

heat equation, 179, 186

Hessenberg matrix, 227

heterogeneous computing, 151

hg, see mercurial

High Performance Fortran (HPF), 109, 112

High-Performance Computing (HPC), 114, 143

HITS, 321–322

Horner's rule, 280, 323

host process, 133

host processor, 151

Householder reflectors, 343
in LU factorization, 197

HPC, see High-Performance Computing

HPF, see High Performance Fortran

hybrid computing, 108–109

hyper-threading, 89

Hyperbolic PDEs, 179

hypercube, 124

hyperthreading, 35, 91, 96

hypervisor, 83

I/O subsystem, 20

IBM, 114, 120, 155–157, 427
BlueGene, 83, 130, 150
BlueGene Q, 33
Power 5, 17
Power series, 80
Power6, 156
Roadrunner, 63, 133

IBVP, see Initial Boundary Value Problem
ICL
DAP, 78, 133

idle, 138

idle time, 106

ILP, see Instruction Level Parallelism

ILU, see Incomplete LU
imbalance
load, 116

incidence matrix, 321

Incomplete LU (ILU), 222

Incomplete LU (ILU)
parallel, 264–265

independent sets, 352

indirect addressing, 208
Inf, 156

Initial Boundary Value Problem (IBVP), 186–190

Initial Boundary Value Problems (IBVP), 179

Initial Value Problem (IVP), 172–178

inner products, 259–260

input redirection, see redirection

instruction

handling

in-order, 15

out-of-order, 15

issue, 79

pipeline, 21

Instruction Level Parallelism (ILP), 14, 20, 62, 85

Intel, 20, 35, 36, 79, 149

80287 co-processor, 170

Haswell, 17

i860, 80, 133

Itanium, 24, 149

MIC, 151

Paragon, 133

Sandy Bridge, 14, 17, 26

Woodcrest, 32

Xeon Phi, 15, 78, 91, 96, 109, 134, 137

bandwidth, 134

inter-node communication, 108

interior nodes, 352

interoperability

C to Fortran, 470–476

interrupt, 157

intra-node communication, 108

irreducible, 354, see reducible

Ising model, 333–334

ispc, 118

iterative methods, 216–233

floating point performance of, 261

iterative refinement, 221

IVP, see Initial Value Problem

Jacobi method, 219

kernel

CUDA, 135

Krylov methods, 216–233

ksh, 367

LAN, see Local Area Network

language

context-free, 360

context-sensitive, 360

recursively enumerable, 360

regular, 360

language interoperability, see interoperability

Lapack, 148, 440

routines, 441–442

Laplace equation, 345, 357

latency, 22, 130, 235

hiding, 23, 39

latency hiding, 118

LATEX, see also TEX, 477–490

Lattice Boltzmann Methods (LBM), 160

LBM, see Lattice Boltzmann Method

leaf nodes, 352

level sets, 215

lexicographic ordering, 183, 184, 265

libraries

creating and using, 390–391

numerical, 430–446, see Lapack, see PETSc

Linda, 114

linear algebra, 337–343

software

dense, 440–446

sparse, 430–440

linear array (of processors), 122

linker, 390

LINPACK

benchmark, 40

Linpack, 148, 440

benchmark, 46, 148, 440

Linux

distributions, 367

kernel, 83

Lisp, 148

Little's law, 33

little-endian, 420, 427

load

balancing, 138–142

balancing, dynamic, 139

balancing,static, 139
rebalancing, 140, 141, 328
redistributing, 140, 141
unbalance, 68, 80, 87, 98, 138
Local Area Network (LAN), 143
local solve, 263
locality, 14
core, 46
in parallel computing, 131–132
spatial, 41, 43
temporal, 41
lock, 90, 92
loop unrolling, 47, 59, 168
LRU, see Least Recently Used
ls, 368
LU factorization, 196–204
computation in parallel, 249–250
solution in parallel, 248–249
M-matrix, 182, 223
machine epsilon, see machine precision
machine precision, 158
Make, 392–405
and LATEX, 404–405
automatic variables, 397
debugging, 403
template rules, 398
malloc, 88
Mandelbrot set, 86, 139
mantissa, 155–159, 162, 163
manycore, 135
MapReduce, 146–148
Markov
chains, 355
matrix, 356
MasPar, 78
master-worker paradigm, 85
matrix
adjacency, 316, see adjacency matrix
nonnegative, 321
norms, 338
associated, 338
permutation, 353
reducible, 322
stochastic, 322

storage, dense, 249–250
storage, sparse, 205–209
strictly upper triangular, 279
structurally symmetric, 140, 210, 258
tesselation, 252
transposition, 57
unit upper triangular, 279
Matrix Market, 208
matrix format, 208
matrix-matrix product, 59, 250
Cannon’s algorithm for, 251
Goto implementation, 59–61, 444
reuse in, 40
matrix-vector product, 58
sparse, see sparse, matrix-vector product
memory
access pattern, 27
banks, 34
distributed, 81
distributed shared, 83
hierarchy, 21
leak, 434, 455
model, 93
pages, 34
shared, 81
stall, 22
violations, 454
virtual shared, 83
wall, 21, 498
memory model
relaxed, 93
memory-bound, 14
mercurial, 415–419
Message Passing Interface (MPI), 102–108
Metropolis algorithm, 334
MGS, see Modified Gram-Schmidt
MIC, see Many Integrated Cores
MIMD, see Multiple Instruction Multiple Data
minimum
spanning tree, 235, 319
minimum degree ordering, 216
MIPS, 80
mkdir, 370
MKL, 137, 440, 444

modified Gramm-Schmidt, see Gram-Schmidt, modified

Moore's law, 61

MPI, see Message Passing Interface
MPI 3.0 draft, 105

MSI, see Modified-Shared-Invalid

MTA, see Multi-Threaded Architecture

multi-colouring, 274

Multi-Threaded Architecture (MTA), 132

multi-threading, 96, 132

multicore, 21, 27, **35–39**, 62, 82, 203, 284
motivated by power, 62

multigrid, 227

Multiple Instruction Multiple Data (MIMD), 80

MUMPS, 497

$n_{1/2}$, 17

N-body problems, 325–329

NaN, 156

natural ordering, 184

NBODY6, 327

nearest neighbour, 123

nested dissection, 215, 265–272, 320

netlib, 443

Neumann boundary condition, 181

Newton method, 223

nm, 390, 470

node, 108, 123

nodes, 350

non-blocking communication, 102, 105, 106, 119

non-local operation, 106

Non-Uniform Memory Access (NUMA), 82, 132

norm, 337

NP-complete, 146

null termination, 475

NUMA, see Non-Uniform Memory Access

NVidia, 134

Tesla, 160

object file, 390, 470

octtree, 327

hashed, 328

odd-even transposition sort, 309, 309

ODE, see Ordinary Differential Equation

offloading, 137

one-sided communication, 107–108, 114–116

OpenMP, 96–99

version 4, 80

Operating System (OS), 95, 367

option pricing, 331

order, 344

Ordinary Differential Equation (ODE), 172–178

OS, see Operating System

out-of-order, see instruction, handling, out-of-order

output redirection, see redirection

overdecomposition, 87, 116, 139, 249

overflow, 154, 156, 450

overflow bit, 154

overlays, 34

owner computes, 138

page table, 35

PageRank, 256, 320, **322–323**

pages

memory, see memory, pages

Parabolic PDEs, 179

parallel fraction, 72

parallel prefix, see prefix operation

Parallel Random Access Machine (PRAM), 71

parallelism

average, 69

data, 78, 84, 111

dynamic, 97

fine-grained, 84

instruction-level, 64

irregular, 86

parallelization

incremental, 97

parameter sweep, 86

partial derivates, 345

Partial Differential Equation (PDE), 179–190

partial differential equations, 345–346

partial ordering, 285

partial pivoting, 195

Partitioned Global Address Space (PGAS), 83, 109–115

path (graph theory), 351
patsubst, 399

PCI bus, 134

PCI-X, 137

PCI-X bus, 136

PDE, see Partial Differential Equation

peak performance, 19, 40, 41, 46

penta-diagonal, 184

permutation, 215

Perron vector, 321, 323

Perron-Frobenius theorem, 342

PETSc, 430–440

PGAS, see Partitioned Global Address Space
pipeline, 16–19, 46–48

depth, 21, 62

flush, 21

length, 21

processor, 79

stall, 21

pivoting, 193–195, 200–201

diagonal, 195

full, 195

partial, 195

pivots, 193

PLapack, 440

point-to-point communication, 102

Poisson equation, 180, 345

posix_memalign, 500

power

consumption, 61–63

efficiency, 78

wall, 62

power method, 321, 322, 340

PowerPC, 79

PRAM, see Parallel Random Access Machine
precision

double, see double precision

extended, see extended precision

machine, see machine precision

of the intermediate result, 162, 168

single, see single precision

preconditioner, 221–223

predicate, 316

prefetch data stream, 32

prefix operation, 279, 311, 361–363

prefix operations

sparse matrix vector product, 362–363

Pregel, 117

Prim's algorithm, 319

prime number

finding, 313–314

priority queue, 320

probability vector, 322

process, 88

affinity, see affinity

program counter, 13, 88, 96

protein interactions, 315

ps, 378

pseudo-random numbers, see random numbers

pthread, 89–91

purify, 463

PVM, 107

QR factorization, 339

Quicksort, 308, 310–312

race condition, 92, 93

radix point, 155

random numbers, 332–333

generator, 332–333

lagged Fibonacci, 332

linear congruential, 332

random placement, 116, 140

rate

computational, 76

RDMA, see Remote Direct Memory Access

re-association, 168

real numbers

representation of, 152

recursive doubling, 19, 29, 275, 277–279

Red Hat, 367

red-black ordering, 184, 272–274

redirection, 375, 377

reduce-scatter, 238, 245, 246, 258

reducible, 210

reduction, 236–237

and thread safety, 92

under multi-threading, 92
reduction operations, 98

redundancy, 121

redundant computation, 282

refinement
    adaptive, 141

region of influence, 179

register, 13, 21, **23–25**
    file, 23
    resident in, 24, 25
    spill, 24, 48
    variable, 25
    vector, 79, 79

Remote Direct Memory Access (RDMA), 107

remote method invocation, 115

remote procedure call, 143

repository, 406

representation error
    absolute, 158
    relative, 158

reproducibility, 168
bitwise, 168

residual, 217

resource contention, 95

revision control systems, see source code control
Riemann sums, 331
ring network, 122
roofline model, 41
round-off error analysis, 161–167
    in parallel computing, 167

round-robin
    storage assignment, 111
    task scheduling, 98

rounding, 157
    correct, 159

row-major, 443, 473

SAN, see Storage Area Network

SAS, see Software As a Service

satisfyability, 67

scalability, 74–77
    strong, 74, 243
    weak, 74, 243

Scalapack, 148, 440

ScaleMP, 83

scaling, 74–77

scheduling
    dynamic, 87, 98
    fair-share, 146
    job, 146
    static, 98

Scons, 392

search direction, 229, 233
sed, 375

segmentation fault, 461

segmentation violation, 454

semaphore, 92

separable problem, 222

separator, 265, 320

Sequent Symmetry, 37

sequential complexity, 309

sequential consistency, 93, 96, 285

sequential fraction, 72

serialized execution, 100

SFC, see Space-Filling Curve

SGE, 146

SGI, 133
UV, 83

sh, 367

shared libraries, 390

shell, 367
    command history, 404
    startup files, **384–385**

shortest path
    all-pairs, 315
    single source, 315

side effects, 148

side-effects, 451

sign bit, 153, 155

significant, 155

significant digits, 163

SIMD, see Single Instruction Multiple Data
    lanes, 79
    width, 41, 118

SIMD Streaming Extensions (SSE), 79, 133, 328

SIMT, see Single Instruction Multiple Thread

Single Instruction Multiple Data (SIMD), 118

Single Instruction Multiple Thread (SIMT), 135
single precision, 136
Single Program Multiple Data (SPMD), 80, 118

skyline storage, 213

Slurm, 146

SM, see Streaming Multiprocessor
small world, 321

smoothers, 227

SMP, see Symmetric Multi Processing
snooping, 38, 122

SOA, see Structure-Of-Arrays
socket, 27, 36, 108, 150

Software As a Service (SAS), 144

SOR, see Successive Over-Relaxation
sorting, 308–313
network, 309, 313

source code control, 406
distributed, 406, 407, 415

source-to-source transformations, 48

SP, see Streaming Processor
space-filling curve, 138, 140–142
spanning tree, 319

Spark, 147

sparse
linear algebra, 204–216

matrix, 182, 205, 353
from PDEs, 183

matrix-vector product
implementation, 208

in graph theory, 315

locality in, 209

parallel, 253–258

parallel setup, 258

performance of, 261

spatial locality, see locality, spatial

SPD, see symmetric positive definite
speculative execution, 21

speedup, 68, 80

SPMD, see Single Program Multiple Data

SRAM, see Static Random-Access Memory

SSE, see SIMD Streaming Extensions
stall, 33

static libraries, 390

Static Random-Access Memory (SRAM), 25

stationary iteration, 218

statistical mechanics, 334

steady state, 177, 179, 180, 187, 346, 496

stochastic matrix, 356

Storage Area Network (SAN), 143

storage by diagonals, 205

stride, 28, 43, 50

strip mining, 59

structurally symmetric, 353

Structure-Of-Arrays (SOA), 117

subdomain, 266

substructuring, 266

Subversion, 406, 408–414

Successive Over-Relaxation (SOR), 221

Sun
Ray, 144

superlinear speedup, 68

superscalar, 14, 20, 85

superstep, 66

supersteps, 108, 116

surface-to-volume, 87

svn, see Subversion

swapped, 34

switch, 126

leaf, 128, 132

symbol table, 459

Symmetric Multi Processing (SMP), 81, 122

symmetric positive definite (SPD), 201, 227, 232,
233

TACC, 128

ranger cluster, 129

Stampede cluster, 129, 134

task

parallelism, 85–86

queue, 86, 95, 317

task parallelism, 67

Taylor series, 174, 347–349

tcsh, 367

temporal locality, see locality, temporal

Tera Computer

MTA, 133

ternary arithmetic, 156

TEX, 477
environment variables, 484
thin client, 144
thread, 88–96, 136
  affinity, see affinity
  blocks, 135
  master, 88
  private data, 91
  safe, 92, 97
  shared data, 91
  spawning, 88
  team, 88
  use in OpenMP, 97
throughput computing, 134
Tianhe-1A, 134
Tianhe-2, 134
Titanium, 112
TLB, see Translation Look-aside Buffer
  miss, 35
top, 386
top 500, 148–150
topology, 119
torus, 123
  clusters based on, 130
TotalView, 458
touch, 370
transactional memory, 93
Translation Look-aside Buffer (TLB), 35, 506
tree graph, 352
tridiagonal, 184
tridiagonal matrix, 182
truncation, 157
truncation error, 175, 181
tuple space, 114
Turing machine, 360
Ubuntu, 367
UMA, see Uniform Memory Access
  unconditionally stable, 177
underflow, 156
  gradual, 156, 157, 169

Unified Parallel C (UPC), 110–112
Uniform Memory Access (UMA), 81, 122, 126
unitary basis transformations, 215
Unix
  user account, 385
  unsigned, 153
UPC, see Unified Parallel C
uptime, 386
utility computing, 143
valgrind, 463–464
value safety, 168
vector
  instructions, 79, 277
  norms, 337
  pipeline, see pipeline, processor
  register, see register, vector
  vector processor, 79, 277
  verbatim mode, 482
version control systems, see source code control
vertices, 350
Virtual memory, 34
virtualization, 143
von Neumann architectures, 12
WAN, see Wide Area Network
wavefront, 276–277
weak scaling, 72
weighted graph, 351
which, 369, 376
who, 385
whoami, 385
Wide Area Network (WAN), 143
wildcard, 373, 399
work pool, 139
World Wide Web, 315
X10, 114
x86, 149, 150
zsh, 367
ISBN 978-1-257-99254-6

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

9 781257 992546

90000

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