Parallelising random reads seems to work well — why?

Question

Consider the following very simple computer program:

for i = 1 to n:
    y[i] = x[p[i]]

Here $x$ and $y$ are $n$-element arrays of bytes, and $p$ is an $n$-element array of words. Here $n$ is large, e.g., $n = 2^{31}$ (so that only a negligible fraction of the data fits in any kind of cache memory).

Assume that $p$ consists of random numbers, uniformly distributed between $1$ and $n$.

From the perspective of modern hardware, this should mean the following:

• reading $p[i]$ is cheap (sequential read)
• reading $x[p[i]]$ is very expensive (random reads; almost all reads are cache misses; we will have to fetch each individual byte from the main memory)
• writing $y[i]$ is cheap (sequential write).

And this is indeed what I am observing. The program is very slow in comparison with a program that does only sequential reads and writes. Great.

Now comes the question: how well does this program parallelise on modern multi-core platforms?

---

My hypothesis was that this program does not parallelise well. After all, the bottleneck is the main memory. A single core is already wasting most of its time just waiting for some data from the main memory.

However, this was not what I observed when I started experimenting with some algorithms where the bottleneck was this kind of operation!

I simply replaced the naive for-loop with an OpenMP parallel for-loop (in essence, it will just split the range $[1,n]$ to smaller parts and run these parts on different CPU cores in parallel).

On low-end computers, speedups were indeed minor. But on higher-end platforms I was surprised that I was getting excellent near-linear speedups. Some concrete examples (the exact timings may be a bit off, there is a lot of random variation; these were just quick experiments):

• 2 x 4-core Xeon (in total 8 cores): factor 5-8 speedups in comparison to single-threaded version.

• 2 x 6-core Xeon (in total 12 cores): factor 8-14 speedups in comparison to single-threaded version.

Now this was totally unexpected. Questions:

1. Precisely why does this kind of program parallelise so well? What happens in the hardware? (My current guess is something along these lines: the random reads from different thread are "pipelined" and the average rate of getting answers to these is much higher than in the case of a single thread.)

2. Is it necessary to use multiple threads and multiple cores to obtain any speedups? If some kind of pipelining indeed takes place in the interface between the main memory and the CPU, couldn't a single-threaded application let the main memory know that it will soon need $x[p[i]]$, $x[p[i+1]]$, ... and the computer could start fetching the relevant cache lines from the main memory? If this is possible in principle, how do I achieve it in practice?

3. What is the right theoretical model that we could use to analyse this kind of programs (and make correct predictions of the performance)?

---

Edit: There is now some source code and benchmark results available here: https://github.com/suomela/parallel-random-read

Some examples of ballpark figures ($n = 2^{32}$):

• approx. 42 ns per iteration (random read) with a single thread
• approx. 5 ns per iteration (random read) with 12 cores.

Answer 1

Forget for a moment all of the issues related to the access to main memory and level 3 cache. From a parallel perspective, ignoring these issues, the program parallelize perfectly when using $p$ processors (or cores), owing to the fact that, once you partition the work to be done through domain decomposition, each core must process either $\left\lfloor {\frac{n}{p}} \right\rfloor$ or $\left\lceil {\frac{n}{p}} \right\rceil$ elements and there is no communication and/or synchronization overhead since there is no functional dependence among the processors. Therefore, ignoring memory issues you expect a speedup equal to $p$.

Now, let's take into account the memory issues. The super-linear speedup you actually observed on your high-end Xeon based node is justified as follows.

A parallel system might exhibit such behaviour if its memory is hierarchical and if access time increases (in discrete steps) with the memory used by the program. In this case the effective computation speed could be slower on a serial processor than on a parallel computer using similar processors. This is because a sequential algorithm using $n$ bytes of memory will use only $n/p$ bytes on each processor of a $p$ processor parallel system, while cache and virtual memory effects could instead reduce the serial processor’s effective computation rate.

For $n = 2^{31}$ bytes, we need 2048 Mbytes of memory. But, when using 12 cores as in your last example, every core needs to deal with only 2048/12 Mbytes of data, which is about 170 Mbytes. High-end Xeon processors are equipped with a cache level 3 whose size ranges from 15 to 30 Mbytes of size. Clearly, with this huge cache size the cache hit ratio is high, and this explains the good or even super-linear speedup-observed.

Regarding your second question, current architectures already prefetch data by evicting cache lines and replacing them as required to exploit temporal and spacial locality of data. But this will not be enough for a single core processing 2048 Mbytes of data. If you restrict $n$ to about 170 Mbytes, then you should see on a single core more or less the same performances, since you are now running under (more or less, not exactly) the same conditions.

Finally, besides QSM (Queueing Shared Memory), I am not aware of any other theoretical parallel model taking into account at the same level the contention for access to shared memory (in your case, when using OpenMP the main memory is shared among the cores, and the cache is always shared as well among the cores). Anyway, even though the model is interesting, it did not obtain great success.

Answer 2

I decided to try out __builtin_prefetch() myself. I'm posting it here as answer in case others want to test it on their machines. The results are close to what Jukka describes: About a 20% decrease in running time when prefetching 20 elements ahead versus prefetching 0 elements ahead.

Results:

prefetch =   0, time = 1.58000
prefetch =   1, time = 1.47000
prefetch =   2, time = 1.39000
prefetch =   3, time = 1.34000
prefetch =   4, time = 1.31000
prefetch =   5, time = 1.30000
prefetch =   6, time = 1.27000
prefetch =   7, time = 1.28000
prefetch =   8, time = 1.26000
prefetch =   9, time = 1.27000
prefetch =  10, time = 1.27000
prefetch =  11, time = 1.27000
prefetch =  12, time = 1.30000
prefetch =  13, time = 1.29000
prefetch =  14, time = 1.30000
prefetch =  15, time = 1.28000
prefetch =  16, time = 1.24000
prefetch =  17, time = 1.28000
prefetch =  18, time = 1.29000
prefetch =  19, time = 1.25000
prefetch =  20, time = 1.24000
prefetch =  19, time = 1.26000
prefetch =  18, time = 1.27000
prefetch =  17, time = 1.26000
prefetch =  16, time = 1.27000
prefetch =  15, time = 1.28000
prefetch =  14, time = 1.29000
prefetch =  13, time = 1.26000
prefetch =  12, time = 1.28000
prefetch =  11, time = 1.30000
prefetch =  10, time = 1.31000
prefetch =   9, time = 1.27000
prefetch =   8, time = 1.32000
prefetch =   7, time = 1.31000
prefetch =   6, time = 1.30000
prefetch =   5, time = 1.27000
prefetch =   4, time = 1.33000
prefetch =   3, time = 1.38000
prefetch =   2, time = 1.41000
prefetch =   1, time = 1.41000
prefetch =   0, time = 1.59000

Code:

#include <stdlib.h>
#include <time.h>
#include <stdio.h>

void cracker(int *y, int *x, int *p, int n, int pf) {
    int i;
    int saved = pf;  /* let compiler optimize address computations */

    for (i = 0; i < n; i++) {
        __builtin_prefetch(&x[p[i+saved]]);
        y[i] += x[p[i]];
    }
}

int main(void) {
    int n = 50000000;
    int *x, *y, *p, i, pf, k;
    clock_t start, stop;
    double elapsed;

    /* set up arrays */
    x = malloc(sizeof(int)*n);
    y = malloc(sizeof(int)*n);
    p = malloc(sizeof(int)*n);
    for (i = 0; i < n; i++)
        p[i] = rand()%n;

    /* warm-up exercise */
    cracker(y, x, p, n, pf);

    k = 20;
    for (pf = 0; pf < k; pf++) {
        start = clock();
        cracker(y, x, p, n, pf);
        stop = clock();
        elapsed = ((double)(stop-start))/CLOCKS_PER_SEC;
        printf("prefetch = %3d, time = %.5lf\n", pf, elapsed);
    }
    for (pf = k; pf >= 0; pf--) {
        start = clock();
        cracker(y, x, p, n, pf);
        stop = clock();
        elapsed = ((double)(stop-start))/CLOCKS_PER_SEC;
        printf("prefetch = %3d, time = %.5lf\n", pf, elapsed);
    }

    return 0;
}

Answer 3

1. DDR3 access is indeed pipelined. http://www.eng.utah.edu/~cs7810/pres/dram-cs7810-protocolx2.pdf slides 20 and 24 show what happens in the memory bus during pipelined read operations.

2. (partially wrong, see below) Multiple threads are not necessary if the CPU architecture supports cache prefetch. Modern x86 and ARM as well as many other architectures have an explicit prefetch instruction. Many additionally attempt to detect patterns in memory accesses and do the prefetching automatically. The software support is compiler-specific, for example GCC and Clang have __builtin_prefech() intrinsic for explicit prefetching.

Intel-style hyperthreading seems to work very well for programs that spend most of their time waiting for cache misses. In my experience, in computation intensive workload the speedup goes very little above the number of physical cores.

EDIT: I was wrong in point 2. It seems that while prefetching can optimize memory access for single core, the combined memory bandwidth of multiple cores is greater than bandwidth of single core. How much greater, depends on the CPU.

The hardware prefetcher and other optimizations together make benchmarking very tricky. It is possible to construct cases where explicit prefetching has a very visible or non-existent effect on performance, this benchmark being one of the latter.