I am trying to solve numerically a set of partial differential equations in three dimensions. In each of the equations the next value of the unknown in a point depends on the current value of each unknown in the closest points.
To write an efficient code I need to keep the points close in the three dimensions close in the (one-dimensional) memory space, so that each value is called from memory just once.
I was thinking of using octtrees, but I was wondering if someone knows a better method.
Octtrees are the way to go. You subdivide the array into 8 octants:
1 2 3 4 --- 5 6 7 8
And then lay them out in memory in the order 1, 2, 3, 4, 5, 6, 7, 8 as above. You repeat this recursively within each octant until you get down to some base size, probably around 128 bytes or so (this is just a guess -- make sure to profile to determine the optimal cutoff point). This has much, much better cache coherency and locality of reference than the naive layout.
One alternative to the tree-method: Use the Morton-Order to encode your data.
In three dimension it goes like this: Take the coordinate components and interleave each bit two zero bits. Here shown in binary: 11111b becomes 1001001001b
A C-function to do this looks like this (shown for clarity and only for 11 bits):
int morton3 (int a)
{
int result = 0;
int i;
for (i=0; i<11; i++)
{
// check if the i'th bit is set.
int bit = a&(1<<i);
if (bit)
{
// if so set the 3*i'th bit in the result:
result |= 1<<(i*3);
}
}
return result;
}
You can use this function to combine your positions like this:
index = morton3 (position.x) +
morton3 (position.y)*2 +
morton3 (position.z)*4;
This turns your three dimensional index into a one dimensional one. Best part of it: Values that are close in 3D space are close in 1D space as well. If you access values close to each other frequently you will also get a very nice speed-up because the morton-order encoding is optimal in terms of cache locality.
For morton3 you better not use the code above. Use a small table to look up 4 or 8 bits at a time and combine them together.
Hope it helps, Nils
The book Foundations of Multidimensional and Metric Data Structures can help you decide which data structure is fastest for range queries: octrees, kd-trees, R-trees, ... It also describes data layouts for keeping points together in memory.
