Lecture 15: Cache-Oblivious Algorithms

• Download English-US transcript (PDF)

ANNOUNCER: The following content is provided under a Creative Commons license. Your support will help MIT Open Courseware continue to offer high quality educational resources for free. To make a donation or to view additional materials from hundreds of MIT courses, visit MIT OpenCourseWare at ocw.mit.edu.

JULIAN SHUN: All right. So today we're going to talk about cache oblivious algorithms. Who remembers what a cache oblivious algorithm is? So what is a cache oblivious algorithm oblivious to? Cache size.

So a cache oblivious algorithm is an algorithm that automatically tunes for the cache size on your machine. So it achieves good cache efficiency. And the code doesn't need to have any knowledge of the cache parameters of your machine. In contrast, a cache-aware algorithm would actually know the parameters of the cache sizes on your machine. And the code would actually put the size of the cache inside.

So today we're going to talk a lot more about cache oblivious algorithms. Last time we talked about one cache oblivious algorithm that was for matrix multiplication. And today we're going to talk about some other ones. So first example I want to talk about is simulation of heat diffusion. So here's a famous equation known as the heat equation.

And this equation is in two dimensions. And we want to simulate this function u that has three parameters t, x, and y. t is the time step. X and y are the x, y coordinates of the 2D space. And we want to know the temperature for each x, y coordinate for any point in time t.

And the 2D heat equation can be modeled using a differential equation. So how many of you have seen differential equations before? OK. So good, most of you. So here I'm showing the equation in two dimensions. But you can get similarly equations for higher dimensions.

So here it says the partial derivative of u with respect to t is equal to alpha. Alpha is what's called the thermo diffusivity. It's a constant times the sum of the second partial derivative of u with respect to x, and the second partial derivative of u with respect to y.

So this is a pretty famous equation. And you can see that we have a single derivative on the left side and a double derivative on the right side. And how do we actually write code to simulate this 2D heat process? So oftentimes in scientific computing, people will come up with these differential equations just to describe physical processes. And then they want to come up with efficient code to actually simulate the physical process.

OK. So here's an example of a 2D heat diffusion. So let's say we started out with a configuration on the left. And here the color corresponds to a temperature. So a brighter color means it's hotter. Yellow is the hottest and blue is the coldest.

In on the left we just have 6172, which is the course number. So if you didn't know that, you're probably in the wrong class. And then afterwards we're going to run it for a couple time steps. And then the heat is going to diffuse to the neighboring regions of the 2D space. So after you run it for a couple of steps, you might get the configuration on the right where the heat is more spread out now.

And oftentimes, you want to run this simulation for a number of time steps until the distribution of heat converges so it becomes stable and that doesn't change by much anymore. And then we stop the simulation. So this is the 1D heat equation. I showed you a 2D one earlier.

But we're actually going to generate code for the 1D heat equation since it's simpler. But all the ideas generalize to higher dimensions. And here's the range of colors corresponding to temperature, so the hottest colors on the left and the coldest colors on the right. And if you had a heat source that's on the left hand side of this bar, then this might possibly be a stable distribution. So if you keep running the simulation, you might get a stable distribution of heat that looks like this.

OK. So how do we actually write code to simulate this differential equation? So one commonly used method is known as finite difference approximation. So we're going to approximate the partial derivative of u with respect to each of its coordinates. So the partial derivative of u with respect to t is approximately equal to u of t plus delta t where delta t is some small value, and x, and then minus u of tx. And then that's all divided by delta t.

So how many of you have seen this approximation method before from your calculus class? OK, good. So as you bring the value of delta t down to 0, then the thing on the right hand side approach is the true partial derivative. So that's a partial derivative with respect to t. We also need to get the partial derivative with respect to x.

And here I'm saying the partial derivative with respect to x is approximately equal to ut of x plus delta x over 2 minus ut of x minus delta x over 2, all divided by delta x. So notice here that instead of adding delta x in the first term and not adding anything in the second term, I'm actually adding delta x over 2 in the first term and subtracting text over 2 in the second term. And it turns out that I can do this with the approximation method.

And it still turns out to be valid as long as the two terms that I'm putting in, their difference is delta x. So here the difference is delta x. And I can basically decide how to split up this delta x term among the two things in the numerator.

And the reason why I chose delta x over 2 here is because the math is just going to work out nicely. And it's going to give us cleaner code. This is just the first partial derivative with respect to x. Actually need the second partial derivative since the right hand side of this equation has the second partial derivative.

So this is what the second partial derivative looks like. So I just take the partial derivative with respect to x of each of the terms in my numerator from the equation above. And then now I can actually plug in the value of this partial derivative by applying the equation above using the arguments t and x plus delta x over 2, and similarly for the second term.

So for the first term when I plug it into the equation for the partial derivative with respect to x, I'm just going to get ut of x plus delta x minus utx. And then for the second term, I'm going to get ut of x minus delta x. And then I subtract another factor of utx.

So that's why I'm subtracting 2 times utx in the numerator here. And then the partial derivative of each of the things in a numerator also have to divide by this delta x term. So on the denominator, I get delta x squared.

So now I have the second partial derivative with respect to x. And I also have the first partial derivative with respect to t. So I can just plug them into my equation above.

So on the left hand side I just have this term here. And I'm multiplying by this alpha constant. And then on this term just comes from here. So this is what the 1d heat equation reduces to using the finite difference approximation method. So any questions on this

So how do we actually write code to simulate this equation here? So we're going to use what's called a stencil computation. And here I'm going to set delta at x and delta t equal to 1, just for simplicity. But in general you can set them to whatever you want. You can make them smaller to have a more fine grained simulation.

So my set delta x and delta eat t equal to 1. Then the denominators of these two terms just become one and I don't need to worry about them. And then I'm going to represent my 2D space using a 2D matrix where the horizontal axis represents values of x, and the vertical axis represents values of t. And I want to fill in all these entries that have a black dot in it. The ones with the orange dot, those are my boundaries.

So those actually fixed throughout the computation. So I'm not going to do any computation on those. Those are just given to me as input. So they could be heat sources if we're doing the heat simulation.

And then now I can actually write code to simulate this equation. So if I want to compute u of t plus 1x, I can just go up here. And I see that that's equal to this thing over here. And then I bring the negative utx term to the right.

So I get ut of x plus alpha times ut x plus 1 minus 2 times utx plus ut x minus 1. As I said before, we just want to keep iterating this until the temperatures becomes stable. So I'm going to proceed in time, which in time is going up here.

And to compute one of these points-- so let's say this is ut plus 1x-- I need to know the value of utx, which is just a thing below me in the matrix. And I also need to know utx plus 1 and utx minus 1. And those are just the things below me and diagonal to either the left or the right side.

So each value here just depends on three other values. And this is called a three point stencil. This is the pattern that this equation is representing. And in general, a stencil computation is going to update each point in an array using a fixed pattern. This is called a stencil.

So I'm going to do the same thing for all of the other points. And here, I'm going to compute all the values of x for a given time step. And then I move on to the next time step.

And then I keep doing this until the distribution of temperatures becomes stable. And then I'm done. OK. So these stencil computations are widely used in scientific computing.

They're used for weather simulations, stock market simulations, fluid dynamics, image processing probability, and so on. So they're used all over the place in science. So this is a very important concept to know.

So let's say I just ran the code as I showed you in the animation. So I completed one row at a time before I moved on to the next row. How would this code perform with respect to caching? Yes?

AUDIENCE: I think if x is less than a third of the cache size, [INAUDIBLE]

JULIAN SHUN: Yeah. So if x is small, this would do pretty well. But what if x is much larger than the size of your cache?

AUDIENCE: [INAUDIBLE].

JULIAN SHUN: Yeah, you. Would do badly and why is that?

AUDIENCE: Because the whole [INAUDIBLE] second row, [INAUDIBLE]

JULIAN SHUN: Yeah. So it turns out that there could be some reuse here, because when I compute the second row, I'm actually using some values that are computed in the first row. But if the row is much larger than the cache size, by the time I get to the second row, the values that I loaded into cache from the first row would have been evicted. And therefore, I'm going to suffer a cache miss again for the second row, for the values I need to load in, even though they could have been used if we made our code have good locality.

Another question I have is if we only cared about the values of x at the most recent time step, do we actually have to keep around this whole 2D matrix? Or can we get by with less storage? Yeah. So how many rows would I have to keep around if I only cared about the most recent time step? Yeah?

AUDIENCE: Two.

JULIAN SHUN: Two. And why is that?

AUDIENCE: So one for the previous time step. One for the current time step. [INAUDIBLE]

JULIAN SHUN: Right. So I need to keep around two rows because I need the row from the previous time step in order to compute the values in the current time step. And after the current time step, I can just swap the roles of the two rows that I'm keeping around, and then reuse the previous row for the next one. Yes?

AUDIENCE: Would you only need one and then a constant amount of extra space, like if you had three extra things, you could probably do it with one.

JULIAN SHUN: So I need to know-- when I'm computing the second row, I need to keep around all of these values that I computed in the first row, because these values get fed to one of the computations in the second row. So I need to actually keep all of them around.

AUDIENCE: I think if you iterate to the right, then you have three that are this one and the next one. Just three. Then you can--

JULIAN SHUN: Oh, I see I see what you're saying. Yeah. So that's actually a good observation. So you only need to keep a constant amount more storage, because you'll just be overwriting the values as you go through the row. So if you keep around one row, some of the values would be for the current time step, and some of them would be from the previous time step. So that's a good observation, yes.

OK. So that code, as we saw, it wasn't very cache efficient. You could make a cache efficient using tiling. But we're going to go straight to the cache oblivious algorithm because it's much cleaner. So let's recall the ideal cache model. We talked about this in the previous lecture.

So here we have a two level hierarchy. We have a cache. And then we have the main memory. The cache has size of M bytes. And a cache line is B bytes.

So you can keep around M over B cache lines in your cache. And if something's in cache and you operate on it, then it doesn't incur any cache misses. But if you have to go to main memory to load the cache line in, then you incur one cache miss.

The ideal cache model assumes that the cache is fully associative. So any cache line can go anywhere in the cache. And it also assumes either an optimal omniscient replacement policy or the LRU policy.

So the optimal omniscient replacement policy knows the sequence of all future requests to memory. And when it needs to evict something, it's going to pick the thing that leads to the fewest cache misses overall to evict. The LRU policy just evict the thing that was least recently used.

But we saw from the previous lecture, that in terms of asymptotic costs, these two replacement policies will give you cache misses within a constant fact of each other. So you can use either one, depending on what's convenient.

And two performance measures that we care about when we're analyzing an algorithm and the ideal cache model are the work and the number of cache misses. So the work is just the total number of operations that the algorithm incurs. And serially, this is just the ordinary running time. And the number of cache misses is the number of cache lines you have to transfer between the cache and your main memory.

So let's assume that we're running an algorithm or analyzing an algorithm in the ideal cache model, and it runs serially. What kinds of cache misses does the ideal cache model not capture? So remember, we talked about several types of cache misses. And there's one type of cache miss that this model doesn't capture when we're running serially.

So let's assume we're running this serially without any parallelism here. So the sharing misses has only come about when you have parallelism. Yes?

AUDIENCE: Conflictness?

JULIAN SHUN: Yes. So the answer is conflictnesses. And why is that? Why does this model not capture it?

AUDIENCE: There's not a specific sets that could get replaced then since it's fully associated.

JULIAN SHUN: Yes. So this is a fully associative cache. So any cache line can go anywhere in the cache. And you can only get conflict misses for set associated schemes where each cache line can only be mapped to a particular set.

And if you have too many cache lines that map to that particular set, then you're going to keep evicting each other even though the rest of the cache could have space. And that's what's called a conflict miss. The ideal cash model does capture capacity misses.

So therefore, it is still a very good model to use at a high level when you're designing efficient algorithms, because it encourages you to optimize for spatial and temporal locality. And once you have a good algorithm in the ideal cache model then you can start dealing with conflict misses using some of the strategies that we talked about last time such as padding or using temporary memory. So any questions on this?

OK. So this is the code that does the heat simulation that we saw earlier. So it's just two for loops, a nested for loop. In the outer loop, we're looping over the time dimension. In the inner loop, we're looping over the space dimension. So we're computing all the values of x before we move on to the next time step.

And then we're storing two rows here and we're using this trick called a even odd trick. And here's how it works. So to access the next row that we want to compute, that we just do a t plus 1 mod 2. And then to access the current row, it's just t mod 2. So this is implicitly going to swap the roles of the two rows that we're keeping around as we progress through time.

And then we're going to set u of t plus 1 mod 2x equal to kernel of u-- pointer to ut mod 2x. And this kernel function is defined up here. And recall, that when we're actually passing a pointer to this kernel function, we can actually treat a pointer as the beginning of an array. So we're using array notation up here inside the kernel function. So the array W is passed as input.

And then we need to return W of 0. That's just the element at the current pointer that we passed to kernel. And then we add alpha times w of negative 1. That's one element before the thing that we're pointing to, minus 2 times W of 0 plus W of 1. W of 1 is the next element that we're pointing to.

OK. So let's look at the caching behavior of this code. So we're going to analyze the cache complexity. And we're going to assume the LRU replacement policy here, because we can.

And as we said before, we're going to loop through one entire row at a time before we go onto the next row. So the number of cache misses I get, assuming that n is greater than M, so that the row size is greater than the cache size, the number of cache misses is theta of NT over B. So how do I get this cache complexity around here?

So how many cache misses do I have to incur for each row of this 2D space that I'm computing? Yes?

AUDIENCE: N over B.

JULIAN SHUN: Right. So I need N over B cache misses for each row. And this is because I can load in B bytes at a time. So I benefit from spatial locality there. And then I have N elements I need to compute. So it's theta of N over B per row.

And as we said before, when we get to the next row, the stuff that we need from the previous row have already been evicted from cache. So I basically have to incur theta of N over B cache misses for every row. And the number of rows I'm going to compute as t. So it's just theta of NT over B. Any questions on this analysis?

So how many of you think we can do better than this? OK. So one person. Two, three. OK. So turns out that we can do better than this. You can actually do better with tiling, but I'm not going to do the tiling version.

I want to do the cache oblivious version. And the cache oblivious version is going to work on trapezoidal regions in the 2D space. And recall that a trapezoid has a top base and a bottom base. And here the top base is at t1, the bottom base is at t0, and the height is just t1 minus t0. And the width of a trapezoid is just the width of it at the midpoint between t1 and t0, so at t1 plus t0 over 2.

So we're going to compute all of the points inside this trapezoid that satisfy these inequalities here. So t has to be greater than or equal to t0 and less than t1. And then x is greater than or equal to x0 plus dx0 times t minus t0. So dx0 is actually the inverse slope here. And then it also has to be less than x1 plus dx1 times t minus t0.

So dx1 is the inverse slope on the other side. And dx0 and dx1 have to be either negative 1, 0, or 1. So negative 1 just corresponds to inverse slope of negative 1, which is also a slope of negative 1. If it's 1, then it's just a slope or inverse of 1. And then if it's 0, then we just have a vertical line.

OK. So the nice property of this trapezoid is that we can actually compute everything inside the trapezoid without looking outside the trapezoid. So we can compute everything here independently of any other trapezoids we might be generating. And we're going to come up with a divide and conquer approach to execute this code. So the divide and conquer algorithm has a base case. So our base case is going to be when the height of the trapezoid is 1.

And when the height is 1, then we're just going to compute all of the values using a simple loop. And any order if the computation inside this loop is valid, since we have all the values in the base of the trapezoid and we can compute the values in the top of the trapezoid in whatever order. They don't depend on each other. So that's a base case. Any questions so far?

So here's one of the recursive cases. It turns out that we're going to have two different types of cuts. The first cut is called a space cut. So I'm going to do a space cut if the width of the trapezoid is greater than or equal to twice the height.

So this means that the trapezoid is too wide. And I'm going to cut it vertically. More specifically, I'm going to cut it with a line, with slope negative 1 going through the center of the trapezoid. And then I'm going to traverse the trapezoid on the left side first. And then after I'm done with that, traverse the trapezoid on the right side.

So can I actually switch the order of this? Can I compute the stuff on the right side before I do this stuff on the left side? No. Why is that?

AUDIENCE: [INAUDIBLE].

JULIAN SHUN: Yeah. So there some points in the right trapezoid that depend on the values from the left trapezoid. And so for the left trapezoid, every point we want to compute, we already have all of its points, assuming that we get all the values of the base points. But for the right hand side, some of the values depend on values in the left trapezoid.

So we can't execute the right trapezoid until we're done with the left trapezoid. And this is the reason why I cut this trapezoid with a slope of negative 1 instead of using a vertical cut. Because if I did a vertical cut then inside both of the trapezoids, I would have points that depend on the other trapezoid.

So this is one of the two cuts. This is called a space cut. And it happens when the trapezoid is too wide. The other cut is the time cut I'm going to cut with respect to the time dimension. And this happens when the trapezoid is too tall, so when the width is less than twice the height of the trapezoid.

Then what I'm going to do is I'm just going to cut it with a horizontal line through the center. And then I'm going to traverse the bottom trapezoid first. And after I'm done with that, I can traverse a top trapezoid. And again, the top trapezoid depends on some points from the bottom trapezoid. So it's I can't switch the order of those. Any questions?

OK. So let's now look at the code that implements this recursive divide and conquer algorithm. So here's the C code. It takes as input t0 and t1. These are the coordinates of the top and the bottom up the trapezoid, or bottom and top of the trapezoid, then x. 0 is the left side of the trapezoid-- of the base of the trapezoid.

dx0 is the inverse slope, and the x1 is the right side of the bottom base of the trapezoid, and dx1 is the inverse slope on the right side. So we're first going to compute the height of our trapezoid. And we're going to let LT be the height. And that's just t1 minus t0. And if the height is 1, then we're just going to use a for loop over all the points in that height 1 trapezoid.

We're going to call this kernel function that we defined before. And otherwise, the height is greater than 1. And we're going to check whether we should do a space cut or time cut. So we do a space cut if the trapezoid is too wide. And this condition inside the if clause is checking if the trapezoid is too wide.

And if so, then we're going to make two recursive calls to trapezoid. And we're going to cut it in the middle using this slope of negative 1. So you see the negative ones here in the recursive calls. And otherwise, we'll do a time cut.

And for the time cut we just cut it in the middle. So we cut it at the value of t that's equal to LT divided by, 2 or t0 plus LT divided by 2. And again, we have two recursive calls to trapezoid.

OK. So even though I'm only generating slopes of negative 1 in this recursive call, this code is going to work even if I have slopes of 1 and 0, because I could start out with slopes of 1 and 0. For example, if I had a rectangular region, then the slopes are just going to be 0, and this code is still going to work. But most of the slopes that I'm dealing with are going to be a negative 1, because those are the new slopes and I'm generating. Any questions?

So this code is very concise. It turns out that even, odd tricks still works here. You can still keep around just two rows, because you're guaranteed that you're not going to overwrite any value until all the things that depend on it are computed. But when you're using just two rows, the values in a particular row might not all correspond to the same time step, because we're not finishing an entire row before we move on to the next one. We're actually partially computing rows.

OK. So let's analyze the cash complexity of this algorithm. Again, we're going to use the recursion tree approach that we talked about last time. And our code is going to split itself into two cell problems at every level until it gets to a leaf. And each leaf represents theta of hw points where h is theta of w. So h is the height. w is the width.

And we have the property that h is equal to theta w because of the nature of the algorithm. When the trapezoid becomes too wide, we're going to make it less wide by doing a space cut. And if it becomes too tall, we're going to do a horizontal cut. So we're guaranteed that the height and the width are going to be within a constant factor of each other when we get to the base case. And each leaf in the base case is going to incur theta of w over B misses because we have to load in the base of the trapezoid from main memory.

And once that's in cache, we can compute all of the other points in the trapezoid without incurring any more cache misses. So the cache misses per leaf is just theta w over B. And we're going to set w equal to theta of M in the analysis, because that's the point when the trapezoid fits into cache.

The algorithm doesn't actually have any knowledge of this M parameter. So it's still going to keep divide and conquering until it gets the base case of size 1. But just for the analysis, we're going to use a base case when w is theta of M.

So the number of leaves we have is theta of NT divided by w because each leaf is a size theta hw. The number of internal notes we have is equal to a number of leaves minus because we have a tree here. But the internal notes don't contribute substantially to the cache complexity, because each of them is just doing a constant number of cache misses to set up the two recursive calls. And it's not doing anything more expensive than that.

So we just need to compute the number of cache misses at the leaves. We have theta of NT over hw leaves, each one of which takes theta of w over B cache misses. And we're going to multiply that out.

So the w term cancels out. We just have NT over hB and we set h and w to be theta of M. So we're just left with NT over MB as our cache complexity. Yes?

AUDIENCE: Can you explain why hB [INAUDIBLE]?

JULIAN SHUN: Sure. So each leaf only incurs theta w over B cache misses because we need to compute the values of the base of the trapezoid. And that's going to incur theta of w over B cache misses because it's w wide, and therefore, everything else is going to fit into cache. So when we compute them, we already have all of our previous values that we want in cache. So that's why it's not going to incur any more cache misses. So does that makes sense?

AUDIENCE: Yeah, I forgot that it was [INAUDIBLE].

JULIAN SHUN: Yeah.

OK. So this was just analysis for one dimension. But it actually generalizes to more than one dimension. So in general, if we have d dimensions, then the cache complexity is going to be theta of NT divided by M to the 1 over d times B.

So if d is 1, then that just reduces to NT over MB. If d is 2, then it's going to be NT over B times square root of M and so on. And it turns out that this bound is also optimal.

So any questions? So compared to the looping code, this code actually has much better cache complexity. It's saving by a factor of M. The looping code had a cache complexity that was theta of NT over B. It didn't have an M in the denominator. OK.

So we're actually going to do a simulation now. We're going to compare how the looping code in a trapezoid code runs. And in this simulation, the green points correspond to a cache hit, the purple points correspond to a cache miss. And we're going assume a fully associative cache using the LRU replacement policy where the cache line size is 4 points and cache size is 32 points. And we're going to set the cache hit latency to be one cycle, and the cache miss latency to be 10 cycles.

So an order of magnitude slower for cache misses. And we're doing this for a rectangular region where N is 95. And we're doing it for 87 time steps. So when we pull out the simulation now.

OK. So on the left hand side, we're going to have the looping code. On the right hand side, we're going to have the trapezoidal code. So let's start this. So you can see that the looping code is going over one row at a time whereas the trapezoidal code is not doing that. It's partially computing one row and then moving on to the next row.

I can also show the cuts that are generated by the trapezoidal algorithm. I have to remember how to do this. So C--

So there there's the cuts that are generated by the trapezoidal algorithm. And I can speed this up. So you can see that the trapezoidal algorithm is incurring much fewer cache misses than the looping code.

So I just make this go faster. And the trapezoid code is going to finish, while the looping code is still slowly making its way up the top. OK. So it's finally done. So any questions on this? Yeah?

AUDIENCE: Why doesn't the trapezoid [INAUDIBLE]?

JULIAN SHUN: In which of the regions? So it's loading-- you get a cache miss for the purple dots here. And then the cache line size is 4.

So you get a cache miss for the first point, and then you hit on the next three points. Then you get another cache miss for the fifth point. And then you hit on the 6, 7, and 8 points.

AUDIENCE: I was indicating the one above it.

JULIAN SHUN: So for the one above it-- so we're assuming that the two arrays fitting cache already. So we don't actually need to load them from memory. The thing above it just depends on the values that we have already computed.

And that fits in cache. Those are ready in cache. This base of the trapezoid is already in cache. And the row right above it, we just need to look at those values. Does that makes sense?

AUDIENCE: OK. Because of the even odd?

JULIAN SHUN: Yeah. Yeah. Any other questions on this? Does anyone want to see this again? Yeah?

So I could let this run for the rest of the lecture, but I have more interesting material that I want to talk about. So let's just stop after this finishes. And as you see again, the looping code is slowly making its way to the top. It's much slower than the trapezoid code.

OK. So that was only for one dimensions. Now let's look at what happens in two dimensions. And here, I'm going to show another demo. And this demo, I'm actually going to run the code for real. The previous demo was just a simulation.

So this is going to happen in real time. And I'm going to simulate the heat in a 2D space. And recall that the colors correspond to the temperature. So a brighter color means it's hotter. So let me pull out the other demo.

OK. So here, my mouse cursor is the source of heat. So you see that it's making the points around my mouse cursor hot. And then it's slowly diffusing its way to the other points. Now I can actually move this so that my heat source changes, and then the heat I generated earlier slowly goes away. OK.

So in the lower left hand corner, I'm showing the number of iterations per second of the code. And we can see that the looping code is taking-- it's doing about 1,560 iterations per second. Let's see what happens when we switch to the trapezoid code. So the trapezoid code is doing about 1,830 iterations per second. So it's a little bit faster, but not by too much.

Does anyone have an idea why we're seeing this behavior? So we said that the trapezoid code incurs many fewer cache misses than the looping code, so we would expect it to be significantly faster. But here it's only a little bit faster. Yeah?

AUDIENCE: [INAUDIBLE].

JULIAN SHUN: Right. So that's a good point. So in 2D you're only saving a factor of square root of M instead of M. But square root of M is still pretty big compared to the speed up we're getting here. So any other ideas? Yeah.

So there is a constant factor in the trapezoidal code. But even after accounting for the constant factor, you should still see a better speed up than this. So even accounting for the constant factors, what other problem might be going on here? Yeah?

AUDIENCE: Is it that we don't actually have an ideal cache?

JULIAN SHUN: Yeah. So that's another good observation. But the caches that we're using, they still should get pretty good cache. I mean, they should still have cache complexly that's pretty close to the ideal cache model. I mean, you might be off by small constant factor, but not by too much. Yeah?

AUDIENCE: Maybe because [INAUDIBLE]

JULIAN SHUN: Sorry. Could you repeat that?

AUDIENCE: There are [INAUDIBLE] this time like [INAUDIBLE]

JULIAN SHUN: Yeah. So OK. So if I move the cursor, it's probably going to be a little bit slower, go slower by a little bit. But that doesn't really affect the performance. I can just leave my cursor there and this is what the iterations per second is. Yes?

AUDIENCE: Maybe there's like, a lot of similar programs doing this [INAUDIBLE].

JULIAN SHUN: Yeah. So there is some other factor dominating. Does anyone have an idea what that factor might be?

AUDIENCE: Rendering?

JULIAN SHUN: No. It's not the rendering. I ran the code without showing the graphics, and performance was similar. Yes?

AUDIENCE: Maybe similar to what she said there could be other things using cache [INAUDIBLE].

JULIAN SHUN: Yes. Yeah. So there could be other things using the cache. But that would be true for both of the programs. And I don't actually have anything that's intensive running, except for PowerPoint. I don't think that uses much of my cache. All right.

So let's look at why this is the case. So it turns out that the hardware is actually helping the looping code. So the question is how come the cache oblivious trapezoidal code can have so many fewer cache misses, but the advantage gained over the looping version is so marginal? Turns out that for the looping code, the hardware is actually helping it by doing hardware prefetching.

And hardware prefetching for the looping code is actually pretty good, because the access pattern is very regular. It's just going one row at a time. So the hardware can predict the memory access pattern of the looping code, and therefore, he can bring in the cache lines that the looping code would need before it actually gets to that part of the computation. So prefetching is helping the looping code. And it's not helping the trapezoid code that much because the access pattern is less regular there.

And prefetching does use memory bandwidth. But when you're using just a single core, you have more than enough bandwidth to take advantage of the hardware prefetching capabilities of the machine. But later on, we'll see that when we're running in parallel the memory bandwidth does become more of an issue when you have multiple processors all using the memory. Yeah? Question?

AUDIENCE: Is there a way of touching a cache [INAUDIBLE] or touching a piece of memory before you need it so that you don't need [INAUDIBLE]

JULIAN SHUN: You can do software prefetching. There are instructions to do software prefetching. Hardware prefetching is usually more efficient, but it's like less flexible than the software prefetching. But here we're not actually doing that. We're just taking advantage of hardware prefetching.

AUDIENCE: [INAUDIBLE]?

JULIAN SHUN: Yeah. So we didn't actually try that. It could benefit a little bit if we used a little bit of software prefetching. Although, I think it would benefit the looping code probably as well if we did that. Yes?

AUDIENCE: Is hardware prefetching [INAUDIBLE]?

JULIAN SHUN: Sorry? Sorry?

AUDIENCE: Is hardware prefetching always enabled?

JULIAN SHUN: Yeah. So hardware prefetching is enabled. It's always done by the hardware on the machines that we're using today. This was a pretty surprising result. But we'll actually see the fact of the memory bandwidth later on when we look at the parallel code. Any other questions before I continue?

OK. So let's now look at the interplay between caching and parallelism. So this was a theorem that we proved in the previous lecture. So let's recall what it says. It says let Q sub p be the number of cache misses in a deterministic Cilk computation when run on p processors, each with a private cache, and let s sub p be the number of successful steals during the computation.

In an ideal cache model with a cache size of M and a block size of B, the number of cache misses on p processors equal to the number of cache misses on one processor plus order number of successful steals times M over B. And last time we also said that the number of successful steals is upper bounded by the span of the computation times the number of processors. If you minimize the span of your computation, then you can also minimize the number of successful steals. And then for low span algorithms, the first term is usually going to dominate the Q1 term.

So I'm not going to go over the proof. We did that last time. The moral of the story is that minimizing cache misses in the serial elision essentially minimizes them into parallel execution, assuming that you have a low span algorithm.

So let's see whether our trapezoidal algorithm works in parallel. So does the space cut work in parallel? Recall that the space cut, I'm cutting it with a slope of negative 1 through the center, because it's too wide. So can I execute the two trapezoids in parallel here? No.

The reason is that the right trapezoid depends on the result of the left trapezoid. So I can't execute them at the same time. But there is a way that I can execute trapezoids in parallel. So instead of just doing the cut through the center, I'm actually going to do a V-cut. So now I have three trapezoids.

The two trapezoid in black-- I can actually execute those in parallel, because they're independent. And everything in those two trapezoids just depends on the base of that trapezoid. And after I'm done with the two trapezoids labeled 1, then I can compute the trapezoid label 2. And this is known as a parallel space cut. It produces two black trapezoids as well as a gray trapezoid and two black trapezoids executed in parallel.

And afterwards, the gray trapezoid executes. And this is done recursively as well. Any questions? Yeah? No. OK.

We also have the time cut. Oh, sorry. So if the trapezoid is inverted, then we're going to do this upside down V-cut. And in this case, we're going to execute the middle trapezoid before we execute the two trapezoids on the side.

For the time cut, it turns out that we're just going to use the same cut as before. And we're just going to execute the two trapezoids serially. So we do get a little bit of parallelism here from the parallel space cut. Let's look at how the parallel codes perform now.

So, OK. So this was a serial looping code. Here's the parallel looping code. So we had 1,450 before.

About 3,700 now. So little over twice the speed up. And this is on a four core machine. It's just on my laptop.

AUDIENCE: [INAUDIBLE]?

JULIAN SHUN: Sorry?

AUDIENCE: [INAUDIBLE]?

JULIAN SHUN: Oh yeah, sure. Yeah, it's slowing down a little bit, but not by too much. OK. Let's look at the trapezoidal code now. So as we saw before, the trapezoid code does about 1,840 iterations per second. And we can paralyze this.

So now it's doing about 5,350 iterations per second. So it's getting about a factor of three speed up. I can move it around a little bit more if you want to see it. So serial trapezoid and parallel trapezoid. Is everyone happy? OK.

Because I had to do this in real time, the input size wasn't actually that big. So I ran the experiment offline without the rendering on a much larger input. So this input is a 3,000 by 3,000 grid. And I did this for 1,000 time steps using four processor cores. And my cache size is 8 megabytes.

So last level cache size. So the input size here is much larger than my last level cache size. And here are the times that I got. So the serial looping code took about 129 seconds. And when we did it in parallel, it was about 66 seconds. So it got about a factor of two speed up, which is consistent with what we saw.

For the trapezoidal code, it actually got a better speed up when we increased the input size. So we got about a factor of four speed up. And this is because for larger input size, the cache efficiency plays a much larger role, because the cache is so small compared to our input size. So here we see that the parallel looping code achieves less than half of the potential speed up, even though the parallel looping code has much more potential parallelism than the trapezoidal code. So trapezoidal code only had a little bit of parallelism only for the space cuts, whereas the trapezoidal code is actually getting pretty good speed up. So this is near linear speed up, since I'm using four cores and it's getting 3.96 x speed up.

So what could be going on here? Another thing to look at is to compare the serial trapezoid code with the serial looping code, as well as the parallel trapezoid code with the parallel looping code. So if you look at the serial trapezoid code, you see that it's about twice as fast as the serial looping code.

But the parallel trapezoid or code is about four times faster than the parallel looping code. And the reason here is that the harbor prefetching can't help the parallel looping code that much. Because when you're running in parallel, all of the cores are using memory. And there's a memory bandwidth bottleneck here. And prefetching actually needs to use memory bandwidth.

So in the serial case, we had plenty of memory bandwidth we could use for a prefetching, but in the parallel case, we don't actually have much parallel-- but much memory bandwidth we can use here. So that's why in a parallel case, the trapezoid code gets a better speed up over the parallel looping code, compared to the serial case. And the trapezoid code also gets better speed up because it does things more locally, so it needs to use less-- fewer memory operations.

And there's a scalability bottleneck at the memory interconnect. But because the trapezoidal code is cache oblivious, it does a lot of work in cache, whereas the looping code does more accesses to the main memory. Any questions on this?

So how do we know when we have a parallel speed up bottleneck, how do we know what's causing it? So there are several main things that we should look at. So we should see if our algorithm has insufficient parallelism, whether the scheduling overhead is dominating, whether we have a lack of memory bandwidth, or whether there is contention going on. And contention can refer to either locking or true and false sharing, which we talked about in the last lecture.

So the first two are usually quite easy to diagnose. You can compute the work in the span of your algorithm, and from that you can get the parallelism. You can also use Cilkscale to help you diagnose the first two problems, because Cilkscale can tell you how much parallelism you're getting in your code. And it can also tell you the burden of parallelism which includes the scheduler overhead.

What about for memory bandwidth? How can we diagnose that? So does anyone have any ideas? So I can tell you one way to do it. I can open up my hardware and take out all of my memory chips except for one of them, and run my serial code.

And if it slows down, then that means it was probably memory bandwidth bound when I did it in parallel. But that's a pretty heavy handed way to diagnose this problem. Is there anything we can do was just software? Yes?

AUDIENCE: Can we simulate like Valgrind would do it and count how memory accesses [INAUDIBLE]

JULIAN SHUN: Yeah, so you could use a tool like Valgrind to count the number of memory accesses. Yes?

AUDIENCE: It's like toolset or something where you can make sure that only one processor is being use for this.

JULIAN SHUN: So you can make sure only one processor is being used, but you can't-- but it might be using like, more memory than just the memory from one chip. There's actually a simpler way to do this. The idea is that we'll just run p identical copies of the serial code.

And then they will all executing in parallel. And if the execution slows down, then that means they were probably contending for memory bandwidth. Does that make sense?

One caveat of this is you can only do this if you have enough physical memory, because when you're running p identical copies, you have to use more DRAM than if you just ran one copy. So you have to have enough physical memory. But oftentimes, you can usually isolate some part of the code that you think has a performance bottleneck, and just execute that part of the program with p copies in parallel.

And hopefully that will take less memory. There are also hardware counters you can check if you have root access to your machine that can measure how much memory bandwidth your program is using. But this is a pretty simple way to do this.

So there are ways to diagnose lack of memory bandwidth. Turns out that contention is much harder to diagnose. There are tools that exist that detect lock contention in an execution, but usually they only detect a contention when the contention actually happens, but the contention doesn't have to happen every time you run your code. So these tools don't detect a potential for lock contention. And potential for true and false sharing is even harder to detect, especially false sharing, because if you're using a bunch of variables in your code, you don't know which of those map to the same cache line.

So this is much harder to detect. Usually when you're trying to debug the speed up bottleneck in your code, you would first look at the first three things here. And then once you eliminated those first few things, then you can start looking at whether contention is causing the problem. Any questions?

OK. So I talked about stencil computation. I want to now talk about another problem, sorting. And we want to do this cache efficiently. OK. So let's first analyze the cache complexity of a standard merge sort. So we first need to analyze the complexity of merging, because this is used as a subroutine in merge sort.

And as you recall in merging, we're given two sorted input arrays. And we want to generate an output array that's also sorted containing all the elements from the two input arrays. And the algorithm is going to maintain a pointer to the head of each of our input arrays.

And then it's going to compare the two elements and take the smaller one and put it into the output, and then increment the pointer for that array. And then we keep doing this, taking the smaller of the two elements until the two input arrays become empty, at which point we're done with the algorithm. We have one sorted output array.

OK. So to merge n elements, the time to do this is just theta of n. Here n is the sum of the sizes of my two input arrays. And this is because I'm only doing a constant amount of work for each of my input elements.

What about the number of cache misses? How many cache misses will incur when I'm merging n elements? Yes?

AUDIENCE: [INAUDIBLE].

JULIAN SHUN: Yeah. So I'm going to incur theta of n over B cache misses because my two input arrays are stored contiguously in memory so I can read them at B bytes at a time with just one cache miss. And then my output array is also stored contiguously so I can write things out B bytes at a time with just one cache miss. I might waste to cache line at the beginning and end of each of my three arrays, but that only affects the bound by a constant factor. So it's theta of n over B.

So now let's look at merge sort. So recall that merge sort has two recursive calls to itself on inputs of half the size. And then it doesn't merge at the end to merge the two sorted outputs of its recursive calls. So if you look at how the recursion precedes, its first going to divide the input array into two halves. It's going to divide it into two halves again again, until you get to the base case of just one element, at which point you return.

And then now we start merging pairs of these elements together. So now I have these arrays of size 2 in sorted order. And I merged pairs of those arrays. And I get subarrays of size 4. And then finally, I do this one more time to get my sorted output.

OK. So let's review the work of merge sort. So what's the recurrence for merge sort if we're computing the work? Yes?

AUDIENCE: [INAUDIBLE].

JULIAN SHUN: Yeah. So that's correct. So I have two subproblems of size N over 2. And then I need to do theta n work to do the merge. And this is case two of master theorem.

So I'm computing log base b of a, which is log base 2 of 2. And that's just 1. And that's the same as the exponent in the term that I'm adding in. So since they're the same, I add in an additional log factor. And my overall work is just theta of n log n.

OK. So now I'm going to solve this recurrence again using the recursion tree method. I'm still going to get theta of n log n. But I'm doing this because it's going to be useful when we analyze the cache complexity. So at the top level I have a problem of size n. And I'm going to branch into two problems of size n over 2.

And when I'm done with them, I have to do a merge, which takes theta of n work. And I'm just putting n here. I'm ignoring the constants. And I'm going to branch again. And each one of these is going to do n over 2 work to merge.

And I'm going to get all the way down to my base case after log base 2 of n levels. The top level is doing n work. Second level is also doing n work. And it's going to be the same for all levels down to the leaves.

Leaves is also doing a linear amount of work. So the overall work is just the work per level times the number of levels. So it's just theta of n log n.

OK. So now let's analyze this with caching. OK. So we said earlier that the cache complexity of the merging subroutine is theta of n over B. And here's the recurrence for the cache complexity of merge sort. So my base case here is when n is less than or equal to cM, for some sufficiently small constant c. And this is because at this point, my problem size fits into cache.

And everything else I do for that problem is still going to be in cache. And to load it into cache, the base case, I need to incur theta and over B cache misses. And otherwise, I'm going to have to recursive calls of size n over 2. And then I need to do theta of n over B cache misses to do the merge of my two results.

So here, my base case is larger than what I did for the work. The algorithm actually is still recursing down to a constant size base case. But just for analysis, I'm stopping the recursion when n is less than or equal to CM.

So let's analyze this. So again, I'm going to have the problems of size n at the beginning. And then I'm going to split into two problems of size n over 2. And then I'm going to have to pay n over B cache misses to merge the results together. Similarly for the next level, now I'm paying n over 2B cache misses for each of my two problems here to do the merge.

And I keep going down until I get to a subproblem of size theta of cM. At that point, it's going to fit in cache. And I don't need to recurse anymore in my analysis.

So number of levels of this recursion tree is just log base 2 of n over cM. So I'm basically chopping off the bottom up this recursion tree. The number of levels I had below this is log base 2 of cM. So I'm taking a log base 2 of n and subtracting log base 2 of cM. And that's equivalent to log base 2 of n divided by cM.

The number of leaves I have is n over cM since each leaf is state of cM large. And the number of cache misses I need to incur-- to process a leaf is just theta of m over B, because I just need to load the input for that subproblem into cache. And then everything else fits in cache.

So for the top level, I'm incurring n over B cache misses. The next level, I'm also incurring n over B cache misses. Same with the third level.

And then for the leaves, I'm incurring m over B times n over cM cache misses. The n over cM is the number of leaves I have and theta of m over B is the number of cache misses per leaf. And that also equals theta of n over B.

So overall, I multiply theta of n over B by the number of levels I have. So the number of levels I have is log base 2 of n over cM. And I just got rid of the constant here, since doesn't affect the asymptotic bound. So the number of cache misses I have is theta of n over B times log base 2 of-- or any base for the log of n over M. So any questions on this analysis?

So I am saving a factor of B here in the first term. So that does reduce. That just gives me a better cache complexity than just a work bound. But for the M term, it's actually inside the denominator of the log. And that doesn't actually help me that much.

So let's look at how much we actually save. So one n is much greater than M. Then log base 2 of n over M is approximately equal to log base 2 of n. So the M term doesn't contribute much to the asymptotic costs, and therefore, compared to the work, we're only saving a factor of B. When n is approximately equal to M, then log base 2 of n over m is constant, and we're saving a factor of B log n.

So we save more when the memory size-- or when the problem size is small. But for large enough problem sizes, we're only saving a factor of B. So does anyone think that we can do better than this?

So I've asked this question several times before, and the answer is always the same. Yes? It's a good answer. So we're going to do this using multiway merging. So instead of just merging two sorted subarrays, we're going to merge together R sorted subarrays.

So we're going to have R subarrays, each of size n over R. And these are sorted. And I'm going to merge them together using what's called a tournament tree. So how the tournament tree works is I'm going to compare the heads of each pair of these subarrays and store it in the node of the tournament tree.

And then after I do that, I compare these two nodes. And I get the minimum of those two. Eventually, after I compare all of the elements, I'm just left with the minimum element at the root of the tournament tree. And then I can place that into my output.

So the first time I want to fill this tournament tree, it's going to theta of R work because there are R nodes in my tournament tree. So when I compare these two elements, the smaller one is 6. For these two, the smaller one is 2. And then I compare 2 and 6, take the smaller one. And then on the other side, I have a 7 here.

So I compare 2 and 7. 2 is smaller, so it appears at the root. And then I know that that's going to be my minimum element among the heads of all of the R subarrays that I'm passing it.

So the first time to generate this tournament tree takes theta of R work because I have to fill in R nodes. But once I generated this tournament tree, for all subsequent rounds, I only need to fill in the path from the element that one to the output or right, or to the root of the tournament tree, because those are the only values that would have changed. So now I'm going to fill in this path here. And this only takes theta of log R work to do it, because the height of this tournament tree is log base 2 of R.

So I'm going to fill this in. Now 14 goes here. 6 is a smaller of the two. And then 6 is a smaller of the two again. So my next element is 6.

And I keep doing this until all of the elements for my R subarrays get put into the output. The total work for merging is going to be theta of R for the first round, plus n log R for all the remaining rounds. And that's just equal to theta of n log R, because we're assuming that n is bigger than R here. Any questions on how the multiway merge works? No?

OK. So let's analyze the work of this multiway merge when used inside merge sort. So the recurrence is going to be as follows. So if n is 1, then we just do a constant amount of work. Otherwise, we have R subproblems of size n over R. And then we're paying theta of n log R to do the multiway merge.

So here's the recursion tree. At the top level, we have a problem of size n. And then we're going to split into R subproblems of size n/R. And then we have to pay n log R work to merge the results of the recursive call together. And then we keep doing this.

Turns out that the work at each level sums up to n log base 2 of R. And the number of levels we have here is log base R of n, because each time we're branching by a factor of R. For the leaves, we have n leaves. And we just pay linear work for that, because we don't have to pay for the merge. We're not doing anymore recursive calls.

So the overall work is going to be theta of n log R times log base R of n, plus n for the leaves, but that's just a lower order term. And if you work out the math, some of these terms are going to cancel out, and you just get theta of log n for the work. So the work is the same as the binary merge sort.

Let's now analyze the cash complexity. So let's assume that we have R less than cM over B for a sufficiently small constant C less than or equal to 1. We're going to consider the R way merging of contiguous arrays of total size n.

And if R is less than cM over B, then we can fit the entire tournament tree into cache. And we can also fit one block from each of the R subarrays into cache. And in that case, the total number of cache misses to do the multiway merge is just theta of n over B, because we just have to go over the n elements in our input arrays.

So the recurrence for the R way merge sort is as follows. So if n is less than cM, then it fits in cache. So the number of cache misses we pay is just theta of n over B.

And otherwise, we have R subproblems of size n over R. And that we add theta of n over B to do the merge of the results of the subproblems. Any questions on the recurrence here? Yes?

AUDIENCE: So how do we pick up the value of R? Does it make it cache oblivious?

JULIAN SHUN: Good question. So we didn't pick the value of R. So this is not a cache oblivious algorithm. And we'll see what to choose for R in a couple of slides.

So let's analyze the cache complexity of this algorithm again, using the recursion tree analysis. So at the top level, we're going to have R subproblems of size n over R. And we have to pay n over B cache misses to merge them. And it turns out that at each level, the number of cache misses we have to pay is n over B, if you work out the math. And the number of levels of this recursion tree is going to be log base R of n over cM, because we're going to stop recurring when our problem size is cM.

And on every level of recursion, we're branching by a factor of R. So our leaf size is cM, therefore the number of leaves we have is n over cM. And for each leaf, it's going to take theta of m over B cache misses to load it into cache. And afterwards, we can do everything in cache.

And multiplying the number of leaves by the cost per leaf, we get theta of n over B cache misses. And therefore, the number of cache misses is n over B times the number of levels. Number of levels is log base R of n over M.

So compared to the binary merge sort algorithm, here we actually have a factor of R in the base of the log. Before, the base of the log was just 2. So now the question is what we're going to set R to be. So again, we have a voodoo parameter. This is not a cache oblivious algorithm.

And we said that R has to be less than or equal to cM over B in order for the analysis to work out. So let's just make it as large as possible. Let's just set R equal to theta of M over B. And now we can see what this complexity works out to be. So we have the total cache assumption.

We also have the fact that log base M of n over M is equal to theta of log n over log M. So the cache complexity is theta of n over B times log base M over B of n over M. But if we have the tall cache assumption that log base M over B is the same as log base of M asymptotically.

So that's how we get to the second line. And then you can rearrange some terms and cancel some terms out. And we'll end up with theta of n log n divided by B log M.

So we're saving a factor of theta of b log M compared to the work of the algorithm, whereas for the binary version of merge sort, we were only saving a factor of B for large enough inputs. So here we get another factor of log M in our savings. So as I said, the binary one cache misses is n log n over B, whereas the multiway one is n log n over B log M. And as longest as n is much greater than M, then we're actually saving much more than the binary version. So we're saving a factor of log M in cache misses.

And let's just ignore the constants here and look at what log M can be in practice. So here are some typical cache sizes. The L1 cache size is 32 kilobytes, so that's 2 to the 15th.

And log base 2 of that is 15. So we get a 15x savings. And then for the larger cache sizes, we get even larger saving. So any questions on this?

So the problem with this algorithm is that it's not cache oblivious. We have to tune the value of R for a particular machine. And even when we're running on the same machine, there could be other jobs running that contend for the cache. Turns out that there are several cache oblivious sorting algorithms. The first one that was developed was by paper by Charles Leiserson, and it's called funnel sort.

The idea here is to recursively sort n to the 1/3 groups of n to the 2/3 elements, and then merge the sorted groups with an n to the 1/3 funnel. And this funnel is called a k funnel, more generally, and it merges together k cubed elements in a k sorted list. And the costs for doing this merge is shown here. And if you plug this into the recurrence, you'll get that, the asymptotic number of cache misses is the same as that of the multiway merge sort algorithm while being cache oblivious. And this bound is actually optimal.

So I'm not going to have time to talk about the details of the funnel sort algorithm. But I do have time to just show you a pretty picture of what the funnel looks like. So this is what a k funnel looks like. It's recursively constructed. We have a bunch of square root of k funnels inside. They're all connected to some buffers.

So they feed elements the buffer, and then the buffers feed element to this output square root of k funnel, which becomes the output for the k funnel. So this whole blue thing is the k funnel. And the small green triangles are square root of k funnels. And the number of cache misses incurred by doing this merge is shown here. And it uses the tall cache assumption for analysis.

So a pretty cool algorithm. And we've posted a paper online that describes this algorithm. So if you're interested in the details, I encourage you to read that paper. There are also many other cache oblivious algorithms out there. So there's been hundreds of papers on cache oblivious algorithms.

Here are some of them. Some of these are also described in the paper. In fact, I think all of these are described in the paper we posted.

There are also some cool cache oblivious data structures that have been developed, such as for B-trees, ordered-file maintenance, and priority queues. So it's a lot of literature on cache oblivious algorithms. And there's also a lot of material online if you're interested in learning more.