Finding Morton-Like Layouts for Multi-Dimensional Arrays Using Evolutionary Algorithms

Stephen Nicholas Swatman

University of Amsterdam
The Netherlands

Ana-Lucia Varbanescu

University of Twente
The Netherlands

Andy D. Pimentel

University of Amsterdam
The Netherlands

Andreas Salzburger

Attila Krasznahorkay

<concept_desc>Software and its engineering Software performance</concept_desc>
<concept_desc>Mathematics of computing Combinatorial optimization</concept_desc>
<concept_desc>Information systems Data layout</concept_desc>

[500]Software and its engineering Software performance
[500]Mathematics of computing Combinatorial optimization
[500]Information systems Data layout

The layout of multi-dimensional data can have a significant impact on the efficacy of hardware caches and, by extension, the performance of applications. Common multi-dimensional layouts include the canonical row-major and column-major layouts as well as the Morton curve layout. In this paper, we describe how the Morton layout can be generalized to a very large family of multi-dimensional data layouts with widely varying performance characteristics. We posit that this design space can be efficiently explored using a combinatorial evolutionary methodology based on genetic algorithms. To this end, we propose a chromosomal representation for such layouts as well as a methodology for estimating the fitness of array layouts using cache simulation. We show that our fitness function correlates to kernel running time in real hardware, and that our evolutionary strategy allows us to find candidates with favorable simulated cache properties in four out of the eight real-world applications under consideration in a small number of generations. Finally, we demonstrate that the array layouts found using our evolutionary method perform well not only in simulated environments but that they can effect significant performance gains—up to a factor ten in extreme cases—in real hardware.

§ INTRODUCTION

Structured multi-dimensional data are ubiquitous in high-performance computing <cit.>: three-dimensional fluid simulations, dense linear algebra operations, and stencil kernels are just a few examples of applications which rely fundamentally on multi-dimensional arrays. In spite of the importance of such applications, however, most modern computer systems have one-dimensional memories: from the perspective of the programmer, memory is nothing more than a very large one-dimensional array of bytes. This discrepancy between application requirements and hardware design requires programmers to carefully consider array layouts: injective functions which translate multi-dimensional indices into one-dimensional memory addresses.

Although array layouts do not impact the functional properties of programs, choosing a suitable layout can significantly impact application performance in modern processors with complex cache hierarchies <cit.>. Exploiting these caches is of critical importance to achieving high performance in all but purely compute-bound applications, but doing so requires locality of access—both temporal and spatial—in memory. Kernels often exhibit locality in multiple dimensions, and a well-chosen array layout maximizes the degree to which this application-level locality is translated to the address-level locality that caches are designed to exploit; as a result, that layout increases the efficacy of hardware caching and—by extension—the performance of an application.

Data in two-dimensions is commonly laid out in row-major order (shown in <ref> for an $8 \times 8$ array) or column-major order (<ref>) which provide good locality of access in a single dimension, but poor locality in all others. Thankfully, the design space for data orderings—in two dimensions or more—extends far beyond these canonical layouts: the Morton layout (<ref>), for example, is a layout based on a space-filling curve which provides balanced locality between multiple dimensions <cit.>. Our work explores a family of data layouts which generalize the Morton order, and allow us to carefully tune the cache behavior in multiple dimensions to match a given application.

The design space of the aforementioned family of data layouts is dauntingly large; indeed, the number of possible layout for arrays at scales applicable to real-world problems is so large that it renders exhaustive search infeasible. In order to find suitable array layouts in tractable amounts of time, we propose to employ genetic algorithms—heuristics known to be able to efficiently find high-quality solutions in large search spaces <cit.>. To this end, we design a chromosomal representation of Morton-like array layouts, as well as a fitness function that uses cache simulation to estimate the performance of individual array layouts. Finally, we evaluate our evolutionary strategy and the array layouts it discovers.

In short, our paper makes the following contributions:

* We characterize the design space given by a generalization of the Morton array layout, and we show that that the size of this design space renders exhaustive search infeasible (<ref>);

* We propose an evolutionary methodology based on genetic algorithms for exploring the aforementioned design space based on the simulated cache-friendliness of layouts (<ref>);

* We design and execute a series of experiments to assess the accuracy of our fitness function, the efficacy of our evolutionary process, and the performance of the discovered array layouts, showing that our method can improve performance up to a factor ten (<ref>).

§ BACKGROUND AND RELATED WORK

In this section, we provide a brief overview of the basic concepts and notations which are essential to the remainder of this paper, and highlight relevant related work.

§.§ Indexing Functions and Canonical Layouts

in 0,...,15
=0 1=5 1=10 1=15 10=1
(p) at (,0) [circle,fill=red,scale=0.25] ;
(p) at (,0) [circle,fill=black,scale=0.25] ;

in 0,...,14

[gray, thick, shorten >=0.2mm, shorten <=0.2mm, ->] (p) -> (p);

[red, thick, shorten >=0.2mm, shorten <=0.2mm, ->] (p0.north) to [out=20, in=160] (p5.north);
[red, thick, shorten >=0.2mm, shorten <=0.2mm, ->] (p5.north) to [out=20, in=160] (p10.north);
[red, thick, shorten >=0.2mm, shorten <=0.2mm, ->] (p10.north) to [out=20, in=160] (p15.north);
(current bounding box.south west) +(-1mm,-1mm) (current bounding box.north east) +(1mm,1mm);

Row-major layout

in 0,...,15
=0 1=3 1=12 1=15 10=1
(p) at (,0) [circle,fill=red,scale=0.25] ;
(p) at (,0) [circle,fill=black,scale=0.25] ;

in 0,...,14

[gray, thick, shorten >=0.2mm, shorten <=0.2mm, ->] (p) -> (p);

[red, thick, shorten >=0.2mm, shorten <=0.2mm, ->] (p0.north) to [out=20, in=160] (p3.north);
[red, thick, shorten >=0.2mm, shorten <=0.2mm, ->] (p3.north) to [out=20, in=160] (p12.north);
[red, thick, shorten >=0.2mm, shorten <=0.2mm, ->] (p12.north) to [out=20, in=160] (p15.north);
(current bounding box.south west) +(-1mm,-1mm) (current bounding box.north east) +(1mm,1mm);

Morton layout

Two-dimensional arrays laid out in memory along the gray arrows. An application accesses the array diagonally along the red arrows. Application locality is shown above, memory locality is shown below.

Dense $n$-dimensional arrays can be imagined as structured grids in which each element is assigned to exactly one point in $\mathbb{N}^n$. In most modern processors, multi-dimensional arrays are a software-level abstraction over the one-dimensional memory of the machine; in order to actually access multi-dimensional data, we need to define a function that converts indices in $n$ dimensions to memory addresses[In reality, address calculations must also consider array offsets (the address of the first element) and scales (the size of each element). We skip over these complications as they are handled transparently by address generation units in modern hardware, and they affect all array layouts in the same manner.]. We refer to the class of such functions as indexing functions, and they are isomorphic to array layouts. In short, an $n$-dimensional indexing functions is an injective (often bijective) function of the following type, where $N_i$ represents the size of the array in the $i$th dimension:

\begin{equation}
f : \bigtimes_{i=0}^{n-1} \llbracket 0,N_i - 1 \rrbracket \to \left\llbracket0,\left(\prod_{i=0}^{n-1} N_i\right) - 1 \right\rrbracket
\end{equation}

In a multi-dimensional grid, we denote the elements along a given axis—that is to say, the sequence of elements for which all indices except one are fixed—as fibers <cit.>. In a two-dimensional case, fibers along the $x$-axis are known as rows, and fibers along the $y$-axis as columns. In order to facilitate the description of arrays in three or more dimensions, we use the term mode-$m$ fibers to describe fibers along the $m$th dimension, such that mode-0 fibers are synonymous with rows, mode-1 fibers refer to columns, and so forth.

The most common group of multi-dimensional indexing functions are the canonical layouts, sometimes known as the lexicographic layouts or, in the two-dimensional case, the row- and column-major layouts. In a canonical layout, one-dimensional array indices are calculated according to the following equation, in which $x_0, \ldots, x_{n-1}$ are components of the $n$-dimensional index, and $N_0, \ldots, N_{n-1}$ represent the size of the array in each dimension:

\begin{equation}
\label{eq:lexicographic}
f(x_0, \ldots, x_{n-1}; N_0, \ldots, N_{n-1}) = \sum_{i=0}^{n-1}\left(\prod_{j=0}^{i-1}N_j\right)x_i
\end{equation}

An important corollary of <ref> is that the mode-0 fibers are contiguous in memory, i.e. $f(x_0 + 1, x_1, \ldots, x_{n-1}) = f(x_0, x_1, \ldots, x_{n-1}) + 1$. It is worth noting that the calculation of addresses in column-major layout—in which the mode-1 fibers are contiguous—is also given by <ref>, with the order of the indices and sizes swapped. The canonical array layouts achieve perfect spatial locality in one dimension: if a kernel accesses memory along mode-$m$ fibers, then a canonical layout where the $m$th dimension is major will provide the optimal translation between locality in the multi-dimensional space to locality in memory. Many real world applications, however, exhibit locality in multiple dimensions; a kernel might, for example, iterate diagonally over an array; an example of this—and the resulting locality in memory—is given in <ref>.

The performance of canonical storage layouts has been studied extensively.
discuss methods for compensating for the weaknesses of canonical layouts using tiling and recursive layouts <cit.>.
Similarly,  propose a variety of strategies that mitigate cache misses in canconical storage layouts for numerical applications <cit.>.
propose a metric for the locality of array layouts <cit.>.
analyze the performance of access patterns in multi-dimensional data in graphics processing units (GPUs) <cit.>.  propose a method for automatically optimizing storage layouts <cit.>.

§.§ Morton Layouts

The Morton order is a notable example of a non-canonical array layout that provides balanced locality in multiple dimensions. It is conceptually simple to understand, efficient to implement in commodity hardware (as we will show in <ref>), and it has been shown to positively affect the efficacy of hardware caches:  show the efficacy of the Morton layout in molecular dynamics applications <cit.>,   describe its benefits in matrix decomposition <cit.>, and  provide an in-depth performance analysis of this array layout in a range of kernels <cit.>.  show the applicability of Morton layouts—as well as other non-canonical layouts—in matrix multiplication <cit.>, and this work is expanded upon in <cit.>. Applications of the Morton order in more than two dimensions have been studied by  <cit.>.  show the applicability of array layouts based on space-filling curves—like the Morton layout—for irregular applications <cit.>. The practical applicability of the Morton layout is further evidenced by the Opie compiler, which employs Morton array layouts natively <cit.>.

The performance benefits of the Morton layout stem from its spatial structure: an example—which justifies why this layout is sometimes known as the Z-order layout—is given in <ref>; note the difference in locality in the address space compared to the canonical layout (<ref>). Beyond low-level memory access patterns, the Morton order layout has also been applied to data movement in distributed systems by  <cit.>, and  have applied the layout to workload distribution in parallel processes <cit.>.  explores a variety of applications of space-filling curves in scientific programs <cit.>.

In the Morton order, multi-dimensional indices can be converted to one-dimensional addresses in a variety of ways. The Moser–de Bruijn sequence <cit.> is commonly used as it allows efficient conversions in two dimensions, but this method requires us to store the Moser–de Bruijn sequence in memory, and accessing this sequence causes additional overhead. Therefore, we prefer a different method based on the interleaving of the (unsigned) binary representation of multi-dimensional indices. As an example, the two-dimensional index $(3, 5)$ can be bijectively mapped into one-dimensional memory by finding the binary expansions of the indices, i.e. $({\color{ForestGreen}011}_2, {\color{Fuchsia}101}_2)$, and interleaving the bits yielding ${\color{Fuchsia}1}{\color{ForestGreen}0}{\color{Fuchsia}0}{\color{ForestGreen}1}{\color{Fuchsia}1}{\color{ForestGreen}1}_2 = 39_{10}$.  This is equivalent to first dilating and shifting the binary expansions of the numbers, and then taking their bitwise disjunction (OR): the first index is dilated yielding ${\color{gray}0}{\color{ForestGreen}0}{\color{gray}0}{\color{ForestGreen}1}{\color{gray}0}{\color{ForestGreen}1}_2$ while the second index is dilated and shifted left yielding ${\color{Fuchsia}1}{\color{gray}0}{\color{Fuchsia}0}{\color{gray}0}{\color{Fuchsia}1}{\color{gray}0}_2$. Taking the bitwise disjunction of these numbers yields the same address as using the interleaving strategy. The computation of Morton indices through bit manipulation can be extended to an arbitrary number of dimensions; the three-dimensional index $(3, 5, 4)$ expands to $({\color{ForestGreen}011}_2, {\color{Fuchsia}101}_2, {\color{RedOrange}100}_2)$, and the resulting memory address is ${\color{RedOrange}1}{\color{Fuchsia}1}{\color{ForestGreen}0}{\color{RedOrange}0}{\color{Fuchsia}0}{\color{ForestGreen}1}{\color{RedOrange}0}{\color{Fuchsia}1}{\color{ForestGreen}1}_2 = 395_{10}$. Note that the relative significance of bits in each of the input indices is preserved in the output address.  present the idea that the Morton layout can be generalized by allowing arbitrary bit-interleaving orders <cit.>, which is foundational to our work. This idea is further expanded on by  <cit.>.

§.§ Genetic Algorithms

Genetic algorithms are a class of heuristics introduced by  which are designed to solve optimization and search problems by emulating the process of evolution as it happens in the natural world <cit.>. In genetic algorithms, generations of individuals, i.e. sets of candidate solutions, iteratively explore a design space through genetic operators. In particular, crossover operators model the combination of the genetic material of two or more individuals <cit.>, and mutation operators model random changes to the gene pool <cit.>. In genetic algorithms, individuals are removed from the population based on their fitness, i.e. the quality of the solution they represent to the problem posed <cit.>. Genetic algorithms have seen successful application in an extremely broad range of fields, ranging from drug discovery <cit.> to music composition <cit.>. Genetic algorithms have also proven useful for design space exploration in computer systems;  shows that they can be applied in the design of embedded systems <cit.>.  show that a broader class of evolutionary approaches can be used in the design of neural networks <cit.>. The optimization problem we consider in this paper is combinatorial in nature, and the application of genetic algorithms to such problems has also been extensively studied and proven across a variety of domains <cit.>

§ GENERALIZED MORTON LAYOUTS

All 20 layouts for $8\times 8$ arrays generated by the family of indexing schemes described in <ref>. Note that <ref> corresponds to a row-major layout, while <ref> corresponds to a column-major layout.

The Morton layout functions by interleaving the bits of the input indices in a fixed pattern: bits are drawn from each of the inputs in a round-robin manner. In this section, we generalize this idea, allowing bits to be interleaved in arbitrary order. This gives rise to more specialized layouts with different structure and, as a result, different extra-functional properties <cit.>. <ref> shows all 20 layouts that are given by different bit interleaving orders for an $8 \times 8$ array. As with the standard Morton layout, the generalized Morton layout can be applied to any number of dimensions. As an example, the following three-dimensional layout selects two bits from the second index, one bit from the third index, then two bits from the first index, etc.:

\begin{equation}
\arraycolsep=1pt
\label{ex:permute:4}
f({\color{ForestGreen}011}_2, {\color{Fuchsia}101}_2, {\color{RedOrange}100}_2) = \begin{array}{r r}&{\color{gray}0}{\color{ForestGreen}0}{\color{gray}00}{\color{ForestGreen}1}{\color{ForestGreen}1}{\color{gray}000}_2\\\vee&{\color{gray}000}{\color{Fuchsia}1}{\color{gray}000}{\color{Fuchsia}01}_2\\&{\color{RedOrange}1}{\color{gray}0}{\color{RedOrange}0}{\color{gray}000}{\color{RedOrange}0}{\color{gray}00}_2\\\hline&{\color{RedOrange}1}{\color{ForestGreen}0}{\color{RedOrange}0}{\color{Fuchsia}1}{\color{ForestGreen}11}{\color{RedOrange}0}{\color{Fuchsia}01}_2\end{array} = 313_{10}
\end{equation}

Our goal in this paper is to find Morton-like layouts, i.e. bit-interleaving patterns, that improve application performance through an increase in cache efficacy. In this section, we will show that the design space for such layouts is very large, motivating the use of genetic algorithms. This necessitates a chromosomal representation of layouts, which we also present in this section. In addition, we describe how the canonical layouts can be described using the same representation, and we delve into practical considerations such as the computational cost of computing indices and support for same-instruction multiple-data (SIMD) processing.

§.§ Enumerating Layouts

We can characterize Morton-like layouts by the bit scattering pattern applied to each of the inputs (e.g., for <ref>, the first index is scattered to the fourth, fifth, and eighth bits). However, such a characterization is unsound in the sense that is allows us to describe invalid layouts: if two bits from any of the input indices are mapped onto the same bit in the output, the bitwise disjunction becomes an information-destroying operation and the layout becomes non-injective—that is, it would cause multiple multi-dimensional indices to map onto the same location in memory, making the layout unusable.

We can instead characterize layouts in a manner that is both complete and sound by enumerating the source of each bit in the output index. In the remainder of this work we shall denote array layouts using sequences of the form $[i_0, \ldots, i_{n-1}]$, indicating the source indices in order of increasing bit significance: the least significant bit in the output index is drawn from the $i_0$th input index, the second-least significant bit is drawn from the $i_1$th input, and the most significant bit is drawn from the $i_{n-1}$th input. Note that each input bit must be used once and only once: whenever a bit is to be drawn from a given input index, we implicitly use the least significant bit for that input which has not yet been consumed. For the layout shown in <ref>, the two least significant bits are drawn from the second input, the third-least significant bit is drawn from the third input, and so forth: the resulting array layout is denoted using the sequence $[1,1,2,0,0,1,2,0,2]$.

The aforementioned characterization of multi-dimensional layouts gives rise to families of layouts. The family of layouts over $n$ inputs, where each input has $b_0,\ldots,b_{n-1}$ bits, is isomorphic to the set of permutations of the multiset $S = \{0:b_0,\ldots,n-1:b_{n-1}\}$. We denote this set of permutations as $\mathfrak{S}(S)$. For convenience, we obviate the intermediate multiset such that $\mathfrak{S}'(b_0,\ldots,b_{n-1}) = \mathfrak{S}(\{0:b_0,\ldots,n-1:b_{n-1}\})$. We can then determine the total number of possible layouts as the number of multiset permutations of $\mathfrak{S}'$ <cit.>:

\begin{equation}
\label{eq:permutation_set_size}
|\mathfrak{S}'(b_0, \ldots, b_{n-1})| = \binom{\sum_{i = 0}^{n-1} b_i}{b_0,\ldots,b_{n-1}} = \frac{\big(\sum_{i = 0}^{n-1} b_i \big)!}{\prod_{i = 0}^{n-1} (b_i!)}
\end{equation}

§.§ Including Canonical Layouts

It is worth noting that canonical layouts over arrays for which the size in each dimension is a power of two are, in fact, members of the family of Morton-like layouts. In order to sketch an informal argument for this, we recall that the indexing function for an $n$-dimensional canonical layout given array sizes $N_0, \ldots, N_{n-1}$ is defined as in <ref>. If we assume that all sizes are powers of two, then the product of these sizes is guaranteed to be itself a power of two. Because multiplication by powers of two can be interpreted as a left-ward shift, the canonical layouts shift each input index $x_0, \ldots, x_n$ to a specific location in the binary expansion of the output index. Furthermore, because we assume $\forall i : x_i < N_i$, each bit in the output is determined by exactly one of the input indices; this allows us to interpret the summation as a series of bit-wise disjunctions, exactly like the definition of our Morton-like layouts. In general, a mode-0-major canonical layout of a $2^{b_0} \times \ldots \times 2^{b_{n-1}}$ array can be characterized—in the the scheme defined in <ref>—by contiguous subsequences of bits, each drawn from the same index, i.e. a sequence of the following form:

\begin{equation}
\label{eq:cacnonical_repr}
[\underbrace{0,\ldots,0}_{b_0\text{ times}}, \underbrace{1,\ldots,1}_{b_1\text{ times}}, \ldots, \underbrace{n-1,\ldots,n-1}_{b_{n-1}\text{ times}}]
\end{equation}

Canonical layouts with different major axes can be constructed by changing the order of the contiguous subsequences. The fact that the canonical layouts are members of the Morton-like family of array layouts allows us to evaluate the performance of these layouts in the exact same framework as the rest of the Morton-like layouts, and we will exploit this in <ref>.

§.§ Hardware-Accelerated Indexing

It is tempting to extend the aforementioned ideas to even more exotic indexing functions, like the Hilbert array layout <cit.>. The computational cost of many such functions renders them impractical, however: if the cost of computing memory addresses is too large, any performance gained by improving the cache-friendliness of a program will be negated. The Morton-like layouts we consider in this work allow efficient index calculations on modern commodity hardware, which we demonstrate in this section.

Under canonical array layouts, indices are calculated either iteratively through repeated addition and multiplication, or in parallel through parallel multiplication followed by reduction through addition. In $n$-dimensional cases both approaches require $n-1$ additions and $n-1$ multiplications, operations which can be efficiently performed on virtually all processors. Specifically, the Intel Haswell and AMD Zen 3 microarchitectures—on which we focus in this work—can perform 64-bit register addition () with a latency 1 cycle and a reciprocal throughput of 0.25 cycles, while they can execute multiplication () with a latency of 3 cycles and a reciprocal throughput of 1 cycle <cit.>.

Our bit-interleaving array layouts rely, in $n$-dimensional cases, on $n-1$ bitwise disjunctions and $n$ bit-scatter operations. Such disjunctions () can be performed with a latency of 1 cycle and a reciprocal throughput of 0.25 cycles—the same as the  instruction—on both of the aforementioned microarchitectures. We perform the bit-scattering operation using the parallel bit deposition () instruction, which is included in the BMI2 extension to the x86-64 instruction set <cit.>. The Intel Haswell and AMD Zen 3 microarchitectures both perform bit deposition with a latency of 3 cycles and a reciprocal throughput of 1 cycle, identical to the  instruction. It follows that Morton-like indexing requires—in theory—only a single additional instruction over canonical index calculation.

The hardware extension required to perform bit deposition is widely supported: BMI2 has been included in Intel processors starting with the Haswell microarchitecture (2013) <cit.>, and in AMD processors starting with the Excavator microarchitecture (2015), albeit in a limited fashion; AMD processors gained full hardware support for these instructions starting with the Zen 3 microarchitecture (2020) <cit.>[Pre-Zen 3 processors supported parts of the BMI2 instruction set—the  and  instructions in particular—through emulation in microcode rather than in hardware, making them very slow.].

Throughput of a kernel calculating array indices using canonical layouts as well as Morton-like layouts on the Intel Haswell microarchitecture as given by OSACA.

In order to further evaluate the competitiveness of Morton-like layouts compared to canonical layouts, we analyze implementations of both indexing schemes over a range of dimensionalities as compiled by gcc 12.3 and clang 15.0 using OSACA 0.5.2 <cit.>. All code was compiled using the  optimization flag. The results of this analysis are shown in <ref>. Over the range of dimensionalities considered, the canonical layouts are consistently faster, i.e. require fewer cycles to compute, than the Morton-like layouts. However, the difference in performance—approximately one cycle—is relatively small and overshadowed by the number of cycled saved due to a reduction in cache misses. Furthermore, we focus primarily on memory-bound applications, in which a small increase in index calculation time is unlikely to affect performance. We conclude, therefore, that Morton-like layouts are competitive with canonical layouts strictly in terms of address computation costs.

§.§ Support for SIMD

An important consideration in the design of array layouts is the ability to vectorize kernels through single-instruction multiple-data (SIMD) operations. Canonical layouts guarantee the contiguity of fibers in the array, which facilitates the (automated) vectorization (e.g. the application of SIMD) of many operations, and this benefit is lost when applying the array layouts discussed in this paper. However, we posit that there remains ample opportunity to accelerate computation on Morton-like arrays using SIMD, and we argue this by distinguishing two classes of computation patterns.

The first class consists of unstructured patterns in which data is operated on element-wise without spatial context, i.e. without consideration of nearby elements; a prominent example of such an operation is matrix addition. In such applications, SIMD can be trivially applied to the underlying one-dimensional memory, regardless of the layout of the data: since elements can be added point-wise in any order, doing so in the order in which the data is laid out in memory is both feasible and enables SIMD.

The second class of problems consists of structured patterns in which operations must be performed in a specific order. A prime example of such an operation is matrix multiplication where the inner product of fibers must be computed. In such cases, it is imperative that fibers can be accessed in contiguous blocks. The size of these blocks depends on the vectorization technology used as well as the size of the data type: in the x86 instruction set, SSE vectorisation requires two consecutive double-precision numbers or four consecutive single-precision numbers <cit.>; the much wider ARM SVE instruction set extension <cit.> may require up to thirty-two consecutive double-precision numbers or sixty-four single-precision numbers.

In order to facilitate vectorization for structured patterns of computation, we can impose certain constraints on the array layouts we consider. Indeed, if the $n$ least-significant bits of an interleaving pattern are all drawn from the $m$th input index, then the layout guarantees that the mode-$m$ fibers in the array are contiguous in blocks of $2^n$ elements. This requirement can be incorporated into the selection of array layouts; for example, we can enable efficient AVX2 vectorisation (with a vector width of 256 bits) using single-precision (32-bit) floating point numbers by ensuring that the three least significant bits in an array layout are drawn from the same source. In other words, we can easily constrain our search space to include only array layouts with properties that favor vectorization, and we believe that doing so will enable SIMD-accelerated computation arrays laid out in Morton-like orders.

§ EXPLORATION THROUGH EVOLUTION

The canonical set of indexing bijections for laying out multi-dimensional memory is small: for two-dimensional data, there are two possible layouts, and the performance of these layouts can be evaluated using exhaustive benchmarks <cit.>. Exhaustively exploring the family of indexing function outlined in <ref>, however, is impractical owing to the sheer number of permissible permutations. Importantly, the number of canonical layouts increases only with the number of dimensions, while the number of Morton-like layouts increases with both the number of dimensions and the size of the array in each of those dimensions. By <ref>, a small $4 \times 4$ array (indexed by two bits in each dimension) can be laid out in $\nicefrac{(2+2)!}{2!2!} = 6$ ways. A larger array of size $\num{4096} \times \num{4096}$ (twelve bits in each dimension) can be laid out in $\nicefrac{(12+12)!}{12!12!} = \num{2704156}$ ways. A three-dimensional array of size $256 \times 256 \times 256$ has the same number of elements as the aforementioned $\num{4096} \times \num{4096}$ array, but permits $\nicefrac{(8+8+8)!}{8!8!8!} = \num{9465511770}$ permutations. As these examples indicate, the number of possible permutations quickly scales beyond what can be feasibly explored through exhaustive search; in order to tackle the explosive growth in the design space for Morton-like layouts, we propose the use of genetic algorithms (<ref>).

§.§ Genetic Algorithm Configuration

In this work, we employ a relatively simple $(\lambda, \mu)$-ES genetic algorithm <cit.>. The chromosomal representations of array layouts is identical to the characterization given in <ref>, and this gives rise to a combinatorial optimization problem. We facilitate the recombination of array layouts into novel layouts using the ordered crossover (OX) operator <cit.>, and we employ inversion-based mutation <cit.>. Our approach differs from classical genetic algorithms in only one significant way: our initial population is not chosen randomly from the solution space. Instead, the initial populations for our evolutionary experiments always consist of two individuals, depicting two canonical layouts for a given array size, as described in <ref>. We choose to do this to ensure that our initial populations are unbiased and deterministic, allowing us to more easily assess the efficacy of our genetic strategy.

§.§ Fitness Function Design

There are two general strategies for evaluating the performance, i.e. fitness, of a given array layout under a given cache hierarchy and access pattern: measurement and simulation. In order to assess fitness through measurement, we execute a program on actual hardware and measuring the running time of the process. Although such a fitness function is conceptually simple, it suffers from two primary flaws:

* measurements are noisy and may suffer from run-to-run variance, which may hinder the performance of genetic algorithms <cit.>—in particular, our genetic algorithm is vulnerable to noise stemming from cache pollution effects; and

* measurements require access to the target hardware, which may be inconvenient or even impossible—for example, in hardware-software co-design scenarios, where the hardware under consideration does not (yet) exist.

For these reasons, we choose not to base our fitness function on measurements.

Instead, we employ simulation for which we need a simulator that can accurately compare the cache performance for different access-patterns on the same cache hierarchy.
For this, we selected pycachesim, a component of the Kerncraft toolkit <cit.>
. We use pycachesim by simulating an access pattern such as matrix multiplication and registering the relevant trace of load and store operations. After all accesses have been recorded, we force a write-back of the caches and collect the number of hits and misses in each cache level. We combine the number of hits in every cache level as well as in main memory with the latency of retrieving data from each of these levels to compute the total number of cycles spent retrieving data from the cache hierarchy. Given an array layout $I$, an access pattern $A$ and a simulated cache hierarchy $H$, we calculate the total number of cycles using the following equation, in which $\mathrm{L}i_\mathrm{hit}$, $\mathrm{L}i_\mathrm{miss}$, and $\mathrm{L}i_\mathrm{lat}$ represent the number of hits, the number of misses, and the latency of the $i$th cache level

, and $M$ represents main memory:

\begin{equation}
\label{eq:fitness-pre}
C(I; A, H) = \mathrm{M}_\mathrm{hit}(I; A, H)\mathrm{M}_\mathrm{lat}(H) + \sum_i \mathrm{L}i_\mathrm{hit}(I; A, H)\mathrm{L}i_\mathrm{lat}(H)
\end{equation}

From this, we compute an approximation of the number of accesses performed per cycle, giving rise to a higher-is-better fitness function defined as follows:

\begin{equation}
\label{eq:fitness}
F(I; A, H) = \frac{\mathrm{L1}_\mathrm{hit}(I; A, H) + \mathrm{L1}_\mathrm{miss}(I; A, H)}{\mathrm{L1}_\mathrm{lat}(H)\cdot C(I; A, H)}
\end{equation}

Intuitively, the numerator in <ref> counts the total number of memory accesses, as all accesses either hit or miss in L1.
The denominator, then, estimates the total number of cycles spent retrieving data from the various cache levels. The denominator is multiplied by a normalizing factor equal to the latency of the L1 cache; it follows from <ref> that the achievable performance is softly bound by the reciprocal of the L1 access latency. Indeed, this performance is achieved if and only if all accesses hit the L1 cache. Normalizing the fitness function using the L1 cache latency improves our ability to compare results between different cache hierarchies.

§ EVALUATION

We evaluate the efficacy of the methods hitherto discussed by demonstrating that

* our fitness function is well-chosen, i.e. that is correlates with performance measurements in real hardware; that

* our evolutionary process is capable of finding novel array layouts with favorable cache properties; and that

* the layouts which are found by our evolutionary process actually lead to relevant performance gains in real hardware.

Our validation is based on eight distinct access patterns and two processors with distinct cache hierarchies.

§.§ Experimental Setup

template <concepts::array<2> M>
void mm_ijk(const M     A, const M     B, M     C)
const auto m = C.get_size();
for (std::size_t i = 0; i < m; ++i)
for (std::size_t j = 0; j < m; ++j)
typename M::value_type acc = 0.;
for (std::size_t k = 0; k < m; ++k)
acc += A.load(i, k) * B.load(k, j);
C.store(acc, i, j);

Example of how an access pattern (MMijk) is described in C++. Metaprogramming allows the same source to be used for both simulation and execution on real hardware.

We consider a set of eight access patterns loosely based on the selection of algorithms used by  <cit.>. The access patterns were picked to represent common real-world applications (dense linear algebra and fluid dynamics), to represent both two-dimensional and three-dimensional applications, and to differ in critical properties such as memory size and number of loads and stores. A description of the access patterns we consider in this paper is given in <ref>.

All our access patterns are described using C++ code—see the example in <ref>—which ensures high performance as opposed to the Python code used for our evolutionary processes; the interaction between the C++ and Python components of our project is managed using pybind11 <cit.>. We use template meta-programming to generalize our access patterns in such a way that a single definition can be used for both simulation and benchmarking without loss of performance due to run-time polymorphism; this eliminates any possible discrepancies between the code used for simulation and the code used for measurement.

Overview of the access patterns used for evaluation, including the use of memory and the number of loads and stores.

Access pattern     Description     Mem. size     Loads     Stores
$\textsc{MMijk}(m; s)$     Multiplication of two $2^m \times 2^m$ matrices, both of $s$-byte real numbers.     $3\cdot s \cdot 2^{2m}\,\si{\byte}$     $2 \cdot 2^{3m}$     $2^{2m}$

$\textsc{MMTijk}(m, n; s)$     Multiplication of a $2^m \times 2^n$ matrix by a transposed $2^m \times 2^n$ matrix.

$\textsc{MMikj}(m; s)$     Same as $\textsc{MMijk}(m;s)$ with the order of the inner loops switched.             $2^{3m}$

$\textsc{MMTikj}(m, n; s)$     Same as $\textsc{MMTijk}(m, n; s)$ with the order of the inner loops switched.

$\textsc{Jacobi2D}(m, n; s)$     Four-point stencil kernel over a $2^m \times 2^n$ array of $s$-byte real numbers.     $2 \cdot s\cdot 2^{m+n}\,\si{\byte}$     $\sim4 \cdot 2^{m + n}$     $2^{m + n}$

$\textsc{Cholesky}(m; s)$     Cholesky–Banachiewicz decomposition of a $2^m \times 2^m$ matrix.     $2 \cdot s \cdot 2^{2m}\,\si{\byte}$     $2 \cdot 2^{2m}$     $\sim\frac{1}{2} \cdot 2^{2m}$

$\textsc{Crout}(m; s)$     Crout decomposition of a $2^m \times 2^m$ matrix of $s$-byte real numbers.         $\frac{7}{2} \cdot 2^{2m}$     $2^{2m}$

$\textsc{Himeno}(m, n, p; s)$     Nineteen-point Himeno stencil <cit.> over $2^m \times 2^n \times 2^p$ arrays.     $12\cdot s \cdot 2^{m + n + p}\,\si{\byte}$     $24 \cdot 2^{m+n+p}$     $2^{m+n+p}$

sets: 64
ways: 8
line: 64
replacement: LRU
write_back: true
store_to: L2
load_from: L2
latency: 4
sets: 512
ways: 8
line: 64
replacement: LRU
write_back: true
store_to: L3
load_from: L3
victim_to: L3
latency: 12
sets: 25600
ways: 16
line: 64
replacement: LRU
write_back: true
latency: 36
first: L1
last: L3
latency: 200

Intel Xeon E5-2660 v3

sets: 64
ways: 8
line: 64
replacement: LRU
write_back: true
store_to: L2
load_from: L2
latency: 7
sets: 1024
ways: 8
line: 64
replacement: LRU
write_back: true
store_to: L3
load_from: L3
victim_to: L3
latency: 12
sets: 32768
ways: 16
line: 64
replacement: LRU
write_back: true
latency: 46
first: L1
last: L3
latency: 200

AMD EPYC 7413

Two examples of cache specifications for different CPU models. Note that these configurations are approximations of the true cache hierarchies.

We conduct our experiments on two different CPUs: the Intel Xeon E5-2660 v3 <cit.> based on the Haswell microarchitecture <cit.>, and the AMD EPYC 7413 <cit.> based on the Zen 3 microarchitecture <cit.>. When we perform experiments on non-simulated Haswell processors we use the the DAS-6 cluster <cit.>, whereas we use a machine located at CERN for experiments on Zen 3 processors. When we perform experiments based on simulation, we use the the DAS-6 cluster <cit.> and configure our cache simulator according the cache configurations shown in <ref> for the Haswell processor, and <ref> for the Zen 3 processor. Note that the cache configurations are based on the accessibility of caches from a single core. This is especially relevant for the L3 cache on the Zen 3 chip, which is shared across groups of cores rather than the entire CPU: in the case of the AMD EPYC 7413, the CPU comes equipped with 128 of L3 cache, but only 32 is accessible from any single core <cit.>. We simplify the cache replacement policies of the actual hardware by assuming LRU caches; in reality, the Haswell caches employ policies consistent with tree-PLRU for the L1 and L2 caches <cit.>, while the L3 cache is consistent with a set-dueling-controlled adaptive insertion policy <cit.>. Cache sizes were gathered from specification documents <cit.>, while cache latencies were obtained optimistically from sources on the fastest load-to-use latencies <cit.>. The Zen 3 L1 cache has a fastest load-to-use latency of four cycles for integers and seven cycles for floating point values <cit.>—we use the latter in our simulations. Finally, we assume a constant 200 cycle access latency for main memory in both systems.

§.§ Fitness Function Validation

The fitness function we use in our evolutionary process (<ref>) is based on simulation results because
simulation yields significant benefits over empirical measurements, primarily in terms of determinism and in the ability to simulate future hardware. However, this strategy is not without risk: the simulation we perform is based on a non-cycle-accurate simulator, uses simplified cache hierarchies, and ignores computation entirely.
Consequently, we must evaluate the usefulness of our fitness function by establishing its correlation with execution time in real hardware.

Ideally, the running time of a kernel using a given array layout would correlate inversely linearly with our fitness function, therefore ensuring two important properties.
Firstly and most importantly, it guarantees that running time decreases monotonically with the value of the fitness function, such that an array layout with a higher fitness value is guaranteed to run more quickly; this allows us to establish a ranking of layouts and enables us to reliably select the best-performing array layout.
Secondly, linear correlation guarantees proportionality between fitness and running time, which facilitates the weighted selection of individuals.

To evaluate the degree to which the aforementioned criteria are met, we randomly select one hundred array layouts for each of the eight access patterns given in <ref>. We then evaluate the simulated fitness and measure the running time in real hardware of each pair of array layout and access pattern. The fitness functions of the pairs are calculated in parallel, as they are designed to be deterministic and impervious to cache pollution or resource contention. The empirical benchmarks are performed sequentially, ensuring that the benchmark is the sole user of the processor caches. All measurements are repeated ten times, and we report the mean and standard deviation of the running time.

Scatter plot of the fitness and measured running time on an Intel Xeon E5-2660 v3 CPU and AMD EPYC 7413 for randomly chosen array layouts.

The results of this experiment are shown in <ref>. The coefficient of variation of the measurements never exceeded a value of $c_{\mathrm{v}} = 0.0801$. Accordingly, we have opted to omit error bars from the figure. Upon visual inspection, it is clear that the correlation between our fitness function and running time is not linear, although the two do appear correlated. We confirm our suspicions of correlation by computing Pearson's coefficient of correlation ($\rho_p$) and Spearman's coefficient of rank correlation ($\rho_s$); the resulting statistics are given in <ref>. We observe that our fitness function and running time correlate moderately to strongly with running time for the Intel Xeon E5 2660 v3 processor, although the correlation is weaker for the AMD EPYC 7413 processor. Although it is clear that there is space for the fitness function to be improved, we believe that it correlates sufficiently with running time to enable its use in genetic algorithms.

Pearson's coefficient of correlation ($\rho_p$) and Spearman's coefficient of rank correlation ($\rho_s$) between our simulation-based fitness function and true running time.

2cIntel E5-2660 v3     2cAMD EPYC 7413
Access pattern     $\rho_{p}$     $\rho_{s}$     $\rho_{p}$     $\rho_{s}$
$\textsc{MMijk}(9; 4)$     -0.672     -0.480     -0.648     -0.489

$\textsc{MMTijk}(9, 9; 4)$     -0.810     -0.896     -0.863     -0.823

$\textsc{MMikj}(9; 4)$     -0.845     -0.815     -0.800     -0.838

$\textsc{MMTikj}(9, 9; 4)$     -0.777     -0.744     -0.291     -0.405

$\textsc{Jacobi2D}(13, 13; 4)$     -0.760     -0.769     -0.390     -0.428

$\textsc{Cholesky}(10; 4)$     -0.827     -0.953     -0.725     -0.892

$\textsc{Crout}(9; 4)$     -0.846     -0.663     -0.213     -0.704

$\textsc{Himeno}(8, 7, 7; 4)$     -0.607     -0.475     -0.561     -0.496

§.§ Genetic Algorithm Performance

To evaluate our evolutionary process (<ref>) as a whole,

we intend to verify that it can, indeed, find Morton-like array layouts that have a higher simulated fitness than the canonical layouts. To this end, we perform the evolutionary process for each combination of our two simulated processors and eight access patterns, giving rise to a total of sixteen experiments. For all of these experiments, we configure our genetic algorithm to use $\mu = 20$, $\lambda = 20$, and a mutation rate of $25\%$. We simulate a total of 20 generations in each case.

<ref> shows a violin plot of the fitness distribution of all individuals considered during the evolutionary process. <ref> shows the evolution of population fitness over the course of our experiments. We notice that for the MMTijk, MMikj, Jacobi2D, and Himeno access patterns, our method does not manage to discover any layouts with higher fitness than the initial population of canonical layouts. In the experiment on the MMijk access pattern, we discover layouts with a fitness $149.8\%$ higher than the canonical layouts on the Intel Xeon E5-2660 v3 processor, and we improve on the fitness of canonical layouts by $187.5\%$ for the AMD EPYC 7413. We also find layouts with improved fitness for the MMTikj ($109.6\%$ and $141.1\%$ for the Intel and AMD processors, respectively), Cholesky ($26.4\%$ and $36.8\%$), and Crout ($545.9\%$ and $541.1\%$) access patterns. It is notable that we are able to find layouts with high fitness in few generations.

Distribution of the fitness for all individuals in eight evolution experiments.

Range of fitness values across eight experiments for the Intel Xeon E5-2660 v3 (blue) and AMD EPYC 7413 (red). Mean fitness values are given by the dashed lines.

§.§ Real-World Performance

In order to evaluate whether the layouts identified by our evolutionary algorithms as superior to canonical layouts are indeed better, we evaluate them on real hardware. We collect the fittest individual from each of the successful evolution experiments, i.e. experiments in which our method improved upon canonical layouts, and evaluate the performance of those layouts compared to the canonical layouts on real hardware. Given that our genetic algorithm discovered superior layouts for four access patters—MMijk, MMTikj, Cholesky, and Crout–and that we evaluate a discovered layout and two canonical layouts for each access pattern, this gives rise to twenty-four experiments. We repeat each experiment ten times to compensate for run-to-run variance.

The results of our experiments are shown in <ref>; they show that some access patterns—the Cholesky pattern in particular—benefit very little from our method, with speed-ups ranging from small on the Haswell processor to insignificant on the Zen 3 processor. The matrix multiplication access patterns benefit more, and performance for these access patterns is improved significantly. The Crout access pattern stands out as achieving very large speedup—up to a factor ten—from our method. It is worth noting that, in most cases, the Zen 3 processor benefits more from our evolutionary methodology than the Haswell processor; we do not currently have a satisfactory explanation for this behavior.

Comparison of running time between the best-performing canonical layout and the best-performing layout found by our evolutionary process for four access patterns.

2lAccess pattern     Best can.     Best evo.     Speedup
2lIntel Xeon E5-2660 v3

$\textsc{MMijk}(11; 4)$     17.84      10.94      $63.1\%$

$\textsc{MMTikj}(11, 11; 4)$     18.13      13.96      $29.9\%$

$\textsc{Cholesky}(12; 4)$     11.84      11.43      $3.6\%$

$\textsc{Crout}(12; 4)$     158.54      43.72      $262.6\%$

2lAMD EPYC 7413

$\textsc{MMijk}(11; 4)$     37.71      9.58      $293.8\%$

$\textsc{MMTikj}(11, 11; 4)$     32.35      15.21      $112.6\%$

$\textsc{Cholesky}(12; 4)$     9.72      9.55      $1.0\%$

$\textsc{Crout}(12; 4)$     232.84      21.03      $1007.0\%$

It is important to note that we do not claim to have discovered a novel way of performing matrix multiplication or matrix decomposition that outperforms existing implementations. Indeed, our experiments are based on relatively naive implementations of these algorithms; high-performance implementations of matrix multiplication commonly rely on tiling to significantly improve the cache behavior of the application <cit.>, and the performance of tiled matrix multiplication surpasses what we achieve in this paper. The purpose of the methodology described in this paper, rather, is to provide an alternative way of improving the cache behavior of an application in a manner which is fully agnostic of the application: unlike tiling and other application-specific optimizations, our methodology of altering the array layouts can be applied to any multi-dimensional problem without the need for application-specific knowledge. In addition, our approach requires few code changes, making it easy to implement.

§ LIMITATIONS AND THREATS TO VALIDITY

Throughout this work, we evaluate cache efficacy through a simplified lens which may reduce the applicability of our methods in more complex, real-world applications. Indeed, we consider accesses to memory in isolation, decoupled from computation and cache-polluting effects. We assume single-threaded execution without scheduling, which means that our caches will not be polluted by processes sharing (parts of) the cache hierarchy, nor will the application have its cached data evicted due to context switching. We also assume scalar, in-order execution of memory accesses. Finally, we take an optimistic view of cache latencies, using the fastest load-to-use latencies provided by hardware manufacturers; in real-world scenarios, cache latencies may be both more pessimistic and less stable than we assume. The results shown in <ref> indicate, however, that our fitness function is sufficiently accurate to be effective in real hardware.

In addition, the family of array layouts described in this work requires array sizes to be powers of two in each dimension. In applications where this is not the case, arrays must be over-allocated. For $n$-dimensional applications, using the layouts described in this paper requires over-allocation by a factor of $\mathcal{O}(2^n)$. Furthermore, applications using such layouts must consider the use of SIMD vectorization: it remains an open question which operations on arrays laid out in non-standard ways can be (automatically) vectorized. We have argued for the feasibility of SIMD in Morton-like arrays in <ref>.

Finally, our work considers only multiset permutations, in which the rank significance of bits in the input indices is preserved. This decision is based on current commodity hardware, which is capable of efficiently permuting bits only under this condition. There exists an even larger family of layouts in which rank bit significance is not preserved[That is to say, the layout $[0_0, 0_1, 1_0, 1_1]$ (which draws its least significant bit from the least significant bit of the first index) is distinct from the layout $[0_1, 0_0, 1_0, 1_1]$ (which instead draws its least significant bit from the second-least significant bit of the first index).]; such layouts could be of practical use in theoretical future processors with more advanced bit manipulation instructions, or in current FPGA and ASIC devices which permit the implementation of custom bit manipulation operations. Although we have not tested our approach on this further generalization, we are confident that an evolutionary approach like the one presented in this paper could be beneficial in exploring this (even larger) design space.

§ CONCLUSIONS AND FUTURE WORK

In this paper, we have discussed a generalization of the Morton layout for multi-dimensional data and we have shown that there exist families of array layouts with strongly varying cache behavior which, in turn, impact the performance of applications. We have shown how these layouts can be systematically described, and that the number of possible layouts quickly exceed the limits of what can be feasibly explored using exhaustive search. We have proposed a method based on evolutionary algorithms for the exploration of the design space of such layouts. We have evaluated the fitness of different array layouts using cache simulation and we have presented results indicating that our fitness function correlates with real world performance. Furthermore, we have shown that the methodology described in this paper can be used to improve the performance of applications on real hardware by up to ten times.

In the future, we intend to investigate the use of multi-objective optimization using NSGA-II <cit.> in order to find array layouts that provide favorable cache behavior across multiple applications. We also intend to explore more advanced genetic algorithms which are known to perform well in combinatorial problems, such as RKGA <cit.> and BRKGA <cit.>. It is our belief that exploring more evolutionary strategies will give us more insight into the convergence properties of various methods, and allow us to select the most efficient one. Although our fitness function correlates with real-world performance, the correlation is not perfect; we believe that the efficacy of our method could be improved through the development of more advanced fitness function, perhaps through the use of machine learning methods. In particular, we believe that the field of metric learning may enable us to develop more accurate fitness functions, and we aim to explore this avenue of research in the future. Finally, we aim to expand our research to a broader range of access patterns and hardware, including graphics processing units (GPUs).

The work presented in this paper was done in the context of the CERN Doctoral Student Programme. Many of the experimental results shown in this paper were gathered on the Advanced School for Computing and Imaging (ASCI) DAS-6 compute cluster <cit.>.

#1     #1#1#1     #1     #1     #1     #1#1          #1#1

[Abel and Reineke(2019)]
authorpersonAndreas Abel and personJan
Reineke. year2019.
uops.info: Characterizing Latency, Throughput,
and Port Usage of Instructions on Intel Microarchitectures. In
booktitleProceedings of the Twenty-Fourth International
Conference on Architectural Support for Programming Languages and Operating
Systems (Providence, Rhode Island, USA) (seriesASPLOS
'19). publisherAssociation for Computing Machinery,
addressNew York, New York, USA, pages673–686.

[Advanced Micro Devices, Inc.(2014)]
authorpersonAdvanced Micro Devices, Inc.
booktitleAMD EPYC 7413.
Advanced Micro Devices, Inc.

[Advanced Micro Devices, Inc.(2020)]
authorpersonAdvanced Micro Devices, Inc.
titleSoftware optimization guide for the AMD family 19h

[Akshintala et al(2019)]
authorpersonAmogh Akshintala, personBhushan
Jain, personChia-Che Tsai, personMichael Ferdman,
and personDonald E. Porter. year2019.
x86-64 Instruction Usage among C/C++
Applications. In booktitleProceedings of the 12th ACM
International Conference on Systems and Storage (Haifa, Israel)
(seriesSYSTOR '19). publisherAssociation for
Computing Machinery, addressNew York, New York, USA,

[Al-Kharusi and Walker(2019)]
authorpersonIbrahim Al-Kharusi and
personDavid W. Walker. year2019.
Locality properties of 3D data orderings with
application to parallel molecular dynamics simulations.
journalThe International Journal of High Performance
Computing Applications volume33, number5
(year2019), pages998–1018.

[Alves et al(2022)]
authorpersonJoão Nuno Ferreira Alves,
personLuís Manuel Silveira Russo, and
personAlexandre Francisco. year2022.
Cache-Oblivious Hilbert Curve-Based Blocking Scheme
for Matrix Transposition.
journalACM Trans. Math. Softw.
volume48, number4, Article articleno37
(datedec year2022), numpages28 pages.

[Anderson and Ferris(1994)]
authorpersonEdward J. Anderson and
personMichael C. Ferris. year1994.
Genetic Algorithms for Combinatorial Optimization:
The Assemble Line Balancing Problem.
journalORSA Journal on Computing
volume6, number2 (year1994),

[Asanovic et al(2009)]
authorpersonKrste Asanovic, personRastislav
Bodik, personJames Demmel, personTony Keaveny,
personKurt Keutzer, personJohn Kubiatowicz,
personNelson Morgan, personDavid Patterson,
personKoushik Sen, personJohn Wawrzynek,
personDavid Wessel, and personKatherine Yelick.
A View of the Parallel Computing Landscape.
journalCommun. ACM volume52,
number10 (dateoct year2009),

authorpersonMichael Bader.
booktitleSpace-filling curves: an introduction with
applications in scientific computing. Vol. volume9.
publisherSpringer Science & Business Media,
addressBerlin, Heidelberg, Germany.

[Bal et al(2016)]
authorpersonHenri Bal, personDick Epema,
personCees de Laat, personRob van Nieuwpoort,
personJohn Romein, personFrank Seinstra,
personCees Snoek, and personHarry Wijshoff.
A Medium-Scale Distributed System for Computer
Science Research: Infrastructure for the Long Term.
journalComputer volume49,
number05 (date5 year2016),

authorpersonJames Carl Bean.
Genetic Algorithms and Random Keys for Sequencing
and Optimization.
journalORSA Journal on Computing
volume6, number2 (year1994),

authorpersonRichard A. Brualdi.
booktitleIntroductory combinatorics
(edition5 ed.).
publisherPearson Education, addressLondon,
United Kingdom.

[Chatterjee et al(1999a)]
authorpersonSiddhartha Chatterjee,
personVibhor V. Jain, personAlvin R. Lebeck,
personShyam Mundhra, and personMithuna
Thottethodi. year1999a.
Nonlinear Array Layouts for Hierarchical Memory
Systems. In booktitleProceedings of the 13th International
Conference on Supercomputing (Rhodes, Greece) (seriesICS
'99). publisherAssociation for Computing Machinery,
addressNew York, New York, USA, pages444–453.

[Chatterjee et al(1999b)]
authorpersonSiddhartha Chatterjee,
personAlvin R. Lebeck, personPraveen K. Patnala,
and personMithuna Thottethodi.
Recursive Array Layouts and Fast Parallel Matrix
Multiplication. In booktitleProceedings of the Eleventh
Annual ACM Symposium on Parallel Algorithms and Architectures (Saint Malo,
France) (seriesSPAA '99). publisherAssociation
for Computing Machinery, addressNew York, New York, USA,

[Che et al(2011)]
authorpersonShuai Che, personJeremy W.
Sheaffer, and personKevin Skadron.
Dymaxion: Optimizing Memory Access Patterns for
Heterogeneous Systems. In booktitleProceedings of 2011
International Conference for High Performance Computing, Networking, Storage
and Analysis (Seattle, Washington) (seriesSC '11).
publisherAssociation for Computing Machinery,
addressNew York, New York, USA, Article articleno13,
numpages11 pages.

[Codreanu et al(2017)]
authorpersonVali Codreanu, personJoerg
Hertzer, personCristian Morales, personJorge
Rodriguez, personOle Widar Saastad, and
personMartin Stachon. year2017.
titleBest practice guide Haswell / Broadwell.

authorpersonLawrence Davis.
Applying Adaptive Algorithms to Epistatic Domains.
In booktitleProceedings of the 9th International Joint
Conference on Artificial Intelligence - Volume 1 (Los Angeles, California,
USA) (seriesIJCAI'85). publisherMorgan Kaufmann
Publishers Inc., addressSan Francisco, California, USA,

[Deb et al(2002)]
authorpersonKalyanmoy Deb, personAmrit
Pratap, personSameer Agarwal, and
personThirunavukarasu Meyarivan.
A fast and elitist multiobjective genetic
algorithm: NSGA-II.
journalIEEE Transactions on Evolutionary
Computation volume6, number2
(year2002), pages182–197.

[Deford and Kalyanaraman(2013)]
authorpersonDaryl Deford and personAnanth
Kalyanaraman. year2013.
Empirical Analysis of Space-Filling Curves for
Scientific Computing Applications. In booktitle2013 42nd
International Conference on Parallel Processing (Lyon, France).
publisherIEEE, addressNew York, New York, USA,

[Eiben and Smith(2015)]
authorpersonAgoston E. Eiben and
personJim E. Smith. year2015.
booktitleRepresentation, Mutation, and
publisherSpringer Berlin Heidelberg,
addressBerlin, Heidelberg, Germany, pages49–78.

[Evers et al(2022)]
authorpersonMark Evers, personLeslie
Barnes, and personMike Clark.
The AMD Next-Generation “Zen 3” Core.
journalIEEE Micro volume42,
number3 (year2022), pages7–12.

[Frigo et al(1999)]
authorpersonMatteo Frigo, personCharles E.
Leiserson, personHarald Prokop, and personSridhar
Ramachandran. year1999.
Cache-oblivious algorithms. In
booktitle40th Annual Symposium on Foundations of Computer
Science (New York, New York, USA) (seriesFOCS '99).
publisherIEEE, addressNew York, New York, USA,

[Gabriel and Wise(2004)]
authorpersonSteven T. Gabriel and
personDavid S. Wise. year2004.
The Opie Compiler from Row-Major Source to
Morton-Ordered Matrices. In booktitleProceedings of the 3rd
Workshop on Memory Performance Issues: In Conjunction with the 31st
International Symposium on Computer Architecture (Munich, Germany)
(seriesWMPI '04). publisherAssociation for
Computing Machinery, addressNew York, New York, USA,

[Gartland-Jones and Copley(2003)]
authorpersonAndrew Gartland-Jones and
personPeter Copley. year2003.
The Suitability of Genetic Algorithms for Musical
journalContemporary Music Review
volume22, number3 (year2003),

[Gonçalves and Resende(2011)]
authorpersonJosé Fernando Gonçalves and
personMauricio G. C. Resende. year2011.
Biased random-key genetic algorithms for
combinatorial optimization.
journalJournal of Heuristics volume17,
number5 (date01 Oct year2011),

[Gottschling et al(2007)]
authorpersonPeter Gottschling,
personDavid S. Wise, and personMichael D. Adams.
Representation-Transparent Matrix Algorithms with
Scalable Performance. In booktitleProceedings of the 21st
Annual International Conference on Supercomputing (Seattle, Washington,
USA) (seriesICS '07). publisherAssociation for
Computing Machinery, addressNew York, New York, USA,

[Gottschling et al(2009)]
authorpersonPeter Gottschling,
personDavid S. Wise, and personAdwait Joshi.
Generic support of algorithmic and structural
recursion for scientific computing.
journalInternational Journal of Parallel, Emergent
and Distributed Systems volume24, number6
(year2009), pages479–503.

[Hammarlund et al(2014)]
authorpersonPer Hammarlund, personAlberto J.
Martinez, personAtiq A. Bajwa, personDavid L. Hill,
personErik Hallnor, personHong Jiang,
personMartin Dixon, personMichael Derr,
personMikal Hunsaker, personRajesh Kumar,
personRandy B. Osborne, personRavi Rajwar,
personRonak Singhal, personReynold D'Sa,
personRobert Chappell, personShiv Kaushik,
personSrinivas Chennupaty, personStephan Jourdan,
personSteve Gunther, personTom Piazza, and
personTed Burton. year2014.
Haswell: The Fourth-Generation Intel Core
journalIEEE Micro volume34,
number2 (year2014), pages6–20.

[Hammer et al(2017)]
authorpersonJulian Hammer, personJan
Eitzinger, personGeorg Hager, and personGerhard
Wellein. year2017.
Kerncraft: A Tool for Analytic Performance Modeling
of Loop Kernels. In booktitleTools for High Performance
Computing 2016 (Stuttgart, Germany),
editorpersonChristoph Niethammer,
personJosé Gracia, personTobias Hilbrich,
personAndreas Knüpfer, personMichael M. Resch,
and personWolfgang E. Nagel (Eds.).
publisherSpringer International Publishing,
addressCham, Switzerland, pages1–22.

[Hegerty et al(2009)]
authorpersonBrian Hegerty, personChih-Cheng
Hung, and personKristen Kasprak.
A comparative study on differential evolution and
genetic algorithms for some combinatorial problems. In
booktitleProceedings of the 8th Mexican International
Conference on Artificial Intelligence (Guanajuato, Mexico),
Vol. volume9. publisherSpringer Verlag,
addressBerlin, Heidelberg, Germany, pages13.

authorpersonDavid Hilbert.
Ueber die stetige Abbildung einer Linie auf ein
journalMath. Ann. volume38,
number3 (year1891), pages459–460.

authorpersonRyutaro Himeno.
titleThe Riken Himeno CFD Benchmark.

authorpersonJohn Henry Holland.
booktitleAdaptation in natural and artificial
systems: an introductory analysis with applications to biology, control, and
artificial intelligence.
publisherMIT Press, addressCambridge,
Massachusetts, USA.

authorpersonJohn Henry Holland.
Genetic algorithms.
journalScientific American volume267,
number1 (year1992), pages66–73.

authorpersonOEIS Foundation Inc.
titleMoser–de Bruijn sequence, entry A000695 in The
On-Line Encyclopedia of Integer Sequences.

[Intel Corporation(2014)]
authorpersonIntel Corporation.
booktitleIntel Xeon Processor E5-2660 v3.
Intel Corporation.

[Intel Corporation(2023)]
authorpersonIntel Corporation.
titleIntel 64 and IA-32 Architectures Software
Developer's Manual.

[Jakob et al(2017)]
authorpersonWenzel Jakob, personJason
Rhinelander, and personDean Moldovan.
titlepybind11 – Seamless operability between C++11 and

[Jang et al(2011)]
authorpersonByunghyun Jang, personDana
Schaa, personPerhaad Mistry, and personDavid
Kaeli. year2011.
Exploiting Memory Access Patterns to Improve Memory
Performance in Data-Parallel Architectures.
journalIEEE Transactions on Parallel and Distributed
Systems volume22, number1 (year2011),

[Kolda and Bader(2009)]
authorpersonTamara G. Kolda and
personBrett W. Bader. year2009.
Tensor Decompositions and Applications.
journalSIAM Rev. volume51,
number3 (year2009), pages455–500.

[Kowarschik and Weiß(2003)]
authorpersonMarkus Kowarschik and
personChristian Weiß. year2003.
booktitleAn Overview of Cache Optimization
Techniques and Cache-Aware Numerical Algorithms.
publisherSpringer Berlin Heidelberg,
addressBerlin, Heidelberg, Germany, pages213–232.

[Laukemann et al(2018)]
authorpersonJan Laukemann, personJulian
Hammer, personJohannes Hofmann, personGeorg Hager,
and personGerhard Wellein. year2018.
Automated Instruction Stream Throughput Prediction
for Intel and AMD Microarchitectures. In booktitle2018
IEEE/ACM Performance Modeling, Benchmarking and Simulation of High
Performance Computer Systems (PMBS) (Dallas, Texas, USA).
publisherIEEE, addressNew York, New York, USA,

[Mellor-Crummey et al(2001)]
authorpersonJohn Mellor-Crummey, personDavid
Whalley, and personKen Kennedy.
Improving Memory Hierarchy Performance for
Irregular Applications Using Data and Computation Reorderings.
journalInternational Journal of Parallel
Programming volume29, number3 (date01
Jun year2001), pages217–247.

[Miller and Goldberg(1995)]
authorpersonBrad L. Miller and
personDavid E. Goldberg. year1995.
Genetic algorithms, tournament selection, and the
effects of noise.
journalComplex systems volume9,
number3 (year1995), pages193–212.

authorpersonGuy Macdonald Morton.
booktitleA Computer Oriented Geodetic Data Base and
a New Technique in File Sequencing.
typeTechnical Report.
institutionInternational Business Machines Company.

authorpersonHeinz Mühlenbein.
Parallel genetic algorithms, population genetics
and combinatorial optimization. In booktitleParallelism,
Learning, Evolution (Wildbad Kreuth, Germany). publisherSpringer
Berlin Heidelberg, addressBerlin, Heidelberg, Germany,

[Park et al(2003)]
authorpersonNeungsoo Park, personBo Hong,
and personViktor K. Prasanna. year2003.
Tiling, block data layout, and memory hierarchy
journalIEEE Transactions on Parallel and Distributed
Systems volume14, number7 (year2003),

[Pavai and Geetha(2016)]
authorpersonG. Pavai and personT. V.
Geetha. year2016.
A Survey on Crossover Operators.
journalACM Comput. Surv. volume49,
number4, Article articleno72 (datedec
year2016), numpages43 pages.

[Pawłowski et al(2019)]
authorpersonFilip Pawłowski, personBora
Uçar, and personAlbert-Jan Nicholas Yzelman.
A multi-dimensional Morton-ordered block storage
for mode-oblivious tensor computations.
journalJournal of Computational Science
volume33 (year2019), pages34–44.

[Perdacher et al(2020)]
authorpersonMartin Perdacher, personClaudia
Plant, and personChristian Böhm.
Improved Data Locality Using Morton-order Curve
on the Example of LU Decomposition. In booktitle2020 IEEE
International Conference on Big Data (Big Data) (Atlanta, Georgia, USA).
publisherIEEE, addressNew York, New York, USA,

authorpersonAndy D. Pimentel.
Exploring Exploration: A Tutorial Introduction to
Embedded Systems Design Space Exploration.
journalIEEE Design & Test volume34,
number1 (year2017), pages77–90.

[Qureshi et al(2008)]
authorpersonMoinuddin K. Qureshi,
personAamer Jaleel, personYale N. Patt,
personSimon C. Steely, and personJoel Emer.
Set-Dueling-Controlled Adaptive Insertion for
High-Performance Caching.
journalIEEE Micro volume28,
number1 (year2008), pages91–98.

[Sapra and Pimentel(2020)]
authorpersonDolly Sapra and personAndy D.
Pimentel. year2020.
Constrained Evolutionary Piecemeal Training to
Design Convolutional Neural Networks. In booktitleTrends in
Artificial Intelligence Theory and Applications. Artificial Intelligence
Practices. publisherSpringer International Publishing,
addressCham, Switzerland, pages709–721.

[Sivanandam and Deepa(2008)]
authorpersonS.N. Sivanandam and personS.N.
Deepa. year2008.
booktitleGenetic Algorithms.
publisherSpringer Berlin Heidelberg,
addressBerlin, Heidelberg, Germany, pages15–37.

[Slowik and Kwasnicka(2020)]
authorpersonAdam Slowik and personHalina
Kwasnicka. year2020.
Evolutionary algorithms and their applications to
engineering problems.
journalNeural Computing and Applications
volume32, number16 (date01 Aug
year2020), pages12363–12379.

[Srinivas and Patnaik(1994)]
authorpersonM. Srinivas and personLalit M.
Patnaik. year1994.
Genetic algorithms: a survey.
journalComputer volume27,
number6 (year1994), pages17–26.

[Stephens et al(2017)]
authorpersonNigel Stephens, personStuart
Biles, personMatthias Boettcher, personJacob Eapen,
personMbou Eyole, personGiacomo Gabrielli,
personMatt Horsnell, personGrigorios Magklis,
personAlejandro Martinez, personNathanael Premillieu,
personAlastair Reid, personAlejandro Rico, and
personPaul Walker. year2017.
The ARM Scalable Vector Extension.
journalIEEE Micro volume37,
number2 (year2017), pages26–39.

[Swatman et al(2023)]
authorpersonStephen Nicholas Swatman,
personAna-Lucia Varbanescu, personAndy Pimentel,
personAndreas Salzburger, and personAttila
Krasznahorkay. year2023.
Systematically Exploring High-Performance
Representations of Vector Fields Through Compile-Time Composition. In
booktitleProceedings of the 2023 ACM/SPEC International
Conference on Performance Engineering (Coimbra, Portugal)
(seriesICPE '23). publisherAssociation for
Computing Machinery, addressNew York, New York, USA,

[Terfloth and Gasteiger(2001)]
authorpersonLothar Terfloth and
personJohann Gasteiger. year2001.
Neural networks and genetic algorithms in drug
journalDrug Discovery Today  volume6
(year2001), pages102–108.

[Thiyagalingam et al(2006)]
authorpersonJeyarajan Thiyagalingam,
personOlav Beckmann, and personPaul H. J. Kelly.
Is Morton layout competitive for large
two-dimensional arrays yet?
journalConcurrency and Computation: Practice and
Experience volume18, number11
(year2006), pages1509–1539.

[Vila et al(2020)]
authorpersonPepe Vila, personPierre Ganty,
personMarco Guarnieri, and personBoris Köpf.
CacheQuery: Learning Replacement Policies from
Hardware Caches. In booktitleProceedings of the 41st ACM
SIGPLAN Conference on Programming Language Design and Implementation
(London, United Kingdom) (seriesPLDI 2020).
publisherAssociation for Computing Machinery,
addressNew York, New York, USA, pages519–532.

authorpersonDavid W. Walker.
Morton ordering of 2D arrays for efficient
access to hierarchical memory.
journalThe International Journal of High Performance
Computing Applications volume32, number1
(year2018), pages189–203.

[Walker and Skjellum(2023)]
authorpersonDavid W. Walker and
personAnthony Skjellum. year2023.
titleThe Impact of Space-Filling Curves on Data Movement
in Parallel Systems.

[arxiv]2307.07828 [cs.DC]

[Weinberg et al(2005)]
authorpersonJonathan Weinberg,
personMichael O. McCracken, personErich Strohmaier,
and personAllan Snavely. year2005.
Quantifying Locality In The Memory Access Patterns
of HPC Applications. In booktitleSC '05: Proceedings of the
2005 ACM/IEEE Conference on Supercomputing (Seattle, Washington, USA).
publisherIEEE, addressNew York, New York, USA,

[Yzelman and Bisseling(2012)]
authorpersonAlbert-Jan Nicholas Yzelman and
personRob Hendrik Bisseling. year2012.
A Cache-Oblivious Sparse Matrix–Vector
Multiplication Scheme Based on the Hilbert Curve. In
booktitleProgress in Industrial Mathematics at ECMI 2010.
publisherSpringer Berlin Heidelberg, addressBerlin,
Heidelberg, Germany, pages627–633.