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,