Computer Science

m._frigo.pdf

Home >Computer Science homework help >Computer Science

Cache-Oblivious Algorithms EXTENDED ABSTRACT

Matteo Frigo Charles E. Leiserson Harald Prokop Sridhar Ramachandran MIT Laboratory for Computer Science, 545 Technology Square, Cambridge, MA 02139�� ! ��"#��#�$� %�'&(�� #��&*)+�%�,&-��"

Abstract This paper presents asymptotically optimal algo- rithms for rectangular matrix transpose, FFT, and sorting on computers with multiple levels of caching. Unlike previous optimal algorithms, these algorithms are cache oblivious: no variables dependent on hardware parameters, such as cache size and cache-line length, need to be tuned to achieve opti- mality. Nevertheless, these algorithms use an optimal amount of work and move data optimally among multiple levels of cache. For a cache with size Z and cache-line length L where Z . Ω / L2 0 the number of cache misses for an m 1 n ma- trix transpose is Θ / 1 2 mn 3 L 0 . The number of cache misses for either an n-point FFT or the sorting of n numbers is Θ / 1 24/ n 3 L 0 / 1 2 logZ n 050 . We also give an Θ / mn p 0 -work al- gorithm to multiply an m 1 n matrix by an n 1 p matrix that incurs Θ / 1 26/ mn 2 n p 2 m p 0 3 L 2 mn p 3 L 7 Z 0 cache faults.

We introduce an “ideal-cache” model to analyze our algo- rithms. We prove that an optimal cache-oblivious algorithm designed for two levels of memory is also optimal for multi- ple levels and that the assumption of optimal replacement in the ideal-cache model can be simulated efficiently by LRU re- placement. We also provide preliminary empirical results on the effectiveness of cache-oblivious algorithms in practice.

1. Introduction Resource-oblivious algorithms that nevertheless use re- sources efficiently offer advantages of simplicity and portability over resource-aware algorithms whose re- source usage must be programmed explicitly. In this paper, we study cache resources, specifically, the hier- archy of memories in modern computers. We exhibit several “cache-oblivious” algorithms that use cache as effectively as “cache-aware” algorithms.

Before discussing the notion of cache obliviousness, we first introduce the 8 Z 9 L : ideal-cache model to study the cache complexity of algorithms. This model, which is illustrated in Figure 1, consists of a computer with a two-level memory hierarchy consisting of an ideal (data) cache of Z words and an arbitrarily large main mem- ory. Because the actual size of words in a computer is typically a small, fixed size (4 bytes, 8 bytes, etc.), we

This research was supported in part by the Defense Advanced Research Projects Agency (DARPA) under Grant F30602-97-1-0270. Matteo Frigo was supported in part by a Digital Equipment Corpora- tion fellowship.

Q cache misses

organized by optimal replacement

strategy

Main Memory

Cache

Z 3 L Cache lines Lines

of length L

CPU

W work

Figure 1: The ideal-cache model

shall assume that word size is constant; the particular constant does not affect our asymptotic analyses. The cache is partitioned into cache lines, each consisting of L consecutive words which are always moved together between cache and main memory. Cache designers typ- ically use L ; 1, banking on spatial locality to amortize the overhead of moving the cache line. We shall gener- ally assume in this paper that the cache is tall:

Z < Ω 8 L2 :�9 (1) which is usually true in practice.

The processor can only reference words that reside in the cache. If the referenced word belongs to a line already in cache, a cache hit occurs, and the word is delivered to the processor. Otherwise, a cache miss oc- curs, and the line is fetched into the cache. The ideal cache is fully associative [20, Ch. 5]: cache lines can be stored anywhere in the cache. If the cache is full, a cache line must be evicted. The ideal cache uses the optimal off-line strategy of replacing the cache line whose next access is furthest in the future [7], and thus it exploits temporal locality perfectly.

Unlike various other hierarchical-memory models [1, 2, 5, 8] in which algorithms are analyzed in terms of a single measure, the ideal-cache model uses two mea- sures. An algorithm with an input of size n is measured by its work complexity W 8 n : —its conventional running time in a RAM model [4]—and its cache complexity Q 8 n; Z 9 L : —the number of cache misses it incurs as a

function of the size Z and line length L of the ideal cache. When Z and L are clear from context, we denote the cache complexity simply as Q 8 n : to ease notation.

We define an algorithm to be cache aware if it con- tains parameters (set at either compile-time or runtime) that can be tuned to optimize the cache complexity for the particular cache size and line length. Otherwise, the algorithm is cache oblivious. Historically, good perfor- mance has been obtained using cache-aware algorithms, but we shall exhibit several optimal1 cache-oblivious al- gorithms.

To illustrate the notion of cache awareness, consider the problem of multiplying two n � n matrices A and B to produce their n � n product C. We assume that the three matrices are stored in row-major order, as shown in Figure 2(a). We further assume that n is “big,” i.e., n ; L, in order to simplify the analysis. The conventional way to multiply matrices on a computer with caches is to use a blocked algorithm [19, p. 45]. The idea is to view each matrix M as consisting of 8 n � s : � 8 n � s : submatrices Mi j (the blocks), each of which has size s � s, where s is a tuning parame- ter. The following algorithm implements this strategy:

BLOCK-MULT / A � B � C � n 0 1 for i � 1 to n 3 s 2 do for j � 1 to n 3 s 3 do for k � 1 to n 3 s 4 do ORD-MULT / Aik � Bk j � Ci j � s 0

The ORD-MULT 8 A 9 B 9 C 9 s : subroutine computes C � C � AB on s � s matrices using the ordinary O 8 s3 : algorithm. (This algorithm assumes for simplicity that s evenly divides n, but in practice s and n need have no special relationship, yielding more complicated code in the same spirit.)

Depending on the cache size of the machine on which BLOCK-MULT is run, the parameter s can be tuned to make the algorithm run fast, and thus BLOCK-MULT is a cache-aware algorithm. To minimize the cache com- plexity, we choose s to be the largest value such that the three s � s submatrices simultaneously fit in cache. An s � s submatrix is stored on Θ 8 s � s2 � L : cache lines. From the tall-cache assumption (1), we can see that s < Θ 8�� Z : . Thus, each of the calls to ORD-MULT runs with at most Z

� L < Θ 8 s2 � L : cache misses needed to

bring the three matrices into the cache. Consequently, the cache complexity of the entire algorithm is Θ 8 1 � n2 � L � 8 n � � Z : 3 8 Z � L : :�< Θ 8 1 � n2 � L � n3 � L � Z : , since

the algorithm has to read n2 elements, which reside on� n2 � L cache lines.

The same bound can be achieved using a simple

1 For simplicity in this paper, we use the term “optimal” as a syn- onym for “asymptotically optimal,” since all our analyses are asymp- totic.

cache-oblivious algorithm that requires no tuning pa- rameters such as the s in BLOCK-MULT. We present such an algorithm, which works on general rectangular matrices, in Section 2. The problems of computing a matrix transpose and of performing an FFT also suc- cumb to remarkably simple algorithms, which are de- scribed in Section 3. Cache-oblivious sorting poses a more formidable challenge. In Sections 4 and 5, we present two sorting algorithms, one based on mergesort and the other on distribution sort, both of which are op- timal in both work and cache misses.

The ideal-cache model makes the perhaps- questionable assumptions that there are only two levels in the memory hierarchy, that memory is man- aged automatically by an optimal cache-replacement strategy, and that the cache is fully associative. We address these assumptions in Section 6, showing that to a certain extent, these assumptions entail no loss of generality. Section 7 discusses related work, and Section 8 offers some concluding remarks, including some preliminary empirical results.

2. Matrix multiplication This section describes and analyzes a cache-oblivious al- gorithm for multiplying an m � n matrix by an n � p ma- trix cache-obliviously using Θ 8 mnp : work and incurring Θ 8 m � n � p � 8 mn � np � mp : � L � mnp � L � Z : cache misses. These results require the tall-cache assumption (1) for matrices stored in row-major layout format, but the assumption can be relaxed for certain other layouts. We also show that Strassen’s algorithm [31] for multi- plying n � n matrices, which uses Θ 8 nlg 7 : work2, incurs Θ 8 1 � n2 � L � nlg 7 � L � Z : cache misses.

In [9] with others, two of the present authors analyzed an optimal divide-and-conquer algorithm for n � n ma- trix multiplication that contained no tuning parameters, but we did not study cache-obliviousness per se. That algorithm can be extended to multiply rectangular matri- ces. To multiply a m � n matrix A and a n � p matrix B, the REC-MULT algorithm halves the largest of the three dimensions and recurs according to one of the following three cases:

A1 A2 � B <

A1B A2B � 9 (2)�

A1 A2 B1B2 � < A1B1 � A2B2 9 (3) A � B1 B2 < � AB1 AB2 �� (4)

In case (2), we have m � max � n 9 p � . Matrix A is split horizontally, and both halves are multiplied by matrix B. In case (3), we have n � max � m 9 p � . Both matrices are

2 We use the notation lg to denote log2.

split, and the two halves are multiplied. In case (4), we have p � max � m 9 n � . Matrix B is split vertically, and each half is multiplied by A. For square matrices, these three cases together are equivalent to the recursive mul- tiplication algorithm described in [9]. The base case oc- curs when m < n < p < 1, in which case the two ele- ments are multiplied and added into the result matrix.

Although this straightforward divide-and-conquer al- gorithm contains no tuning parameters, it uses cache op- timally. To analyze the REC-MULT algorithm, we as- sume that the three matrices are stored in row-major or- der, as shown in Figure 2(a). Intuitively, REC-MULT uses the cache effectively, because once a subproblem fits into the cache, its smaller subproblems can be solved in cache with no further cache misses.

Theorem 1 The REC-MULT algorithm uses Θ 8 mnp : work and incurs Θ 8 m � n � p � 8 mn � np � mp : � L � mnp

� L � Z : cache misses when multiplying an m � n ma-

trix by an n � p matrix. Proof. It can be shown by induction that the work of REC-MULT is Θ 8 mnp : . To analyze the cache misses, let α ; 0 be the largest constant sufficiently small that three submatrices of sizes m

� � n � , n � � p � , and m � � p � , where max � m � 9 n � 9 p � � � α � Z, all fit completely in the cache. We distinguish four cases depending on the initial size of the matrices. Case I: m 9 n 9 p ; α � Z. This case is the most intuitive. The matrices do not fit in cache, since all dimensions are “big enough.” The cache complexity can be described by the recurrence

Q / m � n � p 0�� (5)�� Θ / / mn 2 n p 2 m p 0 3 L 0 if m � n � p � α 7 Z 3 2 � α 7 Z � �2Q / m 3 2 � n � p 0 2 O / 1 0 ow. if m � n and m � p � 2Q / m � n 3 2 � p 0 2 O / 1 0 ow. if n m and n � p � 2Q / m � n � p 3 2 0 2 O / 1 0 otherwise �

The base case arises as soon as all three submatrices fit in cache. The total number of lines used by the three submatrices is Θ 8 8 mn � np � mp : � L : . The only cache misses that occur during the remainder of the recursion are the Θ 8*8 mn � np � mp : � L : cache misses required to bring the matrices into cache. In the recursive cases, when the matrices do not fit in cache, we pay for the cache misses of the recursive calls, which depend on the dimensions of the matrices, plus O 8 1 : cache misses for the overhead of manipulating submatrices. The solution to this recurrence is Q 8 m 9 n 9 p : < Θ 8 mnp � L � Z : . Case II: (m

� α � Z and n 9 p ; α � Z) or (n � α � Z and

m 9 p ; α � Z) or (p � α � Z and m 9 n ; α � Z). Here, we shall present the case where m

� α � Z and n 9 p ; α � Z.

The proofs for the other cases are only small variations of this proof. The REC-MULT algorithm always divides n or p by 2 according to cases (3) and (4). At some point

00 11 22 33 44 55 66 77

88 99

1010 1111 1212 1313 1414 1515

1616 1717 1818 1919 2020 2121 2222 2323

2424 2525 2626 2727 2828 2929 3030 3131

3232 3333 3434 3535 3636 3737 3838 3939

4040 4141 4242 4343 4444 4545 4646 4747

4848 4949 5050 5151 5252 5353 5454 5555

5656 5757 5858 5959 6060 6161 6262 6363

00 11 22 33 44 55 66 77 88 99 1010 1111 1212 1313 1414 1515 1616 1717 1818 1919 2020 2121 2222 2323 2424 2525 2626 2727 2828 2929 3030 3131 3232 3333 3434 3535 3636 3737 3838 3939 4040 4141 4242 4343 4444 4545 4646 4747 4848 4949 5050 5151 5252 5353 5454 5555 5656 5757 5858 5959 6060 6161 6262 6363

(c)

(a)

(d)

(b)

00 11 22 33 44 55 66 77 88 99 1010 1111 1212 1313 1414 1515

1616 1717 1818 1919 2020 2121 2222 2323 2424 2525 2626 2727 2828 2929 3030 3131

3232 3333 3434 3535 3636 3737 3838 3939 4040 4141 4242 4343 4444 4545 4646 4747

4848 4949 5050 5151 5252 5353 5454 5555 5656 5757 5858 5959 6060 6161 6262 6363

00 11 22 33

44 55 66 77

88 99 1010 1111

1212 1313 1414 1515

1616 1717 1818 1919

2020 2121 2222 2323

2424 2525 2626 2727

2828 2929 3030 3131

3232 3333 3434 3535

3636 3737 3838 3939

4040 4141 4242 4343

4444 4545 4646 4747

4848 4949 5050 5151

5252 5353 5454 5555

5656 5757 5858 5959

6060 6161 6262 6363

Figure 2: Layout of a 16 � 16 matrix in (a) row ma- jor, (b) column major, (c) 4 � 4-blocked, and (d) bit- interleaved layouts.

in the recursion, both are small enough that the whole problem fits into cache. The number of cache misses can be described by the recurrence

Q / m � n � p 0�� (6)�� Θ / 1 2 n 2 n p 3 L 2 m 0 if n � p �� α 7 Z 3 2 � α 7 Z � � 2Q / m � n 3 2 � p 0 2 O / 1 0 otherwise if n � p � 2Q / m � n � p 3 2 0 2 O / 1 0 otherwise ;

whose solution is Q 8 m 9 n 9 p : < Θ 8 np � L � mnp � L � Z : . Case III: (n 9 p � α � Z and m ; α � Z) or (m 9 p � α � Z and n ; α � Z) or (m 9 n � α � Z and p ; α � Z). In each of these cases, one of the matrices fits into cache, and the others do not. Here, we shall present the case where n 9 p � α � Z and m ; α � Z. The other cases can be proven similarly. The REC-MULT algorithm always di- vides m by 2 according to case (2). At some point in the recursion, m falls into the range α � Z � 2 � m � α � Z, and the whole problem fits in cache. The number cache misses can be described by the recurrence

Q / m � n 0�� (7)� Θ / 1 2 m 0 if m �� α 7 Z 3 2 � α 7 Z � � 2Q / m 3 2 � n � p 0 2 O / 1 0 otherwise ;

whose solution is Q 8 m 9 n 9 p : < Θ 8 m � mnp � L � Z : . Case IV: m 9 n 9 p � α � Z. From the choice of α, all three matrices fit into cache. The matrices are stored on Θ 8 1 � mn � L � np � L � mp � L : cache lines. Therefore, we have Q 8 m 9 n 9 p : < Θ 8 1 � 8 mn � np � mp : � L : .

We require the tall-cache assumption (1) in these analyses, because the matrices are stored in row-major order. Tall caches are also needed if matrices are stored

in column-major order (Figure 2(b)), but the assumption that Z < Ω 8 L2 : can be relaxed for certain other matrix layouts. The s � s-blocked layout (Figure 2(c)), for some tuning parameter s, can be used to achieve the same bounds with the weaker assumption that the cache holds at least some sufficiently large constant number of lines. The cache-oblivious bit-interleaved layout (Figure 2(d)) has the same advantage as the blocked layout, but no tuning parameter need be set, since submatrices of size O 8 � L : � O 8 � L : are cache-obliviously stored on O 8 1 : cache lines. The advantages of bit-interleaved and re- lated layouts have been studied in [11, 12, 16]. One of the practical disadvantages of bit-interleaved layouts is that index calculations on conventional microprocessors can be costly, a deficiency we hope that processor archi- tects will remedy.

For square matrices, the cache complexity Q 8 n : < Θ 8 n � n2 � L � n3 � L � Z : of the REC-MULT algorithm is the same as the cache complexity of the cache-aware BLOCK-MULT algorithm and also matches the lower bound by Hong and Kung [21]. This lower bound holds for all algorithms that execute the Θ 8 n3 : opera- tions given by the definition of matrix multiplication

ci j < n

∑ k � 1 aikbk j �

No tight lower bounds for the general problem of matrix multiplication are known.

By using an asymptotically faster algorithm, such as Strassen’s algorithm [31] or one of its variants [37], both the work and cache complexity can be reduced. When multiplying n � n matrices, Strassen’s algorithm, which is cache oblivious, requires only 7 recursive multiplica- tions of n

� 2 � n � 2 matrices and a constant number of

matrix additions, yielding the recurrence

Q 8 n : � � Θ 8 1 � n � n2 � L : if n2 � αZ 9 7Q 8 n � 2 : � O 8 n2 � L : otherwise ; (8)

where α is a sufficiently small constant. The solution to this recurrence is Θ 8 n � n2 � L � nlg 7 � L � Z : . 3. Matrix transposition and FFT This section describes a recursive cache-oblivious al- gorithm for transposing an m � n matrix which uses O 8 mn : work and incurs O 8 1 � mn � L : cache misses, which is optimal. Using matrix transposition as a sub- routine, we convert a variant [36] of the “six-step” fast Fourier transform (FFT) algorithm [6] into an optimal cache-oblivious algorithm. This FFT algorithm uses O 8 n lg n : work and incurs O 8 1 � 8 n � L :�8 1 � logZ n : : cache misses.

The problem of matrix transposition is defined as fol- lows. Given an m � n matrix stored in a row-major lay- out, compute and store AT into an n � m matrix B also

stored in a row-major layout. The straightforward algo- rithm for transposition that employs doubly nested loops incurs Θ 8 mn : cache misses on one of the matrices when m � Z � L and n � Z � L, which is suboptimal.

Optimal work and cache complexities can be ob- tained with a divide-and-conquer strategy, however. If n � m, the REC-TRANSPOSE algorithm partitions

A < 8 A1 A2 : 9 B <

B1 B2 �

and recursively executes REC-TRANSPOSE 8 A1 9 B1 : and REC-TRANSPOSE 8 A2 9 B2 : . Otherwise, it divides matrix A horizontally and matrix B vertically and likewise per- forms two transpositions recursively. The next two lem- mas provide upper and lower bounds on the performance of this algorithm.

Lemma 2 The REC-TRANSPOSE algorithm involves O 8 mn : work and incurs O 8 1 � mn � L : cache misses for an m � n matrix. Proof. That the algorithm does O 8 mn : work is straight- forward. For the cache analysis, let Q 8 m 9 n : be the cache complexity of transposing an m � n matrix. We as- sume that the matrices are stored in row-major order, the column-major layout having a similar analysis.

Let α be a constant sufficiently small such that two submatrices of size m � n and n � m, where max � m 9 n � � αL, fit completely in the cache even if each row is stored in a different cache line. We distinguish the three cases.

Case I: max � m 9 n � � αL. Both the matrices fit in O 8 1 : � 2mn � L lines. From the choice of α, the number of lines required is at most Z

� L. Therefore Q 8 m 9 n :�<

O 8 1 � mn � L : . Case II: m

� αL � n or n � αL � m. Suppose first that

m �

αL � n. The REC-TRANSPOSE algorithm divides the greater dimension n by 2 and performs divide and conquer. At some point in the recursion, n falls into the range αL

� 2 �

n �

αL, and the whole problem fits in cache. Because the layout is row-major, at this point the input array has n rows and m columns, and it is laid out in contiguous locations, requiring at most O 8 1 � nm � L : cache misses to be read. The output array consists of nm elements in m rows, where in the worst case every row lies on a different cache line. Consequently, we incur at most O 8 m � nm � L : for writing the output array. Since n � αL � 2, the total cache complexity for this base case is O 8 1 � m : . These observations yield the recurrence Q 8 m 9 n : �

� O 8 1 � m : if n �� αL � 2 9 αL ��9 2Q 8 m 9 n � 2 : � O 8 1 : otherwise ;

whose solution is Q 8 m 9 n : < O 8 1 � mn � L : . The case n

� αL � m is analogous.

Case III: m 9 n ; αL. As in Case II, at some point in the recursion both n and m fall into the range � αL � 2 9 αL � . The whole problem fits into cache and can be solved with at most O 8 m � n � mn � L : cache misses. The cache complexity thus satisfies the recurrence

Q 8 m 9 n : �� O 8 m � n � mn � L : if m 9 n � � αL � 2 9 αL ��9 2Q 8 m � 2 9 n : � O 8 1 : if m � n 9 2Q 8 m 9 n � 2 : � O 8 1 : otherwise;

whose solution is Q 8 m 9 n : < O 8 1 � mn � L : . Theorem 3 The REC-TRANSPOSE algorithm has opti- mal cache complexity.

Proof. For an m � n matrix, the algorithm must write to mn distinct elements, which occupy at least � mn � L � < Ω 8 1 � mn � L : cache lines.

As an example of an application of this cache- oblivious transposition algorithm, in the rest of this sec- tion we describe and analyze a cache-oblivious algo- rithm for computing the discrete Fourier transform of a complex array of n elements, where n is an exact power of 2. The basic algorithm is the well-known “six-step” variant [6, 36] of the Cooley-Tukey FFT algorithm [13]. Using the cache-oblivious transposition algorithm, how- ever, the FFT becomes cache-oblivious, and its perfor- mance matches the lower bound by Hong and Kung [21].

Recall that the discrete Fourier transform (DFT) of an array X of n complex numbers is the array Y given by

Y � i ��< n � 1∑ j � 0 X � j � ω

� i j n 9 (9)

where ωn < e2π � � 1 � n is a primitive nth root of unity, and 0

� i � n. Many algorithms evaluate Equation (9) in

O 8 n lg n : time for all integers n [15]. In this paper, how- ever, we assume that n is an exact power of 2, and we compute Equation (9) according to the Cooley-Tukey al- gorithm, which works recursively as follows. In the base case where n < O 8 1 : , we compute Equation (9) directly. Otherwise, for any factorization n < n1n2 of n, we have Y � i1 � i2n1 � < (10)

n2 � 1 ∑

j2 � 0

� n1 � 1 ∑ j1 � 0 X � j1n2 � j2 � ω � i1 j1n1 � ω � i1 j2n � ω � i2 j2n2 �

Observe that both the inner and outer summations in Equation (10) are DFT’s. Operationally, the computa- tion specified by Equation (10) can be performed by computing n2 transforms of size n1 (the inner sum), mul- tiplying the result by the factors ω � i1 j2n (called the twid- dle factors [15]), and finally computing n1 transforms of size n2 (the outer sum).

We choose n1 to be 2 lg n � 2 � and n2 to be 2 � lg n � 2 � . The recursive step then operates as follows:

1. Pretend that input is a row-major n1 � n2 matrix A.

Transpose A in place, i.e., use the cache-oblivious REC-TRANSPOSE algorithm to transpose A onto an auxiliary array B, and copy B back onto A. Notice that if n1 < 2n2, we can consider the matrix to be made up of records containing two elements.

2. At this stage, the inner sum corresponds to a DFT of the n2 rows of the transposed matrix. Compute these n2 DFT’s of size n1 recursively. Observe that, because of the previous transposition, we are trans- forming a contiguous array of elements.

3. Multiply A by the twiddle factors, which can be computed on the fly with no extra cache misses.

4. Transpose A in place, so that the inputs to the next stage are arranged in contiguous locations.

5. Compute n1 DFT’s of the rows of the matrix recur- sively.

6. Transpose A in place so as to produce the correct output order.

It can be proven by induction that the work com- plexity of this FFT algorithm is O 8 n lg n : . We now an- alyze its cache complexity. The algorithm always op- erates on contiguous data, by construction. Thus, by the tall-cache assumption (1), the transposition oper- ations and the twiddle-factor multiplication require at most O 8 1 � n � L : cache misses. Thus, the cache com- plexity satisfies the recurrence

Q 8 n : � �� O 8 1 � n � L : 9 if n � αZ 9

n1Q 8 n2 : � n2Q 8 n1 : otherwise ;� O 8 1 � n � L : (11) where α ; 0 is a constant sufficiently small that a sub- problem of size αZ fits in cache. This recurrence has solution

Q 8 n :#< O 8 1 � 8 n � L :%8 1 � logZ n :*: 9 which is optimal for a Cooley-Tukey algorithm, match- ing the lower bound by Hong and Kung [21] when n is an exact power of 2. As with matrix multiplication, no tight lower bounds for cache complexity are known for the general DFT problem.

4. Funnelsort Cache-oblivious algorithms, like the familiar two-way merge sort, are not optimal with respect to cache misses. The Z-way mergesort suggested by Aggarwal and Vit- ter [3] has optimal cache complexity, but although it ap- parently works well in practice [23], it is cache aware. This section describes a cache-oblivious sorting algo- rithm called “funnelsort.” This algorithm has optimal

k-merger

buffers

L � k

Figure 3: Illustration of a k-merger. A k-merger is built recursively out of � k “left” � k-mergers L1, L2, � � � , L � k , a series of buffers, and one “right” � k-merger R. O 8 n lg n : work complexity, and optimal O 8 1 � 8 n � L : 8 1 � logZ n : : cache complexity.

Funnelsort is similar to mergesort. In order to sort a (contiguous) array of n elements, funnelsort performs the following two steps:

1. Split the input into n1 � 3 contiguous arrays of size n2 � 3, and sort these arrays recursively.

2. Merge the n1 � 3 sorted sequences using a n1 � 3- merger, which is described below.

Funnelsort differs from mergesort in the way the merge operation works. Merging is performed by a de- vice called a k-merger, which inputs k sorted sequences and merges them. A k-merger operates by recursively merging sorted sequences which become progressively longer as the algorithm proceeds. Unlike mergesort, however, a k-merger suspends work on a merging sub- problem when the merged output sequence becomes “long enough” and resumes work on another merging subproblem.

This complicated flow of control makes a k-merger a bit tricky to describe. Figure 3 shows a representa- tion of a k-merger, which has k sorted sequences as in- puts. Throughout its execution, the k-merger maintains the following invariant.

Invariant Each invocation of a k-merger outputs the next k3 elements of the sorted sequence obtained by merging the k input sequences.

A k-merger is built recursively out of � k-mergers in the following way. The k inputs are partitioned into � k sets of � k elements, which form the input to the � k� k-mergers L1 9 L2 9 � � � 9 L � k in the left part of the figure. The outputs of these mergers are connected to the inputs

of � k buffers. Each buffer is a FIFO queue that can hold 2k3 � 2 elements. Finally, the outputs of the buffers are connected to the � k inputs of the � k-merger R in the right part of the figure. The output of this final � k- merger becomes the output of the whole k-merger. The intermediate buffers are overdimensioned, since each can hold 2k3 � 2 elements, which is twice the number k3 � 2 of elements output by a � k-merger. This additional buffer space is necessary for the correct behavior of the algorithm, as will be explained below. The base case of the recursion is a k-merger with k < 2, which produces k3 < 8 elements whenever invoked.

A k-merger operates recursively in the following way. In order to output k3 elements, the k-merger invokes R k3 � 2 times. Before each invocation, however, the k- merger fills all buffers that are less than half full, i.e., all buffers that contain less than k3 � 2 elements. In order to fill buffer i, the algorithm invokes the corresponding left merger Li once. Since Li outputs k

3 � 2 elements, the buffer contains at least k3 � 2 elements after Li finishes.

It can be proven by induction that the work com- plexity of funnelsort is O 8 n lg n : . We will now analyze the cache complexity. The goal of the analysis is to show that funnelsort on n elements requires at most Q 8 n : cache misses, where

Q 8 n :#< O 8 1 � 8 n � L : 8 1 � logZ n :*: � In order to prove this result, we need three auxiliary lem- mas. The first lemma bounds the space required by a k-merger.

Lemma 4 A k-merger can be laid out in O 8 k2 : contigu- ous memory locations.

Proof. A k-merger requires O 8 k2 : memory locations for the buffers, plus the space required by the � k- mergers. The space S 8 k : thus satisfies the recurrence

S 8 k : � 8 � k � 1 : S 8 � k : � O 8 k2 :�9 whose solution is S 8 k :�< O 8 k2 : .

In order to achieve the bound on Q 8 n : , the buffers in a k-merger must be maintained as circular queues of size k. This requirement guarantees that we can man- age the queue cache-efficiently, in the sense stated by the next lemma.

Lemma 5 Performing r insert and remove operations on a circular queue causes in O 8 1 � r � L : cache misses as long as two cache lines are available for the buffer.

Proof. Associate the two cache lines with the head and tail of the circular queue. If a new cache line is read during a insert (delete) operation, the next L � 1 insert (delete) operations do not cause a cache miss.

The next lemma bounds the cache complexity of a k-merger.

Lemma 6 If Z < Ω 8 L2 : , then a k-merger operates with at most

QM 8 k : < O 8 1 � k � k3 � L � k3 logZ k � L : cache misses.

Proof. There are two cases: either k � α � Z or k ; α � Z, where α is a sufficiently small constant. Case I: k � α � Z. By Lemma 4, the data structure associated with the k-merger requires at most O 8 k2 :�< O 8 Z : contiguous memory locations, and therefore it fits into cache. The k-merger has k input queues from which it loads O 8 k3 : elements. Let ri be the number of elements extracted from the ith input queue. Since k � α � Z and the tall-cache assumption (1) implies that L < O 8�� Z : , there are at least Z � L < Ω 8 k : cache lines available for the input buffers. Lemma 5 applies, whence the total number of cache misses for accessing the input queues is

∑ i � 1 O 8 1 � ri � L : < O 8 k � k3 � L : �

Similarly, Lemma 4 implies that the cache complexity of writing the output queue is O 8 1 � k3 � L : . Finally, the algorithm incurs O 8 1 � k2 � L : cache misses for touching its internal data structures. The total cache complexity is therefore QM 8 k : < O 8 1 � k � k3 � L : . Case I: k ; α � Z. We prove by induction on k that whenever k ; α � Z, we have

QM 8 k : � ck3 logZ k � L � A 8 k :�9 (12) where A 8 k :�< k 8 1 � 2c logZ k � L :�< o 8 k3 : . This particular value of A 8 k : will be justified at the end of the analysis.

The base case of the induction consists of values of k such that αZ1 � 4 � k � α � Z. (It is not sufficient only to consider k < Θ 8�� Z : , since k can become as small as Θ 8 Z1 � 4 : in the recursive calls.) The analysis of the first case applies, yielding QM 8 k : < O 8 1 � k � k3 � L : . Be- cause k2 ; α � Z < Ω 8 L : and k < Ω 8 1 : , the last term dominates, which implies QM 8 k :+< O 8 k3 � L : . Conse- quently, a big enough value of c can be found that satis- fies Inequality (12).

For the inductive case, suppose that k ; α � Z. The k-merger invokes the � k-mergers recursively. Since αZ1 � 4 � � k � k, the inductive hypothesis can be used to bound the number QM 8 � k : of cache misses incurred by the submergers. The “right” merger R is invoked exactly k3 � 2 times. The total number l of invocations of “left” mergers is bounded by l � k3 � 2 � 2 � k. To see why, con- sider that every invocation of a left merger puts k3 � 2 el- ements into some buffer. Since k3 elements are output and the buffer space is 2k2, the bound l � k3 � 2 � 2 � k follows.

Before invoking R, the algorithm must check every buffer to see whether it is empty. One such check re- quires at most � k cache misses, since there are � k buffers. This check is repeated exactly k3 � 2 times, lead- ing to at most k2 cache misses for all checks. These considerations lead to the recurrence

QM 8 k : � � 2k3 � 2 � 2 � k � QM 8 � k : � k2 � Application of the inductive hypothesis and the choice A 8 k :�< k 8 1 � 2c logZ k � L : yields Inequality (12) as fol- lows:

QM 8 k : � � 2k3 � 2 � 2 � k � QM 8 � k : � k2 �

2 � k3 � 2 � � k � � ck3 � 2 logZ k

2L � A 8 � k : � � k2

� ck3 logZ k

� L � k2 8 1 � c logZ k � L :

�

� 2k3 � 2 � 2 � k � A 8 � k :�

ck3 logZ k � L � A 8 k : �

Theorem 7 To sort n elements, funnelsort incurs O 8 1 � 8 n � L :-8 1 � logZ n : : cache misses. Proof. If n � αZ for a small enough constant α, then the algorithm fits into cache. To see why, observe that only one k-merger is active at any time. The biggest k-merger is the top-level n1 � 3-merger, which requires O 8 n2 � 3 : � O 8 n : space. The algorithm thus can operate in O 8 1 � n � L : cache misses.

If N ; αZ, we have the recurrence Q 8 n : < n1 � 3Q 8 n2 � 3 : � QM 8 n1 � 3 : �

By Lemma 6, we have QM 8 n1 � 3 :�< O 8 1 � n1 � 3 � n � L � n logZ n

� L : .

By the tall-cache assumption (1), we have n � L <

Ω 8 n1 � 3 : . Moreover, we also have n1 � 3 < Ω 8 1 : and lg n < Ω 8 lg Z : . Consequently, QM 8 n1 � 3 : < O 8 n logZ n � L : holds, and the recurrence simplifies to

Q 8 n :#< n1 � 3Q 8 n2 � 3 : � O 8 n logZ n � L : � The result follows by induction on n.

This upper bound matches the lower bound stated by the next theorem, proving that funnelsort is cache- optimal.

Theorem 8 The cache complexity of any sorting algo- rithm is Q 8 n :#< Ω 8 1 � 8 n � L :-8 1 � logZ n :*: . Proof. Aggarwal and Vitter [3] show that there is an Ω 8 8 n � L : logZ � L 8 n � Z :*: bound on the number of cache misses made by any sorting algorithm on their “out-of- core” memory model, a bound that extends to the ideal- cache model. The theorem can be proved by apply- ing the tall-cache assumption Z < Ω 8 L2 : and the trivial lower bounds of Q 8 n : < Ω 8 1 : and Q 8 n :#< Ω 8 n � L : .

5. Distribution sort

In this section, we describe another cache-oblivious op- timal sorting algorithm based on distribution sort. Like the funnelsort algorithm from Section 4, the distribution- sorting algorithm uses O 8 n lg n : work to sort n elements, and it incurs O 8 1 � 8 n � L :%8 1 � logZ n : : cache misses. Unlike previous cache-efficient distribution-sorting al- gorithms [1, 3, 25, 34, 36], which use sampling or other techniques to find the partitioning elements before the distribution step, our algorithm uses a “bucket splitting” technique to select pivots incrementally during the dis- tribution step.

Given an array A (stored in contiguous locations) of length n, the cache-oblivious distribution sort operates as follows:

1. Partition A into � n contiguous subarrays of size� n. Recursively sort each subarray. 2. Distribute the sorted subarrays into q buckets

B1 9 � � � 9 Bq of size n1 9 � � � 9 nq, respectively, such that 1. max � x � x � Bi � � min � x � x � Bi � 1 � for i <

1 9 2 9 � � � 9 q � 1. 2. ni

� 2 � n for i < 1 9 2 9 � � � 9 q.

(See below for details.)

3. Recursively sort each bucket.

4. Copy the sorted buckets to array A.

A stack-based memory allocator is used to exploit spatial locality.

The goal of Step 2 is to distribute the sorted subarrays of A into q buckets B1 9 B2 9 � � � 9 Bq. The algorithm main- tains two invariants. First, at any time each bucket holds at most 2 � n elements, and any element in bucket Bi is smaller than any element in bucket Bi � 1. Second, every bucket has an associated pivot. Initially, only one empty bucket exists with pivot ∞.

The idea is to copy all elements from the subarrays into the buckets while maintaining the invariants. We keep state information for each subarray and bucket. The state of a subarray consists of the index next of the next element to be read from the subarray and the bucket number bnum where this element should be copied. By convention, bnum < ∞ if all elements in a subarray have been copied. The state of a bucket consists of the pivot and the number of elements currently in the bucket.

We would like to copy the element at position next of a subarray to bucket bnum. If this element is greater than the pivot of bucket bnum, we would increment bnum un- til we find a bucket for which the element is smaller than the pivot. Unfortunately, this basic strategy has poor caching behavior, which calls for a more complicated procedure.

The distribution step is accomplished by the recur- sive procedure DISTRIBUTE 8 i 9 j 9 m : which distributes elements from the ith through 8 i � m � 1 : th subarrays into buckets starting from B j . Given the precondition that each subarray i 9 i � 1 9 � � � 9 i � m � 1 has its bnum � j, the execution of DISTRIBUTE 8 i 9 j 9 m : enforces the post- condition that subarrays i 9 i � 1 9 � � � 9 i � m � 1 have their bnum � j � m. Step 2 of the distribution sort invokes DISTRIBUTE 8 1 9 1 9 � n : . The following is a recursive im- plementation of DISTRIBUTE:

DISTRIBUTE 8 i 9 j 9 m : 1 if m < 1 2 then COPYELEMS 8 i 9 j : 3 else DISTRIBUTE 8 i 9 j 9 m � 2 : 4 DISTRIBUTE 8 i � m � 2 9 j 9 m � 2 : 5 DISTRIBUTE 8 i 9 j � m � 2 9 m � 2 : 6 DISTRIBUTE 8 i � m � 2 9 j � m � 2 9 m � 2 :

In the base case, the procedure COPYELEMS 8 i 9 j : copies all elements from subarray i that belong to bucket j. If bucket j has more than 2 � n elements af- ter the insertion, it can be split into two buckets of size at least � n. For the splitting operation, we use the deter- ministic median-finding algorithm [14, p. 189] followed by a partition.

Lemma 9 The median of n elements can be found cache-obliviously using O 8 n : work and incurring O 8 1 � n � L : cache misses.

Proof. See [14, p. 189] for the linear-time median find- ing algorithm and the work analysis. The cache com- plexity is given by the same recurrence as the work com- plexity with a different base case.

Q 8 m :#< �� O 8 1 � m � L : if m � αZ 9

Q 8 � m � 5 ��: � Q 8 7m � 10 � 6 : otherwise ;� O 8 1 � m � L : where α is a sufficiently small constant. The result fol- lows.

In our case, we have buckets of size 2 � n � 1. In ad- dition, when a bucket splits, all subarrays whose bnum is greater than the bnum of the split bucket must have their bnum’s incremented. The analysis of DISTRIBUTE is given by the following lemma.

Lemma 10 The distribution step involves O 8 n : work, incurs O 8 1 � n � L : cache misses, and uses O 8 n : stack space to distribute n elements.

Proof. In order to simplify the analysis of the work used by DISTRIBUTE, assume that COPYELEMS uses O 8 1 : work for procedural overhead. We will account for the work due to copying elements and splitting of buck- ets separately. The work of DISTRIBUTE is described by

the recurrence

T 8 c : < 4T 8 c � 2 : � O 8 1 : � It follows that T 8 c :�< O 8 c2 : , where c < � n initially. The work due to copying elements is also O 8 n : .

The total number of bucket splits is at most � n. To see why, observe that there are at most � n buckets at the end of the distribution step, since each bucket contains at least � n elements. Each split operation involves O 8 � n : work and so the net contribution to the work is O 8 n : . Thus, the total work used by DISTRIBUTE is W 8 n : < O 8 T 8 � n :*: � O 8 n : � O 8 n :�< O 8 n : .

For the cache analysis, we distinguish two cases. Let α be a sufficiently small constant such that the stack space used fits into cache. Case I, n

� αZ: The input and the auxiliary space of

size O 8 n : fit into cache using O 8 1 � n � L : cache lines. Consequently, the cache complexity is O 8 1 � n � L : . Case II, n ; αZ: Let R 8 c 9 m : denote the cache misses incurred by an invocation of DISTRIBUTE 8 a 9 b 9 c : that copies m elements from subarrays to buckets. We first prove that R 8 c 9 m : < O 8 L � c2 � L � m � L : , ignoring the cost splitting of buckets, which we shall account for sep- arately. We argue that R 8 c 9 m : satisfies the recurrence

R 8 c 9 m : � �� O 8 L � m � L : if c � αL 9

∑ i � 1 R 8 c

� 2 9 mi : otherwise ; (13)

where ∑4i � 1 mi < m, whose solution is R 8 c 9 m :�< O 8 L � c2 � L � m � L : . The recursive case c ; αL follows im-

mediately from the algorithm. The base case c �

αL can be justified as follows. An invocation of DISTRIBUTE 8 a 9 b 9 c : operates with c subarrays and c buckets. Since there are Ω 8 L : cache lines, the cache can hold all the auxiliary storage involved and the currently accessed element in each subarray and bucket. In this case there are O 8 L � m � L : cache misses. The initial ac- cess to each subarray and bucket causes O 8 c : < O 8 L : cache misses. Copying the m elements to and from con- tiguous locations causes O 8 1 � m � L : cache misses.

We still need to account for the cache misses caused by the splitting of buckets. Each split causes O 8 1 �� n � L : cache misses due to median finding (Lemma 9) and partitioning of � n contiguous elements. An addi- tional O 8 1 � � n � L : misses are incurred by restoring the cache. As proven in the work analysis, there are at most� n split operations. By adding R 8 � n 9 n : to the split complexity, we conclude that the total cache complexity of the distribution step is O 8 L � n � L � � n 8 1 � � n � L : :�< O 8 n � L : .

The analysis of distribution sort is given in the next theorem. The work and cache complexity match lower bounds specified in Theorem 8.

Theorem 11 Distribution sort uses O 8 n lg n : work and incurs O 8 1 � 8 n � L :�8 1 � logZ n : : cache misses to sort n elements.

Proof. The work done by the algorithm is given by

W 8 n :#< � nW 8 � n : � q∑ i � 1W 8 ni : � O 8 n :�9

where each ni �

2 � n and ∑ ni < n. The solution to this recurrence is W 8 n :�< O 8 n lg n : .

The space complexity of the algorithm is given by

S 8 n : � S 8 2 � n : � O 8 n :�9 where the O 8 n : term comes from Step 2. The solution to this recurrence is S 8 n :�< O 8 n : .

The cache complexity of distribution sort is described by the recurrence

Q 8 n : � �� O 8 1 � n � L : if n � αZ 9� nQ 8 � n : � ∑qi � 1 Q 8 ni : otherwise ;� O 8 1 � n � L :

where α is a sufficiently small constant such that the stack space used by a sorting problem of size αZ, in- cluding the input array, fits completely in cache. The base case n

� αZ arises when both the input array A

and the contiguous stack space of size S 8 n : < O 8 n : fit in O 8 1 � n � L : cache lines of the cache. In this case, the algorithm incurs O 8 1 � n � L : cache misses to touch all involved memory locations once. In the case where n ; αZ, the recursive calls in Steps 1 and 3 cause Q 8 � n : � ∑qi � 1 Q 8 ni : cache misses and O 8 1 � n � L : is the cache complexity of Steps 2 and 4, as shown by Lemma 10. The theorem follows by solving the recur- rence.

6. Theoretical justifications for the ideal- cache model

How reasonable is the ideal-cache model for algorithm design? The model incorporates four major assumptions that deserve scrutiny:

� optimal replacement, � exactly two levels of memory, � automatic replacement, � full associativity.

Designing algorithms in the ideal-cache model is easier than in models lacking these properties, but are these assumptions too strong? In this section we show that cache-oblivious algorithms designed in the ideal-cache model can be efficiently simulated by weaker models.

The first assumption that we shall eliminate is that of optimal replacement. Our strategy for the simula- tion is to use an LRU (least-recently used) replacement strategy [20, p. 378] in place of the optimal and om- niscient replacement strategy. We start by proving a

lemma that bounds the effectiveness of the LRU simu- lation. We then show that algorithms whose complex- ity bounds satisfy a simple regularity condition (includ- ing all algorithms heretofore presented) can be ported to caches incorporating an LRU replacement policy.

Lemma 12 Consider an algorithm that causes Q

� 8 n; Z 9 L : cache misses on a problem of size n using a 8 Z 9 L : ideal cache. Then, the same algorithm incurs Q 8 n; Z 9 L : � 2Q � 8 n; Z � 2 9 L : cache misses on a 8 Z 9 L : cache that uses LRU replacement.

Proof. Sleator and Tarjan [30] have shown that the cache misses on a 8 Z 9 L : cache using LRU replacement are 8 Z � L : � 8 8 Z � Z � : � L � 1 : -competitive with optimal replacement on a 8 Z � 9 L : ideal cache if both caches start empty. It follows that the number of misses on a 8 Z 9 L : LRU-cache is at most twice the number of misses on a 8 Z � 2 9 L : ideal-cache. Corollary 13 For any algorithm whose cache- complexity bound Q 8 n; Z 9 L : in the ideal-cache model satisfies the regularity condition

Q 8 n; Z 9 L :�< O 8 Q 8 n; 2Z 9 L : : 9 (14) the number of cache misses with LRU replacement is Θ 8 Q 8 n; Z 9 L : : . Proof. Follows directly from (14) and Lemma 12.

The second assumption we shall eliminate is the as- sumption of only two levels of memory. Although mod- els incorporating multiple levels of caches may be nec- essary to analyze some algorithms, for cache-oblivious algorithms, analysis in the two-level ideal-cache model suffices. Specifically, optimal cache-oblivious algo- rithms also perform optimally in computers with mul- tiple levels of LRU caches. We assume that the caches satisfy the inclusion property [20, p. 723], which says that the values stored in cache i are also stored in cache i � 1 (where cache 1 is the cache closest to the proces- sor). We also assume that if two elements belong to the same cache line at level i, then they belong to the same line at level i � 1. Moreover, we assume that cache i � 1 has strictly more cache lines than cache i. These as- sumptions ensure that cache i � 1 includes the contents of cache i plus at least one more cache line.

The multilevel LRU cache operates as follows. A hit on an element in cache i is served by cache i and is not seen by higher-level caches. We consider a line in cache i � 1 to be marked if any element stored on the line be- longs to cache i. When cache i misses on an access, it recursively fetches the needed line from cache i � 1, re- placing the least-recently accessed unmarked cache line. The replaced cache line is then brought to the front of cache 8 i � 1 : ’s LRU list. Because marked cache lines are

never replaced, the multilevel cache maintains the inclu- sion property. The next lemma, whose proof is omitted, asserts that even though a cache in a multilevel model does not see accesses that hit at lower levels, it neverthe- less behaves like the first-level cache of a simple two- level model, which sees all the memory accesses.

Lemma 14 A 8 Zi 9 Li : -cache at a given level i of a mul- tilevel LRU model always contains the same cache lines as a simple 8 Zi 9 Li : -cache managed by LRU that serves the same sequence of memory accesses.

Lemma 15 An optimal cache-oblivious algorithm whose cache complexity satisifies the regularity condi- tion (14) incurs an optimal number of cache misses on each level3 of a multilevel cache with LRU replacement.

Proof. Let cache i in the multilevel LRU model be a8 Zi 9 Li : cache. Lemma 14 says that the cache holds ex- actly the same elements as a 8 Zi 9 Li : cache in a two-level LRU model. From Corollary 13, the cache complex- ity of a cache-oblivious algorithm working on a 8 Zi 9 Li : LRU cache lower-bounds that of any cache-aware algo- rithm for a 8 Zi 9 Li : ideal cache. A 8 Zi 9 Li : level in a mul- tilevel cache incurs at least as many cache misses as a 8 Zi 9 Li : ideal cache when the same algorithm is executed.

Finally, we remove the two assumptions of automatic replacement and full associativity. Specifically, we shall show that a fully associative LRU cache can be main- tained in ordinary memory with no asymptotic loss in expected performance.

Lemma 16 A 8 Z 9 L : LRU-cache can be maintained us- ing O 8 Z : memory locations such that every access to a cache line in memory takes O 8 1 : expected time. Proof. Given the address of the memory location to be accessed, we use a 2-universal hash function [24, p. 216] to maintain a hash table of cache lines present in the memory. The Z

� L entries in the hash table

point to linked lists in a heap of memory that contains Z � L records corresponding to the cache lines. The 2-

universal hash function guarantees that the expected size of a chain is O 8 1 : . All records in the heap are organized as a doubly linked list in the LRU order. Thus, the LRU policy can be implemented in O 8 1 : expected time using O 8 Z � L : records of O 8 L : words each.

3 Alpern, Carter and Feig [5] show that optimality on each level of memory in the UMH model does not necessarily imply global optimal- ity. The UMH model incorporates a single cost measure that combines the costs of work and cache faults at each of the levels of memory. By analyzing the levels independently, our multilevel ideal-cache model remains agnostic about the various schemes by which work and cache faults might be combined.

Theorem 17 An optimal cache-oblivious algorithm whose cache-complexity bound satisfies the regularity condition (14) can be implemented optimally in expec- tation in multilevel models with explicit memory man- agement.

Proof. Combine Lemma 15 and Lemma 16.

Corollary 18 The recursive cache-oblivious algorithms for matrix multiplication, matrix transpose, FFT, and sorting are optimal in multilevel models with explicit memory management.

Proof. Their complexity bounds satisfy the regularity condition (14).

It can also be shown [26] that cache-oblivous algo- rithms satisfying (14) are also optimal (in expectation) in the previously studied SUMH [5, 34] and HMM [1] models. Thus, all the algorithmic results in this paper apply to these models, matching the best bounds previ- ously achieved.

Other simulation results can be shown. For example, by using the copying technique of [22], cache-oblivious algorithms for matrix multiplication and other problems can be designed that are provably optimal on direct- mapped caches.

7. Related work In this section, we discuss the origin of the notion of cache-obliviousness. We also give an overview of other hierarchical memory models.

Our research group at MIT noticed as far back as 1994 that divide-and-conquer matrix multiplication was a cache-optimal algorithm that required no tuning, but we did not adopt the term “cache-oblivious” until 1997. This matrix-multiplication algorithm, as well as a cache- oblivious algorithm for LU-decomposition without piv- oting, eventually appeared in [9]. Shortly after leaving our research group, Toledo [32] independently proposed a cache-oblivious algorithm for LU-decomposition with pivoting. For n � n matrices, Toledo’s algorithm uses Θ 8 n3 : work and incurs Θ 8 1 � n2 � L � n3 � L � Z : cache misses. More recently, our group has produced an FFT library called FFTW [18], which in its most recent incar- nation [17], employs a register-allocation and schedul- ing algorithm inspired by our cache-oblivious FFT al- gorithm. The general idea that divide-and-conquer en- hances memory locality has been known for a long time [29].

Previous theoretical work on understanding hierar- chical memories and the I/O-complexity of algorithms has been studied in cache-aware models lacking an auto- matic replacement strategy, although [10, 28] are recent

0.05

0.1

0.15

0.2

0.25

0 200 400 600 800 1000 1200

T im

e (

m ic

ro se

co n

d s)

iterative recursive

Figure 4: Average time to transpose an N � N matrix, divided by N 2.

exceptions. Hong and Kung [21] use the red-blue peb- ble game to prove lower bounds on the I/O-complexity of matrix multiplication, FFT, and other problems. The red-blue pebble game models temporal locality using two levels of memory. The model was extended by Savage [27] for deeper memory hierarchies. Aggarwal and Vitter [3] introduced spatial locality and investigated a two-level memory in which a block of P contiguous items can be transferred in one step. They obtained tight bounds for matrix multiplication, FFT, sorting, and other problems. The hierarchical memory model (HMM) by Aggarwal et al. [1] treats memory as a linear array, where the cost of an access to element at location x is given by a cost function f 8 x : . The BT model [2] extends HMM to support block transfers. The UMH model by Alpern et al. [5] is a multilevel model that allows I/O at different levels to proceed in parallel. Vitter and Shriver introduce parallelism, and they give algorithms for ma- trix multiplication, FFT, sorting, and other problems in both a two-level model [35] and several parallel hierar- chical memory models [36]. Vitter [33] provides a com- prehensive survey of external-memory algorithms.

8. Conclusion The theoretical work presented in this paper opens two important avenues for future research. The first is to determine the range of practicality of cache-oblivious algorithms, or indeed, of any algorithms developed in the ideal-cache model. The second is to resolve, from a complexity-theoretic point of view, the relative strengths of cache-oblivious and cache-aware algorithms. To con- clude, we discuss each of these avenues in turn.

Figure 4 compares per-element time to transpose a matrix using the naive iterative algorithm employing a doubly nested loop with the recursive cache-oblivious REC-TRANSPOSE algorithm from Section 3. The two algorithms were evaluated on a 450 megahertz AMD K6III processor with a 32-kilobyte 2-way set-associative L1 cache, a 64-kilobyte 4-way set-associative L2 cache, and a 1-megabyte L3 cache of unknown associativ- ity, all with 32-byte cache lines. The code for REC- TRANSPOSE was the same as presented in Section 3, ex-

0.02

0.04

0.06

0.08

0.1

0.12

0 100 200 300 400 500 600

T im

e (

m ic

ro se

co n

d s)

iterative recursive

Figure 5: Average time taken to multiply two N � N matrices, divided by N 3.

cept that the divide-and-conquer structure was modified to produce exact powers of 2 as submatrix sizes wher- ever possible. In addition, the base cases were “coars- ened” by inlining the recursion near the leaves to in- crease their size and overcome the overhead of proce- dure calls. (A good research problem is to determine an effective compiler strategy for coarsening base cases automatically.)

Although these results must be considered prelimi- nary, Figure 4 strongly indicates that the recursive al- gorithm outperforms the iterative algorithm throughout the range of matrix sizes. Moreover, the iterative al- gorithm behaves erratically, apparently due to so-called “conflict” misses [20, p. 390], where limited cache asso- ciativity interacts with the regular addressing of the ma- trix to cause systematic interference. Blocking the itera- tive algorithm should help with conflict misses [22], but it would make the algorithm cache aware. For large ma- trices, the recursive algorithm executes in less than 70% of the time used by the iterative algorithm, even though the transpose problem exhibits no temporal locality.

Figure 5 makes a similar comparison between the naive iterative matrix-multiplication algorithm, which uses three nested loops, with the O 8 n3 : -work recur- sive REC-MULT algorithm described in Section 2. This problem exhibits a high degree of temporal locality, which REC-MULT exploits effectively. As the figure shows, the average time used per integer multiplication in the recursive algorithm is almost constant, which for large matrices, is less than 50% of the time used by the iterative variant. A similar study for Jacobi multipass filters can be found in [26].

Several researchers [12, 16] have also observed that recursive algorithms exhibit performance advantages over iterative algorithms for computers with caches. A comprehensive empirical study has yet to be done, how- ever. Do cache-oblivious algorithms perform nearly as well as cache-aware algorithms in practice, where con- stant factors matter? Does the ideal-cache model cap- ture the substantial caching concerns for an algorithms designer?

An anecdotal affirmative answer to these questions is exhibited by the popular FFTW library [17, 18], which

uses a recursive strategy to exploit caches in Fourier transform calculations. FFTW’s code generator pro- duces straight-line “codelets,” which are coarsened base cases for the FFT algorithm. Because these codelets are cache oblivious, a C compiler can perform its register allocation efficiently, and yet the codelets can be gen- erated without knowing the number of registers on the target architecture.

To close, we mention two theoretical avenues that should be explored to determine the complexity- theoretic relationship between cache-oblivious algo- rithms and cache-aware algorithms.

Separation: Is there a gap in asymptotic complexity between cache-aware and cache-oblivious algorithms? It appears that cache-aware algorithms should be able to use caches better than cache-oblivious algorithms, since they have more knowledge about the system on which they are running. Do there exist problems for which this advantage is asymptotically significant, for example an Ω 8 lg Z : advantage? Bilardi and Peserico [8] have re- cently taken some steps in proving a separation.

Simulation: Is there a limit as to how much better a cache-aware algorithm can be than a cache-oblivious algorithm for the same problem? That is, given a class of optimal cache-aware algorithms to solve a single problem, can we construct a good cache-oblivious al- gorithm that solves the same problem with only, for example, O 8 lg Z : loss of efficiency? Perhaps simula- tion techniques can be used to convert a class of effi- cient cache-aware algorithms into a comparably efficient cache-oblivious algorithm.

Acknowledgments

Thanks to Bobby Blumofe, now of the University of Texas at Austin, who sparked early discussions at MIT about what we now call cache obliviousness. Thanks to Gianfranco Bilardi of University of Padova, Sid Chat- terjee of University of North Carolina, Chris Joerg of Compaq CRL, Martin Rinard of MIT, Bin Song of MIT, Sivan Toledo of Tel Aviv University, and David Wise of Indiana University for helpful discussions. Thanks also to our anonymous reviewers.

References [1] A. Aggarwal, B. Alpern, A. K. Chandra, and M. Snir.

A model for hierarchical memory. In Proceedings of the 19th Annual ACM Symposium on Theory of Computing (STOC), pages 305–314, May 1987.

[2] A. Aggarwal, A. K. Chandra, and M. Snir. Hierarchi- cal memory with block transfer. In 28th Annual Sym- posium on Foundations of Computer Science (FOCS), pages 204–216, Los Angeles, California, 12–14 Oct. 1987. IEEE.

[3] A. Aggarwal and J. S. Vitter. The input/output complex- ity of sorting and related problems. Communications of the ACM, 31(9):1116–1127, Sept. 1988.

[4] A. V. Aho, J. E. Hopcroft, and J. D. Ullman. The Design and Analysis of Computer Algorithms. Addison-Wesley Publishing Company, 1974.

[5] B. Alpern, L. Carter, and E. Feig. Uniform memory hier- archies. In Proceedings of the 31st Annual IEEE Sym- posium on Foundations of Computer Science (FOCS), pages 600–608, Oct. 1990.

[6] D. H. Bailey. FFTs in external or hierarchical memory. Journal of Supercomputing, 4(1):23–35, May 1990.

[7] L. A. Belady. A study of replacement algorithms for vir- tual storage computers. IBM Systems Journal, 5(2):78– 101, 1966.

[8] G. Bilardi and E. Peserico. Efficient portability across memory hierarchies. Unpublished manuscript, 1999.

[9] R. D. Blumofe, M. Frigo, C. F. Joerg, C. E. Leiserson, and K. H. Randall. An analysis of dag-consistent dis- tributed shared-memory algorithms. In Proceedings of the Eighth Annual ACM Symposium on Parallel Algo- rithms and Architectures (SPAA), pages 297–308, Padua, Italy, June 1996.

[10] L. Carter and K. S. Gatlin. Towards an optimal bit- reversal permutation program. In Proceedings of the 39th Annual Symposium on Foundations of Computer Science, pages 544–555. IEEE Computer Society Press, 1998.

[11] S. Chatterjee, V. V. Jain, A. R. Lebeck, and S. Mundhra. Nonlinear array layouts for hierarchical memory sys- tems. In Proceedings of the ACM International Confer- ence on Supercomputing, Rhodes, Greece, June 1999.

[12] S. Chatterjee, A. R. Lebeck, P. K. Patnala, and M. Thot- tethodi. Recursive array layouts and fast parallel matrix multiplication. In Proceedings of the Eleventh ACM Sym- posium on Parallel Algorithms and Architectures (SPAA), June 1999.

[13] J. W. Cooley and J. W. Tukey. An algorithm for the ma- chine computation of the complex Fourier series. Mathe- matics of Computation, 19:297–301, Apr. 1965.

[14] T. H. Cormen, C. E. Leiserson, and R. L. Rivest. In- troduction to Algorithms. MIT Press and McGraw Hill, 1990.

[15] P. Duhamel and M. Vetterli. Fast Fourier transforms: a tutorial review and a state of the art. Signal Processing, 19:259–299, Apr. 1990.

[16] J. D. Frens and D. S. Wise. Auto-blocking matrix- multiplication or tracking BLAS3 performance from source code. In Proceedings of the Sixth ACM SIG- PLAN Symposium on Principles and Practice of Parallel Programming (PPoPP), pages 206–216, Las Vegas, NV, June 1997.

[17] M. Frigo. A fast Fourier transform compiler. In Proceed- ings of the ACM SIGPLAN’99 Conference on Program- ming Language Design and Implementation (PLDI), At- lanta, Georgia, May 1999.

[18] M. Frigo and S. G. Johnson. FFTW: An adaptive soft- ware architecture for the FFT. In Proceedings of the In- ternational Conference on Acoustics, Speech, and Signal Processing, Seattle, Washington, May 1998.

[19] G. H. Golub and C. F. van Loan. Matrix Computations. Johns Hopkins University Press, 1989.

[20] J. L. Hennessy and D. A. Patterson. Computer Architec- ture: A Quantitative Approach. Morgan Kaufmann, 2nd edition, 1996.

[21] J.-W. Hong and H. T. Kung. I/O complexity: the red-blue pebbling game. In Proceedings of the 13th Annual ACM Symposium on Theory of Computing (STOC), pages 326– 333, Milwaukee, 1981.

[22] M. S. Lam, E. Rothberg, and M. E. Wolf. The cache per- formance and optimizations of blocked algortihms. In Fourth International Conference on Architectural Sup- port for Programming Languages and Operating Systems (ASPLOS), pages 63–74, Santa Clara, CA, Apr. 1991. ACM SIGPLAN Notices 26:4.

[23] A. LaMarca and R. E. Ladner. The influence of caches on the performance of sorting. Proceedings of the Eighth Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), pages 370–377, 1997.

[24] R. Motwani and P. Raghavan. Randomized Algorithms. Cambridge University Press, 1995.

[25] M. H. Nodine and J. S. Vitter. Deterministic distribution sort in shared and distributed memory multiprocessors. In Proceedings of the Fifth ACM Symposium on Paral- lel Algorithms and Architectures (SPAA), pages 120–129, Velen, Germany, 1993.

[26] H. Prokop. Cache-oblivious algorithms. Master’s thesis, Massachusetts Institute of Technology, June 1999.

[27] J. E. Savage. Extending the Hong-Kung model to mem- ory hierarchies. In D.-Z. Du and M. Li, editors, Com- puting and Combinatorics, volume 959 of Lecture Notes in Computer Science, pages 270–281. Springer Verlag, 1995.

[28] S. Sen and S. Chatterjee. Towards a theory of cache- efficient algorithms. Unpublished manuscript, 1999.

[29] R. C. Singleton. An algorithm for computing the mixed radix fast Fourier transform. IEEE Transactions on Audio and Electroacoustics, AU-17(2):93–103, June 1969.

[30] D. D. Sleator and R. E. Tarjan. Amortized efficiency of list update and paging rules. Communications of the ACM, 28(2):202–208, Feb. 1985.

[31] V. Strassen. Gaussian elimination is not optimal. Nu- merische Mathematik, 13:354–356, 1969.

[32] S. Toledo. Locality of reference in LU decomposition with partial pivoting. SIAM Journal on Matrix Analysis and Applications, 18(4):1065–1081, Oct. 1997.

[33] J. S. Vitter. External memory algorithms and data struc- tures. In J. Abello and J. S. Vitter, editors, External Memory Algorithms and Visualization, DIMACS Series in Discrete Mathematics and Theoretical Computer Sci- ence. American Mathematical Society Press, Providence, RI, 1999.

[34] J. S. Vitter and M. H. Nodine. Large-scale sorting in uniform memory hierarchies. Journal of Parallel and Distributed Computing, 17(1–2):107–114, January and February 1993.

[35] J. S. Vitter and E. A. M. Shriver. Algorithms for par- allel memory I: Two-level memories. Algorithmica, 12(2/3):110–147, August and September 1994.

[36] J. S. Vitter and E. A. M. Shriver. Algorithms for parallel memory II: Hierarchical multilevel memories. Algorith- mica, 12(2/3):148–169, August and September 1994.

[37] S. Winograd. On the algebraic complexity of functions. Actes du Congrès International des Mathématiciens, 3:283–288, 1970.