Access optimized data structure for representing integer lattice

Question

Consider the integer lattice in $2d$, namely the set $\mathbb{Z}^2 = \{(x,y): x,y\in \mathbb{Z}\}$, and let $u:\mathbb{Z}^2 \to \mathbb{R} $ be a function defined on some bounded subset of $\mathbb{Z}^2$. The discrete (graph) Laplace of $u$, denote it by $\Delta^1$, at point $(x,y)$ is defined as $$ \tag{1} \Delta^1 u(x,y) = u(x + 1, y) + u(x-1,y) + u(x, y+1) + u(x, y-1) - 4u(x,y), $$ i.e. we sum over all lattice neighbors of $(x,y)$ and subtract the value of the function at the given point multiplied by $4$ (which is the $2d$ since dimension $d=2$).

To compute the discrete Laplace numerically, I represent the lattice (part of it, of course) as a $2d$-array, and thus the values of $u$ are being stored in a $2d$ matrix replicating the grid. Given the nature of $\Delta^1$, when I need to compute it at some index $(i,j)$ of my array, then I need to jump over a row of the array to access the values of neighbors of $(i,j)$. This is obviously not an efficient way of accessing the array, as we force the pointer to repeatedly jump over large chunks of memory locations. Indeed, for certain problems when computing the Laplace of a given function excessively (in some iterative process, for example), even if the size of the matrix is moderately small (say $512\times512$ ) the computation time is enormous compared to the one with the same number of operations but without array access.

The question: is there an efficient way (a data structure) to represent a $2d$ grid to optimize for array access operation of the form described above? Say, specifically for access pattern as in $(1)$ above.

I am aware of spacial and temporal optimization of arrays, as well as that repeatedly accessing the same memory locations the compiler might promote those pointers into registry. However, the access pattern in my case is not very regular.

The question is rather naive of course, but as a $2d$ grid is a ubiquitous object, I'd much appreciate if you can share your insights on efficient data structures concerning it.

Answer 1

In essence, you are asking whether you can enumerate the integer lattice sites within your domain from $1$ to $N$ in such a way that accessing the east/west/north/south neighbors of a location $n$ requires accessing positions $n_e,n_w,n_n,n_s$ so that the distance between index $n$ and these four indices is minimal.

This problem is equivalent to shuffling the rows and columns of a matrix so that its bandwidth becomes minimal. Algorithms such as Cuthill-McKee are used to compute permutations of rows and columns to achieve this.

However, it can be shown that if your domain is, say, a square part of the integer lattice with maximal distance $d$ between nodes in both $x$ and $y$ direction (i.e., you have a total of $N=(d+1)^2$ lattice sites), then there is no ordering that reduces the bandwidth (i.e., $\max\{|n-n_E|,|n-n_W|,|n-n_N|,|n-n_S|\}$) to a value below $d$. In other words, you cannot avoid accessing elements that are stored far away in memory.

Answer 2

You can probably speed things up a little bit by storing the array in 4x4 square subarrays so that each of them fit in cache line (64 bytes = 4x4 32-bit integers). This changes the probability distribution of number of memory accesses from [0,0,1] (mean number of accesses = 3) (for 1,2,or 3 accesses) to [1/4,1/2,1/4] (mean number of accesses = 2).

Maybe some space-filling curve can also help a bit https://en.wikipedia.org/wiki/Hilbert_curve