Computationally efficient three dimensional arrays in C

Question

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.

Answer 1

Have you tried looking in Numerical recipes?

Answer 2

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.

Answer 3

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.

Answer 4

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