# Access optimized data structure for representing integer lattice

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.

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.

• thanks for your reply, the idea with shuffling was nice (+1) The domain I'm dealing with is indeed a square. I understand that it's not reasonable to expect complete elimination of the need for access to far located elements of the array. But imagine we know in advance (by some heuristics) roughly how many times we'll need to access index $(i,j)$ and compute $\Delta^1$ there. So the need for access is weighted. In this case, are you aware of any algorithms that shuffle the matrix that minimize the access attempts given the initial weights ? Or perhaps one can just solve it on their own ?! – Hayk Oct 17 '18 at 20:52
• Every time you access one lattice site, you will also have to access all of its neighbors. There's just no way around it. So the weights will all be equal -- nothing you can do about that. – Wolfgang Bangerth Oct 18 '18 at 3:29
• Suppose I know that my domain is a square of size $512\times512$ (say) and that during the course of my computation I'll need to access the site $(0,0)$ $N_1$ times, the site $(0,5)$ (e.g.) $N_2$ times, ... and the rest of the sites a lot more less. Now if $N_1$ and $N_2$ are huge it makes sense to do away the plain array structure and put those sites and their neighbors closer in memory. In general case the weight for each site shows the N of access required, it is the access order that is unknown. Perhaps I'm missing something... Thanks for your comments, I'll think about it a bit more. – Hayk Oct 18 '18 at 5:20
• I see. Yes, you could try to use this. The Cuthill-McKee algorithm generally puts the starting vertices and its neighbors close to each other. You could start the Cuthill-McKee algorithm with not one lattice site, but several -- exactly the ones that are visited most often. I will say that I expect that not to change very much in practice. – Wolfgang Bangerth Oct 18 '18 at 13:11

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

• thanks for the reply. I'm not sure if the space-filling curve might be of help (unless you don't mean accessing the array along the curve), as I have no way to control the access order. But your suggestion of "block-storage" is interesting, and might be helpful (+1). – Hayk Oct 18 '18 at 5:27
• That space-filling curve is similar to the block storage. Normal, the data in a 2D array is stored in 1D memory line by line. If you store the data in a way the Hilbert curve does, locations that are close together in 2D might on average become closer together in 1D. I'm not sure whether it's any better than the block storage though. – Ark-kun Oct 18 '18 at 7:53