# Introduction to High Performance Scientific Computing

Victor Eijkhout

► To cite this version:

Victor Eijkhout. Introduction to High Performance Scientific Computing. 2015. [hal-01257189](https://hal.science/hal-01257189)

HAL Id: hal-01257189

https://hal.science/hal-01257189

Submitted on 19 Jan 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Introduction to High Performance Scientific Computing

Evolving Copy - open for comments

Victor Eijkhout

with

Edmond Chow, Robert van de Geijn

2nd edition, revision 2016
<figure><img src="image_2.png" /></figure>

Introduction to High-Performance Scientific Computing © Victor Eijkhout, distributed under a Creative Commons Attribution 3.0 Unported (CC BY 3.0) license and made possible by funding from The Saylor Foundation http://www.saylor.org.
## Preface

The field of high performance scientific computing lies at the crossroads of a number of disciplines and skill sets, and correspondingly, for someone to be successful at using high performance computing in science requires at least elementary knowledge of and skills in all these areas. Computations stem from an application context, so some acquaintance with physics and engineering sciences is desirable. Then, problems in these application areas are typically translated into linear algebraic, and sometimes combinatorial, problems, so a computational scientist needs knowledge of several aspects of numerical analysis, linear algebra, and discrete mathematics. An efficient implementation of the practical formulations of the application problems requires some understanding of computer architecture, both on the CPU level and on the level of parallel computing. Finally, in addition to mastering all these sciences, a computational scientist needs some specific skills of software management.

While good texts exist on numerical modeling, numerical linear algebra, computer architecture, parallel computing, performance optimization, no book brings together these strands in a unified manner. The need for a book such as the present became apparent to the author working at a computing center: users are domain experts who not necessarily have mastery of all the background that would make them efficient computational scientists. This book, then, teaches those topics that seem indispensable for scientists engaging in large-scale computations.

The contents of this book are a combination of theoretical material and self-guided tutorials on various practical skills. The theory chapters have exercises that can be assigned in a classroom, however, their placement in the text is such that a reader not inclined to do exercises can simply take them as statement of fact.

The tutorials should be done while sitting at a computer. Given the practice of scientific computing, they have a clear Unix bias.

**Public draft** This book is open for comments. What is missing or incomplete or unclear? Is material presented in the wrong sequence? Kindly mail me with any comments you may have.

You may have found this book in any of a number of places; the authoritative download location is http://www.tacc.utexas.edu/~eijkhout/istc/istc.html. That page also links to lulu.com where you can get a nicely printed copy.

Victor Eijkhout eijkhout@tacc.utexas.edu
Research Scientist
Texas Advanced Computing Center
The University of Texas at Austin

**Acknowledgement** Helpful discussions with Kazushige Goto and John McCalpin are gratefully acknowledged. Thanks to Dan Stanzione for his notes on cloud computing, Ernie Chan for his notes on scheduling of block algorithms, and John McCalpin for his analysis of the top500. Thanks to Elie de Brauwer, Susan Lindsey, and Lorenzo Pesce for proofreading and many comments. Edmond Chow wrote the chapter on Molecular Dynamics. Robert van de Geijn contributed several sections on dense linear algebra.
# Introduction

Scientific computing is the cross-disciplinary field at the intersection of modeling scientific processes, and the use of computers to produce quantitative results from these models. It is what takes a domain science and turns it into a computational activity. As a definition, we may posit

> The efficient computation of constructive methods in applied mathematics.

This clearly indicates the three branches of science that scientific computing touches on:

* Applied mathematics: the mathematical modeling of real-world phenomena. Such modeling often leads to implicit descriptions, for instance in the form of partial differential equations. In order to obtain actual tangible results we need a constructive approach.

* Numerical analysis provides algorithmic thinking about scientific models. It offers a constructive approach to solving the implicit models, with an analysis of cost and stability.

* Computing takes numerical algorithms and analyzes the efficacy of implementing them on actually existing, rather than hypothetical, computing engines.

One might say that ‘computing’ became a scientific field in its own right, when the mathematics of real-world phenomena was asked to be constructive, that is, to go from proving the existence of solutions to actually obtaining them. At this point, algorithms become an object of study themselves, rather than a mere tool.

The study of algorithms became especially important when computers were invented. Since mathematical operations now were endowed with a definable time cost, complexity of algorithms became a field of study; since computing was no longer performed in ‘real’ numbers but in representations in finite bitstrings, the accuracy of algorithms needed to be studied. Some of these considerations in fact predate the existence of computers, having been inspired by computing with mechanical calculators.

A prime concern in scientific computing is efficiency. While to some scientists the abstract fact of the existence of a solution is enough, in computing we actually want that solution, and preferably yesterday. For this reason, in this book we will be quite specific about the efficiency of both algorithms and hardware. It is important not to limit the concept of efficiency to that of efficient use of hardware. While this is important, the difference between two algorithmic approaches can make optimization for specific hardware a secondary concern.

This book aims to cover the basics of this gamut of knowledge that a successful computational scientist needs to master. It is set up as a textbook for graduate students or advanced undergraduate students; others can use it as a reference text, reading the exercises for their information content.
# Contents

I Theory 11

1 Single-processor Computing 12

1.1 The Von Neumann architecture 12

1.2 Modern processors 14

1.3 Memory Hierarchies 21

1.4 Multicore architectures 35

1.5 Locality and data reuse 39

1.6 Programming strategies for high performance 46

1.7 Power consumption 61

1.8 Review questions 63

2 Parallel Computing 64

2.1 Introduction 64

2.2 Functional parallelism versus data parallelism 67

2.3 Quantifying parallelism 68

2.4 Parallel Computers Architectures 77

2.5 Different types of memory access 81

2.6 Granularity of parallelism 83

2.7 Parallel programming 87

2.8 Topologies 119

2.9 Multi-threaded architectures 132

2.10 Co-processors 133

2.11 Remaining topics 138

3 Computer Arithmetic 152

3.1 Integers 152

3.2 Real numbers 154

3.3 Round-off error analysis 161

3.4 Compilers and round-off 168

3.5 More about floating point arithmetic 169

3.6 Conclusions 171

4 Numerical treatment of differential equations 172

4.1 Initial value problems 172

4.2 Boundary value problems 179

4.3 Initial boundary value problem 186
<table><tbody><tr><td>5</td><td>Numerical linear algebra</td><td>191</td></tr><tr><td>5.1</td><td>Elimination of unknowns</td><td>191</td></tr><tr><td>5.2</td><td>Linear algebra in computer arithmetic</td><td>194</td></tr><tr><td>5.3</td><td>LU factorization</td><td>196</td></tr><tr><td>5.4</td><td>Sparse matrices</td><td>204</td></tr><tr><td>5.5</td><td>Iterative methods</td><td>216</td></tr><tr><td>5.6</td><td>Further Reading</td><td>234</td></tr><tr><td>6</td><td>High performance linear algebra</td><td>235</td></tr><tr><td>6.1</td><td>Collective operations</td><td>235</td></tr><tr><td>6.2</td><td>Parallel dense matrix-vector product</td><td>238</td></tr><tr><td>6.3</td><td>LU factorization in parallel</td><td>248</td></tr><tr><td>6.4</td><td>Matrix-matrix product</td><td>250</td></tr><tr><td>6.5</td><td>Sparse matrix-vector product</td><td>252</td></tr><tr><td>6.6</td><td>Computational aspects of iterative methods</td><td>259</td></tr><tr><td>6.7</td><td>Parallel preconditioners</td><td>262</td></tr><tr><td>6.8</td><td>Ordering strategies and parallelism</td><td>265</td></tr><tr><td>6.9</td><td>Operator splitting</td><td>275</td></tr><tr><td>6.10</td><td>Parallelism and implicit operations</td><td>276</td></tr><tr><td>6.11</td><td>Grid updates</td><td>281</td></tr><tr><td>6.12</td><td>Block algorithms on multicore architectures</td><td>284</td></tr></tbody></table>

## II Applications 289

<table><tbody><tr><td>7</td><td>Molecular dynamics</td><td>290</td></tr><tr><td>7.1</td><td>Force Computation</td><td>291</td></tr><tr><td>7.2</td><td>Parallel Decompositions</td><td>295</td></tr><tr><td>7.3</td><td>Parallel Fast Fourier Transform</td><td>301</td></tr><tr><td>7.4</td><td>Integration for Molecular Dynamics</td><td>304</td></tr><tr><td>8</td><td>Combinatorial algorithms</td><td>308</td></tr><tr><td>8.1</td><td>Sorting</td><td>308</td></tr><tr><td>8.2</td><td>Prime number finding</td><td>313</td></tr><tr><td>9</td><td>Graph analytics</td><td>315</td></tr><tr><td>9.1</td><td>Traditional graph algorithms</td><td>315</td></tr><tr><td>9.2</td><td>'Real world' graphs</td><td>320</td></tr><tr><td>9.3</td><td>Hypertext algorithms</td><td>321</td></tr><tr><td>9.4</td><td>Large-scale computational graph theory</td><td>323</td></tr><tr><td>10</td><td>N-body problems</td><td>325</td></tr><tr><td>10.1</td><td>The Barnes-Hut algorithm</td><td>326</td></tr><tr><td>10.2</td><td>The Fast Multipole Method</td><td>327</td></tr><tr><td>10.3</td><td>Full computation</td><td>327</td></tr><tr><td>10.4</td><td>Implementation</td><td>328</td></tr><tr><td>11</td><td>Monte Carlo Methods</td><td>330</td></tr><tr><td>11.1</td><td>Motivation</td><td>330</td></tr></tbody></table>
11.2 *Parallel Random Number Generation* 332

11.3 *Examples* 333

### III Appendices 335

12 *Linear algebra* 337

12.1 *Norms* 337

12.2 *Gram-Schmidt orthogonalization* 338

12.3 *The power method* 340

12.4 *Nonnegative matrices; Perron vectors* 342

12.5 *The Gershgorin theorem* 342

12.6 *Householder reflectors* 343

13 *Complexity* 344

14 *Partial Differential Equations* 345

14.1 *Partial derivatives* 345

14.2 *Poisson or Laplace Equation* 345

14.3 *Heat Equation* 346

14.4 *Steady state* 346

15 *Taylor series* 347

16 *Graph theory* 350

16.1 *Definitions* 350

16.2 *Common types of graphs* 351

16.3 *Graph colouring and independent sets* 352

16.4 *Graphs and matrices* 352

16.5 *Spectral graph theory* 356

17 *Automata theory* 359

17.1 *Finite State Automata* 359

17.2 *General discussion* 359

18 *Parallel Prefix* 361

18.1 *Parallel prefix* 361

18.2 *Sparse matrix vector product as parallel prefix* 362

### IV Tutorials 365

19 *Unix intro* 367

19.1 *Files and such* 367

19.2 *Text searching and regular expressions* 373

19.3 *Command execution* 375

19.4 *Scripting* 379

19.5 *Expansion* 382

19.6 *Startup files* 384

19.7 *Shell interaction* 385

19.8 *The system and other users* 385

19.9 *The sed and awk tools* 386
<table>
  <tbody>
    <tr>
      <td>19.10</td>
      <td>Review questions</td>
      <td>387</td>
    </tr>
    <tr>
      <td>20</td>
      <td>Compilers and libraries</td>
      <td>388</td>
    </tr>
    <tr>
      <td>20.1</td>
      <td>An introduction to binary files</td>
      <td>388</td>
    </tr>
    <tr>
      <td>20.2</td>
      <td>Simple compilation</td>
      <td>388</td>
    </tr>
    <tr>
      <td>20.3</td>
      <td>Libraries</td>
      <td>390</td>
    </tr>
    <tr>
      <td>21</td>
      <td>Managing projects with Make</td>
      <td>392</td>
    </tr>
    <tr>
      <td>21.1</td>
      <td>A simple example</td>
      <td>392</td>
    </tr>
    <tr>
      <td>21.2</td>
      <td>Variables and template rules</td>
      <td>397</td>
    </tr>
    <tr>
      <td>21.3</td>
      <td>Miscellania</td>
      <td>400</td>
    </tr>
    <tr>
      <td>21.4</td>
      <td>Shell scripting in a Makefile</td>
      <td>402</td>
    </tr>
    <tr>
      <td>21.5</td>
      <td>Practical tips for using Make</td>
      <td>403</td>
    </tr>
    <tr>
      <td>21.6</td>
      <td>A Makefile for LATEX</td>
      <td>404</td>
    </tr>
    <tr>
      <td>22</td>
      <td>Source code control</td>
      <td>406</td>
    </tr>
    <tr>
      <td>22.1</td>
      <td>Workflow in source code control systems</td>
      <td>406</td>
    </tr>
    <tr>
      <td>22.2</td>
      <td>Subversion or SVN</td>
      <td>408</td>
    </tr>
    <tr>
      <td>22.3</td>
      <td>Mercurial or hg</td>
      <td>415</td>
    </tr>
    <tr>
      <td>23</td>
      <td>Scientific Data Storage</td>
      <td>420</td>
    </tr>
    <tr>
      <td>23.1</td>
      <td>Introduction to HDF5</td>
      <td>420</td>
    </tr>
    <tr>
      <td>23.2</td>
      <td>Creating a file</td>
      <td>421</td>
    </tr>
    <tr>
      <td>23.3</td>
      <td>Datasets</td>
      <td>422</td>
    </tr>
    <tr>
      <td>23.4</td>
      <td>Writing the data</td>
      <td>426</td>
    </tr>
    <tr>
      <td>23.5</td>
      <td>Reading</td>
      <td>428</td>
    </tr>
    <tr>
      <td>24</td>
      <td>Scientific Libraries</td>
      <td>430</td>
    </tr>
    <tr>
      <td>24.1</td>
      <td>The Portable Extendable Toolkit for Scientific Computing</td>
      <td>430</td>
    </tr>
    <tr>
      <td>24.2</td>
      <td>Libraries for dense linear algebra: BLAS, Lapack, Scalapack</td>
      <td>440</td>
    </tr>
    <tr>
      <td>25</td>
      <td>Plotting with GNUplot</td>
      <td>447</td>
    </tr>
    <tr>
      <td>25.1</td>
      <td>Usage modes</td>
      <td>447</td>
    </tr>
    <tr>
      <td>25.2</td>
      <td>Plotting</td>
      <td>448</td>
    </tr>
    <tr>
      <td>25.3</td>
      <td>Workflow</td>
      <td>449</td>
    </tr>
    <tr>
      <td>26</td>
      <td>Good coding practices</td>
      <td>450</td>
    </tr>
    <tr>
      <td>26.1</td>
      <td>Defensive programming</td>
      <td>450</td>
    </tr>
    <tr>
      <td>26.2</td>
      <td>Guarding against memory errors</td>
      <td>454</td>
    </tr>
    <tr>
      <td>26.3</td>
      <td>Testing</td>
      <td>457</td>
    </tr>
    <tr>
      <td>27</td>
      <td>Debugging</td>
      <td>458</td>
    </tr>
    <tr>
      <td>27.1</td>
      <td>Invoking gdb</td>
      <td>458</td>
    </tr>
    <tr>
      <td>27.2</td>
      <td>Finding errors</td>
      <td>460</td>
    </tr>
    <tr>
      <td>27.3</td>
      <td>Memory debugging with Valgrind</td>
      <td>463</td>
    </tr>
    <tr>
      <td>27.4</td>
      <td>Stepping through a program</td>
      <td>464</td>
    </tr>
    <tr>
      <td>27.5</td>
      <td>Inspecting values</td>
      <td>466</td>
    </tr>
    <tr>
      <td>27.6</td>
      <td>Breakpoints</td>
      <td>466</td>
    </tr>
    <tr>
      <td>27.7</td>
      <td>Parallel debugging</td>
      <td>467</td>
    </tr>
    <tr>
      <td>27.8</td>
      <td>Further reading</td>
      <td>469</td>
    </tr>
    <tr>
      <td>28</td>
      <td>C/Fortran interoperability</td>
      <td>470</td>
    </tr>
  </tbody>
</table>
28.1 Linker conventions 470

28.2 Arrays 473

28.3 Strings 474

28.4 Subprogram arguments 476

28.5 Input/output 476

28.6 Fortran/C interoperability in Fortran2003 476

29 LaTeX for scientific documentation 477

29.1 The idea behind LATEX, some history of TTEX 477

29.2 A gentle introduction to LaTeX 478

29.3 A worked out example 484

29.4 Where to take it from here 490

V Projects, codes, index 491

30 Class projects 492

30.1 Cache simulation and analysis 492

30.2 Evaluation of Bulk Synchronous Programming 494

30.3 Heat equation 495

30.4 The memory wall 498

31 Codes 499

31.1 Preliminaries 499

31.2 Cache size 500

31.3 Cachelines 502

31.4 Cache associativity 504

31.5 TLB 506

31.6 Unrepresentible numbers 509

32 Bibliography, index, and list of acronyms 510

32.1 Bibliography 510

32.2 List of acronyms 519

32.3 Index 520
Contents
# PART I

## THEORY
# Chapter 1

## Single-processor Computing

In order to write efficient scientific codes, it is important to understand computer architecture. The difference in speed between two codes that compute the same result can range from a few percent to orders of magnitude, depending only on factors relating to how well the algorithms are coded for the processor architecture. Clearly, it is not enough to have an algorithm and 'put it on the computer': some knowledge of computer architecture is advisable, sometimes crucial.

Some problems can be solved on a single CPU, others need a parallel computer that comprises more than one processor. We will go into detail on parallel computers in the next chapter, but even for parallel processing, it is necessary to understand the individual CPUs.

In this chapter, we will focus on what goes on inside a CPU and its memory system. We start with a brief general discussion of how instructions are handled, then we will look into the arithmetic processing in the processor core; last but not least, we will devote much attention to the movement of data between memory and the processor, and inside the processor. This latter point is, maybe unexpectedly, very important, since memory access is typically much slower than executing the processor's instructions, making it the determining factor in a program's performance; the days when 'flop<fnref n="1" /> counting' was the key to predicting a code's performance are long gone. This discrepancy is in fact a growing trend, so the issue of dealing with memory traffic has been becoming more important over time, rather than going away.

This chapter will give you a basic understanding of the issues involved in CPU design, how it affects performance, and how you can code for optimal performance. For much more detail, see an online book about PC architecture [95], and the standard work about computer architecture, Hennesey and Patterson [84].

### 1.1 The Von Neumann architecture

While computers, and most relevantly for this chapter, their processors, can differ in any number of details, they also have many aspects in common. On a very high level of abstraction, many architectures can be described as *von Neumann architectures*. This describes a design with an undivided memory that stores both program and data ('stored program'), and a processing unit that executes the instructions, operating on the data in 'fetch, execute, store' cycle<fnref n="2" />.

<footnote n="1">1. Floating Point Operation.</footnote>

<footnote n="2">2. This model with a prescribed sequence of instructions is also referred to as *control flow*. This is in contrast to *data flow*, which we will see in section 6.12.</footnote>
This setup distinguishes modern processors for the very earliest, and some special purpose contemporary, designs where the program was hard-wired. It also allows programs to modify themselves or generate other programs, since instructions and data are in the same storage. This allows us to have editors and compilers: the computer treats program code as data to operate on<fnref n="3" />. In this book we will not explicitly discuss compilers, the programs that translate high level languages to machine instructions. However, on occasion we will discuss how a program at high level can be written to ensure efficiency at the low level.

In scientific computing, however, we typically do not pay much attention to program code, focusing almost exclusively on data and how it is moved about during program execution. For most practical purposes it is as if program and data are stored separately. The little that is essential about instruction handling can be described as follows.

The machine instructions that a processor executes, as opposed to the higher level languages users write in, typically specify the name of an operation, as well as of the locations of the operands and the result. These locations are not expressed as memory locations, but as *registers*: a small number of named memory locations that are part of the CPU<fnref n="4" />. As an example, here is a simple C routine

```c
void store(double *a, double *b, double *c) {
    *c = *a + *b;
}
```

and its X86 assembler output, obtained by<fnref n="5" /> `gcc -O2 -S -o -store.c`:

```
.text
.p2align 4,,15
.globl store
.type store, @function
store:
    movsd (%rdi), %xmm0 # Load *a to %mma0
    addsd (%rsi), %mma0 # Load *b and add to %mma0
    movsd %mma0, (%rdx) # Store to *c
    ret
```

The instructions here are:

* A load from memory to register;

* Another load, combined with an addition;

* Writing back the result to memory.

Each instruction is processed as follows:

* Instruction fetch: the next instruction according to the *program counter* is loaded into the processor. We will ignore the questions of how and from where this happens.

<footnote n="3">3. At one time, the stored program concept was include as an essential component the ability for a running program to modify its own source. However, it was quickly recognized that this leads to unmaintainable code, and is rarely done in practice [41].</footnote>

<footnote n="4">4. Direct-to-memory architectures are rare, though they have existed. The Cyber 205 supercomputer in the 1980s could have 3 data streams, two from memory to the processor, and one back from the processor to memory, going on at the same time. Such an architecture is only feasible if memory can keep up with the processor speed, which is no longer the case these days.</footnote>

<footnote n="5">5. This is 64-bit output; add the option `-m64` on 32-bit systems.</footnote>
* Instruction decode: the processor inspects the instruction to determine the operation and the operands.

* Memory fetch: if necessary, data is brought from memory into a register.

* Execution: the operation is executed, reading data from registers and writing it back to a register.

* Write-back: for store operations, the register contents is written back to memory.

The case of array data is a little more complicated: the element loaded (or stored) is then determined as the base address of the array plus an offset.

In a way, then, the modern CPU looks to the programmer like a von Neumann machine. There are various ways in which this is not so. For one, while memory looks randomly addressable<fnref n="6" />, in practice there is a concept of *locality*: once a data item has been loaded, nearby items are more efficient to load, and reloading the initial item is also faster.

Another complication to this story of simple loading of data is that contemporary CPUs operate on several instructions simultaneously, which are said to be ‘in flight’, meaning that they are in various stages of completion. Of course, together with these simultaneous instructions, their inputs and outputs are also being moved between memory and processor in an overlapping manner. This is the basic idea of the *superscalar* CPU architecture, and is also referred to as *Instruction Level Parallelism* (ILP). Thus, while each instruction can take several clock cycles to complete, a processor can complete one instruction per cycle in favourable circumstances; in some cases more than one instruction can be finished per cycle.

The main statistic that is quoted about CPUs is their Gigahertz rating, implying that the speed of the processor is the main determining factor of a computer's performance. While speed obviously correlates with performance, the story is more complicated. Some algorithms are *cpu-bound*, and the speed of the processor is indeed the most important factor; other algorithms are *memory-bound*, and aspects such as bus speed and cache size, to be discussed later, become important.

In scientific computing, this second category is in fact quite prominent, so in this chapter we will devote plenty of attention to the process that moves data from memory to the processor, and we will devote relatively little attention to the actual processor.

## 1.2 Modern processors

Modern processors are quite complicated, and in this section we will give a short tour of what their constituent parts. Figure 1.1 is a picture of the *die* of an *Intel Sandy Bridge* processor. This chip is about an inch in size and contains close to a billion transistors.

### 1.2.1 The processing cores

In the Von Neumann model there is a single entity that executes instructions. This has not been the case in increasing measure since the early 2000s. The Sandy Bridge pictured in figure 1.1 has eight cores, each of which is an independent unit executing a stream of instructions. In this chapter we will mostly discuss aspects of a single core; section 1.4 will discuss the integration aspects of the multiple cores.

<footnote n="6">6. There is in fact a theoretical model for computation called the ‘Random Access Machine’; we will briefly see its parallel generalization in section 2.3.2.</footnote>
1.2.1.1 Instruction handling

The Von Neumann model is also unrealistic in
that it assumes that all instructions are exe-
cuted strictly in sequence. Increasingly, over
the last twenty years, processor have used out-
of-order instruction handling, where instruc-
tions can be processed in a different order than
the user program specifies. Of course the pro-
cessor is only allowed to re-order instructions
if that leaves the result of the execution intact!

In the block diagram (figure 1.2) you see var-
ious units that are concerned with instrunc-
tion handling: This cleverness actually costs
considerable energy, as well as sheer amount
of transistors. For this reason, processors such
as the Intel Xeon Phi use in-order instruction
handling.

1.2.1.2 Floating point units

In scientific computing we are mostly inter-
ested in what a processor does with float-
ing point data. Computing with integers or
booleans is typically of less interest. For this
reason, cores have considerable sophistication
for dealing with numerical data.

For instance, while past processors had just a
single Floating Point Unit (FPU), these days
they will have multiple, capable of executing
simultaneously.

For instance, often there are separate addition and multiplication units; if the compiler can find addition and multiplication operations that are independent, it can schedule them so as to be executed simultaneously, thereby doubling the performance of the processor. In some cases, a processor will have multiple addition or multiplication units.

Another way to increase performance is to have a Fused Multiply-Add (FMA) unit, which can execute the instruction $x \leftarrow ax + b$ in the same amount of time as a separate addition or multiplication. Together with pipelining (see below), this means that a processor has an asymptotic speed of several floating point operations per clock cycle.

Incidentally, there are few algorithms in which division operations are a limiting factor. Correspondingly, the division operation is not nearly as much optimized in a modern CPU as the additions and multiplications

<figure><img src="image_3.png" /><figcaption>Figure 1.1: The Intel Sandybridge processor die</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 2-1. Intel microarchitecture code name Sandy Bridge Pipeline Functionality</figcaption></figure>

<caption>Figure 1.2: Block diagram of the Intel Sandy Bridge core</caption>

are. Division operations can take 10 or 20 clock cycles, while a CPU can have multiple addition and/or multiplication units that (asymptotically) can produce a result per cycle.

### 1.2.1.3 Pipelining

The floating point add and multiply units of a processor are pipelined, which has the effect that a stream of independent operations can be performed at an asymptotic speed of one result per clock cycle.

The idea behind a pipeline is as follows. Assume that an operation consists of multiple simpler operations, and that for each suboperation there is separate hardware in the processor. For instance, an addition instruction can have the following components:

* Decoding the instruction, including finding the locations of the operands.

* Copying the operands into registers ('data fetch').

* Aligning the exponents; the addition $.35 \times 10^{-1} + .6 \times 10^{-2}$ becomes $.35 \times 10^{-1} + .06 \times 10^{-1}$.

* Executing the addition of the mantissas, in this case giving .41.
<table><caption>Table 1.1: Floating point capabilities (per core) of several processor architectures, and DAXPY cycle number for 8 operands</caption><thead><tr><th>Processor</th><th>year</th><th>add/mult/fma units (count×width)</th><th>daxpy cycles (arith vs load/store)</th></tr></thead><tbody><tr><td>MIPS R10000</td><td>1996</td><td>1 × 1 + 1 × 1 + 0</td><td>8/24</td></tr><tr><td>Alpha EV5</td><td>1996</td><td>1 × 1 + 1 × 1 + 0</td><td>8/12</td></tr><tr><td>IBM Power5</td><td>2004</td><td>0 + 0 + 2 × 1</td><td>4/12</td></tr><tr><td>AMD Bulldozer</td><td>2011</td><td>2 × 2 + 2 × 2 + 0</td><td>2/4</td></tr><tr><td>Intel Sandy Bridge</td><td>2012</td><td>1 × 4 + 1 × 4 + 0</td><td>2/4</td></tr><tr><td>Intel Haswell</td><td>2014</td><td>0 + 0 + 2 × 4</td><td>1/2</td></tr></tbody></table>

* Normalizing the result, in this example to .41 × 10⁻¹. (Normalization in this example does not do anything. Check for yourself that in .3 × 10⁰ + .8 × 10⁰ and .35 × 10⁻³ + (−.34) × 10⁻³ there is a non-trivial adjustment.)

* Storing the result.

These parts are often called the ‘stages’ or ‘segments’ of the pipeline.

If every component is designed to finish in 1 clock cycle, the whole instruction takes 6 cycles. However, if each has its own hardware, we can execute two operations in less than 12 cycles:

* Execute the decode stage for the first operation;

* Do the data fetch for the first operation, and at the same time the decode for the second.

* Execute the third stage for the first operation and the second stage of the second operation simultaneously.

* Et cetera.

You see that the first operation still takes 6 clock cycles, but the second one is finished a mere 1 cycle later.

Let us make a formal analysis of the speedup you can get from a pipeline. On a traditional FPU, producing *n* results takes $t(n) = n\ell\tau$ where *ℓ* is the number of stages, and *τ* the clock cycle time. The rate at which results are produced is the reciprocal of $t(n)/n$: $r_{\text{serial}} = (\ell\tau)^{-1}$.

On the other hand, for a pipelined FPU the time is $t(n) = [s + \ell + n - 1]\tau$ where *s* is a setup cost: the first operation still has to go through the same stages as before, but after that one more result will be produced each cycle. We can also write this formula as

$$t(n) = [n + n_{1/2}]\tau.$$

**Exercise 1.1.** Let us compare the speed of a classical FPU, and a pipelined one. Show that the result rate is now dependent on *n*: give a formula for $r(n)$, and for $r_∞ = \lim_{n→∞} r(n)$. What is the asymptotic improvement in *r* over the non-pipelined case? Next you can wonder how long it takes to get close to the asymptotic behaviour. Show that for *n* = *n*₁/₂ you get $r(n) = r_∞/2$. This is often used as the definition of *n*₁/₂.

Since a vector processor works on a number of instructions simultaneously, these instructions have to be independent. The operation $\forall_i: a_i \leftarrow b_i + c_i$ has independent additions; the operation $\forall_i: a_{i+1} \leftarrow$
<figure><img src="image_5.png" /><figcaption>Figure 1.3: Schematic depiction of a pipelined operation</figcaption></figure>

$a_i b_i + c_i$ feeds the result of one iteration ($a_i$) to the input of the next ($a_{i+1} = ...$), so the operations are not independent.

A pipelined processor can speed up operations by a factor of 4, 5, 6 with respect to earlier CPUs. Such numbers were typical in the 1980s when the first successful vector computers came on the market. These days, CPUs can have 20-stage pipelines. Does that mean they are incredibly fast? This question is a bit complicated. Chip designers continue to increase the clock rate, and the pipeline segments can no longer finish their work in one cycle, so they are further split up. Sometimes there are even segments in which nothing happens: that time is needed to make sure data can travel to a different part of the chip in time.

The amount of improvement you can get from a pipelined CPU is limited, so in a quest for ever higher performance several variations on the pipeline design have been tried. For instance, the Cyber 205 had separate addition and multiplication pipelines, and it was possible to feed one pipe into the next without data going back to memory first. Operations like $\forall_i: a_i \leftarrow b_i + c_i \cdot d_i$ were called ‘linked triads’ (because of the number of paths to memory, one input operand had to be scalar).

**Exercise 1.2.** Analyse the speedup and $n_{1/2}$ of linked triads.

Another way to increase performance is to have multiple identical pipes. This design was perfected by the NEC SX series. With, for instance, 4 pipes, the operation $\forall_i: a_i \leftarrow b_i + c_i$ would be split module 4, so that the first pipe operated on indices $i = 4 \cdot j$, the second on $i = 4 \cdot j + 1$, etcetera.

**Exercise 1.3.** Analyze the speedup and $n_{1/2}$ of a processor with multiple pipelines that operate in parallel. That is, suppose that there are $p$ independent pipelines, executing the same instruction, that can each handle a stream of operands.
(You may wonder why we are mentioning some fairly old computers here: true pipeline supercomputers hardly exist anymore. In the US, the Cray X1 was the last of that line, and in Japan only NEC still makes them. However, the functional units of a CPU these days are pipelined, so the notion is still important.)

Exercise 1.4. The operation

```c
for (i) {
    x[i+1] = a[i] * x[i] + b[i];
}
```

can not be handled by a pipeline because there is a dependency between input of one iteration of the operation and the output of the previous. However, you can transform the loop into one that is mathematically equivalent, and potentially more efficient to compute. Derive an expression that computes $x[i+2]$ from $x[i]$ without involving $x[i+1]$. This is known as recursive doubling. Assume you have plenty of temporary storage. You can now perform the calculation by

* Doing some preliminary calculations;

* computing $x[i], x[i+2], x[i+4], \dots$, and from these,

* compute the missing terms $x[i+1], x[i+3], \dots$.

Analyze the efficiency of this scheme by giving formulas for $T_0(n)$ and $T_s(n)$. Can you think of an argument why the preliminary calculations may be of lesser importance in some circumstances?

### 1.2.1.4 Peak performance

Thanks to pipelining, for modern CPUs there is a simple relation between the clock speed and the peak performance. Since each FPU can produce one result per cycle asymptotically, the peak performance is the clock speed times the number of independent FPUs. The measure of floating point performance is 'floating point operations per second', abbreviated flops. Considering the speed of computers these days, you will mostly hear floating point performance being expressed in 'gigaflops': multiples of $10^9$ flops.

## 1.2.2 8-bit, 16-bit, 32-bit, 64-bit

Processors are often characterized in terms of how big a chunk of data they can process as a unit. This can relate to

* The width of the path between processor and memory: can a 64-bit floating point number be loaded in one cycle, or does it arrive in pieces at the processor.

* The way memory is addressed: if addresses are limited to 16 bits, only 64,000 bytes can be identified. Early PCs had a complicated scheme with segments to get around this limitation: an address was specified with a segment number and an offset inside the segment.

* The number of bits in a register, in particular the size of the integer registers which manipulate data address; see the previous point. (Floating point register are often larger, for instance 80 bits in the x86 architecture.) This also corresponds to the size of a chunk of data that a processor can operate on simultaneously.
* The size of a floating point number. If the arithmetic unit of a CPU is designed to multiply 8-byte numbers efficiently ('double precision'; see section 3.2.2) then numbers half that size ('single precision') can sometimes be processed at higher efficiency, and for larger numbers ('quadruple precision') some complicated scheme is needed. For instance, a quad precision number could be emulated by two double precision numbers with a fixed difference between the exponents.

These measurements are not necessarily identical. For instance, the original Pentium processor had 64-bit data busses, but a 32-bit processor. On the other hand, the Motorola 68000 processor (of the original Apple Macintosh) had a 32-bit CPU, but 16-bit data busses.

The first Intel microprocessor, the 4004, was a 4-bit processor in the sense that it processed 4 bit chunks. These days, 64 bit processors are becoming the norm.

### 1.2.3 Caches: on-chip memory

The bulk of computer memory is in chips that are separate from the processor. However, there is usually a small amount (typically a few megabytes) of on-chip memory, called the *cache*. This will be explained in detail in section 1.3.4.

### 1.2.4 Graphics, controllers, special purpose hardware

One difference between 'consumer' and 'server' type processors is that the consumer chips devote considerable real-estate on the processor chip to graphics. Processors for cell phones and tablets can even have dedicated circuitry for security or mp3 playback. Other parts of the processor are dedicated to communicating with memory or the *I/O subsystem*. We will not discuss those aspects in this book.

### 1.2.5 Superscalar processing and instruction-level parallelism

In the von Neumann model processors operate through *control flow*: instructions follow each other linearly or with branches without regard for what data they involve. As processors became more powerful and capable of executing more than one instruction at a time, it became necessary to switch to the *data flow* model. Such *superscalar* processors analyze several instructions to find data dependencies, and execute instructions in parallel that do not depend on each other.

This concept is also known as *Instruction Level Parallelism* (ILP), and it is facilitated by various mechanisms:

* multiple-issue: instructions that are independent can be started at the same time;

* pipelining: already mentioned, arithmetic units can deal with multiple operations in various stages of completion;

* branch prediction and speculative execution: a compiler can 'guess' whether a conditional instruction will evaluate to true, and execute those instructions accordingly;

* out-of-order execution: instructions can be rearranged if they are not dependent on each other, and if the resulting execution will be more efficient;

* prefetching: data can be speculatively requested before any instruction needing it is actually encountered (this is discussed further in section 1.3.5).
Above, you saw pipelining in the context of floating point operations. Nowadays, the whole CPU is pipelined. Not only floating point operations, but any sort of instruction will be put in the *instruction pipeline* as soon as possible. Note that this pipeline is no longer limited to identical instructions: the notion of pipeline is now generalized to any stream of partially executed instructions that are simultaneously “in flight”.

As clock frequency has gone up, the processor pipeline has grown in length to make the segments executable in less time. You have already seen that longer pipelines have a larger $n_{1/2}$, so more independent instructions are needed to make the pipeline run at full efficiency. As the limits to instruction-level parallelism are reached, making pipelines longer (sometimes called ‘deeper’) no longer pays off. This is generally seen as the reason that chip designers have moved to *multicore* architectures as a way of more efficiently using the transistors on a chip; section 1.4.

There is a second problem with these longer pipelines: if the code comes to a branch point (a conditional or the test in a loop), it is not clear what the next instruction to execute is. At that point the pipeline can *stall*. CPUs have taken to *speculative execution* for instance, by always assuming that the test will turn out true. If the code then takes the other branch (this is called a *branch misprediction*), the pipeline has to be *flushed* and restarted. The resulting delay in the execution stream is called the *branch penalty*.

## 1.3 Memory Hierarchies

We will now refine the picture of the Von Neuman architecture, in which data is loaded immediately from memory to the processors, where it is operated on. This picture is unrealistic because of the so-called *memory wall* [166]: the memory is too slow to load data into the process at the rate the processor can absorb it. Specifically, a single load can take 1000 cycles, while a processor can perform several operations per cycle. (After this long wait for a load, the next load can come faster, but still too slow for the processor. This matter of wait time versus throughput will be addressed below in section 1.3.2.)

In reality, there will be various memory levels in between the FPU and the main memory: the registers and the caches, together called the *memory hierarchy*. These try to alleviate the memory wall problem by making recently used data available quicker than it would be from main memory. Of course, this presupposes that the algorithm and its implementation allow for data to be used multiple times. Such questions of *data reuse* will be discussed in more detail in section 1.5.1.

Both registers and caches are faster than main memory to various degrees; unfortunately, the faster the memory on a certain level, the smaller it will be. These differences in size and access speed lead to interesting programming problems, which we will discuss later in this chapter, and particularly section 1.6.

We will now discuss the various components of the memory hierarchy and the theoretical concepts needed to analyze their behaviour.

### 1.3.1 Busses

The wires that move data around in a computer, from memory to cpu or to a disc controller or screen, are called *busses*. The most important one for us is the *Front-Side Bus (FSB)* which connects the processor
to memory. In one popular architecture, this is called the ‘north bridge’, as opposed to the ‘south bridge’
which connects to external devices, with the exception of the graphics controller.

<figure><img src="image_6.png" /><figcaption>Figure 1.4: Bus structure of a processor</figcaption></figure>

The bus is typically slower than the processor, operating with clock frequencies slightly in excess of 1GHz, which is a fraction of the CPU clock frequency. This is one reason that caches are needed; the fact that a processors can consume many data items per clock tick contributes to this. Apart from the frequency, the bandwidth of a bus is also determined by the number of bits that can be moved per clock cycle. This is typically 64 or 128 in current architectures. We will now discuss this in some more detail.

**1.3.2 Latency and Bandwidth**

Above, we mentioned in very general terms that accessing data in registers is almost instantaneous, whereas
loading data from memory into the registers, a necessary step before any operation, incurs a substantial
delay. We will now make this story slightly more precise.

There are two important concepts to describe the movement of data: *latency* and *bandwidth*. The assumption here is that requesting an item of data incurs an initial delay; if this item was the first in a stream of data, usually a consecutive range of memory addresses, the remainder of the stream will arrive with no further delay at a regular amount per time period.

**Latency** is the delay between the processor issuing a request for a memory item, and the item actually arriving. We can distinguish between various latencies, such as the transfer from memory to cache, cache to register, or summarize them all into the latency between memory and processor. Latency is measured in (nano) seconds, or clock periods.

If a processor executes instructions in the order they are found in the assembly code, then execu-
tion will often *stall* while data is being fetched from memory; this is also called *memory stall*.
For this reason, a low latency is very important. In practice, many processors have ‘out-of-order
execution’ of instructions, allowing them to perform other operations while waiting for the re-
quested data. Programmers can take this into account, and code in a way that achieves *latency*
hiding; see also section 1.5.1. Graphics Processing Units (GPUs) (see section 2.10.3) can switch very quickly between threads in order to achieve latency hiding.

**Bandwidth** is the rate at which data arrives at its destination, after the initial latency is overcome. Bandwidth is measured in bytes (kilobytes, megabytes, gigabytes) per second or per clock cycle. The bandwidth between two memory levels is usually the product of the cycle speed of the channel (the bus speed) and the bus width: the number of bits that can be sent simultaneously in every cycle of the bus clock.

The concepts of latency and bandwidth are often combined in a formula for the time that a message takes from start to finish:

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

where $\alpha$ is the latency and $\beta$ is the inverse of the bandwidth: the time per byte.

Typically, the further away from the processor one gets, the longer the latency is, and the lower the bandwidth. These two factors make it important to program in such a way that, if at all possible, the processor uses data from cache or register, rather than from main memory. To illustrate that this is a serious matter, consider a vector addition

$$\begin{array}{l} \texttt{for } (i) \\ a[i] = b[i] + c[i] \end{array}$$

Each iteration performs one floating point operation, which modern CPUs can do in one clock cycle by using pipelines. However, each iteration needs two numbers loaded and one written, for a total of 24 bytes<fnref n="7" /> of memory traffic. Typical memory bandwidth figures (see for instance figure 1.5) are nowhere near 24 (or 32) bytes per cycle. This means that, without caches, algorithm performance can be bounded by memory performance. Of course, caches will not speed up every operations, and in fact will have no effect on the above example. Strategies for programming that lead to significant cache use are discussed in section 1.6.

The concepts of latency and bandwidth will also appear in parallel computers, when we talk about sending data from one processor to the next.

### 1.3.3 Registers

Every processor has a small amount of memory that is internal to the processor: the *registers*, or together the *register file*. The registers are what the processor actually operates on: an operation such as

$$a := b + c$$

is actually implemented as

* load the value of $b$ from memory into a register,

* load the value of $c$ from memory into another register,

* compute the sum and write that into yet another register, and

* write the sum value back to the memory location of $a$.

<footnote n="7">7. Actually, $a[i]$ is loaded before it can be written, so there are 4 memory access, with a total of 32 bytes, per iteration.</footnote>
Looking at assembly code (for instance the output of a compiler), you see the explicit load, compute, and store instructions.

Compute instructions such as add or multiply only operate on registers. For instance, in assembly language you will see instructions such as

```
addl %eax, %edx
```

which adds the content of one register to another. As you see in this sample instruction, registers are not numbered, as opposed to memory addresses, but have distinct names that are referred to in the assembly instruction. Typically, a processor has 16 or 32 floating point registers; the *Intel Itanium* was exceptional with 128 floating point registers.

Registers have a high bandwidth and low latency because they are part of the processor. You can consider data movement to and from registers as essentially instantaneous.

In this chapter you will see stressed that moving data from memory is relatively expensive. Therefore, it would be a simple optimization to leave data in register when possible. For instance, if the above computation is followed by a statement

```
a := b + c
```

```
d := a + e
```

the computed value of a could be left in register. This optimization is typically performed as a *compiler optimization*: the compiler will simply not generate the instructions for storing and reloading a. We say that a stays *resident in register*.

Keeping values in register is often done to avoid recomputing a quantity. For instance, in

```
t1 = sin(alpha) * x + cos(alpha) * y;
```

```
t2 = -cos(alsph) * x + sin(alpha) * y;
```

the sine and cosine quantity will probably be kept in register. You can help the compiler by explicitly introducing temporary quantities:

```
s = sin(alpha); c = cos(alpha);
```

```
t1 = s * x + c * y;
```

```
t2 = -c * x + s * y
```

Of course, there is a limit to how many quantities can be kept in register; trying to keep too many quantities in register is called *register spill* and lowers the performance of a code.

Keeping a variable in register is especially important if that variable appears in an inner loop. In the computation

```
for i=1,length
```

```
a[i] = b[i] * c
```

the quantity c will probably be kept in register by the compiler, but in
for k=1,nvectors
for i=1,length
a[i,k] = b[i,k] * c[k]

it is a good idea to introduce explicitly a temporary variable to hold `c[k]`. In C, you can give a hint to the compiler to keep a variable in register by declaring it as a *register variable*:

register double t;

### 1.3.4 Caches

In between the registers, which contain the immediate input and output data for instructions, and the main memory where lots of data can reside for a long time, are various levels of *cache memory*, that have lower latency and higher bandwidth than main memory and where data are kept for an intermediate amount of time. Caches typically consist of *Static Random-Access Memory (SRAM)*, which is faster than then *Dynamic Random-Access Memory (DRAM)* used for the main memory, but is also more expensive.

Data from memory travels through the caches to wind up in registers. The advantage to having cache memory is that if a data item is reused shortly after it was first needed, it will still be in cache, and therefore it can be accessed much faster than if it would have to be brought in from memory.

#### 1.3.4.1 A motivating example

As an example, let's suppose a variable `x` is used twice, and its uses are too far apart that it would stay resident in register:

... = ... x ..... // instruction using x
......... // several instructions not involving x
... = ... x ..... // instruction using x

The assembly code would then be

* load `x` from memory into register; operate on it;

* do the intervening instructions;

* load `x` from memory into register; operate on it;

With a cache, the assembly code stays the same, but the actual behaviour of the memory system now becomes:

* load `x` from memory into cache, and from cache into register; operate on it;

* do the intervening instructions;

* request `x` from memory, but since it is still in the cache, load it from the cache into register; operate on it.

Since loading from cache is faster than loading from main memory, the computation will now be faster. Caches are fairly small, so values can not be kept there indefinitely. We will see the implications of this in the following discussion.
There is an important difference between cache memory and registers: while data is moved into register by explicit assembly instructions, the move from main memory to cache is entirely done by hardware. Thus cache use and reuse is outside of direct programmer control. Later, especially in sections 1.5.2 and 1.6, you will see how it is possible to influence cache use indirectly.

1.3.4.2 Cache levels, speed and size

The caches are called ‘level 1’ and ‘level 2’ (or, for short, L1 and L2) cache; some processors can have an L3 cache. The L1 and L2 caches are part of the die, the processor chip, although for the L2 cache that is a relatively recent development; the L3 cache is off-chip. The L1 cache is small, typically around 16Kbyte. Level 2 (and, when present, level 3) cache is more plentiful, up to several megabytes, but it is also slower. Unlike main memory, which is expandable, caches are fixed in size. If a version of a processor chip exists with a larger cache, it is usually considerably more expensive.

Data needed in some operation gets copied into the various caches on its way to the processor. If, some instructions later, a data item is needed again, it is first searched for in the L1 cache; if it is not found there, it is searched for in the L2 cache; if it is not found there, it is loaded from main memory. Finding data in cache is called a *cache hit*, and not finding it a *cache miss*.

Figure 1.5 illustrates the basic facts of the *cache hierarchy*, in this case for the *Intel Sandy Bridge* chip: the closer caches are to the FPUs, the faster, but also the smaller they are. Some points about this figure.

<figure><img src="image_7.png" /><figcaption>Figure 1.5: Memory hierarchy of an Intel Sandy Bridge, characterized by speed and size.</figcaption></figure>

* Loading data from registers is so fast that it does not constitute a limitation on algorithm execution speed. On the other hand, there are few registers. Each core has 16 general purpose registers, and 16 SIMD registers.

* The L1 cache is small, but sustains a bandwidth of 32 bytes, that is 4 double precision number, per cycle. This is enough to load two operands each for two operations, but note that the core can actually perform 4 operations per cycle. Thus, to achieve peak speed, certain operands need to stay in register: typically, L1 bandwidth is enough for about half of peak performance.
* The bandwidth of the L2 and L3 cache is nominally the same as of L1. However, this bandwidth is partly wasted on coherence issues.

* Main memory access has a latency of more than 100 cycles, and a bandwidth of 4.5 bytes per cycle, which is about 1/7th of the L1 bandwidth. However, this bandwidth is shared by the multiple cores of a processor chip, so effectively the bandwidth is a fraction of this number. Most clusters will also have more than one *socket* (processor chip) per node, typically 2 or 4, so some bandwidth is spent on maintaining *cache coherence* (see section 1.4), again reducing the bandwidth available for each chip.

On level 1, there are separate caches for instructions and data; the L2 and L3 cache contain both data and instructions.

You see that the larger caches are increasingly unable to supply data to the processors fast enough. For this reason it is necessary to code in such a way that data is kept as much as possible in the highest cache level possible. We will discuss this issue in detail in the rest of this chapter.

**Exercise 1.5.** The L1 cache is smaller than the L2 cache, and if there is an L3, the L2 is smaller than the L3. Give a practical and a theoretical reason why this is so.

### 1.3.4.3 Types of cache misses

There are three types of cache misses.

As you saw in the example above, the first time you reference data you will always incur a cache miss. This is known as a *compulsory cache miss* since these are unavoidable. Does that mean that you will always be waiting for a data item, the first time you need it? Not necessarily: section 1.3.5 explains how the hardware tries to help you by predicting what data is needed next.

The next type of cache misses is due to the size of your working set: a *capacity cache miss* is caused by data having been overwritten because the cache can simply not contain all your problem data. (Section 1.3.4.5 discusses how the processor decides what data to overwrite.) If you want to avoid this type of misses, you need to partition your problem in chunks that are small enough that data can stay in cache for an appreciable time. Of course, this presumes that data items are operated on multiple times, so that there is actually a point in keeping it in cache; this is discussed in section 1.5.1.

Finally, there are *conflict misses* caused by one data item being mapped to the same cache location as another, while both are still needed for the computation, and there would have been better candidates to evict. This is discussed in section 1.3.4.9.

In a *multicore* context there is a further type of cache miss: the *invalidation miss*. This happens if an item in cache has become invalid because another core changed the value of the corresponding memory address. The core will then have to reload this address.

### 1.3.4.4 Reuse is the name of the game

The presence of one or more caches is not immediately a guarantee for high performance: this largely depends on the *memory access pattern* of the code, and how well this exploits the caches. The first time that an item is referenced, it is copied from memory into cache, and through to the processor registers. The
latency and bandwidth for this are not mitigated in any way by the presence of a cache. When the same item is referenced a second time, it may be found in cache, at a considerably reduced cost in terms of latency and bandwidth: caches have shorter latency and higher bandwidth than main memory.

We conclude that, first, an algorithm has to have an opportunity for data reuse. If every data item is used only once (as in addition of two vectors), there can be no reuse, and the presence of caches is largely irrelevant. A code will only benefit from the increased bandwidth and reduced latency of a cache if items in cache are referenced more than once; see section 1.5.1 for a detailed discussion.. An example would be the matrix-vector multiplication $y = Ax$ where each element of $x$ is used in $n$ operations, where $n$ is the matrix dimension.

Secondly, an algorithm may theoretically have an opportunity for reuse, but it needs to be coded in such a way that the reuse is actually exposed. We will address these points in section 1.5.2. This second point especially is not trivial.

Some problems are small enough that they fit completely in cache, at least in the L3 cache. This is something to watch out for when benchmarking, since it gives a too rosy picture of processor performance.

#### 1.3.4.5 Replacement policies

Data in cache and registers is placed there by the system, outside of programmer control. Likewise, the system decides when to overwrite data in the cache or in registers if it is not referenced in a while, and as other data needs to be placed there. Below, we will go into detail on how caches do this, but as a general principle, a Least Recently Used (LRU) cache replacement policy is used: if a cache is full and new data needs to be placed into it, the data that was least recently used is *flushed*, meaning that it is overwritten with the new item, and therefore no longer accessible. LRU is by far the most common replacement policy; other possibilities are FIFO (first in first out) or random replacement.

**Exercise 1.6.** Sketch a simple scenario, and give some (pseudo) code, to argue that LRU is preferable over FIFO as a replacement strategy.

#### 1.3.4.6 Cache lines

Data movement between memory and cache, or between caches, is not done in single bytes, or even words. Instead, the smallest unit of data moved is called a *cache line*, sometimes called a *cache block*. A typical cache line is 64 or 128 bytes long, which in the context of scientific computing implies 8 or 16 double precision floating point numbers. The cache line size for data moved into L2 cache can be larger than for data moved into L1 cache.

It is important to acknowledge the existence of cache lines in coding, since any memory access costs the transfer of several words (see section 1.6.4 for some examples). An efficient program then tries to use the other items on the cache line, since access to them is effectively free. This phenomenon is visible in code that accesses arrays by *stride*: elements are read or written at regular intervals.

<figure><img src="image_8.png" /><figcaption>Figure 1.6: Accessing 4 elements at stride 1</figcaption></figure>
Stride 1 corresponds to sequential access of an array:

```
for (i=0; i<N; i++)
    ... = ... x[i] ...
```

Let us use as illustration a case with 4 words per cacheline. Requesting the first elements loads the whole cacheline that contains it into cache. A request for the 2nd, 3rd, and 4th element can then be satisfied from cache, meaning with high bandwidth and low latency.

A larger stride

```
for (i=0; i<N; i+=stride)
    ... = ... x[i] ...
```

<figure><img src="image_9.png" /><figcaption>Figure 1.7: Accessing 4 elements at stride 3</figcaption></figure>

implies that in every cache line only certain elements are used. We illustrate that with stride 3: requesting the first elements loads a cacheline, and this cacheline also contains the second element. However, the third element is on the next cacheline, so loading this incurs the latency and bandwidth of main memory. The same holds for the fourth element. Loading four elements now needed loading three cache lines instead of one, meaning that two-thirds of the available bandwidth has been wasted. (This second case would also incur three times the latency of the first, if it weren't for a hardware mechanism that notices the regular access patterns, and pre-emptively loads further cachelines; see section 1.3.5.)

Some applications naturally lead to strides greater than 1, for instance, accessing only the real parts of an array of complex numbers (for some remarks on the practical realization of complex numbers see section 3.5.5). Also, methods that use *recursive doubling* often have a code structure that exhibits non-unit strides

```
for (i=0; i<N/2; i++)
    x[i] = y[2*i];
```

In this discussion of cachelines, we have implicitly assumed the beginning of a cacheline is also the beginning of a word, be that an integer or a floating point number. This need not be true: an 8-byte floating point number can be placed straddling the boundary between two cachelines. You can image that this is not good for performance. Section 31.1.3 discusses ways to address *cacheline boundary alignment* in practice.

### 1.3.4.7 Cache mapping

Caches get faster, but also smaller, the closer to the FPUs they get, yet even the largest cache is considerably smaller than the main memory size. We already noted that this has implications for the cache replacement strategy. Another issue we need to address in this context is that of *cache mapping*, which is the question of ‘if an item is placed in cache, where does it get placed’. This problem is generally addressed by mapping the (main memory) address of the item to an address in cache, leading to the question ‘what if two items get mapped to the same address’.
<figure><img src="image_10.png" /><figcaption>Figure 1.8: Direct mapping of 32-bit addresses into a 64K cache</figcaption></figure>

1.3.4.8 Direct mapped caches

The simplest cache mapping strategy is direct mapping. Suppose that memory addresses are 32 bits long, so that they can address 4G bytes<fnref n="8" />; suppose further that the cache has 8K words, that is, 64K bytes, needing 16 bits to address. Direct mapping then takes from each memory address the last ('least significant') 16 bits, and uses these as the address of the data item in cache; see figure 1.8.

Direct mapping is very efficient because of its address calculations can be performed very quickly, leading to low latency, but it has a problem in practical applications. If two items are addressed that are separated by 8K words, they will be mapped to the same cache location, which will make certain calculations inefficient.

Example:

```c
double A[3][8192];
for (i=0; i<512; i++)
    a[2][i] = ( a[0][i] + a[1][i] ) / 2.;
```

or in Fortran:

```fortran
real*8 A(8192,3);
do i=1,512
    a(i,3) = ( a(i,1)+a(i,2) )/2
end do
```

Here, the locations of `a[0][i]`, `a[1][i]`, and `a[2][i]` (or `a(i,1)`, `a(i,2)`, `a(i,3)`)) are 8K from each other for every `i`, so the last 16 bits of their addresses will be the same, and hence they will be mapped to the same location in cache. The execution of the loop will now go as follows:

* The data at `a[0][0]` is brought into cache and register. This engenders a certain amount of latency. Together with this element, a whole cache line is transferred.

* The data at `a[1][0]` is brought into cache (and register, as we will not remark anymore from now on), together with its whole cache line, at cost of some latency. Since this cache line is mapped to the same location as the first, the first cache line is overwritten.

<footnote n="8">8. We implicitly use the convention that K,M,G suffixes refer to powers of 2 rather than 10: 1K=1024, 1M=1,048,576, 1G=1,073,741,824.</footnote>
* In order to write the output, the cache line containing `a[2][0]` is brought into memory. This is again mapped to the same location, causing flushing of the cache line just loaded for `a[1][0]`.

* In the next iteration, `a[0][1]` is needed, which is on the same cache line as `a[0][0]`. However, this cache line has been flushed, so it needs to be brought in anew from main memory or a deeper cache level. In doing so, it overwrites the cache line that holds `a[2][0]`.

* A similar story hold for `a[1][1]`: it is on the cache line of `a[1][0]`, which unfortunately has been overwritten in the previous step.

If a cache line holds four words, we see that each four iterations of the loop involve eight transfers of elements of `a`, where two would have sufficed, if it were not for the cache conflicts.

**Exercise 1.7.** In the example of direct mapped caches, mapping from memory to cache was done by using the final 16 bits of a 32 bit memory address as cache address. Show that the problems in this example go away if the mapping is done by using the first ('most significant') 16 bits as the cache address. Why is this not a good solution in general?

### 1.3.4.9 Associative caches

The problem of cache conflicts, outlined in the previous section, would be solved if any data item could go to any cache location. In that case there would be no conflicts, other than the cache filling up, in which case a cache replacement policy (section 1.3.4.5) would flush data to make room for the incoming item. Such a cache is called *fully associative*, and while it seems optimal, it is also very costly to build, and much slower in use than a direct mapped cache.

For this reason, the most common solution is to have a *k*-way *associative cache*, where *k* is at least two. In this case, a data item can go to any of *k* cache locations. Code would have to have a *k* + 1-way conflict before data would be flushed prematurely as in the above example. In that example, a value of *k* = 2 would suffice, but in practice higher values are often encountered. Figure 1.9 illustrates the mapping of

<figure><img src="image_11.png" /><figcaption>Figure 1.9: Two caches of 12 elements: direct mapped (left) and 3-way associative (right)</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 1.10: Prefetch stream generated by equally spaced requests</figcaption></figure>

memory addresses to cache locations for a direct mapped and a 3-way associative cache. Both caches has 12 elements, but addressed diferently. The direct mapped cache (left) will have a conflict between memory address 0 and 12, but in the 3-way associative cache these two addresses can be mapped to any of three elements.

As a practical example, the *Intel Woodcrest* processor has an L1 cache of 32K bytes that is 8-way set associative with a 64 byte cache line size, and an L2 cache of 4M bytes that is 8-way set associative with a 64 byte cache line size. On the other hand, the *AMD Barcelona* chip has 2-way associativity for the L1 cache, and 8-way for the L2. A higher associativity ('way-ness') is obviously desirable, but makes a processor slower, since determining whether an address is already in cache becomes more complicated. For this reason, the associativity of the L1 cache, where speed is of the greatest importance, is typically lower than of the L2.

Exercise 1.8. Write a small cache simulator in your favourite language. Assume a *k*-way as- sociative cache of 32 entries and an architecture with 16 bit addresses. Run the fol- lowing experiment for *k* = 1, 2, 4, ...:

1. Let *k* be the associativity of the simulated cache.

2. Write the translation from 16 bit memory addresses to 32/*k* cache addresses.

3. Generate 32 random machine addresses, and simulate storing them in cache.

Since the cache has 32 entries, optimally the 32 addresses can all be stored in cache. The chance of this actually happening is small, and often the data of one address will be evicted from the cache (meaning that it is overwritten) when another address conflicts with it. Record how many addresses, out of 32, are actually stored in the cache at the end of the simulation. Do step 3 100 times, and plot the results; give median and average value, and the standard deviation. Observe that increasing the associativity improves the number of addresses stored. What is the limit behaviour? (For bonus points, do a formal statistical analysis.)

1.3.5 Prefetch streams

In the traditional von Neumann model (section 1.1), each instruction contains the location of its operands,
so a CPU implementing this model would make a separate request for each new operand. In practice, often
subsequent data items are adjacent or regularly spaced in memory. The memory system can try to detect
such data patterns by looking at cache miss points, and request a prefetch data stream (figure 1.10).
In its simplest form, the CPU will detect that consecutive loads come from two consecutive cache lines, and automatically issue a request for the next following cache line. This process can be repeated or extended if the code makes an actual request for that third cache line. Since these cache lines are now brought from memory well before they are needed, prefetch has the possibility of eliminating the latency for all but the first couple of data items.

The concept of cache miss now needs to be revisited a little. From a performance point of view we are only interested in stalls on cache misses, that is, the case where the computation has to wait for the data to be brought in. Data that is not in cache, but can be brought in while other instructions are still being processed, is not a problem. If an 'L1 miss' is understood to be only a 'stall on miss', then the term 'L1 cache refill' is used to describe all cacheline loads, whether the processor is stalling on them or not.

Since prefetch is controlled by the hardware, it is also described as hardware prefetch. Prefetch streams can sometimes be controlled from software, though often it takes assembly code, for instance inline assembly instructions in your C code, to do so. On the IBM BlueGene Q prefetch can be activated by program directives.

### 1.3.6 Concurrency and memory transfer

In the discussion about the memory hierarchy we made the point that memory is slower than the processor. As if that is not bad enough, it is not even trivial to exploit all the bandwidth that memory offers. In other words, if you don't program carefully you will get even less performance than you would expect based on the available bandwidth. Let's analyze this.

The memory system typically has a bandwidth of more than one floating point number per cycle, so you need to issue that many requests per cycle to utilize the available bandwidth. This would be true even with zero latency; since there is latency, it takes a while for data to make it from memory and be processed. Consequently, any data requested based on computations on the first data has to be requested with a delay at least equal the memory latency.

For full utilization of the bandwidth, at all times a volume of data equal to the bandwidth times the latency has to be in flight. Since these data have to be independent, we get a statement of Little's law [116]:

$$ \text{Concurrency} = \text{Bandwidth} \times \text{Latency}. $$

This is illustrated in figure 1.11. The problem with maintaining this concurrency is not that a program does not have it; rather, the program is to get the compiler and runtime system recognize it. For instance, if a loop traverses a long array, the compiler will not issue a large number of memory requests. The prefetch mechanism (section 1.3.5) will issue some memory requests ahead of time, but typically not enough. Thus, in order to use the available bandwidth, multiple streams of data need to be under way simultaneously. Therefore, we can also phrase Little's law as

$$ \text{Effective throughput} = \text{Expressed concurrency}/\text{Latency}. $$

### 1.3.7 Memory banks

Above, we discussed issues relating to bandwidth. You saw that memory, and to a lesser extent caches, have a bandwidth that is less than what a processor can maximally absorb. The situation is actually even worse
<figure><img src="image_13.png" /><figcaption>Figure 1.11: Illustration of Little's Law that states how much independent data needs to be in flight</figcaption></figure>

than the above discussion made it seem. For this reason, memory is often divided into *memory banks* that are interleaved: with four memory banks, words 0, 4, 8, ... are in bank 0, words 1, 5, 9, ... are in bank 1, etcetera.

Suppose we now access memory sequentially, then such 4-way interleaved memory can sustain four times the bandwidth of a single memory bank. Unfortunately, accessing by stride 2 will halve the bandwidth, and larger strides are even worse. In practice the number of memory banks will be higher, so that strided memory access with small strides will still have the full advertised bandwidth.

This concept of banks can also apply to caches. For instance, the cache lines in the L1 cache of the *AMD Barcelona* chip are 16 words long, divided into two interleaved banks of 8 words. This means that sequential access to the elements of a cache line is efficient, but strided access suffers from a deteriorated performance.

### 1.3.8 TLB and virtual memory

All of a program's data may not be in memory simultaneously. This can happen for a number of reasons:

* The computer serves multiple users, so the memory is not dedicated to any one user;

* The computer is running multiple programs, which together need more than the physically available memory;

* One single program can use more data than the available memory.

For this reason, computers use *Virtual memory*: if more memory is needed than is available, certain blocks of memory are written to disc. In effect, the disc acts as an extension of the real memory. This means that a block of data can be anywhere in memory, and in fact, if it is *swapped* in and out, it can be in different locations at different times. Swapping does not act on individual memory locations, but rather on *memory pages*: contiguous blocks of memory, from a few kilobytes to megabytes in size. (In an earlier generation of operating systems, moving memory to disc was a programmer's responsibility. Pages that would replace each other were called *overlays*.)
For this reason, we need a translation mechanism from the memory addresses that the program uses to the actual addresses in memory, and this translation has to be dynamic. A program has a ‘logical data space’ (typically starting from address zero) of the addresses used in the compiled code, and this needs to be translated during program execution to actual memory addresses. For this reason, there is a *page table* that specifies which memory pages contain which logical pages.

However, address translation by lookup in this table is slow, so CPUs have a *Translation Look-aside Buffer (TLB)*. The TLB is a cache of frequently used Page Table Entries: it provides fast address translation for a number of pages. If a program needs a memory location, the TLB is consulted to see whether this location is in fact on a page that is remembered in the TLB. If this is the case, the logical address is translated to a physical one; this is a very fast process. The case where the page is not remembered in the TLB is called a *TLB miss*, and the page lookup table is then consulted, if necessary bringing the needed page into memory. The TLB is (sometimes fully) associative (section 1.3.4.9), using an LRU policy (section 1.3.4.5).

A typical TLB has between 64 and 512 entries. If a program accesses data sequentially, it will typically alternate between just a few pages, and there will be no TLB misses. On the other hand, a program that access many random memory locations can experience a slowdown because of such misses.

Section 1.6.5 and appendix 31.5 discuss some simple code illustrating the behaviour of the TLB.

[There are some complications to this story. For instance, there is usually more than one TLB. The first one is associated with the L2 cache, the second one with the L1. In the *AMD Opteron*, the L1 TLB has 48 entries, and is fully (48-way) associative, while the L2 TLB has 512 entries, but is only 4-way associative. This means that there can actually be TLB conflicts. In the discussion above, we have only talked about the L2 TLB. The reason that this can be associated with the L2 cache, rather than with main memory, is that the translation from memory to L2 cache is deterministic.]

## 1.4 Multicore architectures

In recent years, the limits of performance have been reached for the traditional processor chip design.

* Clock frequency can not be increased further, since it increases energy consumption, heating the chips too much; see section 1.7.

* It is not possible to extract more Instruction Level Parallelism (ILP) from codes, either because of compiler limitations, because of the limited amount of intrinsically available parallelism, or because branch prediction makes it impossible (see section 1.2.5).

One of the ways of getting a higher utilization out of a single processor chip is then to move from a strategy of further sophistication of the single processor, to a division of the chip into multiple processing 'cores'<fnref n="9" />. The separate cores can work on unrelated tasks, or, by introducing what is in effect data parallelism (section 2.4.1), collaborate on a common task at a higher overall efficiency.

This solves the above two problems:

<footnote n="9">9. Another solution is Intel's *hyperthreading*, which lets a processor mix the instructions of several instruction streams. The benefits of this are strongly dependent on the individual case. However, this same mechanism is exploited with great success in GPUs; see section 2.10.3. For a discussion see section 2.7.1.10</footnote>
* Two cores at a lower frequency can have the same throughput as a single processor at a higher frequency; hence, multiple cores are more energy-efficient.

* Discovered ILP is now replaced by explicit task parallelism, managed by the programmer.

While the first multicore CPUs were simply two processors on the same die, later generations incorporated L3 or L2 caches that were shared between the two processor cores; see figure 1.12 This design makes it

<figure><img src="image_14.png" /><figcaption>Figure 1.12: Cache hierarchy in a single-core and dual-core chip</figcaption></figure>

efficient for the cores to work jointly on the same problem. The cores would still have their own L1 cache, and these separate caches lead to a *cache coherence* problem; see section 1.4.1 below.

We note that the term ‘processor’ is now ambiguous: it can refer to either the chip, or the processor core on the chip. For this reason, we mostly talk about a *socket* for the whole chip and *core* for the part containing one arithmetic and logic unit and having its own registers. Currently, CPUs with 4 or 6 cores are common, even in laptops, and Intel and AMD are marketing 12-core chips. The core count is likely to go up in the future: Intel has already shown an 80-core prototype that is developed into the 48 core ‘Single-chip Cloud Computer’, illustrated in fig 1.13. This chip has a structure with 24 dual-core ‘tiles’ that are connected through a 2D mesh network. Only certain tiles are connected to a memory controller, others can not reach memory other than through the on-chip network.

With this mix of shared and private caches, the programming model for multicore processors is becoming a hybrid between shared and distributed memory:

*   **Core** The cores have their own private L1 cache, which is a sort of distributed memory. The above mentioned Intel 80-core prototype has the cores communicating in a distributed memory fashion.

*   **Socket** On one socket, there is often a shared L2 cache, which is shared memory for the cores.

*   **Node** There can be multiple sockets on a single ‘node’ or motherboard, accessing the same shared memory.
<figure><img src="image_15.png" /><figcaption>Figure 1.13: Structure of the Intel Single-chip Cloud Computer chip</figcaption></figure>

*Network* Distributed memory programming (see the next chapter) is needed to let nodes communicate. Historically, multicore architectures have a precedent in multiprocessor shared memory designs (section 2.5.1) such as the *Sequent Symmetry* and the *Alliant FX/8*. Conceptually the program model is the same, but the technology now allows to shrink a multiprocessor board to a multicore chip.

### 1.4.1 Cache coherence

With parallel processing, there is the potential for a conflict if more than one processor has a copy of the same data item. The problem of ensuring that all cached data are an accurate copy of main memory, is referred to as *cache coherence*: if one processor alters its copy, the other copy needs to be updated.

In distributed memory architectures, a dataset is usually partitioned disjointly over the processors, so conflicting copies of data can only arise with knowledge of the user, and it is up to the user to deal with the problem. The case of shared memory is more subtle: since processes access the same main memory, it would seem that conflicts are in fact impossible. However, processor typically have some private cache, which contains copies of data from memory, so conflicting copies can occur. This situation arises in particular in multicore designs.

Suppose that two cores have a copy of the same data item in their (private) L1 cache, and one modifies its copy. Now the other has cached data that is no longer an accurate copy of its counterpart: the processor will *invalidate* that copy of the item, and in fact its whole cacheline. When the process needs access to the item again, it needs to reload that cacheline. The alternative is for any core that alters data to send that cacheline to the other cores. This strategy probably has a higher overhead, since other cores are not likely to have a copy of a cacheline.

This process of updating or invalidating cachelines is known as *maintaining cache coherence*, and it is done on a very low level of the processor, with no programmer involvement needed. (This makes updating memory locations an *atomic operation*; more about this in section 2.7.1.5.) However, it will slow down the computation, and it wastes bandwidth to the core that could otherwise be used for loading or storing operands.
The state of a cache line with respect to a data item in main memory is usually described as one of the following:

* Scratch: the cache line does not contain a copy of the item;

* Valid: the cache line is a correct copy of data in main memory;

* Reserved: the cache line is the *only* copy of that piece of data;

* Dirty: the cache line has been modified, but not yet written back to main memory;

* Invalid: the data on the cache line is also present on other processors (it is not *reserved*), and another process has modified its copy of the data.

A simpler variant of this is the Modified-Shared-Invalid (MSI) coherence protocol, where a cache line can be in the following states on a given core:

* Modified: the cacheline has been modified, and needs to be written to the backing store. This writing can be done when the line is *evicted*, or it is done immediately, depending on the write-back policy.

* Shared: the line is present in at least one cache and is unmodified.

* Invalid: the line is not present in the current cache, or it is present but a copy in another cache has been modified.

These states control the movement of cachelines between memory and the caches. For instance, suppose a core does a read to a cacheline that is invalid on that core. It can then load it from memory or get it from another cache, which may be faster. (Finding whether a line exists (in state M or S) on another cache is called *snooping*; an alternative is to maintain cache directories.) If the line is Shared, it can now simply be copied; if it is in state M in the other cache, that core first needs to write it back to memory.

**Exercise 1.9.** Consider two processors, a data item $x$ in memory, and cachelines $x_1, x_2$ in the private caches of the two processors to which $x$ is mapped. Describe the transitions between the states of $x_1$ and $x_2$ under reads and writes of $x$ on the two processors. Also indicate which actions cause memory bandwidth to be used. (This list of transitions is a *Finite State Automaton* (FSA); see section 17.)

Variants of the MSI protocol add an 'Exclusive' or 'Owned' state for increased efficiency.

### 1.4.1.1 Snooping

The *snooping* mechanism can also be used as follows: if a core 'listens in' on all bus traffic, it can invalidate or update its own cacheline copies when another core modifies its copy. Invalidating is cheaper than updating since it is a bit operation, while updating involves copying the whole cacheline.

**Exercise 1.10.** When would updating pay off? Write a simple cache simulator to evaluate this question.

### 1.4.1.2 False sharing

The cache coherence problem can even appear if the cores access different items. For instance, a declaration

`double x,y;`
will likely allocate x and y next to each other in memory, so there is a high chance they fall on the same cacheline. Now if one core updates x and the other y, this cacheline will continuously be moved between the cores. This is called *false sharing—underline*. The most common case of false sharing happens when threads update consecutive locations of an array:

local_results[thread_num] = ....

## 1.4.2 Computations on multicore chips

Multicore chips offer a form of shared memory parallelism and can be programmed as such, for instance using OpenMP (section 2.7.2). We will go discuss in some detail the scheduling of linear algebra operations on multicore chips; section 6.12.

# 1.5 Locality and data reuse

By now it should be clear that there is more to the execution of an algorithm than counting the operations: the data transfer involved is important, and can in fact dominate the cost. Since we have caches and registers, the amount of data transfer can be minimized by programming in such a way that data stays as close to the processor as possible. Partly this is a matter of programming cleverly, but we can also look at the theoretical question: does the algorithm allow for it to begin with.

It turns out that in scientific computing there is often data often interacts mostly with data that is close by in some sense, which will lead to data locality; section 1.5.2. Often such locality derives from the nature of the application, as in the case of the Partial Differential Equations (PDEs) you will see in chapter 4. In other cases such as molecular dynamics (chapter 7) there is no such intrinsic locality because all particles interact with all others, and considerable programming cleverness is needed to get high performance.

### 1.5.1 Data reuse and arithmetic intensity

In the previous sections you learned that processor design is somewhat unbalanced: loading data is slower than executing the actual operations. This imbalance is large for main memory and less for the various cache levels. Thus we are motivated to keep data in cache and keep the amount of *data reuse* as high as possible.

Of course, we need to determine first if the computation allows for data to be reused. For this we define the *arithmetic intensity* of an algorithm as follows:

If $n$ is the number of data items that an algorithm operates on, and $f(n)$ the number of operations it takes, then the arithmetic intensity is $f(n)/n$.

(We can measure data items in either floating point numbers or bytes. The latter possibility makes it easier to relate arithmetic intensity to hardware specifications of a processor.)

Arithmetic intensity is also related to *latency hiding*: the concept that you can mitigate the negative performance impact of data loading behind computational activity going on. For this to work, you need more computations than data loads to make this hiding effective. And that is the very definition of computational intensity: a high ratio of operations per byte/word/number loaded.
### 1.5.1.1 Examples

Consider for example the vector addition

$$ \forall i: x_i \leftarrow x_i + y_i. $$

This involves three memory accesses (two loads and one store) and one operation per iteration, giving a arithmetic intensity of 1/3. The axpy (for 'a times x plus y) operation

$$ \forall i: x_i \leftarrow x_i + a \cdot y_i $$

has two operations, but the same number of memory access since the one-time load of a is amortized. It is therefore more efficient than the simple addition, with a reuse of 2/3.

The inner product calculation

$$ \forall i: s \leftarrow s + x_i \cdot y_i $$

is similar in structure to the axpy operation, involving one multiplication and addition per iteration, on two vectors and one scalar. However, now there are only two load operations, since s can be kept in register and only written back to memory at the end of the loop. The reuse here is 1.

Next, consider the *matrix-matrix product*:

$$ \forall i,j: c_{ij} = \sum_k a_{ik} b_{kj}. $$

This involves $3n^2$ data items and $2n^3$ operations, which is of a higher order. The arithmetic intensity is $O(n)$, meaning that every data item will be used $O(n)$ times. This has the implication that, with suitable programming, this operation has the potential of overcoming the bandwidth/clock speed gap by keeping data in fast cache memory.

**Exercise 1.11.** The matrix-matrix product, considered as *operation*, clearly has data reuse by the above definition. Argue that this reuse is not trivially attained by a simple implementation. What determines whether the naive implementation has reuse of data that is in cache?

[In this discussion we were only concerned with the number of operations of a given *implementation*, not the mathematical *operation*. For instance, there are ways of performing the matrix-matrix multiplication and Gaussian elimination algorithms in fewer than $O(n^3)$ operations [151, 134]. However, this requires a different implementation, which has its own analysis in terms of memory access and reuse.]

The matrix-matrix product is the heart of the *LINPACK* benchmark [43]; see section 2.11.5. Using this as the solve measure of *benchmarking* a computer may give an optimistic view of its performance: the matrix-matrix product is an operation that has considerable data reuse, so it is relatively insensitive to memory bandwidth and, for parallel computers, properties of the network. Typically, computers will attain 60–90% of their *peak performance* on the Linpack benchmark. Other benchmark may give considerably lower figures.
### 1.5.1.2 The roofline model

There is an elegant way of talking about how arithmetic intensity, which is a statement about the ideal algorithm, not its implementation, interacts with hardware parameters and the actual implementation to determine performance. This is known as the *roofline model* [165], and it expresses the basic fact that performance is bounded by two factors, illustrated in the left graph of figure 1.14.

1. The *peak performance*, indicated by the horizontal line at the top of the graph, is an absolute bound on the performance<fnref n="10" />, achieved only if every aspect of a CPU (pipelines, multiple floating point units) are perfectly used. The calculation of this number is purely based on CPU properties and clock cycle; it is assumed that memory bandwidth is not a limiting factor.

2. The number of operations per second is also limited by the product of the bandwidth, an absolute number, and the arithmetic intensity:

$$ \frac{\text{operations}}{\text{second}} = \frac{\text{operations}}{\text{data item}} \cdot \frac{\text{data items}}{\text{second}} $$

This is depicted by the linearly increasing line in the graph.

The *roofline model* is an elegant way of expressing that various factors lower the ceiling. For instance, if an algorithm fails to use the full *SIMD width*, this imbalance lowers the attainable peak. The second graph in figure 1.14 indicates various factors that lower the ceiling. There are also various factors that lower the available bandwidth, such as imperfect data hiding. This is indicated by a lowering of the sloping roofline in the third graph.

For a given arithmetic intensity, the performance is determined by where its vertical line intersects the roof line. If this is at the horizontal part, the computation is called *compute-bound*: performance is determined by characteristics of the processor, and bandwidth is not an issue. On the other hand, if that vertical line intersects the sloping part of the roof, the computation is called *bandwidth-bound*: performance is determined by the memory subsystem, and the full capacity of the processor is not used.

**Exercise 1.12.** How would you determine whether a given program kernel is bandwidth or compute bound?

## 1.5.2 Locality

Since using data in cache is cheaper than getting data from main memory, a programmer obviously wants to code in such a way that data in cache is reused. While placing data in cache is not under explicit programmer control, even from assembly language, in most CPUs<fnref n="11" />, it is still possible, knowing the behaviour of the caches, to know what data is in cache, and to some extent to control it.

The two crucial concepts here are *temporal locality* and *spatial locality*. Temporal locality is the easiest to explain: this describes the use of a data element within a short time of its last use. Since most caches have a LRU replacement policy (section 1.3.4.5), if in between the two references less data has been referenced than the cache size, the element will still be in cache and therefore be quickly accessible. With other replacement policies, such as random replacement, this guarantee can not be made.

<footnote n="10">10. An old joke states that the peak performance is that number that the manufacturer guarantees you will never exceed</footnote>

<footnote n="11">11. Low level memory access can be controlled by the programmer in the Cell processor and in some GPUs.</footnote>
<figure><img src="image_16.png" /><figcaption>Figure 1.14: Illustration of factors determining performance in the roofline model</figcaption></figure>
1.5.2.1 Temporal locality

As an example of temporal locality, consider the repeated use of a long vector:

```c
for (loop=0; loop<10; loop++) {
    for (i=0; i<N; i++) {
        ... = ... x[i] ...
    }
}
```

Each element of x will be used 10 times, but if the vector (plus other data accessed) exceeds the cache size, each element will be flushed before its next use. Therefore, the use of x[i] does not exhibit temporal locality: subsequent uses are spaced too far apart in time for it to remain in cache.

If the structure of the computation allows us to exchange the loops:

```c
for (i=0; i<N; i++) {
    for (loop=0; loop<10; loop++) {
        ... = ... x[i] ...
    }
}
```

the elements of x are now repeatedly reused, and are therefore more likely to remain in the cache. This rearranged code displays better temporal locality in its use of x[i].

1.5.2.2 Spatial locality

The concept of spatial locality is slightly more involved. A program is said to exhibit spatial locality if it references memory that is 'close' to memory it already referenced. In the classical von Neumann architecture with only a processor and memory, spatial locality should be irrelevant, since one address in memory can be as quickly retrieved as any other. However, in a modern CPU with caches, the story is different. Above, you have seen two examples of spatial locality:

* Since data is moved in cache lines rather than individual words or bytes, there is a great benefit to coding in such a manner that all elements of the cacheline are used. In the loop
  ```c
  for (i=0; i<N*s; i+=s) {
    ... x[i] ...
  }
  ```

spatial locality is a decreasing function of the stride s.

Let S be the cacheline size, then as s ranges from 1...S, the number of elements used of each cacheline goes down from S to 1. Relatively speaking, this increases the cost of memory traffic in the loop: if s = 1, we load 1/S cachelines per element; if s = S, we load one cacheline for each element. This effect is demonstrated in section 1.6.4.

* A second example of spatial locality worth observing involves the TLB (section 1.3.8). If a program references elements that are close together, they are likely on the same memory page, and address translation through the TLB will be fast. On the other hand, if a program references
many widely disparate elements, it will also be referencing many different pages. The resulting TLB misses are very costly; see also section 1.6.5.

Exercise 1.13. Consider the following pseudocode of an algorithm for summing *n* numbers *x*[i] where *n* is a power of 2:

for s=2,4,8,...,n/2,n:
    for i=0 to n-1 with steps s:
        x[i] = x[i] + x[i+s/2]
sum = x[0]

Analyze the spatial and temporal locality of this algorithm, and contrast it with the
standard algorithm

sum = 0
for i=0,1,2,...,n-1
    sum = sum + x[i]

Exercise 1.14. Consider the following code, and assume that `nvectors` is small compared to the cache size, and `length` large.

for (k=0; k<nvectors; k++)
    for (i=0; i<length; i++)
        a[k,i] = b[i] * c[k]

How do the following concepts relate to the performance of this code:

• Reuse

• Cache size

• Associativity

Would the following code where the loops are exchanged perform better or worse,
and why?

for (i=0; i<length; i++)
    for (k=0; k<nvectors; k++)
        a[k,i] = b[i] * c[k]

1.5.2.3 Examples of locality

Let us examine locality issues for a realistic example. The matrix-matrix multiplication *C* ← *A* · *B* can be computed in several ways. We compare two implementations, assuming that all matrices are stored by rows, and that the cache size is insufficient to store a whole row or column.

for i=1..n      for i=1..n
for j=1..n      for k=1..n
for k=1..n      for j=1..n
    c[i,j] += a[i,k] * b[k,j]    c[i,j] += a[i,k] * b[k,j]

These implementations are illustrated in figure 1.15 The first implemenation constructs the (i, j) element
<figure><img src="image_17.png" /><figcaption>Figure 1.15: Two loop orderings for the $C \leftarrow A \cdot B$ matrix-matrix product</figcaption></figure>

of $C$ by the inner product of a row of $A$ and a column of $B$, in the second a row of $C$ is updated by scaling rows of $B$ by elements of $A$.

Our first observation is that both implementations indeed compute $C \leftarrow C + A \cdot B$, and that they both take roughly $2n^3$ operations. However, their memory behaviour, including spatial and temporal locality is very different.

c[i,j] In the first implementation, $c[i,j]$ is invariant in the inner iteration, which constitutes temporal locality, so it can be kept in register. As a result, each element of $C$ will be loaded and stored only once.

In the second implementation, $c[i,j]$ will be loaded and stored in each inner iteration. In particular, this implies that there are now $n^3$ store operations, a factor of $n$ more than the first implementation.

a[i,k] In both implementations, $a[i,k]$ elements are accessed by rows, so there is good spatial locality, as each loaded cacheline will be used entirely. In the second implementation, $a[i,k]$ is invariant in the inner loop, which constitutes temporal locality; it can be kept in register. As a result, in the second case $A$ will be loaded only once, as opposed to $n$ times in the first case.

b[k,j] The two implementations differ greatly in how they access the matrix $B$. First of all, $b[k,j]$ is never invariant so it will not be kept in register, and $B$ engenders $n^3$ memory loads in both cases. However, the access patterns differ.

In second case, $b[k,j]$ is accessed by rows so there is good spatial locality: cachelines will be fully utilized after they are loaded.

In the first implementation, $b[k,j]$ is accessed by columns. Because of the row storage of the matrices, a cacheline contains a part of a row, so for each cacheline loaded, only one element is used in the columnwise traversal. This means that the first implementation has more loads for $B$ by a factor of the cacheline length. There may also be TLB effects.

Note that we are not making any absolute predictions on code performance for these implementations, or even relative comparison of their runtimes. Such predictions are very hard to make. However, the above discussion identifies issues that are relevant for a wide range of classical CPUs.

**Exercise 1.15.** There are more algorithms for computing the product $C \leftarrow A \cdot B$. Consider the following:

```
for k=1..n:
    for i=1..n:
        for j=1..n:
```
$$c[i,j] += a[i,k]*b[k,j]$$

Analyze the memory traffic for the matrix C, and show that it is worse than the two algorithms given above.

### 1.5.2.4 Core locality

The above concepts of spatial and temporal locality were mostly properties of programs, although hardware properties such as cacheline length and cache size play a role in analyzing the amount of locality. There is a third type of locality that is more intimately tied to hardware: core locality.

A code's execution is said to exhibit core locality if write accesses that are spatially or temporally close are performed on the same core or processing unit. The issue here is that of cache coherence (section 1.4.1) where two cores both have a copy of a certain cacheline in their local stores. If they both read from it there is no problem. However, if one of them writes to it, the coherence protocol will copy the cacheline to the other core's local store. This takes up precious memory bandwidth, so it is to be avoided.

Core locality is not just a property of a program, but also to a large extent of how the program is executed in parallel.

## 1.6 Programming strategies for high performance

In this section we will look at how different ways of programming can influence the performance of a code. This will only be an introduction to the topic; for further discussion see the book by Goedeker and Hoisie [66].

The full listings of the codes and explanations of the data graphed here can be found in appendix 31. All performance results were obtained on the AMD Opteron processors of the Ranger computer [139].

### 1.6.1 Peak performance

For marketing purposes, it may be desirable to define a 'top speed' for a CPU. Since a pipelined floating point unit can yield one result per cycle asymptotically, you would calculate the theoretical peak performance as the product of the clock speed (in ticks per second), number of floating point units, and the number of cores; see section 1.4. This top speed is unobtainable in practice, and very few codes come even close to it. The Linpack benchmark is one of the measures how close you can get to it; the parallel version of this benchmark is reported in the 'top 500'; see section 2.11.5.

### 1.6.2 Pipelining

In section 1.2.1.3 you learned that the floating point units in a modern CPU are pipelined, and that pipelines require a number of independent operations to function efficiently. The typical pipelineable operation is a vector addition; an example of an operation that can not be pipelined is the inner product accumulation
```c
for (i=0; i<N; i++)
    s += a[i]*b[i]
```

The fact that `s` gets both read and written halts the addition pipeline. One way to fill the floating point pipeline is to apply loop unrolling:

```c
for (i = 0; i < N/2-1; i++) {
    sum1 += a[2*i] * b[2*i];
    sum2 += a[2*i+1] * b[2*i+1];
}
```

Now there are two independent multiplies in between the accumulations. With a little indexing optimization this becomes:

```c
for (i = 0; i < N/2-1; i++) {
    sum1 += *(a + 0) * *(b + 0);
    sum2 += *(a + 1) * *(b + 1);

    a += 2; b += 2;
}
```

A first observation about this code is that we are implicitly using associativity and commutativity of addition: while the same quantities are added, they are now in effect added in a different order. As you will see in chapter 3, in computer arithmetic this is not guaranteed to give the exact same result.

In a further optimization, we disentangle the addition and multiplication part of each instruction. The hope is that while the accumulation is waiting for the result of the multiplication, the intervening instructions will keep the processor busy, in effect increasing the number of operations per second.

```c
for (i = 0; i < N/2-1; i++) {
    temp1 = *(a + 0) * *(b + 0);
    temp2 = *(a + 1) * *(b + 1);

    sum1 += temp1; sum2 += temp2;

    a += 2; b += 2;
}
```

Finally, we realize that the furthest we can move the addition away from the multiplication, is to put it right in front of the multiplication of the next iteration:

```c
for (i = 0; i < N/2-1; i++) {
    sum1 += temp1;
    temp1 = *(a + 0) * *(b + 0);

    sum2 += temp2;
```
```R
temp2 = *(a + 1) * *(b + 1);

a += 2; b += 2;
}
s = temp1 + temp2;
```

Of course, we can unroll the operation by more than a factor of two. While we expect an increased performance because of the longer sequence of pipelined operations, large unroll factors need large numbers of registers. Asking for more registers than a CPU has is called register spill, and it will decrease performance.

Another thing to keep in mind is that the total number of operations is unlikely to be divisible by the unroll factor. This requires cleanup code after the loop to account for the final iterations. Thus, unrolled code is harder to write than straight code, and people have written tools to perform such source-to-source transformations automatically.

Cycle times for unrolling the inner product operation up to six times are given in table 1.2. Note that the timings do not show a monotone behaviour at the unrolling by four. This sort of variation is due to various memory-related factors.

<table><caption>Table 1.2: Cycle times for the inner product operation, unrolled up to six times</caption><tr><td>1</td><td>2</td><td>3</td><td>4</td><td>5</td><td>6</td></tr><tr><td>6794</td><td>507</td><td>340</td><td>359</td><td>334</td><td>528</td></tr></table>

### 1.6.3 Cache size

Above, you learned that data from L1 can be moved with lower latency and higher bandwidth than from L2, and L2 is again faster than L3 or memory. This is easy to demonstrate with code that repeatedly access the same data:

```R
for (i=0; i<NRUNS; i++)
  for (j=0; j (<size; j++))
    array[j] = 2.3*array[j]+1.2;
```

If the size parameter allows the array to fit in cache, the operation will be relatively fast. As the size of the dataset grows, parts of it will evict other parts from the L1 cache, so the speed of the operation will be determined by the latency and bandwidth of the L2 cache. This can be seen in figure 1.16. The full code is given in section 31.2.

**Exercise 1.16.** Argue that with a large enough problem and an LRU replacement policy (section 1.3.4.5) essentially all data in the L1 will be replaced in every iteration of the outer loop. Can you write an example code that will let some of the L1 data stay resident?

Often, it is possible to arrange the operations to keep data in L1 cache. For instance, in our example, we could write
<figure><img src="image_18.png" /><figcaption>Figure 1.16: Average cycle count per operation as function of the dataset size</figcaption></figure>

```c
for (b=0; b < size / l1size; b++) {
    blockstart = 0;
    for (i=0; i < NRUNS; i++) {
        for (j=0; j < l1size; j++)
            array blo ckstart + j) = 2.3 * array blo ckstart + j) + 1.2;
    }
    blockstart += l1size;
}
```

assuming that the L1 size divides evenly in the dataset size. This strategy is called *cache blocking* or *blocking for cache reuse*.

**Exercise 1.17.** To arrive at the blocked code, the loop over *j* was split into a loop over blocks and an inner loop over the elements of the block; the outer loop over *i* was then exchanged with the loop over the blocks. In this particular example you could also simply exchange the *i* and *j* loops. Why may this not be optimal for performance?

### 1.6.4 Cache lines

Since data is moved from memory to cache in consecutive chunks named cachelines (see section 1.3.4.6), code that does not utilize all data in a cacheline pays a bandwidth penalty. This is born out by a simple code

```c
for (i=0, n=0; i < L1WORDS; i++, n += stride)
    array[n] = 2.3 * array[n] + 1.2;
```
Here, a fixed number of operations is performed, but on elements that are at distance stride. As this stride increases, we expect an increasing runtime, which is born out by the graph in figure 1.17.

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

The graph also shows a decreasing reuse of cachelines, defined as the number of vector elements divided by the number of L1 misses (on stall; see section 1.3.5).

The full code is given in section 31.3.

1.6.5 TLB

As explained in section 1.3.8, the Translation Look-aside Buffer (TLB) maintains a small list of frequently used memory pages and their locations; addressing data that are located on one of these pages is much faster than data that are not. Consequently, one wants to code in such a way that the number of pages accessed is kept low.

Consider code for traversing the elements of a two-dimensional array in two different ways.

```c
#define INDEX(i,j,m,n) i+j*m
array = (double*) malloc(m*n*sizeof(double));

/* traversal #1 */
for (j=0; j<n; j++)
    for (i=0; i<m; i++)
        array[INDEX(i,j,m,n)] = array[INDEX(i,j,m,n)]+1;

/* traversal #2 */
```
<figure><img src="image_20.png" /><figcaption>Figure 1.18: Number of TLB misses per column as function of the number of columns; columnwise traversal of the array.</figcaption></figure>

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

The results (see Appendix 31.5 for the source code) are plotted in figures 1.19 and 1.18.

Using $m = 1000$ means that, on the *AMD Opteron* which has pages of 512 doubles, we need roughly two pages for each column. We run this example, plotting the number ‘TLB misses’, that is, the number of times a page is referenced that is not recorded in the TLB.

1. In the first traversal this is indeed what happens. After we touch an element, and the TLB records the page it is on, all other elements on that page are used subsequently, so no further TLB misses occur. Figure 1.19 shows that, with increasing $n$, the number of TLB misses per column is roughly two.

2. In the second traversal, we touch a new page for every element of the first row. Elements of the second row will be on these pages, so, as long as the number of columns is less than the number of TLB entries, these pages will still be recorded in the TLB. As the number of columns grows, the number of TLB increases, and ultimately there will be one TLB miss for each element access. Figure 1.18 shows that, with a large enough number of columns, the number of TLB misses per column is equal to the number of elements per column.

### 1.6.6 Cache associativity

There are many algorithms that work by recursive division of a problem, for instance the *Fast Fourier Transform (FFT)* algorithm. As a result, code for such algorithms will often operate on vectors whose
<figure><img src="image_21.png" /><figcaption>Figure 1.19: Number of TLB misses per column as function of the number of columns; rowwise traversal of the array.</figcaption></figure>

length is a power of two. Unfortunately, this can cause conflicts with certain architectural features of a CPU, many of which involve powers of two.

In section 1.3.4.8 you saw how the operation of adding a small number of vectors

$$ \forall j: y_j = y_j + \sum_{i=1}^{m} x_{i,j} $$

is a problem for direct mapped caches or set-associative caches with associativity.

As an example we take the AMD Opteron, which has an L1 cache of 64K bytes, and which is two-way set associative. Because of the set associativity, the cache can handle two addresses being mapped to the same cache location, but not three or more. Thus, we let the vectors be of size $n = 4096$ doubles, and we measure the effect in cache misses and cycles of letting $m = 1, 2, ...$.

First of all, we note that we use the vectors sequentially, so, with a cacheline of eight doubles, we should ideally see a cache miss rate of $1/8$ times the number of vectors $m$. Instead, in figure 1.20 we see a rate approximately proportional to $m$, meaning that indeed cache lines are evicted immediately. The exception here is the case $m = 1$, where the two-way associativity allows the cachelines of two vectors to stay in cache.

Compare this to figure 1.21, where we used a slightly longer vector length, so that locations with the same $j$ are no longer mapped to the same cache location. As a result, we see a cache miss rate around $1/8$, and a smaller number of cycles, corresponding to a complete reuse of the cache lines.

Two remarks: the cache miss numbers are in fact lower than the theory predicts, since the processor will use prefetch streams. Secondly, in figure 1.21 we see a decreasing time with increasing $m$; this is probably due
<figure><img src="image_22.png" /><figcaption>Figure 1.20: The number of L1 cache misses and the number of cycles for each *j* column accumulation, vector length 4096</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 1.21: The number of L1 cache misses and the number of cycles for each *j* column accumulation, vector length 4096 + 8</figcaption></figure>
to a progressively more favourable balance between load and store operations. Store operations are more
expensive than loads, for various reasons.

1.6.7 Loop tiling

In some cases performance can be increased by breaking up a loop into two nested loops, an outer one for the blocks in the iteration space, and an inner one that goes through the block:

for (i=0; i<n; i++)
...

becomes

bs = ... /* the blocksize */
nblocks = n/bs /* assume that n is a multiple of bs */
for (b=0; b<nblocks; b++)
    for (i=b*bs, j=0; j<bs; i++, j++)
        ...

For a single loop this may not make any difference, but given the right context it may. For instance, if an
array is repeatedly used, but it is too large to fit into cache:

for (n=0; n<10; n++)
    for (i=0; i<100000; i++)
        ...
    = ...
    x[i] ...

then loop tiling may lead to a situation where the array is divided into blocks that will fit in cache:

bs = ... /* the blocksize */
for (b=0; b<100000/bs; b++)
    for (n=0; n<10; n++)
        for (i=b*bs; i<(b+1)*bs; i++)
            ...
        = ...
        x[i] ...

For this reason, loop tiling is also known as *cache blocking*. The block size depends on how much data is accessed in the loop body; ideally you would try to make data reused in L1 cache, but it is also possible to block for L2 reuse. Of course, L2 reuse will not give as high a performance as L1 reuse.

Exercise 1.18. Analyze this example. When is x brought into cache, when is it reused, and when is it flushed? What is the required cache size in this example? Rewrite this example, using a constant

#define L1SIZE 65536

In section 1.5.2 we looked at the matrix-matrix multiplication, and concluded that little data could be kept
in cache. With loop tiling we can improve this situation. For instance, the standard way of writing this
product
```
for i=1..n
    for j=1..n
        for k=1..n
            c[i,j] += a[i,k]*b[k,j]
```

can only be optimized to keep `c[i, j]` in register:

```
for i=1..n
    for j=1..n
        s = 0
        for k=1..n
            s += a[i,k]*b[k,j]
            c[i,j] += s
```

Using loop tiling we can easily keep parts of `a[i, :]` in cache, assuming that `a` is stored by rows:

```
for kk=1..n/bs
    for i=1..n
        for j=1..n
            s = 0
            for k=(kk-1)*bs+1..kk*bs
                s += a[i,k]*b[k,j]
            c[i,j] += s
```

### 1.6.8 Optimization strategies

<figure><img src="image_24.png" /><figcaption>Figure 1.22: Performance of naive and optimized implementations of the Discrete Fourier Transform</figcaption></figure>

Figures 1.22 and 1.23 show that there can be wide discrepancy between the performance of naive implementations of an operation (sometimes called the ‘reference implementation’), and optimized implementations. Unfortunately, optimized implementations are not simple to find. For one, since they rely on blocking, their
<figure><img src="image_25.png" /><figcaption>Figure 1.23: Performance of naive and optimized implementations of the matrix-matrix product</figcaption></figure>

loop nests are double the normal depth: the matrix-matrix multiplication becomes a six-deep loop. Then, the optimal block size is dependent on factors like the target architecture.

We make the following observations:

* Compilers are not able to extract anywhere close to optimal performance<fnref n="12" />.

* There are *autotuning* projects for automatic generation of implementations that are tuned to the architecture. This approach can be moderately to very successful. Some of the best known of these projects are Atlas [162] for Blas kernels, and Spiral [138] for transforms.

### 1.6.9 Cache aware and cache oblivious programming

Unlike registers and main memory, both of which can be addressed in (assembly) code, use of caches is implicit. There is no way a programmer can load data explicitly to a certain cache, even in assembly language.

However, it is possible to code in a 'cache aware' manner. Suppose a piece of code repeatedly operates on an amount of data that less data than the cache size. We can assume that the first time the data is accessed, it is brought into cache; the next time it is accessed it will already be in cache. On the other hand, if the amount of data is more than the cache size<fnref n="13" />, it will partly or fully be flushed out of cache in the process of accessing it.

We can experimentally demonstrate this phenomenon. With a very accurate counter, the code fragment

```
for (x=0; x<NX; x++)
    for (i=0; i<N; i++)
        a[i] = sqrt(a[i]);
```

<footnote n="12">12. Presenting a compiler with the reference implementation may still lead to high performance, since some compilers are trained to recognize this operation. They will then forego translation and simply replace it by an optimized variant.</footnote>

<footnote n="13">13. We are conveniently ignoring matters of set-associativity here, and basically assuming a fully associative cache.</footnote>
will take time linear in N up to the point where a fills the cache. An easier way to picture this is to compute
a normalized time, essentially a time per execution of the inner loop:

```c
    t = time();
    for (x=0; x<NX; x++)
        for (i=0; i<N; i++)
            a[i] = sqrt(a[i]);
    t = time() - t;
    t_normalized = t / (N * NX);
```

The normalized time will be constant until the array a fills the cache, then increase and eventually level off again. (See section 1.6.3 for an elaborate discussion.)

The explanation is that, as long as `a[0]...a[N-1]` fit in L1 cache, the inner loop will use data from the L1 cache. Speed of access is then determined by the latency and bandwidth of the L1 cache. As the amount of data grows beyond the L1 cache size, some or all of the data will be flushed from the L1, and performance will be determined by the characteristics of the L2 cache. Letting the amount of data grow even further, performance will again drop to a linear behaviour determined by the bandwidth from main memory.

If you know the cache size, it is possible in cases such as above to arrange the algorithm to use the cache optimally. However, the cache size is different per processor, so this makes your code not portable, or at least its high performance is not portable. Also, blocking for multiple levels of cache is complicated. For these reasons, some people advocate *cache oblivious programming* [61].

Cache oblivious programming can be described as a way of programming that automatically uses all levels of the *cache hierarchy*. This is typically done by using a *divide-and-conquer* strategy, that is, recursive subdivision of a problem.

As a simple example of cache oblivious programming is the *matrix transposition* operation $B \leftarrow A^t$. First we observe that each element of either matrix is accessed once, so the only reuse is in the utilization of cache lines. If both matrices are stored by rows and we traverse $B$ by rows, then $A$ is traversed by columns, and for each element accessed one cacheline is loaded. If the number of rows times the number of elements per cacheline is more than the cachesize, lines will be evicted before they can be reused.

<figure><img src="image_26.png" /><figcaption>Figure 1.24: Matrix transpose operation, with simple and recursive traversal of the source matrix</figcaption></figure>
In a cache oblivious implementation we divide A and B as 2 × 2 block matrices, and recursively compute $B_{11} \leftarrow A_{11}^t$, $B_{12} \leftarrow A_{21}^t$, et cetera; see figure 1.24. At some point in the recursion, blocks $A_{ij}$ will now be small enough that they fit in cache, and the cachelines of A will be fully used. Hence, this algorithm improves on the simple one by a factor equal to the cacheline size.

The cache oblivious strategy can often yield improvement, but it is not necessarily optimal. In the matrix-matrix product it improves on the naive algorithm, but it is not as good as an algorithm that is explicitly designed to make optimal use of caches [72].

### 1.6.10 Case study: Matrix-vector product

Let us consider in some detail the *matrix-vector product*

$$ \forall i,j: y_i \leftarrow a_{ij} \cdot x_j $$

This involves $2n^2$ operations on $n^2 + 2n$ data items, so reuse is $O(1)$: memory accesses and operations are of the same order. However, we note that there is a double loop involved, and the $x, y$ vectors have only a single index, so each element in them is used multiple times.

Exploiting this theoretical reuse is not trivial. In

```c
/* variant 1 */
for (i)
    for (j)
        y[i] = y[i] + a[i][j] * x[j];
```

the element $y[i]$ seems to be reused. However, the statement as given here would write $y[i]$ to memory in every inner iteration, and we have to write the loop as

```c
/* variant 2 */
for (i) {
    s = 0;
    for (j)
        s = s + a[i][j] * x[j];
    y[i] = s;
}
```

to ensure reuse. This variant uses $2n^2$ loads and $n$ stores.

This code fragment only exploits the reuse of $y$ explicitly. If the cache is too small to hold the whole vector $x$ plus a column of $a$, each element of $x$ is still repeatedly loaded in every outer iteration.

Reversing the loops as

```c
/* variant 3 */
for (j)
    for (i)
        y[i] = y[i] + a[i][j] * x[j];
```
exposes the reuse of x, especially if we write this as

```c
/* variant 3 */
for (j) {
    t = x[j];
    for (i)
        y[i] = y[i] + a[i][j] * t;
}
```

but now y is no longer reused. Moreover, we now have $2n^2+n$ loads, comparable to variant 2, but $n^2$ stores, which is of a higher order.

It is possible to get reuse both of x and y, but this requires more sophisticated programming. The key here is split the loops into blocks. For instance:

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

This is also called *loop unrolling*, or *strip mining*. The amount by which you unroll loops is determined by the number of available registers.

### 1.6.11 Case study: Goto matrix-matrix product

In section 1.5.1 we already argued that the *matrix-matrix product* (or *dgemm* in Blas terms) has a large amount of possible data reuse: there are $O(n^3)$ operations on $O(n^2)$ data. We will now consider an implementation, due to Kazushige Goto [72], that indeed achieves close to peak performance.

The *matrix-matrix* algorithm has three loops, each of which we can block, giving a six-way nested loop. Since there are no recurrences on the output elements, all resulting loop exchanges are legal. Combine this with the fact that the loop blocking introduces three blocking parameters, and you'll see that the number of potential implementations is enormous. Here we present the global reasoning that underlies the Goto implementation; for a detailed discussion see the paper cited.

We start by writing the product $C \leftarrow A \cdot B$ (or $C \leftarrow C + AB$ according to the Blas standard) as a sequence of low-rank updates:

$$C_{**} = \sum_{k} A_{*k} B_{k*}$$

See figure 1.25. Next we derive the ‘block-panel’ multiplication by multiplying a block of A by a ‘sliver’ of B; see figure 1.26. Finally, the inner algorithm accumulates a small row $C_{i,*}$, typically of small size such as 4, by accumulating:
<figure><img src="image_27.png" /><figcaption>Figure 1.25: Matrix-matrix multiplication as a sequence of low-rank updates</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Figure 1.26: The block-panel multiplication in the matrix-matrix algorithm</figcaption></figure>

```fortran
// compute C[i,*] :
for k:
    C[i,*] = A[i,k] * B[k,*]
```

See figure 1.27. Now this algorithm is tuned.

<figure><img src="image_29.png" /><figcaption>Figure 1.27: The register-resident kernel of the matrix-matrix multiply</figcaption></figure>

* We need enough registers for $C[i,*]$, $A[i,k]$ and $B[k,*]$. On current processors that means that we accumulate four elements of $C$.

* Those elements of $C$ are accumulated, so they stay in register and the only data transfer is loading of the elements of $A$ and $B$; there are no stores!

* The elements $A[i,k]$ and $B[k,*]$ stream from L1.

* Since the same block of $A$ is used for many successive slivers of $B$, we want it to stay resident; we choose the blocksize of $A$ to let it stay in L2 cache.
* In order to prevent TLB problems, *A* is stored by rows. If we start out with a matrix in (Fortran) column-major storage, this means we have to make a copy. Since copying is of a lower order of complexity, this cost is amortizable.

# 1.7 Power consumption

Another important topic in high performance computers is their power consumption. Here we need to distinguish between the power consumption of a single processor chip, and that of a complete cluster.

As the number of components on a chip grows, its power consumption would also grow. Fortunately, in a counter acting trend, miniaturization of the chip features has simultaneously been reducing the necessary power. Suppose that the feature size $\lambda$ (think: thickness of wires) is scaled down to $s\lambda$ with $s < 1$. In order to keep the electric field in the transistor constant, the length and width of the channel, the oxide thickness, substrate concentration density and the operating voltage are all scaled by the same factor.

## 1.7.1 Derivation of scaling properties

The properties of *constant field scaling* or *Dennard scaling* [15, 37] are an ideal-case description of the properties of a circuit as it is miniaturized. One important result is that power density stays constant as chip features get smaller, and the frequency is simultaneously increased.

The basic properties derived from circuit theory are that, if we scale feature size down by $s$:

<table><tr><td>Feature size</td><td>~ s</td></tr><tr><td>Voltage</td><td>~ s</td></tr><tr><td>Current</td><td>~ s</td></tr><tr><td>Frequency</td><td>~ s</td></tr></table>

Then we can derive that

$$ \text{Power} = V \cdot I \sim s^2, $$

and because the total size of the circuit also goes down with $s^2$, the power density stays the same. Thus, it also becomes possible to put more transistors on a circuit, and essentially not change the cooling problem.

This result can be considered the driving force behind *Moore's law*, which states that the number of transistors in a processor doubles every 18 months.

The frequency-dependent part of the power a processor needs comes from charging and discharging the capacitance of the circuit, so

$$ \begin{array}{ll} \text{Charge} & q = CV \\ \text{Work} & W = qV = CV^2 \\ \text{Power} & W/\text{time} = WF = CV^2F \end{array} \qquad (1.1) $$

This analysis can be used to justify the introduction of multicore processors.
1.7.2 Multicore

At the time of this writing (circa 2010), miniaturization of components has almost come to a standstill, because further lowering of the voltage would give prohibitive leakage. Conversely, the frequency can not be scaled up since this would raise the heat production of the chip too far. Figure 1.28 gives a dramatic

<figure><img src="image_30.png" /><figcaption>Figure 1.28: Projected heat dissipation of a CPU if trends had continued – this graph courtesy Pat Helsinger
illustration of the heat that a chip would give off, if single-processor trends had continued.</figcaption></figure>

One conclusion is that computer design is running into a power wall, where the sophistication of a single core can not be increased any further (so we can for instance no longer increase *ILP* and *pipeline depth*) and the only way to increase performance is to increase the amount of explicitly visible parallelism. This development has led to the current generation of multicore processors; see section 1.4. It is also the reason GPUs with their simplified processor design and hence lower energy consumption are attractive; the same holds for Field-Programmable Gate Arrays (FPGAs). One solution to the power wall problem is introduction of multicore processors. Recall equation 1.1, and compare a single processor to two processors at half the frequency. That should have the same computing power, right? Since we lowered the frequency, we can lower the voltage if we stay with the same process technology.

The total electric power for the two processors (cores) is, ideally,

$$
\left.
\begin{array}{l}
C_{\text{multi}} = 2C \\
F_{\text{multi}} = F/2 \\
V_{\text{multi}} = V/2
\end{array}
\right\}
\Rightarrow P_{\text{multi}} = P/4.
$$

In practice the capacitance will go up by a little over 2, and the voltage can not quite be dropped by 2, so
it is more likely that $P_{\text{multi}} \approx 0.4 \times P$ [23]. Of course the integration aspects are a little more complicated
in practice [16]; the important conclusion is that now, in order to lower the power (or, conversely, to allow
further increase in performance while keeping the power constant) we now have to start programming in
parallel.
1.7.3 Total computer power

The total power consumption of a parallel computer is determined by the consumption per processor and the number of processors in the full machine. At present, this is commonly several Megawatts. By the above reasoning, the increase in power needed from increasing the number of processors can no longer be offset by more power-effective processors, so power is becoming the overriding consideration as parallel computers move from the petascale (attained in 2008 by the IBM Roadrunner) to a projected exascale.

In the most recent generations of processors, power is becoming an overriding consideration, with influence in unlikely places. For instance, the Single Instruction Multiple Data (SIMD) design of processors (see section 2.4.1, in particular section 2.4.1.2) is dictated by the power cost of instruction decoding.

**1.8 Review questions**

For the true/false questions, give short explanation if you choose the ‘false’ answer.

Exercise 1.19. True or false. The code

for (i=0; i<N; i++)
a[i] = b[i]+1;

touches every element of a and b once, so there will be a cache miss for each element.

Exercise 1.20. Give an example of a code fragment where a 3-way associative cache will have conflicts, but a 4-way cache will not.

Exercise 1.21. Consider the matrix-vector product with an $N \times N$ matrix. What is the needed cache size to execute this operation with only compulsory cache misses? Your answer depends on how the operation is implemented: answer separately for rowwise and columnwise traversal of the matrix, where you can assume that the matrix is always stored by rows.
# Chapter 2

## Parallel Computing

The largest and most powerful computers are sometimes called ‘supercomputers’. For the last two decades, this has, without exception, referred to parallel computers: machines with more than one CPU that can be set to work on the same problem.

Parallelism is hard to define precisely, since it can appear on several levels. In the previous chapter you already saw how inside a CPU several instructions can be ‘in flight’ simultaneously. This is called *instruction-level parallelism*, and it is outside explicit user control: it derives from the compiler and the CPU deciding which instructions, out of a single instruction stream, can be processed simultaneously. At the other extreme is the sort of parallelism where more than one instruction stream is handled by multiple processors, often each on their own circuit board. This type of parallelism is typically explicitly scheduled by the user.

In this chapter, we will analyze this more explicit type of parallelism, the hardware that supports it, the programming that enables it, and the concepts that analyze it.

### 2.1 Introduction

In scientific codes, there is often a large amount of work to be done, and it is often regular to some extent, with the same operation being performed on many data. The question is then whether this work can be sped up by use of a parallel computer. If there are *n* operations to be done, and they would take time *t* on a single processor, can they be done in time *t*/*p* on *p* processors?

Let us start with a very simple example. Adding two vectors of length *n*

```fortran
for (i=0; i<n; i++)
    a[i] = b[i] + c[i];
```

can be done with up to *n* processors. In the idealized case with *n* processors, each processor has local scalars *a*, *b*, *c* and executes the single instruction *a*=*b*+*c*. This is depicted in figure 2.1.

In the general case, where each processor executes something like

```fortran
for (i=my_low; i<my_high; i++)
    a[i] = b[i] + c[i];
```
<figure><img src="image_31.png" /><figcaption>Figure 2.1: Parallelization of a vector addition</figcaption></figure>

execution time is linearly reduced with the number of processors. If each operation takes a unit time, the
original algorithm takes time $n$, and the parallel execution on $p$ processors $n/p$. The parallel algorithm is
faster by a factor of $p^1$.

Next, let us consider summing the elements of a vector. We again assume that each processor contains just
a single array element. The sequential code:

```fortran
s = 0;
for (i=0; i<n; i++)
    s += x[i]
```

is no longer obviously parallel, but if we recode the loop as

```fortran
for (s=2; s<2*n; s*=2)
    for (i=0; i<n-s/2; i+=s)
        x[i] += x[i+s/2]
```

there is a way to parallelize it: every iteration of the outer loop is now a loop that can be done by $n/s$
processors in parallel. Since the outer loop will go through $\log_2 n$ iterations, we see that the new algorithm
has a reduced runtime of $n/p \cdot \log_2 n$. The parallel algorithm is now faster by a factor of $p/\log_2 n$. This is
depicted in figure 2.2.

Even from these two simple examples we can see some of the characteristics of parallel computing:

* Sometimes algorithms need to be rewritten slightly to make them parallel.

* A parallel algorithm may not show perfect speedup.

There are other things to remark on. In the first case, if each processors has its $x_i$, $y_i$ in a local store the algo-
rithm can be executed without further complications. In the second case, processors need to communicate
data among each other and we haven't assigned a cost to that yet.

<footnote n="1">1. We ignore lower order errors in this result when *p* does not divide perfectly in *n*. We will also, in general, ignore matters of loop overhead.</footnote>
<figure><img src="image_32.png" /><figcaption>Figure 2.2: Parallelization of a vector addition</figcaption></figure>

First let us look systematically at communication. We can take the parallel algorithm in the right half of figure 2.2 and turn it into a tree graph (see Appendix 16) by defining the inputs as leave nodes, all partial sums as interior nodes, and the root as the total sum. There is an edge from one node to another if the first is input to the (partial) sum in the other. This is illustrated in figure 2.3. In this figure nodes are horizontally aligned with other computations that can be performed simultaneously; each level is sometimes called a *superstep* in the computation. Nodes are vertically aligned if they are computed on the same processors, and an arrow corresponds to a communication if it goes from one processor to another. The vertical alignment

<figure><img src="image_33.png" /><figcaption>Figure 2.3: Communication structure of a parallel vector addition</figcaption></figure>

in figure 2.3 is not the only one possible. If nodes are shuffled within a superstep or horizontal level, a different communication pattern arises.

**Exercise 2.1.** Consider placing the nodes within a superstep on random processors. Show that, if no two nodes wind up on the same processor, at most twice the number of communications is performed from the case in figure 2.3.

**Exercise 2.2.** Can you draw the graph of a computation that leaves the sum result on each processor? There is a solution that takes twice the number of supersteps, and there
is one that takes the same number. In both cases the graph is no longer a tree, but a more general Directed Acyclic Graph (DAG).

Processors are often connected through a network, and moving data through this network takes time. This introduces a concept of distance between the processors. This is easily seen in figure 2.3 where the processors are linearly ordered. If the network only connects a processor with its immediate neighbours, each iteration of the outer loop increases the distance over which communication takes place.

**Exercise 2.3.** Assume that an addition takes a certain unit time, and that moving a number from one processor to another takes that same unit time. Show that the communication time equals the computation time.

Now assume that sending a number from processor $p$ to $p \pm k$ takes time $k$. Show that the execution time of the parallel algorithm now is of the same order as the sequential time.

The summing example made the unrealistic assumption that every processor initially stored just one vector element: in practice we will have $P < N$, and every processor stores a number of vector elements. The obvious strategy is to give each processor a consecutive stretch of elements, but sometimes the obvious strategy is not the best.

**Exercise 2.4.** Consider the case of summing 8 elements with 4 processors. Show that some of the edges in the graph of figure 2.3 no longer correspond to actual communications. Now consider summing 16 elements with, again, 4 processors. What is the number of communication edges this time?

These matters of algorithm adaptation, efficiency, and communication, are crucial to all of parallel computing. We will return to these issues in various guises throughout this chapter.

## 2.2 Functional parallelism versus data parallelism

From the above introduction we can describe parallelism as finding independent operations in the execution of a program. In all of the examples these independent operations were in fact identical operations, but applied to different data items. We could call this *data parallelism*: the same operation is applied in parallel to many data elements. This is in fact a common scenario in scientific computing: parallelism often stems from the fact that a data set (vector, matrix, graph,...) is spread over many processors, each working on its part of the data.

The term *data parallelism* is traditionally mostly applied if the operation is a single instruction; in the case of a subprogram it is often called *task parallelism*.

It is also possible to find independence, not based on data elements, but based on the instructions themselves. Traditionally, compilers analyze code in terms of ILP: independent instructions can be given to separate floating point units, or reordered, for instance to optimize register usage (see also section 2.6.2). ILP is one case of *functional parallelism*; on a higher level, functional parallelism can be obtained by considering independent subprograms, often called *task parallelism*; see section 2.6.3.

Some examples of functional parallelism are Monte Carlo simulations, and other algorithms that traverse a parametrized search space, such as boolean *satisfyability* problems.
## 2.3 Quantifying parallelism

There are two important reasons for using a parallel computer: to have access to more memory or to obtain higher performance. It is easy to characterize the gain in memory, as the total memory is the sum of the individual memories. The speed of a parallel computer is harder to characterize. This section will have an extended discussion on theoretical measures for expressing and judging the gain in execution speed from going to a parallel architecture.

### 2.3.1 Definitions

#### 2.3.1.1 Speedup and efficiency

A simple approach to defining speedup is to let the same program run on a single processor, and on a parallel machine with $p$ processors, and to compare runtimes. With $T_1$ the execution time on a single processor and $T_p$ the time on $p$ processors, we define the speedup as $S_p = T_1/T_p$. (Sometimes $T_1$ is defined as 'the best time to solve the problem on a single processor', which allows for using a different algorithm on a single processor than in parallel.) In the ideal case, $T_p = T_1/p$, but in practice we don't expect to attain that, so $S_p \le p$. To measure how far we are from the ideal speedup, we introduce the efficiency $E_p = S_p/p$. Clearly, $0 < E_p \le 1$.

There is a practical problem with this definition: a problem that can be solved on a parallel machine may be too large to fit on any single processor. Conversely, distributing a single processor problem over many processors may give a distorted picture since very little data will wind up on each processor. Below we will discuss more realistic measures of speed-up.

There are various reasons why the actual speed is less than $P$. For one, using more than one processors necessitates communication, which is overhead that was not part of the original computation. Secondly, if the processors do not have exactly the same amount of work to do, they may be idle part of the time (this is known as *load unbalance*), again lowering the actually attained speedup. Finally, code may have sections that are inherently sequential.

Communication between processors is an important source of a loss of efficiency. Clearly, a problem that can be solved without communication will be very efficient. Such problems, in effect consisting of a number of completely independent calculations, is called *embarrassingly parallel*; it will have close to a perfect speedup and efficiency.

**Exercise 2.5.** The case of speedup larger than the number of processors is called *superlinear speedup*. Give a theoretical argument why this can never happen.

In practice, superlinear speedup can happen. For instance, suppose a problem is too large to fit in memory, and a single processor can only solve it by swapping data to disc. If the same problem fits in the memory of two processors, the speedup may well be larger than 2 since disc swapping no longer occurs. Having less, or more localized, data may also improve the cache behaviour of a code.

#### 2.3.1.2 Critical path

The above definitions of speedup and efficiency made an implicit assumption that parallel work can be arbitrarily subdivided. As you saw in the summing example in section 2.1, this may not always be the
case: there can be dependencies between operations, which limits the amount of parallelism that can be employed.

In order to take this into account, we define the following concepts:

### Definition 1

$T_1$ : the time the computation takes on a single processor

$T_p$ : the time the computation takes with $p$ processors

$T_\infty$ : the time the computation takes if unlimited processors are available

$P_\infty$ : the value of $p$ for which $T_p = T_\infty$

With these concepts, we can define the *average parallelism* of an algorithm as $T_1/T_\infty$.

We will now give a few illustrations by showing a graph of tasks and their dependencies. We assume for simplicity that each node is a unit time task.

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

The maximum number of processors that can be used is 2 and the average parallelism is $4/3$:

$$T_1 = 4, \quad T_\infty = 3 \quad \Rightarrow T_1/T_\infty = 4/3$$

$$T_2 = 3, \quad S_2 = 4/3, \quad E_2 = 2/3$$

$$P_\infty = 2$$

The maximum number of processors that can be used is 3 and the average parallelism is $9/5$; efficiency is maximal for $p = 2$:

$$T_1 = 9, \quad T_\infty = 5 \quad \Rightarrow T_1/T_\infty = 9/5$$

$$T_2 = 6, \quad S_2 = 3/2, \quad E_2 = 3/4$$

$$T_3 = 5, \quad S_3 = 9/5, \quad E_3 = 3/5$$

$$P_\infty = 3$$

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

The maximum number of processors that can be used is 4 and that is also the average parallelism; the figure illustrates a parallelization with $P = 3$ that has efficiency $\equiv 1$:

$$T_1 = 12, \quad T_\infty = 4 \quad \Rightarrow T_1/T_\infty = 3$$

$$T_2 = 6, \quad S_2 = 2, \quad E_2 = 1$$

$$T_3 = 4, \quad S_3 = 3, \quad E_3 = 1$$

$$T_4 = 3, \quad S_4 = 4, \quad E_4 = 1$$

$$P_\infty = 4$$

Based on these examples, you probably see that there are two extreme cases:
* If every task depends on precisely on other, you get a chain of dependencies, and $T_p = T_1$ for any $p$.

* On the other hand, if all tasks are independent (and $p$ divides their number) you get $T_p = T_1/p$ for any $p$.

* In a slightly less trivial scenario than the previous, consider the case where the critical path is of length $m$, and in each of these $m$ steps there are $p-1$ independent tasks, or at least: dependent only on tasks in the previous step. There will then be perfect parallelism in each of the $m$ steps, and we can express $T_p = T_1/p$ or $T_p = m + (T_1 - m)/p$.

That last statement actually holds in general. This is known as Brent's theorem:

**Theorem 1** Let $m$ be the total number of tasks, $p$ the number of processors, and $t$ the length of a critical path. Then the computation can be done in

$$T_p = t + \frac{m-t}{p}.$$

*Proof.* Divide the computation in steps, such that tasks in step $i+1$ are independent of each other, and only dependent on step $i$. Let $s_i$ be the number of tasks in step $i$, then the time for that step is $\lceil \frac{s_i}{p} \rceil$. Summing over $i$ gives

$$T_p = \sum_{i}^{t} \lceil \frac{s_i}{p} \rceil \le \sum_{i}^{t} \frac{s_i + p - 1}{p} = t + \sum_{i}^{t} \frac{s_i - 1}{p} = t + \frac{m-t}{p}.$$

### 2.3.2 Asymptotics

If we ignore limitations such as that the number of processors has to be finite, or the physicalities of the interconnect between them, we can derive theoretical results on the limits of parallel computing. This section will give a brief introduction to such results, and discuss their connection to real life high performance computing.

Consider for instance the matrix-matrix multiplication $C = AB$, which takes $2N^3$ operations where $N$ is the matrix size. Since there are no dependencies between the operations for the elements of $C$, we can perform them all in parallel. If we had $N^2$ processors, we could assign each to an $(i, j)$ coordinate in $C$, and have it compute $c_{ij}$ in $2N$ time. Thus, this parallel operation has efficiency 1, which is optimal.

**Exercise 2.6.** Show that this algorithm ignores some serious issues about memory usage:

* If the matrix is kept in shared memory, how many simultaneous reads from each memory locations are performed?

* If the processors keep the input and output to the local computations in local storage, how much duplication is there of the matrix elements?

Adding $N$ numbers $\{x_i\}_{i=1...N}$ can be performed in $\log_2 N$ time with $N/2$ processors. As a simple example, consider the sum of $n$ numbers: $s = \sum_{i=1}^{n} a_i$. If we have $n/2$ processors we could compute:

1. Define $s_i^{(0)} = a_i$.

2. Iterate with $j = 1, ..., \log_2 n$:
3. Compute $n/2^j$ partial sums $s_i^{(j)} = s_{2i}^{(j-1)} + s_{2i+1}^{(j-1)}$

We see that the $n/2$ processors perform a total of $n$ operations (as they should) in $\log_2 n$ time. The efficiency of this parallel scheme is $O(1/\log_2 n)$, a slowly decreasing function of $n$.

**Exercise 2.7.** Show that, with the scheme for parallel addition just outlined, you can multiply two matrices in $\log_2 N$ time with $N^3/2$ processors. What is the resulting efficiency?

It is now a legitimate theoretical question to ask

* If we had infinitely many processors, what is the lowest possible time complexity for matrix-matrix multiplication, or

* Are there faster algorithms that still have $O(1)$ efficiency?

Such questions have been researched (see for instance [82]), but they have little bearing on high performance computing.

A first objection to these kinds of theoretical bounds is that they implicitly assume some form of shared memory. In fact, the formal model for these algorithms is called a *Parallel Random Access Machine (PRAM)*, where the assumption is that every memory location is accessible to any processor. Often an additional assumption is made that multiple access to the same location are in fact possible<fnref n="2" />. These assumptions are unrealistic in practice, especially in the context of scaling up the problem size and the number of processors. A further objection to the PRAM model is that even on a single processor it ignores the memory hierarchy; section 1.3.

But even if we take distributed memory into account, theoretical results can still be unrealistic. The above summation algorithm can indeed work unchanged in distributed memory, except that we have to worry about the distance between active processors increasing as we iterate further. If the processors are connected by a linear array, the number of 'hops' between active processors doubles, and with that, asymptotically, the computation time of the iteration. The total execution time then becomes $n/2$, a disappointing result given that we throw so many processors at the problem.

What if the processors are connected with a hypercube topology? It is not hard to see that the summation algorithm can then indeed work in $\log_2 n$ time. However, as $n \to \infty$, can we build a sequence of hypercubes of $n$ nodes and keep the communication time between two connected constant? Since communication time depends on latency, which partly depends on the length of the wires, we have to worry about the physical distance between nearest neighbours.

The crucial question here is whether the hypercube (an $n$-dimensional object) can be embedded in 3-dimensional space, while keeping the distance (measured in meters) constant between connected neighbours. It is easy to see that a 3-dimensional grid can be scaled up arbitrarily while maintaining a unit wire length, but the question is not clear for a hypercube. There, the length of the wires may have to increase as $n$ grows, which runs afoul of the finite speed of electrons.

We sketch a proof (see [56] for more details) that, in our three dimensional world and with a finite speed of light, speedup is limited to $\sqrt[4]{n}$ for a problem on $n$ processors, no matter the interconnect. The argument goes as follows. Consider an operation that involves collecting a final result on one processor. Assume that each processor takes a unit volume of space, produces one result per unit time, and can send one data

<footnote n="2">2. This notion can be made precise; for instance, one talks of a CREW-PRAM, for Concurrent Read, Exclusive Write PRAM.</footnote>
item per unit time. Then, in an amount of time t, at most the processors in a ball with radius t, that is, $O(t^3)$ processors can contribute to the final result; all others are too far away. In time T, then, the number of operations that can contribute to the final result is $\int_0^T t^3 dt = O(T^4)$. This means that the maximum achievable speedup is the fourth root of the sequential time.

Finally, the question ‘what if we had infinitely many processors’ is not realistic as such, but we will allow it in the sense that we will ask the weak scaling question (section 2.3.4) ‘what if we let the problem size and the number of processors grow proportional to each other’. This question is legitimate, since it corresponds to the very practical deliberation whether buying more processors will allow one to run larger problems, and if so, with what ‘bang for the buck’.

### 2.3.3 Amdahl's law

One reason for less than perfect speedup is that parts of a code can be inherently sequential. This limits the parallel efficiency as follows. Suppose that 5% of a code is sequential, then the time for that part can not be reduced, no matter how many processors are available. Thus, the speedup on that code is limited to a factor of 20. This phenomenon is known as Amdahl's Law [2], which we will now formulate.

Let $F_s$ be the sequential fraction and $F_p$ be the parallel fraction (or more strictly: the ‘parallelizable’ fraction) of a code, respectively. Then $F_p + F_s = 1$. The parallel execution time $T_p$ on p processors is the sum of the part that is sequential $T_1 F_s$ and the part that can be parallelized $T_1 F_p / P$:

$$T_P = T_1(F_s + F_p/P). \qquad (2.1)$$

As the number of processors grows $P \to \infty$, the parallel execution time now approaches that of the sequential fraction of the code: $T_P \downarrow T_1 F_s$. We conclude that speedup is limited by $S_P \le 1/F_s$ and efficiency is a decreasing function $E \sim 1/P$.

The sequential fraction of a code can consist of things such as I/O operations. However, there are also parts of a code that in effect act as sequential. Consider a program that executes a single loop, where all iterations can be computed independently. Clearly, this code is easily parallelized. However, by splitting the loop in a number of parts, one per processor, each processor now has to deal with loop overhead: calculation of bounds, and the test for completion. This overhead is replicated as many times as there are processors. In effect, loop overhead acts as a sequential part of the code.

**Exercise 2.8.** Let's do a specific example. Assume that a code has a setup that takes 1 second and a parallelizable section that takes 1000 seconds on one processor. What are the speedup and efficiency if the code is executed with 100 processors? What are they for 500 processors? Express your answer to at most two significant digits.

**Exercise 2.9.** Investigate the implications of Amdahl's law: if the number of processors P increases, how does the parallel fraction of a code have to increase to maintain a fixed efficiency?

#### 2.3.3.1 Amdahl's law with communication overhead

In a way, Amdahl's law, sobering as it is, is even optimistic. Parallelizing a code will give a certain speedup, but it also introduces communication overhead that will lower the speedup attained. Let us refine our model
of equation (2.1) (see [108, p. 367]):

$$T_p = T_1(F_s + F_p/P) + T_c,$$

where $T_c$ is a fixed communication time.

To assess the influence of this communication overhead, we assume that the code is fully parallelizable, that is, $F_p = 1$. We then find that

$$S_p = \frac{T_1}{T_1/p + T_c} \qquad (2.2)$$

For this to be close to $p$, we need $T_c \ll T_1/p$ or $p \ll T_1/T_c$. In other words, the number of processors should not grow beyond the ratio of scalar execution time and communication overhead.

### 2.3.3.2 Gustafson's law

Amdahl's law was thought to show that large numbers of processors would never pay off. However, the implicit assumption in Amdahl's law is that there is a fixed computation which gets executed on more and more processors. In practice this is not the case: typically there is a way of scaling up a problem (in chapter 4 you will learn the concept of 'discretization'), and one tailors the size of the problem to the number of available processors.

A more realistic assumption would be to say that there is a sequential fraction independent of the problem size, and parallel fraction that can be arbitrarily replicated. To formalize this, instead of starting with the execution time of the sequential program, let us start with the execution time of the parallel program, and say that

$$T_p = T(F_s + F_p) \quad \text{with } F_s + F_p = 1.$$

Now we have two possible definitions of $T_1$. First of all, there is the $T_1$ you get from setting $p = 1$ in $T_p$. (Convince yourself that that is actually the same as $T_p$.) However, what we need is $T_1$ describing the time to do all the operations of the parallel program. This is:

$$T_1 = F_s T + p \cdot F_p T.$$

This gives us a speedup of

$$S_p = \frac{T_1}{T_p} = \frac{F_s + p \cdot F_p}{F_s + F_p} = F_s + p \cdot F_p = p - (p - 1) \cdot F_s.$$

That is, speedup is now a function that decreases from $p$, linearly with $p$.

As with Amdahl's law, we can investigate the behaviour of Gustafson's law if we include communication overhead. Let's go back to equation (2.2) for a perfectly parallelizable problem, and approximate it as

$$S_p = p(1 - \frac{T_c}{T_1} p).$$

Now, under the assumption of a problem that is gradually being scaled up, $T_c, T_1$ become functions of $p$. We see that if $T_1(p) \sim pT_c(p)$, we get linear speedup that is a constant fraction away from 1. In general we can not take this analysis further; in section 6.2.2 you'll see a detailed analysis of an example.
### 2.3.3.3 Amdahl's law and hybrid programming

Above, you learned about hybrid programming, a mix between distributed and shared memory programming. This leads to a new form of Amdahl's law.

Suppose we have $p$ nodes with $c$ cores each, and $F_p$ describes the fraction of the code that uses $c$-way thread parallelism. We assume that the whole code is fully parallel over the $p$ nodes. The ideal speed up would be $pc$, and the ideal parallel running time $T_1/(pc)$, but the actual running time is

$$T_{p,c} = T_1 \left( \frac{F_s}{p} + \frac{F_p}{pc} \right) = \frac{T_1}{pc} (F_s c + F_p) = \frac{T_1}{pc} (1 + F_s(c-1)).$$

**Exercise 2.10.** Show that the speedup $T_1/T_{p,c}$ can be approximated by $p/F_s$.

In the original Amdahl's law, speedup was limited by the sequential portion to a fixed number $1/F_s$, in hybrid programming it is limited by the task parallel portion to $p/F_s$.

### 2.3.4 Scalability

Above, we remarked that splitting a given problem over more and more processors does not make sense: at a certain point there is just not enough work for each processor to operate efficiently. Instead, in practice, users of a parallel code will either choose the number of processors to match the problem size, or they will solve a series of increasingly larger problems on correspondingly growing numbers of processors. In both cases it is hard to talk about speedup. Instead, the concept of *scalability* is used.

We distinguish two types of scalability. So-called *strong scalability* is in effect the same as speedup, discussed above. We say that a program shows strong scalability if, partitioned over more and more processors, it shows perfect or near perfect speedup. Typically, one encounters statements like ‘this problem scales up to 500 processors’, meaning that up to 500 processors the speedup will not noticeably decrease from optimal. It is not necessary for this problem to fit on a single processor: often a smaller number such as 64 processors is used as the baseline from which scalability is judged.

More interestingly, *weak scalability* is a more vaguely defined term. It describes that, as problem size and number of processors grow in such a way that the amount of data per processor stays constant, the speed in operations per second of each processor also stays constant. This measure is somewhat hard to report, since the relation between the number of operations and the amount of data can be complicated. If this relation is linear, one could state that the amount of data per processor is kept constant, and report that parallel execution time is constant as the number of processors grows.

**Exercise 2.11.** We can formulate strong scaling as a runtime that is inversely proportional to the number of processors:

$$t = c/p.$$

Show that on a log-log plot, that is, you plot the logarithm of the runtime against the logarithm of the number of processors, you will get a straight line with slope $-1$. Can you suggest a way of dealing with a non-parallelizable section, that is, with a runtime $t = c_1 + c_2/p$?
Exercise 2.12. Suppose you are investigating the weak scalability of a code. After running it for a couple of sizes and corresponding numbers of processes, you find that in each case the flop rate is roughly the same. Argue that the code is indeed weakly scalable.

Exercise 2.13. In the above discussion we always implicitly compared a sequential algorithm and the parallel form of that same algorithm. However, in section 2.3.1 we noted that sometimes speedup is defined as a comparison of a parallel algorithm with the best sequential algorithm for the same problem. With that in mind, compare a parallel sorting algorithm with runtime (log n)² (for instance, bitonic sort; section 8.1) to the best serial algorithm, which has a running time of n log n.

Show that in the weak scaling case of n = p speedup is p/ log p. Show that in the strong scaling case speedup is a descending function of n.

Although in industry parlance the term ‘scalability’ is sometimes applied to architectures or whole computer systems<fnref n="3" />, in scientific computing scalability is a property of an algorithm and the way it is parallelized on an architecture, in particular noting the way data is distributed. In section 6.2.2 you will find an analysis of the matrix-vector product operation: distributing a matrix by block rows turns out not to be scalable, but a two-dimensional distribution by submatrices is.

2.3.5 Simulation scaling

In most discussions of weak scaling we assume that the amount of work and the amount of storage are lin-
early related. This is not always the case; for instance the operation complexity of a matrix-matrix product
is N<fnref n="3" /> for N² data. If you linearly increase the number of processors, the work will go up with a higher
power.

A similar effect comes into play if you simulate time-dependent PDEs<fnref n="4" />. Here, the total work is a product of the work per time step and the number of time steps. These two numbers are related; in section 4.1.2 you saw that the time step has a certain minimum size as a function of the space discretization. Thus, the number of time steps will go up as the work per time step goes up.

Rather than investigating scalability from the point of the running of an algorithm, in this section we will look at the case where the simulated time *S* and the running time *T* are constant, and we look at how this influences the amount of memory we need. This corresponds to the following real-life scenario: you have a simulation that models a certain amount of real-world time in a certain amount of running time; now you buy a bigger computer, and you wonder what size problem you can solve in the same running time and maintaining the same simulated time.

Let *m* be the memory per processor, and *P* the number of processors, giving:

$M = Pm \quad \text{total memory.}$

<footnote n="3">3. “A scalable computer is a computer designed from a small number of basic components, without a single bottleneck component, so that the computer can be incrementally expanded over its designed scaling range, delivering linear incremental performance for a well-defined set of scalable applications. General-purpose scalable computers provide a wide range of processing, memory size, and I/O resources. Scalability is the degree to which performance increments of a scalable computer are linear” [8].</footnote>

<footnote n="4">4. This uses concepts from chapter 4.</footnote>
If $d$ is the number of space dimensions of the problem, typically 2 or 3, we get

$$ \Delta x = 1/M^{1/d} \quad \text{grid spacing.} $$

For stability this limits the time step $\Delta t$ to

$$ \Delta t = \begin{cases} \Delta x = 1 / M^{1/d} & \text{hyperbolic case} \\ \Delta x^2 = 1 / M^{2/d} & \text{parabolic case} \end{cases} $$

(noting that the hyperbolic case was not discussed in chapter 4.) With a simulated time $S$, we find

$$ k = S/\Delta t \quad \text{time steps.} $$

If we assume that the individual time steps are perfectly parallelizable, that is, we use explicit methods, or implicit methods with optimal solvers, we find a running time

$$ T = kM/P = \frac{S}{\Delta t}m. $$

Setting $T/S = C$, we find

$$ m = C\Delta t, $$

that is, the amount of memory per processor goes down as we increase the processor count. (What is the missing step in that last sentence?)

Further analyzing this result, we find

$$ m = C\Delta t = c \begin{cases} 1 / M^{1/d} & \text{hyperbolic case} \\ 1 / M^{2/d} & \text{parabolic case} \end{cases} $$

Substituting $M = Pm$, we find ultimately

$$ m = C \begin{cases} 1 / P^{1/(d+1)} & \text{hyperbolic} \\ 1 / P^{2/(d+2)} & \text{parabolic} \end{cases} $$

that is, the memory per processor that we can use goes down as a higher power of the number of processors.

### 2.3.6 Other scaling measures

Amdahl's law above was formulated in terms of the execution time on one processor. In many practical situations this is unrealistic, since the problems executed in parallel is too large for any single processor. Some formula manipulation gives us quantities that are to an extent equivalent, but that do not rely on this single-processor number [127].

For starters, applying the definition $S_p(n) = \frac{T_1(n)}{T_p(n)}$ to strong scaling, we observe that $T_1(n)/n$ is the sequential time per operation, and its inverse $n/T_1(n)$ can be called the sequential computational rate, denoted $R_1(n)$. Similarly defining a 'parallel computational rate'

$$ R_p(n) = n/T_p(n) \qquad (2.3) $$
<figure><img src="image_36.png" /><figcaption>Figure 2.4: The four classes of the Flynn’s taxonomy</figcaption></figure>

we find that

$$S_p(n) = R_p(n)/R_1(n)$$

In strong scaling $R_1(n)$ will be a constant, so we make a logarithmic plot of speedup, purely based on measuring $T_p(n)$.

## 2.4 Parallel Computers Architectures

For quite a while now, the top computers have been some sort of parallel computer, that is, an architecture that allows the simultaneous execution of multiple instructions or instruction sequences. One way of characterizing the various forms this can take is due to Flynn [57]. *Flynn’s taxonomy* characterizes architectures by whether the *data flow* and *control flow* are shared or independent. The following four types result (see also figure 2.4):

* **SISD** Single Instruction Single Data: this is the traditional CPU architecture: at any one time only a single instruction is executed, operating on a single data item.

* **SIMD** Single Instruction Multiple Data: in this computer type there can be multiple processors, each operating on its own data item, but they are all executing the same instruction on that data item. Vector computers (section 2.4.1.1) are typically also characterized as SIMD.

* **MISD** Multiple Instruction Single Data. No architectures answering to this description exist; one could argue that redundant computations for safety-critical applications are an example of MISD.

* **MIMD** Multiple Instruction Multiple Data: here multiple CPUs operate on multiple data items, each executing independent instructions. Most current parallel computers are of this type.

We will now discuss SIMD and MIMD architectures in more detail.

### 2.4.1 SIMD

Parallel computers of the SIMD type apply the same operation simultaneously to a number of data items. The design of the CPUs of such a computer can be quite simple, since the arithmetic unit does not need
separate logic and instruction decoding units: all CPUs execute the same operation in lock step. This makes SIMD computers excel at operations on arrays, such as

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

and, for this reason, they are also often called *array processors*. Scientific codes can often be written so that a large fraction of the time is spent in array operations.

On the other hand, there are operations that can not be executed efficiently on an array processor. For instance, evaluating a number of terms of a recurrence $x_{i+1} = ax_i + b_i$ involves that many additions and multiplications, but they alternate, so only one operation of each type can be processed at any one time. There are no arrays of numbers here that are simultaneously the input of an addition or multiplication.

In order to allow for different instruction streams on different parts of the data, the processor would have a 'mask bit' that could be set to prevent execution of instructions. In code, this typically looks like

```c
where (x>0) {
    x[i] = sqrt(x[i])
```

The programming model where identical operations are applied to a number of data items simultaneously, is known as *data parallelism*.

Such array operations can occur in the context of physics simulations, but another important source is graphics applications. For this application, the processors in an array processor can be much weaker than the processor in a PC: often they are in fact bit processors, capable of operating on only a single bit at a time. Along these lines, *ICL* had the 4096 processor *DAP* [92] in the 1980s, and *Goodyear* built a 16K processor *MPP* [7] in the 1970s.

Later, the *Connection Machine* (CM-1, CM-2, CM-5) were quite popular. While the first Connection Machine had bit processors (16 to a chip), the later models had traditional processors capable of floating point arithmetic, and were not true SIMD architectures. All were based on a hyper-cube interconnection network; see section 2.8.5. Another manufacturer that had a commercially successful array processor was *MasPar*; figure 2.5 illustrates the architecture. You clearly see the single control unit for a square array of processors, plus a network for doing global operations.

Supercomputers based on array processing do not exist anymore, but the notion of SIMD lives on in various guises. For instance, GPUs are SIMD-based, enforced through their *CUDA* programming language. Also, the *Intel Xeon Phi* has a strong SIMD component. While early SIMD architectures were motivated by minimizing the number of transistors necessary, these modern co-processors are motivated by power efficiency considerations. Processing instructions (known as *instruction*

<figure><img src="image_37.png" /></figure>
*issue*) is actually expensive compared to a floating point operation, in time, energy, and chip real estate needed. Using SIMD is then a way to economize on the last two measures.

### 2.4.1.1 Pipelining

A number of computers have been based on a *vector processor* or *pipeline processor* design. The first commercially successful supercomputers, the Cray-1 and the Cyber-205 were of this type. In recent times, the Cray-X1 and the NEC SX series have featured vector pipes. The ‘Earth Simulator’ computer [143], which led the TOP500 (section 2.11.5) for 3 years, was based on NEC SX processors. The general idea behind pipelining was described in section 1.2.1.3.

While supercomputers based on pipeline processors are in a distinct minority, pipelining is now mainstream in the superscalar CPUs that are the basis for *clusters*. A typical CPU has pipelined floating point units, often with separate units for addition and multiplication; see section 1.2.1.3.

However, there are some important differences between pipelining in a modern superscalar CPU and in, more old-fashioned, vector units. The pipeline units in these vector computers are not integrated floating point units in the CPU, but can better be considered as attached vector units to a CPU that itself has a floating point unit. The vector unit has *vector registers*<fnref n="5" /> with a typical length of 64 floating point numbers; there is typically no ‘vector cache’. The logic in vector units is also simpler, often addressable by explicit vector instructions. Superscalar CPUs, on the other hand, are fully integrated in the CPU and geared towards exploiting data streams in unstructured code.

### 2.4.1.2 True SIMD in CPUs and GPUs

True SIMD array processing can be found in modern CPUs and GPUs, in both cases inspired by the parallelism that is needed in graphics applications.

Modern CPUs from Intel and AMD, as well as PowerPC chips, have *vector instructions* that can perform multiple instances of an operation simultaneously. On Intel processors this is known as *SIMD Streaming Extensions (SSE)* or *Advanced Vector Extensions (AVX)*. These extensions were originally intended for graphics processing, where often the same operation needs to be performed on a large number of pixels. Often, the data has to be a total of, say, 128 bits, and this can be divided into two 64-bit reals, four 32-bit reals, or a larger number of even smaller chunks such as 4 bits.

The AVX instructions are based on up to 512-bit wide SIMD, that is, eight floating point numbers can be processed simultaneously. Just as single floating point operations operate on data in registers (section 1.3.3), vector operations use *vector registers*. The locations in a vector register are sometimes referred to as *SIMD lanes*.

The use of SIMD is mostly motivated by power considerations. Decoding instructions is actually more power consuming than executing them, so SIMD parallelism is a way to save power.

Current compilers can generate SSE or AVX instructions automatically; sometimes it is also possible for the user to insert pragmas, for instance with the Intel compiler:

<footnote n="5">5. The Cyber205 was an exception, with direct-to-memory architecture.</footnote>
```c
void func(float *restrict c, float *restrict a,
          float *restrict b, int n)
{
    #pragma vector always
    for (int i=0; i<n; i++)
        c[i] = a[i] * b[i];
}
```

Use of these extensions often requires data to be aligned with cache line boundaries (section 1.3.4.6), so there are special `allocate` and `free` calls that return aligned memory.

Version 4 of OpenMP also has directives for indicating SIMD parallelism.

Array processing on a larger scale can be found in GPUs. A GPU contains a large number of simple processors, ordered in groups of 32, typically. Each processor group is limited to executing the same instruction. Thus, this is true example of SIMD processing. For further discussion, see section 2.10.3.

### 2.4.2 MIMD / SPMD computers

By far the most common parallel computer architecture these days is called Multiple Instruction Multiple Data (MIMD): the processors execute multiple, possibly differing instructions, each on their own data. Saying that the instructions differ does not mean that the processors actually run different programs: most of these machines operate in Single Program Multiple Data (SPMD) mode, where the programmer starts up the same executable on the parallel processors. Since the different instances of the executable can take differing paths through conditional statements, or execute differing numbers of iterations of loops, they will in general not be completely in sync as they were on SIMD machines. If this lack of synchronization is due to processors working on different amounts of data, it is called load unbalance, and it is a major source of less than perfect speedup; see section 2.11.1.

There is a great variety in MIMD computers. Some of the aspects concern the way memory is organized, and the network that connects the processors. Apart from these hardware aspects, there are also differing ways of programming these machines. We will see all these aspects below. Many machines these days are called clusters. They can be built out of custom or commodity processors (if they consist of PCs, running Linux, and connected with Ethernet, they are referred to as Beowulf clusters [77]); since the processors are independent they are examples of the MIMD or SPMD model.

### 2.4.3 The commoditization of supercomputers

In the 1980s and 1990s supercomputers were radically different from personal computer and mini or super-mini computers such as the DEC PDP and VAX series. The SIMD vector computers had one (CDC Cyber205 or Cray-1), or at most a few (ETA-10, Cray-2, Cray X/MP, Cray Y/MP), extremely powerful processors, often a vector processor. Around the mid-1990s clusters with thousands of simpler (micro) processors started taking over from the machines with relative small numbers of vector pipes (see http://www.top500.org/lists/1994/11). At first these microprocessors (IBM Power series, Intel i860, MIPS, DEC Alpha) were still much more powerful than 'home computer' processors, but later
<figure><img src="image_38.png" /><figcaption>Figure 2.6: References to identically named variables in the distributed and shared memory case</figcaption></figure>

this distinction also faded to an extent. Currently, many of the most powerful clusters are powered by es-
sentially the same Intel Xeon and AMD Opteron chips that are available on the consumer market. Others
use IBM Power Series or other ‘server’ chips. See section 2.11.5 for illustrations of this history since 1993.

**2.5 Different types of memory access**

In the introduction we defined a parallel computer as a setup where multiple processors work together on the same problem. In all but the simplest cases this means that these processors need access to a joint pool of data. In the previous chapter you saw how, even on a single processor, memory can have a hard time keeping up with processor demands. For parallel machines, where potentially several processors want to access the same memory location, this problem becomes even worse. We can characterize parallel machines by the approach they take to the problem of reconciling multiple accesses, by multiple processes, to a joint pool of data.

The main distinction here is between *distributed memory* and *shared memory*. With distributed memory, each processor has its own physical memory, and more importantly its own *address space*. Thus, if two processors refer to a variable *x*, they access a variable in their own local memory. On the other hand, with shared memory, all processors access the same memory; we also say that they have a *shared address space*. Thus, if two processors both refer to a variable *x*, they access the same memory location.

**2.5.1 Symmetric Multi-Processors: Uniform Memory Access**

Parallel programming is fairly simple if any processor can access any memory location. For this reason,
there is a strong incentive for manufacturers to make architectures where processors see no difference
between one memory location and another: any memory location is accessible to every processor, and the
access times do not differ. This is called *Uniform Memory Access (UMA)*, and the programming model for
architectures on this principle is often called *Symmetric Multi Processing (SMP)*.

There are a few ways to realize an SMP architecture. Current desktop computers can have a few processors accessing a shared memory through a single memory bus; for instance Apple markets a model with 2 six-core processors. Having a memory bus that is shared between processors works only for small numbers of
<figure><img src="image_39.png" /><figcaption>Figure 2.7: Non-uniform memory access in a four-socket motherboard</figcaption></figure>

processors; for larger numbers one can use a *crossbar* that connects multiple processors to multiple memory banks; see section 2.8.6.

On *multicore* processors there is uniform memory access of a different type: the cores typically have a *shared cache*, typically the L3 or L2 cache.

### 2.5.2 Non-Uniform Memory Access

The NUMA approach based on shared memory is obviously limited to a small number of processors. The crossbar networks are expandable, so they would seem the best choice. However, in practice one puts processors with a local memory in a configuration with an exchange network. This leads to a situation where a processor can access its own memory fast, and other processors' memory slower. This is one case of so-called *Non-Uniform Memory Access (NUMA)*: a strategy that uses physically distributed memory, abandoning the uniform access time, but maintaining the logically shared address space: each processor can still access any memory location.

Figure 2.7 illustrates NUMA in the case of the four-socket motherboard of the Ranger supercomputer. Each chip has its own memory (8Gb) but the motherboard acts as if the processors have access to a shared pool of 32Gb. Obviously, accessing the memory of another processor is slower than accessing local memory. In addition, note that each processor has three connections that could be used to access other memory, but the rightmost two chips use one connection to connect to the network. This means that accessing each other's memory can only happen through an intermediate processor, slowing down the transfer, and tying up that processor's connections.

While the NUMA approach is convenient for the programmer, it offers some challenges for the system. Imagine that two different processors each have a copy of a memory location in their local (cache) memory. If one processor alters the content of this location, this change has to be propagated to the other processors. If both processors try to alter the content of the one memory location, the behaviour of the program can become undetermined.
Keeping copies of a memory location synchronized is known as *cache coherence* (see section 1.4.1 for further details); a multi-processor system using it is sometimes called a ‘cache-coherent NUMA’ or *ccNUMA* architecture.

Taking NUMA to its extreme, it is possible to have a software layer that makes network-connected processors appear to operate on shared memory. This is known as *distributed shared memory* or *virtual shared memory*. In this approach a *hypervisor* offers a shared memory API, by translating system calls to distributed memory management. This shared memory API can be utilized by the *Linux kernel*, which can support 4096 threads.

Among current vendors only SGI (the *UV line*) and Cray (the *XE6*) market products with large scale NUMA. Both offer strong support for *Partitioned Global Address Space (PGAS) languages*; see section 2.7.5. There are vendors, such as *ScaleMP*, that offer a software solution to distributed shared memory on regular clusters.

### 2.5.3 Logically and physically distributed memory

The most extreme solution to the memory access problem is to offer memory that is not just physically, but that is also logically distributed: the processors have their own address space, and can not directly see another processor’s memory. This approach is often called ‘distributed memory’, but this term is unfortunate, since we really have to consider the questions separately whether memory *is* distributed and whether *it appears* distributed. Note that NUMA also has physically distributed memory; the distributed nature of it is just not apparent to the programmer.

With *logically and physically distributed memory*, the only way one processor can exchange information with another is through passing information explicitly through the network. You will see more about this in section 2.7.3.3.

This type of architecture has the significant advantage that it can scale up to large numbers of processors: the *IBM BlueGene* has been built with over 200,000 processors. On the other hand, this is also the hardest kind of parallel system to program.

Various kinds of hybrids between the above types exist. In fact, most modern clusters will have NUMA nodes, but a distributed memory network between nodes.

## 2.6 Granularity of parallelism

Let us take a look at the question ‘how much parallelism is there in a program execution’. There is the theoretical question of the absolutely maximum number of actions that can be taken in parallel, but we also need to wonder what kind of actions these are and how hard it is to actually execute them in parallel, as well has how efficient the resulting execution is.

The discussion in this section will be mostly on a conceptual level; in section 2.7 we will go into some detail on how parallelism can actually be programmed.
2.6.1 Data parallelism

It is fairly common for a program that have loops with a simple body, that gets executed for all elements in a large data set:

```
for (i=0; i<1000000; i++)
    a[i] = 2*b[i];
```

Such code is considered an instance of *data parallelism* or *fine-grained parallelism*. If you had as many processors as array elements, this code would look very simple: each processor would execute the statement

```
a = 2*b
```

on its local data.

If your code consists predominantly of such loops over arrays, it can be executed efficiently with all processors in lockstep. Architectures based on this idea, where the processors can in fact *only* work in lockstep, have existed, see section 2.4.1. Such fully parallel operations on arrays appear in computer graphics, where every pixel of an image is processed independently. For this reason, GPUs (section 2.10.3) are strongly based on *data parallelism*.

Continuing the above example for a little bit, consider the operation

```
for 0 <= i < max do
    i left = mod (i - 1, max)
    i right = mod (i + 1, max)
    a_i = (b_{i left} + b_{i right})/2
```

On a data parallel machine, that could be implemented as

```
bleft  left  shiftright(b)
bright right  shiftright(b)
a      left  (bleft + bright)/2
```

where the `shiftright/right` instructions cause a data item to be sent to the processor with a number lower or higher by 1. For this second example to be efficient, it is necessary that each processor can communicate quickly with its immediate neighbours, and the first and last processor with each other.

In various contexts such a 'blur' operations in graphics, it makes sense to have operations on 2D data:

```
for 0 < i < m do
    for 0 < j < n do
        a_{ij} = (b_{ij-1} + b_{ij+1} + b_{i-1j} + b_{i+1j})
```

and consequently processors have be able to move data to neighbours in a 2D grid.
## 2.6.2 Instruction-level parallelism

In ILP, the parallelism is still on the level of individual instructions, but these need not be similar. For instance, in

$$a \leftarrow b + c$$

$$d \leftarrow e * f$$

the two assignments are independent, and can therefore be executed simultaneously. This kind of parallelism is too cumbersome for humans to identify, but compilers are very good at this. In fact, identifying ILP is crucial for getting good performance out of modern superscalar CPUs.

## 2.6.3 Task-level parallelism

At the other extreme from data and instruction-level parallelism, task parallelism is about identifying whole subprograms that can be executed in parallel. As an example, searching in a tree data structure could be implemented as follows:

```
if optimal (root) then
    exit
```

```
else
    parallel: SearchInTree (leftchild),SearchInTree (rightchild)
        Procedure SearchInTree(root)
```

The search tasks in this example are not synchronized, and the number of tasks is not fixed: it can grow arbitrarily. In practice, having too many tasks is not a good idea, since processors are most efficient if they work on just a single task. Tasks can then be scheduled as follows:

```
while there are tasks left do
    wait until a processor becomes inactive;
    spawn a new task on it
```

(There is a subtle distinction between the two previous pseudo-codes. In the first, tasks were self-scheduling: each task spawned off two new ones. The second code is an example of the master-worker paradigm: there is one central task which lives for the duration of the code, and which spawns and assigns the worker tasks.)

Unlike in the data parallel example above, the assignment of data to processor is not determined in advance in such a scheme. Therefore, this mode of parallelism is most suited for thread-programming, for instance through the OpenMP library; section 2.7.2.

Let us consider a more serious example of task-level parallelism.

A finite element mesh is, in the simplest case, a collection of triangles that covers a 2D object. Since angles that are too acute should be avoided, the Delauney mesh refinement process can take certain triangles, and replace them by better shaped ones. This is illustrated in figure 2.8: the black triangles violate some angle condition, so either they themselves get subdivided, or they are joined with some neighbouring ones (rendered in grey) and then jointly redivided.

In pseudo-code, this can be implemented as in figure 2.9. (This figure and code are to be found in [103],
<figure><img src="image_40.png" /><figcaption>Figure 2.8: A mesh before and after refinement</figcaption></figure>

```c
Mesh m = /* read in initial mesh */
WorkList wl;
wl.addeditor(badTriangles());

while (wl.size() != 0) do
    Element e = wl.get(); //get bad triangle
    if (e no longer in mesh) continue;
    Cavity c = new Cavity(e);
    c.expand();
    c.retriangulate();
    mesh.update(c);
    wl.add(badTriangles());
```

<caption>Figure 2.9: Task queue implementation of Delauney refinement</caption>

which also contains a more detailed discussion.)

It is clear that this algorithm is driven by a worklist (or *task queue*) data structure that has to be shared between all processes. Together with the dynamic assignment of data to processes, this implies that this type of *irregular parallelism* is suited to shared memory programming, and is much harder to do with distributed memory.

### 2.6.4 Conveniently parallel computing

In certain contexts, a simple, often single processor, calculation needs to be performed on many different inputs. Since the computations have no data dependencies and need not done in any particular sequence, this is often called *embarassingly parallel* or *conveniently parallel* computing. This sort of parallelism can happen at several levels. In examples such as calculation of the *Mandelbrot set* or evaluating moves in a *chess game*, a subroutine-level computation is invoked for many parameter values. On a coarser level it can be the case that a simple program needs to be run for many inputs. In this case, the overall calculation is referred to as a *parameter sweep*.

### 2.6.5 Medium-grain data parallelism

The above strict realization of data parallelism assumes that there are as many processors as data elements. In practice, processors will have much more memory than that, and the number of data elements is likely to be far larger than the processor count of even the largest computers. Therefore, arrays are grouped onto processors in subarrays. The code then looks like this:
```c
my_lower_bound = // some processor-dependent number
my_upper_bound = // some processor-dependent number
for (i=my_lower_bound; i<my_upper_bound; i++)
    // the loop body goes here
```

This model has some characteristics of data parallelism, since the operation performed is identical on a large number of data items. It can also be viewed as task parallelism, since each processor executes a larger section of code, and does not necessarily operate on equal sized chunks of data.

### 2.6.6 Task granularity

In the previous subsections we considered different level of finding parallel work, or different ways of dividing up work so as to find parallelism. There is another way of looking at this: we define the *granularity* of a parallel scheme as the amount of work (or the task size) that a processing element can perform before having to communicate or synchronize with other processing elements.

In ILP we are dealing with very fine-grained parallelism, on the order of a single instruction or just a few instructions. In true task parallelism the granularity is much coarser.

The interesting case here is data parallelism, where we have the freedom to choose the task sizes. On SIMD machines we can choose a granularity of a single instruction, but, as you saw in section 2.6.5, operations can be grouped into medium-sized tasks. Thus, operations that are data parallel can be executed on distributed memory clusters, given the right balance between the number of processors and total problem size.

Exercise 2.14. Discuss choosing the right granularity for a data parallel operation such as averaging on a two-dimensional grid. Show that there is a *surface-to-volume* effect: the amount of communication is of a lower order than the computation. This means that, even if communication is much slower than computation, increasing the task size will still give a balanced execution.

Unfortunately, choosing a large task size to overcome slow communication may aggregate another problem: aggregating these operations may give tasks with varying running time, causing *load unbalance*. One solution here is to use an *overdecomposition* of the problem: create more tasks then there are processing elements, and assign multiple tasks to a processor (or assign tasks dynamically) to even out irregular running times. This is known as *dynamic scheduling*, and the examples in section 2.6.3 illustrate this; see also section 2.7.2.1. An example of *overdecomposition* in linear algebra is discussed in section 6.3.2.

## 2.7 Parallel programming

Parallel programming is more complicated than sequential programming. While for sequential programming most programming languages operate on similar principles (some exceptions such as functional or logic languages aside), there is a variety of ways of tackling parallelism. Let's explore some of the concepts and practical aspects.

There are various approaches to parallel programming. First of all, there does not seem to be any hope of a *parallelizing compiler* that can automagically transform a sequential program into a parallel one. Apart