Alternatively, we can use recursive halving :

for (s=(n+1)/2; s>1; s/=2)
  for (i=0; i<n; i+=1)
    x[i] += x[i+s]

Cache memory 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, pages, and virtual memory

crumb trail: > sequential > Memory Hierarchies > TLB, pages, 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.

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.

1.3.8.1 Large pages

crumb trail: > sequential > Memory Hierarchies > TLB, pages, and virtual memory > Large pages

In very irregular applications, for instance databases, the page table can get very large as more-or-less random data is brought into memory. However, sometimes these pages show some amount of clustering, meaning that if the page size had been larger, the number of needed pages would be greatly reduced. For this reason, operating systems can have support for large pages  , typically of size around 2$\,$Mb. (Sometimes `huge pages' are used; for instance the Intel Knights Landing has Gigabyte pages.)

The benefits of large pages are application-dependent: if the small pages have insufficient clustering, use of large pages may fill up memory prematurely with the unused parts of the large pages.

1.3.8.2 TLB

crumb trail: > sequential > Memory Hierarchies > TLB, pages, and virtual memory > TLB

However, address translation by lookup in this table is slow, so CPUs have a 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.10  ), using an LRU policy (section  1.3.4.6  ).

A typical TLB has between 64 and 512 entries. If a program accesseses 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 accesses many random memory locations can experience a slowdown because of such misses. The set of pages that is in current use is called the `working set'.

Section  1.7.6 and appendix  sec:tlb-code discuss some simple code illustrating the behavior 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.]

Use of large pages also reduces the number of potential TLB misses, since the working set of pages can be reduced.

1.4 Multicore architectures

crumb trail: > sequential > 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.8.1  .
• It is not possible to extract more 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'. The separate cores can work on unrelated tasks, or, by introducing what is in effect data parallelism (section  2.3.1  ), collaborate on a common task at a higher overall efficiency [Olukotun:1996:single-chip]  .

Remark 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.9.3  . For a discussion see section  2.6.1.9  .
End of remark

This solves the above two problems:

• 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  .

FIGURE 1.12: Cache hierarchy in a single-core and dual-core chip.

This design makes it 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.

QUOTE 1.13: Structure of the Intel Single-chip Cloud Computer chip.

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.
• [ 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.4.1  ) such as the Sequent Symmetry and the \indexterm{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

crumb trail: > sequential > Multicore architectures > 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, processors typically have some private cache that 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 \indexterm{atomic operation}; more about this in section  2.6.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 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; see below.) 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 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 FSA ; see section  app:fsa  .)
End of exercise

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

1.4.1.1 Solutions to cache coherence

crumb trail: > sequential > Multicore architectures > Cache coherence > Solutions to cache coherence

There are two basic mechanisms for realizing cache coherence: snooping and directory-based schemes.

In the snooping mechanism, any request for data is sent to all caches, and the data is returned if it is present anywhere; otherwise it is retrieved from memory. In a variation on this scheme, a core `listens in' on all bus traffic, so that 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 When would updating pay off? Write a simple cache simulator to evaluate this question.
End of exercise

Since snooping often involves broadcast information to all cores, it does not scale beyond a small number of cores. A solution that scales better is using a tag directory : a central directory that contains the information on what data is present in some cache, and what cache it is in specifically. For processors with large numbers of cores (such as the Intel Xeon Phi  ) the directory can be distributed over the cores.

1.4.1.2 Tag directories

crumb trail: > sequential > Multicore architectures > Cache coherence > Tag directories

In multicore processors with distributed, but coherent, caches (such as the Intel Xeon Phi  ) the tag directories can themselves be distributed. This increases the latency of cache lookup.

1.4.2 False sharing

crumb trail: > sequential > Multicore architectures > False sharing

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

  double x,y;

The most common case of false sharing happens when threads update consecutive locations of an array. For instance, in the following OpenMP fragment all threads update their own location in an array of partial results:

  local_results = new double[num_threads];
#pragma omp parallel
{
  int thread_num = omp_get_thread_num();
  for (int i=my_lo; i<my_hi; i++)
    local_results[thread_num] = ... f(i) ...
}
global_result = g(local_results)

That said, false sharing is less of a problem in modern CPUs than it used to be. Let's consider a code that explores the above scheme:

// falsesharing-omp.cxx
  for (int spacing : {16,12,8,4,3,2,1,0} ) {
      int iproc = omp_get_thread_num();
      floattype *write_addr = results.data() + iproc*spacing;
      for (int r=0; r<how_many_repeats; r++) {
	for (int w=0; w<stream_length; w++) {
	  *write_addr += *( read_stream+w );
	}
      }

Executing this on the Intel Core i5 of an Apple Macbook Air shows a modest performance degradation: \lstinputlisting[basicstyle=\ttfamily\scriptsize]{code/hardware/falsesharing-macair.runout}

On the other hand, on the Intel Cascade Lake of the TACC Frontera cluster we see no such thing: \lstinputlisting[basicstyle=\ttfamily\scriptsize]{code/hardware/falsesharing-shared.runout}

The reason is that here the hardware caches the accumulator variable, and does not write to memory until the end of the loop. This obviates all problems with false sharing.

We can force false sharing problems by forcing the writeback to memory, for instance with an OpenMP atomic directive: \lstinputlisting[basicstyle=\ttfamily\scriptsize]{code/hardware/falsesharing-atomic.runout}

(Of course, this has a considerable performance penalty by itself.)

1.4.3 Computations on multicore chips

crumb trail: > sequential > Multicore architectures > Computations on multicore chips

There are various ways that a multicore processor can lead to increased performance. First of all, in a desktop situation, multiple cores can actually run multiple programs. More importantly, we can use the parallelism to speed up the execution of a single code. This can be done in two different ways.

The MPI library (section  2.6.3.3  ) is typically used to communicate between processors that are connected through a network. However, it can also be used in a single multicore processor: the MPI calls then are realized through shared memory copies.

Alternatively, we can use the shared memory and shared caches and program using threaded systems such as OpenMP (section  2.6.2  ). The advantage of this mode is that parallelism can be much more dynamic, since the runtime system can set and change the correspondence between threads and cores during the program run.

We will discuss in some detail the scheduling of linear algebra operations on multicore chips; section  6.11  .

1.4.4 TLB shootdown

crumb trail: > sequential > Multicore architectures > TLB shootdown

Section  1.3.8.2 explained how the TLB is used to cache the translation from logical address, and therefore logical page, to physical page. The TLB is part of the memory unit of the socket  , so in a multi-socket design, it is possible for a process on one socket to change the page mapping, which makes the mapping on the other incorrect.

One solution to this problem is called TLB shoot-down : the process changing the mapping generates an Inter-Processor Interrupt  , which causes the other processors to rebuild their TLB.

1.5 Node architecture and sockets

crumb trail: > sequential > Node architecture and sockets

In the previous sections we have made our way down through the memory hierarchy, visiting registers and various cache levels, and the extent to which they can be private or shared. At the bottom level of the memory hierarchy is the memory that all cores share. This can range from a few Gigabyte on a lowly laptop to a few Terabyte in some supercomputer centers.

While this memory is shared between all cores, there is some structure to it. This derives from the fact that a cluster node can have more than one socket  , that is, processor chip.

\hbox{ }

FIGURE 1.14: Left: a four-socket design. Right: a two-socket design with co-processor.

The shared memory on the node is typically spread over banks that are directly attached to one particular socket. This is for instance illustrated in figure  1.14  , which shows the four-socket node of the TACC Ranger cluster supercomputer (no longer in production) and the two-socket node of the TACC Stampede cluster supercomputer which contains an Intel Xeon Phi co-processor. In both designs you clearly see the memory chips that are directly connected to the sockets.

1.5.1 Design considerations

crumb trail: > sequential > Node architecture and sockets > Design considerations

Consider a supercomputer cluster to be built-up out of $N$ nodes, each of $S$ sockets, each with $C$ cores. We can now wonder, `if $S>1$, why not lower $S$, and increase $N$ or $C$?' Such questions have answers as much motivated by price as by performance.

• Increasing the number of cores per socket may run into limitations of chip fabrication  .
• Increasing the number of cores per socket may also lower the available bandwidth per core. This is a consideration for HPC applications where cores are likely to be engaged in the same activity; for more heterogeneous workloads this may matter less.
• For fixed $C$, lowering $S$ and increasing $N$ mostly raises the price because the node price is a weak function of $S$.
• On the other hand, raising $S$ makes the node architecture more complicated, and may lower performance because of increased complexity of maintaining coherence (section  1.4.1  ). Again, this is less of an issue with heterogeneous workloads, so high values of $S$ are more common in web servers than in HPC installations.

1.5.2 NUMA phenomena

crumb trail: > sequential > Node architecture and sockets > NUMA phenomena

The nodes illustrated above are examples of NUMA design: for a process running on some core, the memory attached to its socket is slightly faster to access than the memory attached to another socket.

One result of this is the first-touch phenomenon. Dynamically allocated memory is not actually allocated until it's first written to. Consider now the following OpenMP (section  2.6.2  ) code:

double *array = (double*)malloc(N*sizeof(double));
for (int i=0; i<N; i++)
   array[i] = 1;
#pragma omp parallel for
for (int i=0; i<N; i++)
   .... lots of work on array[i] ...

The solution here is to also make the initialization loop parallel, even if the amount of work in it may be negligible.

1.6 Locality and data reuse

crumb trail: > sequential > 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 data often interacts mostly with data that is close by in some sense, which will lead to data locality; section  1.6.2  . Often such locality derives from the nature of the application, as in the case of the PDEs you will see in chapter  Numerical treatment of differential equations  . In other cases such as molecular dynamics (chapter  app:md  ) there is no such intrinsic locality because all particles interact with all others, and considerable programming cleverness is needed to get high performance.

1.6.1 Data reuse and arithmetic intensity

crumb trail: > sequential > Locality and data reuse > 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.6.1.1 Example: vector operations

crumb trail: > sequential > Locality and data reuse > Data reuse and arithmetic intensity > Example: vector operations

Consider for example the vector addition \begin{equation} \forall_i\colon x_i\leftarrow x_i+y_i. \end{equation} This involves three memory accesses (two loads and one store) and one operation per iteration, giving an arithmetic intensity of $1/3$.

The BLAS axpy (for ` a  times x plus  y  ) operation \begin{equation} \forall_i\colon x_i\leftarrow a\,x_i+ y_i \end{equation} 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 \begin{equation} \forall_i\colon s\leftarrow s+x_i\cdot y_i \end{equation} 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$.

1.6.1.2 Example: matrix operations

crumb trail: > sequential > Locality and data reuse > Data reuse and arithmetic intensity > Example: matrix operations

Next, consider the matrix-matrix product : \begin{equation} \forall_{i,j}\colon c_{ij} = \sum_k a_{ik}b_{kj}. \end{equation} 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 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?
End of exercise

[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  [St:gaussnotoptimal,Pa:combinations]  . 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 \indextermbus{LINPACK} {benchmark}  [Dongarra1987LinpackBenchmark] ; see section  2.11.6  . Using this as the sole 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.6.1.3 The roofline model

crumb trail: > sequential > Locality and data reuse > Data reuse and arithmetic intensity > 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   [Williams:2009:roofline]  , and it expresses the basic fact that performance is bounded by two factors, illustrated in the first graph of figure  1.15  .

FIGURE 1.15: Illustration of factors determining performance in the roofline model.

1. The peak performance  , indicated by the horizontal line at the top of the graph, is an absolute bound on the performance\footnote {An old joke states that the peak performance is that number that the manufacturer guarantees you will never exceed.}, 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: \begin{equation} \frac{\hbox{\it operations}}{\hbox{\it second}}= \frac{\hbox{\it operations}}{\hbox{\it data item}}\cdot \frac{\hbox{\it data items}}{\hbox{\it second}} \end{equation} This is depicted by the linearly increasing line in the 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 How would you determine whether a given program kernel is bandwidth or compute bound?
End of exercise

1.6.2 Locality

crumb trail: > sequential > Locality and data reuse > 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 (low level memory access can be controlled by the programmer in the Cell processor and in some GPUs.), in most CPUs  , it is still possible, knowing the behavior 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 \index{temporal locality|see{locality, temporal}} and spatial locality \index{spatial locality|see{locality, spatial}}. 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 an LRU replacement policy (section  1.3.4.6  ), 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.

1.6.2.1 Temporal locality

crumb trail: > sequential > Locality and data reuse > Locality > Temporal locality

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

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

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

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

1.6.2.2 Spatial locality

crumb trail: > sequential > Locality and data reuse > Locality > 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 line s 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

  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\ldots\mathtt{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 $\mathtt{s}=1$, we load $1/\mathtt{S}$ cachelines per element; if $\mathtt{s}=\mathtt{S}$, we load one cacheline for each element. This effect is demonstrated in section  1.7.5  .

• A second example of spatial locality worth observing involves the TLB (section  1.3.8.2  ). 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.7.6  .

Exercise 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]

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

Exercise 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]

• 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.6.2.3 Examples of locality

crumb trail: > sequential > Locality and data reuse > Locality > Examples of locality

Let us examine locality issues for a realistic example. The matrix-matrix multiplication $C\leftarrow A\cdot 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 j=1..n
    for k=1..n
      c[i,j] += a[i,k]*b[k,j]

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

These implementations are illustrated in figure  1.6.2.3

\caption{Two loop orderings for the $C\leftarrow A\cdot B$ matrix-matrix product.}

The first implementation constructs the $(i,j)$ element 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 behavior, 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 in 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.

Exercise 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]

1.6.2.4 Core locality

crumb trail: > sequential > Locality and data reuse > Locality > 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  . This comes into play with multicore, multi-threaded programs.

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. Core locality is not just a property of a program, but also to a large extent of how the program is executed in parallel.

The following issues are at play here.

• There are performance implications to cache coherence (section  1.4.1  ) if 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 will be affected if the OS is allowed do thread migration  , thereby invalidating cache contents. Matters of fixing thread affinity are discussed in \PCSEref[chapter]{ch:omp-affinity} and \PCSEref[chapter]{ch:hybrid}.
• In multi- socket systems, there is also the phenomenon of first touch : data gets allocated on the socket where it is first initialized. This means that access from the other socket(s) may be considerably slower. See \PCSEref{sec:first-touch}.

1.7 Programming strategies for high performance

crumb trail: > sequential > 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.

The full listings of the codes and explanations of the data graphed here can be found in chapter  app:codes  .

1.7.1 Peak performance

crumb trail: > sequential > Programming strategies for high performance > 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.6  .

1.7.2 Bandwidth

crumb trail: > sequential > Programming strategies for high performance > Bandwidth

You have seen in section  1.6.1.3 that algorithms can be bandwidth-bound  , meaning that they are more limited by the available bandwidth than by available processing power. In that case we can wonder if the amount of bandwidth per core is a function of the number of cores: it is conceivable that the total bandwidth is less than the product of the number of cores and the bandwidth to a single core.

We test this on TACC's Frontera cluster, with dual-socket Intel Cascade Lake processors, with a total of 56 cores per node.

We let each core execute a simple streaming kernel:

// allocation.cxx
template <typename R>
R Cache<R>::sumstream(int repeats,size_t length,size_t byte_offset /* =0 default */ ) const {
  R s{0};
  size_t loc_offset = byte_offset/sizeof(R);
  const R *start_point = thecache.data()+loc_offset;
  for (int r=0; r<repeats; r++)
    for (int w=0; w<length; w++)
      s += *( start_point+w ) * r;
  return s;
};

// bandwidth.cxx
    vector<double> results(nthreads,0.);
    for ( int t=0; t<nthreads; t++) {
      auto start_point = t*stream_length;
      threads.push_back
	( thread( [=,&results] () {
	    results[t] = memory.sumstream(how_many_repeats,stream_length,start_point);
	  } ) );
    }
    for ( auto &t : threads )
      t.join();

We see that for 40 threads there is essentially no efficiency loss, but with 56 threads there is a 10 percent loss.

OMP_PROC_BIND=true ./bandwidth -t 56 -s 8 -k 64
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                Processing 8192 words
                over up to 56 threads
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Overhead for 1 threads:       66 usec
Threads  1.. ran for   136965 usec (efficiency: 100% )
Overhead for 8 threads:      187 usec
Threads  8.. ran for   137492 usec (efficiency: 99% )
Overhead for 16 threads:      260 usec
Threads 16.. ran for   137642 usec (efficiency: 99% )
Overhead for 24 threads:      406 usec
Threads 24.. ran for   137930 usec (efficiency: 99% )
Overhead for 32 threads:      559 usec
Threads 32.. ran for   138183 usec (efficiency: 99% )
Overhead for 40 threads:      793 usec
Threads 40.. ran for   138405 usec (efficiency: 98% )
Overhead for 48 threads:     1090 usec
Threads 48.. ran for   143611 usec (efficiency: 95% )
Overhead for 56 threads:     1191 usec
Threads 56.. ran for   149732 usec (efficiency: 91% )

1.7.3 Pipelining

crumb trail: > sequential > Programming strategies for high performance > 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

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

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

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

  a += 2; b += 2;
}

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.

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;
}

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

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

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

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.

1.7.3.1 Semantics of unrolling

crumb trail: > sequential > Programming strategies for high performance > Pipelining > Semantics of unrolling

An observation about the unrollling transformation 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  Computer Arithmetic  , in computer arithmetic this is not guaranteed to give the exact same result.

For this reason, a compiler will only apply this transformation if explicitly allowed. For instance, for the Intel compiler  , the option fp-model precise indicates that code transformations should preserve the semantics of floating point arithmetic, and unrolling is not allowed. On the other hand fp-model fast indicates that floating point arithmetic can be sacrificed for speed.

\lstinputlisting{code/hardware/unroll-clx.runout}

1.7.4 Cache size

crumb trail: > sequential > Programming strategies for high performance > 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 accesses the same data:

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

Exercise Argue that with a large enough problem and an LRU replacement policy (section  1.3.4.6  ) 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?
End of exercise

Often, it is possible to arrange the operations to keep data in L1 cache. For instance, in our example, we could write

for (b=0; b<size/l1size; b++) {
  blockstart = 0;
  for (i=0; i<NRUNS; i++) {
    for (j=0; j<l1size; j++)
      array[blockstart+j] = 2.3*array[blockstart+j]+1.2; 
  }
  blockstart += l1size;
}

Remark Like unrolling, blocking code may change the order of evaluation of expressions. Since floating point arithmetic is not associative % , blocking is not a transformation that compilers are allowed to make.
End of remark

1.7.4.1 Measuring cache performance at user level

crumb trail: > sequential > Programming strategies for high performance > Cache size > Measuring cache performance at user level

You can try to write a loop over a small array as above, and execute it many times, hoping to observe performance degradation when the array gets larger than the cache size. The number of iterations should be chosen so that the measured time is well over the resolution of your clock.

This runs into a couple of unforesoon problems. The timing for a simple loop nest

for (int irepeat=0; irepeat<how_many_repeats; irepeat++) {
  for (int iword=0; iword<cachesize_in_words; iword++)
    memory[iword] += 1.1;
}    

To see what the compiler does with this fragment, let your compiler generate an optimization report For the Intel compiler use -qopt-report  . In this report you see that the compiler has decided to exchange the loops: Each element of the array is then loaded only once.

remark #25444: Loopnest Interchanged: ( 1 2 ) --> ( 2 1 )
remark #15542: loop was not vectorized: inner loop was already vectorized

In an attempt to prevent this loop exchange, you can try to make the inner loop more too complicated for the compiler to analyze. You could for instance turn the array into a sort of linked list that you traverse:

// setup
for (int iword=0; iword<cachesize_in_words; iword++)
    memory[iword] = (iword+1) % cachesize_in_words

// use:
ptr = 0
for (int iword=0; iword<cachesize_in_words; iword++)
    ptr = memory[ptr];

To stymie this bit of cleverness you need to make the linked list more random:

for (int iword=0; iword<cachesize_in_words; iword++)
    memory[iword] = random() % cachesize_in_words

The code for this is given in section  sec:cachesize-code  .

Exercise While the strategy just sketched will demonstrate the existence of cache sizes, it will not report the maximal bandwidth that the cache supports. What is the problem and how would you fix it?
End of exercise

1.7.4.2 Detailed timing

crumb trail: > sequential > Programming strategies for high performance > Cache size > Detailed timing

If we have access to a cycle-accurate timer or the hardware counters, we can actually plot the number of cycles per access.

Such a plot is given in figure  1.7.4.2  .

\caption{Average cycle count per operation as function of the dataset size.}

The full code is given in section  sec:cachesize-code  .

1.7.5 Cache lines and striding

crumb trail: > sequential > Programming strategies for high performance > Cache lines and striding

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

for (i=0,n=0; i<L1WORDS; i++,n+=stride)
  array[n] = 2.3*array[n]+1.2;

QUOTE 1.16: Run time in kcycles and L1 reuse as a function of stride.

the graph in figure  1.16  .

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  ).

1.16 only plots up to stride 8, while 1.16  .

The full code is given in section  sec:cacheline-code  .

The effects of striding can be mitigated by the bandwidth and cache behavior of a processor. Consider some run on the Intel Cascade Lake processor of the Frontera cluster at TACC, (28-cores per socket, dual socket, for a total of 56 cores per node). We measure the time-per-operation on a simple streaming kernel, using increasing strides. Table  1.16 reports in the second column indeed a per-operation time that goes up linearly with the stride.

\toprule stride	nsec/word
	56 cores, 3M	56 cores, 0.3M	28 cores, 3M
\midrule 1	7.268	1.368	1.841
2	13.716	1.313	2.051
3	20.597	1.319	2.852
4	27.524	1.316	3.259
5	34.004	1.329	3.895
6	40.582	1.333	4.479
7	47.366	1.331	5.233
8	53.863	1.346	5.773
\bottomrule

\caption{Time per operation in nanoseconds as a function of striding on 56 cores of Frontera, per-core datasize 3.2M.}

However, this is for a dataset that overflows the L2 cache. If we make this run contained in the L2 cache, as reported in the 3rd column, this increase goes away as there is enough bandwdidth available to stream strided data at full speed from L2 cache.

TIKZPICTURE 1.17: Nano second per access for various core and stride counts.

1.7.6 TLB

crumb trail: > sequential > Programming strategies for high performance > TLB

As explained in section  1.3.8.2  , the TLB maintains a small list of frequently used memory pages and their locations; addressing data that are location 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.

#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 */
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  sec:tlb-code for the source code) are plotted in figures 1.7.6 and  1.7.6  .

\caption{Number of TLB misses per column as function of the number of columns; columnwise traversal of the array.}

\caption{Number of TLB misses per column as function of the number of columns; rowwise traversal of the array.}

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.7.6 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.7.6 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.7.7 Cache associativity

crumb trail: > sequential > Programming strategies for high performance > Cache associativity

There are many algorithms that work by recursive division of a problem, for instance the FFT algorithm. As a result, code for such algorithms will often operate on vectors whose 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.9 you saw how the operation of adding a small number of vectors \begin{equation} \forall_j\colon y_j= y_j+\sum_{i=1}^mx_{i,j} \end{equation} 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,\ldots$.

\caption{The number of L1 cache misses and the number of cycles for each $j$ column accumulation, vector length $4096$.}

\caption{The number of L1 cache misses and the number of cycles for each $j$ column accumulation, vector length $4096+8$.}

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.7.7 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.7.7  , 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.7.7 we see a decreasing time with increasing $m$; this is probably due to a progressively more favorable balance between load and store operations. Store operations are more expensive than loads, for various reasons.

1.7.8 Loop nests

crumb trail: > sequential > Programming strategies for high performance > Loop nests

If your code has nested loops  , and the iterations of the outer loop are independent, you have a choice which loop to make outer and which to make inner.

Exercise Give an example of a doubly-nested loop where the loops can be exchanged; give an example where this can not be done. If at all possible, use practical examples from this book.
End of exercise

If you have such choice, there are many factors that can influence your decision.

Programming language: C versus Fortran
If your loop describes the $(i,j)$ indices of a two-dimensional array, it is often best to let the $i$-index be in the inner loop for Fortran, and the $j$-index inner for C.

Exercise Can you come up with at least two reasons why this is possibly better for performance?
End of exercise

However, this is not a hard-and-fast rule. It can depend on the size of the loops, and other factors. For instance, in the matrix-vector product, changing the loop ordering changes how the input and output vectors are used.

Parallelism model
If you want to parallelize your loops with OpenMP  , you generally want the outer loop to be larger than the inner. Having a very short outer loop is definitely bad. A short inner loop can also compiler}. On the other hand, if you are targeting a GPU  , you want the large loop to be the inner one. The unit of parallel work should not have branches or loops.

Other effects of loop ordering in OpenMP are discussed in \PCSEref{sec:omp-row-col-major}.

1.7.9 Loop tiling

crumb trail: > sequential > Programming strategies for high performance > 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. This is known as \indextermbusdef{loop}{tiling}: the (short) inner loop is a tile, many consecutive instances of which form the iteration space.

For instance

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

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 (n=0; n<10; n++)
  for (i=0; i<100000; i++)
    ... = ...x[i] ...

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] ...

Exercise 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

For a less trivial example, let's look at matrix transposition $A\leftarrow B^t$. Ordinarily you would traverse the input and output matrices:

// regular.c
for (int i=0; i<N; i++)
  for (int j=0; j<N; j++)
    A[i][j] = B[j][i];

// blocked.c
for (int ii=0; ii<N; ii+=blocksize)
  for (int jj=0; jj<N; jj+=blocksize)
    for (int i=ii*blocksize; i<MIN(N,(ii+1)*blocksize); i++)
      for (int j=jj*blocksize; j<MIN(N,(jj+1)*blocksize); j++)
        A[i][j] = B[j][i];

FIGURE 1.18: Regular and blocked traversal of a matrix.

Figure  1.18 shows how one of the matrices is traversed in a different order from its storage order, for instance columnwise while it is stored by rows. This has the effect that each element load transfers a cacheline, of which only one element is immediately used. In the regular traversal, this streams of cachelines quickly overflows the cache, and there is no reuse. In the blocked traversal, however, only a small number of cachelines is traversed before the next element of these lines is needed. Thus there is reuse of cachelines, or spatial locality  .

The most important example of attaining performance through blocking is the matrix!matrix product!tiling  . In section  1.6.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]

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

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.7.10 Gradual underflow

crumb trail: > sequential > Programming strategies for high performance > Gradual underflow

The way a processor handles underflow can have an effect on the performance of a code. See section  3.8.1  .

1.7.11 Optimization strategies

crumb trail: > sequential > Programming strategies for high performance > Optimization strategies

QUOTE 1.19: Performance of naive and optimized implementations of the Discrete Fourier Transform.

QUOTE 1.20: Performance of naive and optimized implementations of the matrix-matrix product.

Figures 1.19 and 1.20 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 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

  (Note: {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.}  )

  .

• 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  [atlas-parcomp] for Blas kernels, and Spiral  [spiral] for transforms.

1.7.12 Cache aware and cache oblivious programming

crumb trail: > sequential > Programming strategies for high performance > 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 is less 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

(Note: {We are conveniently ignoring matters of set-associativity here, and basically assuming a fully associative cache.}  )

 , 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]);

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 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 behavior 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 programming}  [Frigo:oblivious]  .

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 \indextermbus{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.

\caption{Matrix transpose operation, with simple and recursive traversal of the source matrix.}

In a cache oblivious implementation we divide $A$ and $B$ as $2\times2$ block matrices, and recursively compute $B_{11}\leftarrow A_{11}^t$, $B_{12}\leftarrow A_{21}^t$, et cetera; see figure  1.7.12  . 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  [GotoGeijn:2008:Anatomy]  .

See section  6.8.4 for a discussion of such techniques in stencil computations.

1.7.13 Case study: Matrix-vector product

crumb trail: > sequential > Programming strategies for high performance > Case study: Matrix-vector product

Let us consider in some detail the matrix-vector product \begin{equation} \forall_{i,j}\colon y_i\leftarrow a_{ij}\cdot x_j \end{equation} 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

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

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

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

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

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

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

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;
}

1.8 Further topics

crumb trail: > sequential > Further topics

1.8.1 Power consumption

crumb trail: > sequential > Further topics > Power consumption

\SetBaseLevel 2

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.8.2 Derivation of scaling properties

crumb trail: > sequential > Further topics > Derivation of scaling properties

The properties of constant field scaling or Dennard scaling   [Bohr:30yearDennard,Dennard:scaling] 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$: \begin{equation} \begin{array}{lc}\toprule \hbox{Feature size}&\sim s\\ \hbox{Voltage}&\sim s\\ \hbox{Current}&\sim s \\ \hbox{Frequency}&\sim s\\ \bottomrule \end{array} \end{equation}

Then we can derive that \begin{equation} \hbox{Power} = V\cdot I \sim s^2, \end{equation} 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.

\toprule Year	\#transistors	Processor
\midrule 1975	$3, 000$	6502
1979	$30, 000$	8088
1985	$300, 000$	386
1989	$1, 000, 000$	486
1995	$6, 000, 000$	Pentium Pro
2000	$40, 000, 000$	Pentium 4
2005	$100, 000, 000$	2-core Pentium D
2008	$700, 000, 000$	8-core Nehalem
2014	$6, 000, 000, 000$	18-core Haswell
2017	$20, 000, 000, 000$	32-core AMD Epyc
2019	$40, 000, 000, 000$	64-core AMD Rome
\bottomrule

Chronology of Moore's Law (courtesy Jukka Suomela, `Programming Parallel Computers')

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{equation} \begin{array}{ll} \toprule \hbox{Charge}&q=CV\\ \hbox{Work}&W=qV=CV^2\\ \hbox{Power}&W/\hbox{time}=WF=CV^2F \\ \bottomrule \end{array} \label{eq:power} \end{equation} This analysis can be used to justify the introduction of multicore processors.

1.8.3 Multicore

crumb trail: > sequential > Further topics > 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.

\caption{Projected heat dissipation of a CPU if trends had continued -- this graph courtesy Pat Helsinger.}

Figure  1.8.3 gives a dramatic illustration of the heat that a chip would give off, if single-processor trends had continued.

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 FPGAs  . One solution to the power wall problem is introduction of multicore processors. Recall equation  1.8.2  , 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, \begin{equation} \left. \begin{array}{c} C_{\mathrm{multi}} = 2C\\ F_{\mathrm{multi}} = F/2\\ V_{\mathrm{multi}} = V/2\\ \end{array} \right\} \Rightarrow P_{\mathrm{multi}} = P/4. \end{equation} 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_{\mathrm{multi}} \approx 0.4\times P$  [Chandrakasa:transformations]  . Of course the integration aspects are a little more complicated in practice  [Bohr:ISSCC2009] ; 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.8.4 Total computer power

crumb trail: > sequential > Further topics > 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 SIMD design of processors (see section  2.3.1  , in particular section  2.3.1.2  ) is dictated by the power cost of instruction decoding.

\SetBaseLevel 1

1.8.5 Operating system effects

crumb trail: > sequential > Further topics > Operating system effects

HPC practitioners typically don't worry much about the \indexacdef{OS}. However, sometimes the presence of the OS can be felt, influencing performance. The reason for this is the periodic interrupt  , where the operating system upwards of 100 times per second interrupts the current process to let another process or a system daemon have a time slice  .

If you are running basically one program, you don't want the overhead and jitter  , the unpredictability of process runtimes, this introduces. Therefore, computers have existed that basically dispensed with having an OS to increase performance.

The periodic interrupt has further negative effects. For instance, it pollutes the cache and TLB  . As a fine-grained effect of jitter, it degrades performance of codes that rely on barriers between threads, such as frequently happens in OpenMP (section  2.6.2  ).

In particular in financial applications  , where very tight synchronization is important, have adopted a Linux kernel mode where the periodic timer ticks only once a second, rather than hundreds of times. This is called a tickless kernel  .

1.9 Review questions

crumb trail: > sequential > Review questions

For the true/false questions, give short explanation if you choose the `false' answer.

Exercise True or false. The code

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

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

Exercise 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.
End of exercise