(De/Re)-Composition of Data-Parallel Computations via Multi-Dimensional Homomorphisms

We define input views as any function that compute an MDA from a collection of (user-defined) BUFs. For example, in the case of MatVec, its input view takes as input two BUFs—a matrix and a vector—and it yields a two-dimensional MDA containing pairs of matrix and vector elements (illustrated in Figure 1). In contrast, the input view of Jacobi1D takes as input a single BUF (representing a vector) only, and it computes an MDA containing triples of BUF elements (Figure 2).

In the following, we introduce higher-order function inp_view which conveniently computes important input views from user-defined index functions \(\mathfrak{idx}_{b,a}:\{0,1,\dotsc\}\to\{0,1,\dotsc\}\), \(b\in[1,B]_{\mathbb{N}}\), \(a\in[1,A_{b}]_{\mathbb{N}}\), in a uniform, structured manner. Here, \(B\in\mathbb{N}\) represents the number of BUFs that the computed input view will take as input, and \(A_{b}\) represents the number of accesses to the \(b\)-th BUF required for computing an individual MDA element.

In the case of MatVec (Figure 1), we use \(B:=2\) because MatVec has two input BUFs: a matrix \(M\) (the first input of MatVec and thus identified by \(b=1\)) and a vector \(v\) (identified by \(b=2\)). For the number of accesses, we use for the matrix \(A_{1}:=1\), as one element is accessed within matrix \(M\) to compute an individual MDA element—matrix element \(M[i,k]\) for computing MDA element at position \((i,k)\). For the vector, we use \(A_{2}:=1\), as the single element \(v[k]\) is accessed within the vector. The index functions of MatVec are \(\mathfrak{idx}_{1,1}(i,k):=(i,k)\) (used to access the matrix) and \(\mathfrak{idx}_{2,1}(i,k):=(k)\) (used for the vector).

In contrast, for Jacobi1D (Figure 2), we use \(B:=1\) because Jacobi1D has vector \(v\) as its only input, and we use \(A_{1}:=3\) because the vector is accessed three times to compute an individual MDA element at arbitrary position \(i\): first access \(v[i+0]\), second access \(v[i+1]\), and third access \(v[i+2]\). The index functions of Jacobi1D are \(\mathfrak{idx}_{1,1}(i):=(i+0)\), \(\mathfrak{idx}_{1,2}(i):=(i+1)\), and \(\mathfrak{idx}_{1,3}(i):=(i+2)\).

Figures 10 and 11 use the examples MatVec and Jacobi1D to informally illustrate how function inp_view uses index functions to compute input views. In the two figures, we use domain-specific identifiers for better clarity: in the case of MatVec, we use for its two input BUFs the identifiers \(M\) and \(v\) instead of \(\mathfrak{b}_{1}\) and \(\mathfrak{b}_{2}\), as well as identifiers \(i\) and \(j\) instead of \(i_{1}\) and \(i_{2}\) for index variables; for Jacobi1D, we use identifier \(v\) instead of \(\mathfrak{b}_{1}\) and \(i\) instead of \(i_{1}\).

Higher-order function inp_view takes as input a collection of index functions of types IDX-FCT, and it computes an input view of type IV (Definition 6) based on the index functions, as illustrated in Figures 10 and 11.

Note that function inp_view is not capable of computing every kind of input view function (Definition 6). For example, inp_view cannot be used for computing MDAs that are required for expressing computations on sparse data formats [Hall, 2020], because such MDAs need dynamically accessing BUFs. This limitation of inp_view can be relaxed by generalizing our index functions toward taking additional, dynamic input arguments, which we consider as future work (as outlined in Section 8).

An output view is the counterpart of an input view: in contrast to an input view which maps BUFs to an MDA, an output view maps an MDA to a collection of BUFs. In the following, we define output views, and we introduce higher-order function out_view which computes output views in a structured manner (analogously to function inp_view for input views).

Figures 12 and 13 illustrate output views informally using the examples transposed Matrix Multiplication and Double Reduction.

In the case of transposed matrix multiplication (Figure 12), the computed output MDA (the computation of matrix multiplication is presented later and not relevant for our following considerations) is stored via an output view as a matrix in a transposed fashion, using index function \((i,j,0)\mapsto(j,i)\). Here, the MDA’s third dimension (accessed via index \(0\)) represents the so-called reduction dimension of matrix multiplication, and it contains only one element after the computation, as all elements in this dimension are combined via addition.

For double reduction (Figures 13), we combine the elements within the vector twice—once using operator \(\oplus\) (e.g., \(\oplus=+\) addition) and once using operator \(\otimes\) (e.g, \(\otimes=*\) multiplication). The final outcome of double reduction is a singleton MDA containing a pair of two elements that represent the combined vector elements (e.g., the elements’ sum and product). We store this MDA via an output view as two individual scalar values, using index functions \((0)\mapsto()\) ¹¹ for both pair elements.

We define higher-order function out_view formally as follows.

We use view functions to transform data from their domain-specific representation (represented in our formalism as BUFs, Definition 5) to our internal, MDA-based representation (via input views) and back (via output views), as also illustrated in Figure 5. In our implementation presented later, we aim to access data uniformly in the form of MDAs, thereby being independent of domain-specific data representations. However, we aim to store the data physically in the domain-specific format, as such format is usually the more efficient representation. For example, we aim to store the input data of MatVec in the domain-specific matrix and vector format, rather than as an MDA, because the input MDA of MatVec contains many redundancies—each vector element once per row of the input matrix (as illustrated in Figure 10).

The following lemma proves that functions inp_view and out_view are invertible and that they are each others inverses. Consequently, the lemma shows how we can arbitrarily switch between the domain-specific and our MDA-based representation, and consequently also that we can implicitly identify MDAs with the domain-specific data representation. For example, for computing MatVec, we will specify the computations via pattern md_hom which operates on MDAs (see Figure 5), but we use the view functions in our implementation to implicitly forward the MDA accesses to the physically stored BUF representation.

The following figure illustrates the lemma using as example the inverse of MatVec’s input view (shown in Figure 10):

We show how our high-level representation is used for expressing linear algebra routines: (1) Dot (Dot Product); (2) MatVec (Matrix-Vector Multiplication); (3) MatMul (Matrix Multiplication); (4) \(\texttt{MatMul}^{\texttt{T}}\) (Transposed Matrix Multiplication) which computes matrix multiplication on transposed input and output matrices; (5) bMatMul (batched Matrix Multiplication) where multiple matrix multiplications are computed using matrices of the same sizes.

We can observe from the subfigure that our high-level expressions for the routines naturally evolve from each other. For example, the md_hom instance for MatVec differs from the md_hom instance for Dot by only containing a further concatenation dimension ++ for its \(i\) dimension. We consider this close relation between the high-level expressions of MatVec and Dot in our approach as natural and favorable, as MatVec can be considered as computing multiple times Dot—one computation of Dot for each value of MatVec’s \(i\) dimension. Similarly, the md_hom instance for MatMul is very similar to the expression of MatVec, by containing the further concatenation dimension \(j\) for MatMul’s \(j\) dimension. The same applies to bMatMul: its md_hom instance is the expression of MatMul augmented with one further concatenation dimension.

Regarding \(\texttt{MatMul}^{\texttt{T}}\), the basic computation part of \(\texttt{MatMul}^{\texttt{T}}\) and MatMul are the same, which is exactly reflected in our formalisms: both \(\texttt{MatMul}^{\texttt{T}}\) and MatMul are expressed using exactly the same md_hom instances. The differences between \(\texttt{MatMul}^{\texttt{T}}\) and MatMul lies only in the data accesses—transposed accesses in the case of \(\texttt{MatMul}^{\texttt{T}}\) and non-transposed accesses in the case of MatMul. Data accesses are expressed in our formalism, in a structured way, via view functions (as discussed in Section 2.3): for example, for \(\texttt{MatMul}^{\texttt{T}}\), we use for its first input matrix \(A\) the index function \((i,j,k)\mapsto(k,i)\) for transposed access, instead of using index function \((i,j,k)\mapsto(i,k)\) as for MatMul’s non-transposed accesses.

Note that all md_hom instances in the subfigure are well defined according to Lemma 1.

We show how convolution-style stencil computations are expressed in our high-level representation: (1) Conv2D expresses a standard convolution that uses a two-dimensional sliding window [Podlozhnyuk, 2007]; (2) MCC expresses a so-called Multi-Channel Convolution [Dumoulin and Visin, 2018]—a generalization of Conv2D that is heavily used in the area of deep learning; (3) MCC_Capsule is a recent generalization of MCC [Hinton et al., 2018] which attracted high attention due to its relevance for advanced deep learning neural networks [Barham and Isard, 2019].

While our md_hom instances for convolutions are quite similar to those of linear algebra routines (they all use multiplication \(*\) as scalar function and a mix of concatenations ++ and point-wise additions \(+\) as combine operators), the index functions used for the view functions of convolutions are notably different from those used for linear algebra routines: the index functions of convolutions contain arithmetic expressions (e.g., p+r and q+s) and thus access neighboring elements in their input—a typical access pattern in stencil computations that requires special optimizations [Hagedorn et al., 2018]. Moreover, convolution-style computations are often high-dimensional (e.g., \(10\) dimensions in the case of MCC_Capsule), whereas linear algebra routines usually rely on less dimensions. Our experiments in Section 5 confirm that respecting the data access patterns and the high dimensionality of convolutions in the optimization process (as in our approach, which we discuss later) often achieves significantly higher performance than using optimizations chosen toward linear algebra routines, as in vendor libraries provided by NVIDIA and Intel for convolutions [Li et al., 2016].

We show how quantum chemistry computation Coupled Cluster (CCSD(T)) [Kim et al., 2019] is expressed in our high-level representation. The computation of CCSD(T) notably differs from those of linear algebra routines and convolution-style stencils, by accessing its high-dimensional input data in sophisticated transposed fashions: for example, the view function of CCSD(T)’s instance one (denoted as I1 in the subfigure) uses indices a and b to access the last two dimensions of its \(A\) input tensor (rather than the first two dimensions of the tensor, as would be the case for non-transposed accesses). For brevity, the subfigure presents only two CCSD(T) instances—in our experiments in Section 5, we present experimental results for nine different real-world CCSD(T) instances.

The subfigures present computations whose scalar functions and combine operators are different from those used in Subfigures 1–3 (which are in Subfigures 1–3 straightforward multiplications \(*\), concatenation ++, and point-wise additions \(+\) only). For example, Jacobi stencils (Subfigure 4) use as scalar function the Jacobi-specific computation \(\texttt{J}_{\texttt{nD}}\) [Cecilia et al., 2012], and Probabilistic Record Linkage (PRL) [Christen, 2012], which is heavily used in data mining to identify duplicate entries in a database, uses a PRL-specific both scalar function wght and combine operator \(\texttt{max}_{\texttt{PRL}}\) (point-wise combination via the PRL-specific binary operator \(\texttt{max}_{\texttt{PRL}}\)) [Rasch et al., 2019b]. Histograms, in their generalized version [Henriksen et al., 2020] (denoted as GenHisto in Subfigure 6), use an arbitrary, user-defined scalar function \(f\) and a user-defined associative and commutative combine operator \(\oplus\); the standard histogram variant Histo is then a particular instance of GenHist, for \(\oplus=+\) (point-wise addition) and \(f=f_{\texttt{Histo}}\), where \(f_{\texttt{Histo}}(e,b)=1\) iff \(e=b\) and \(f_{\texttt{Histo}}(e,b)=0\) otherwise.

We show how typical map and reduce patterns [González-Vélez and Leyton, 2010] are implemented in our high-level representation. Examples map(f) and reduce(\(\oplus\)) (discussed in Examples 3 and 4) are simple and thus straightforwardly expressed in our representation. In contrast, example reduce(\(\oplus,\otimes\)) is more complex and shows how reduce(\(\oplus\)) is extended to combine the input vector simultaneously twice—once combining vector elements via operator \(\oplus\) and once using operator \(\otimes\). The outcome of reduce(\(\oplus,\otimes\)) are two scalars—one representing the result of combination via \(\oplus\) and the other of combination via \(\otimes\)—which we map via the output view to output elements \(\texttt{O}_{\texttt{1}}\) (result of \(\oplus\)) and \(\texttt{O}_{\texttt{2}}\) (result of \(\otimes\)), correspondingly; this is also illustrated in Figure 13.

We present prefix-sum computations [Blelloch, 1990] which differ from the computations in Subfigures 1–7 in terms of their combine operators: the operator used for expressing computations in Subfigure 8 is different from concatenation (Example 1) and point-wise combinations (Example 2). Computation scan(\(\oplus\)) uses as combine operator \(\mbox{++}_{\texttt{prefix-sum}}(\oplus)\) which computes prefix-sum [Gorlatch and Lengauer, 1997] (formally defined by Rasch [2024], Section B.9) according to binary operator \(\oplus\), and Maximum Bottom Box Sum (MBBS) [Farzan and Nicolet, 2019] uses a particular instance of prefix-sum for \(\oplus=+\) (addition).

The de-composition phase (steps 1–\(7\) in Figure 15) partitions input MDA \({{}^{\downarrow}}{\mathfrak{a}}\) (in the top right of Figure 15) to the structure \({{}^{\downarrow}}\mathfrak{a}_{f}^{ \lt \dotsc \gt }\) (bottom right) which we refer to as low-level MDA and define formally in the next subsection. The low-level MDA represents a partitioning of MDA \({{}^{\downarrow}}{\mathfrak{a}}\) (a.k.a hierarchical, multi-dimensional tiling in programming), where each particular choice of indices \(p^{1}_{1}\in[0,2)_{\mathbb{N}_{0}}\), \(p^{1}_{2}\in[0,4)_{\mathbb{N}_{0}}\), \(p^{2}_{1}\in[0,8)_{\mathbb{N}_{0}}\), \(p^{2}_{2}\in[0,16)_{\mathbb{N}_{0}}\), \(p^{3}_{1}\in[0,32)_{\mathbb{N}_{0}}\), \(p^{3}_{2}\in[0,64)_{\mathbb{N}_{0}}\) refers to an MDA that represents an individual part of MDA \({{}^{\downarrow}}{\mathfrak{a}}\) (a.k.a. tile in programming—informally illustrated in Figure 7). The partitions are arranged on multiple layers (indicated by the \(p\)’s superscripts) and in multiple dimensions (indicated by subscripts)—as illustrated in Figure 16—according to the memory/COR of the target architecture and dimensions of the MDH computation: we partition for each of the target architecture’s three layers (HM, L1, COR) and in each of the two dimensions of the MDH (dimensions 1 and 2, as we use example MatVec in Figure 15, which represents a two-dimensional MDH computation). Consequently, our partitioning approach allows efficiently exploiting each particular layer of the target architecture (both memory and core layers), and also optimizing for both dimensions of the target computation (in the case of MatVec, the \(i\)-dimension and also the \(k\)-dimension—see Figure 1), allowing fine-grained optimizations.

We compute the partitionings of MDAs by applying the concatenation operator (Example 1) inversely (indicated by using \(=:\) instead of \(:=\) in the top right part of Figure 15). For example, we partition in Figure 15 MDA \({{}^{\downarrow}}{\mathfrak{a}}\) first via the inverse of \(\mbox{++}^{\texttt{(HM,x)}}_{1}\) in dimension 1 (indicated by the subscript 1 of \(\mbox{++}^{\texttt{(HM,x)}}_{1}\); the superscript (HM,x) is explained later) into two parts, as \(p^{1}_{1}\) iterates over interval \([0,2)_{\mathbb{N}_{0}}=\{0,1\}\) which consists of two elements (\(0\) and 1)—the interval is chosen arbitrarily for this example. Afterward, each of the obtained parts is further partitioned, in the second dimension, via \(\mbox{++}^{\texttt{(HM,y)}}_{2}\) into four parts (\(p^{1}_{2}\) iterates over \([0,4)_{\mathbb{N}_{0}}=\{0,1,2,3\}\) which consists of four elements). The \({(2*4)}\)-many HM parts are then each further partitioned in both dimensions for the COR layer into \((8*16)\) parts, and each individual COR part is again partitioned for the L1 layer into \((32*64)\) parts, resulting in \((2*8*32)*(4*16*64)=512*4096\) parts in total.

We always use a full partitioning in our low-level expressions,¹² i.e., each particular choice of indices \(p^{1}_{1}\), \(p^{1}_{2}\), \(p^{2}_{1}\), \(p^{2}_{2}\), \(p^{3}_{1}\), \(p^{3}_{2}\) points to an MDA that contains a single element only (in Figure 16, the individual elements are denoted via symbol \(\times\), in the bottom part of the figure). By relying on a full partitioning, we can apply scalar function \(f\) to the fully partitioned MDAs later in the scalar phase (described in the next paragraph). This is because function \(f\) is defined on scalar values (Definition 4) to make defining scalar functions more convenient for the user (as discussed in Section2.2).

The superscript of combine operators, e.g., (COR,x) of operator \(\mbox{++}^{\texttt{(COR,x)}}_{1}\), is a so-called operator tag (formal definition given in the next section). Such a tag indicates to our code generator whether its combine operator is assigned to a memory layer (and thus computed sequentially in our generated code) or to a COR (and thus computed in parallel). For example, tag (COR,x) indicates that parts processed by operator \(\mbox{++}^{\texttt{(COR,x)}}_{1}\) should be computed by cores COR, and thus in parallel; the dimension tag x indicates that COR layer’s x dimension should be used for computing the operator (we use dimension x for our example architecture as an analogous concept to CUDA’s thread/block dimensions x,y,z for GPU architectures [NVIDIA, 2022g]), as we also discuss in the next section. In contrast, tag (HM,x) refers to a memory layer (HM) and thus, operator \(\mbox{++}^{\texttt{(HM,x)}}_{1}\) is computed sequentially. Since the current state-of-practice programming approaches (OpenMP, CUDA, OpenCL, \(\dotsc\)) have no explicit notion of memory tiles (e.g., by offering the potential variables tileIdx.x/tileIdx.y/tileIdx.z, as analogous concepts to CUDA variables threadIdx.x/threadIdx.y/threadIdx.z), the dimensions tag x in (HM,x) is currently ignored by our code generator, because HM refers to a memory layer.

Note that the number of parts (e.g., 2 parts on layer 1 in dimension 1, and \(4\) parts on layer 1 in dimension 2 \(\dotsc\)), the combine operators’ tags, and our partition order (e.g., first partitioning in MDA’s dimension 1 and afterwards in dimension 2) are chosen arbitrarily for this example. These choices are critical for performance and should be optimized¹³ for a particular target architecture and characteristics of the input and output data (size, memory layouts, etc.), as we discuss in detail later in this section.

In the scalar phase (step 8 in Figure 15), we apply MDH’s scalar function \(f\) to the individual MDA elementsfor each particular choice of indices \(p^{1}_{1}\), \(p^{1}_{2}\), \(p^{2}_{1}\), \(p^{2}_{2}\), \(p^{3}_{1}\), \(p^{3}_{2}\), which results in

In the figure, \(\vec{f}\) is the slight adaption of function \(f\) that operates on a singleton MDA, rather than a scalar (see Footnote 9).

Annotation \(\rightarrow\) \(\lt\) (1,2), \(\dotsc\)\(\gt\) indicates the application order of applying scalar function (in this example, first iterating over \(p^{1}_{1}\), then over \(p^{1}_{2}\), etc.), and we use annotation \(\rightarrow\) \(\lt\) (HM,x), \(\dotsc\)\(\gt\) to indicate how the scalar computation is assigned to the target architecture (this is described in detail later in this section). Annotations \(\rightarrow\) M: HM, v: L1 and \(\rightarrow\) w: L1 (in the bottom part of Figure 15) indicate the memory regions to be used for reading and writing the input scalar of function \(f\) (also described later in detail).

Finally, the re-composition phase (steps \(9\)–\(15\) in Figure 15) combines the computed parts \({{}^{\uparrow}}\mathfrak{a}_{f}^{ \lt p^{1}_{1},p^{1}_{2}\ |\ p^{2}_{1},p^{2}_{2 }\ |\ p^{3}_{1},p^{3}_{2} \gt }\) (bottom left in the figure) to the final result \({{}^{\uparrow}}{\mathfrak{a}}\) (top left) via MDH’s combine operators, which are in the case of matrix-vector multiplication \(\circledast_{1}:=\mbox{++}\) (concatenation) and \(\circledast_{2}:=+\) (point-wise addition). In this example, we first combine the L1 parts in dimension 2 and then in dimension 1; afterward, we combine the COR parts in both dimensions, and finally the HM parts. Analogously to before, this order of combine operators and their tags are chosen arbitrarily for this example and should be auto-tuned for high performance.

In the de- and re-composition phases, the arrow notation below combine operators allow efficiently exploiting architecture’s memory hierarchy, by indicating the memory region to read from (de-composition phase) or to write to (re-composition phase); the annotations also indicate the memory layouts to use. We exploit these memory and layout information in both (1) our code generation process to assign combine operators’ input and output data to memory regions and to chose memory layouts for the data (row major, column major, etc.); (2) our formalism to specify constraints of programming models, e.g., that in CUDA, results of GPU cores can only be combined in designated memory regions [NVIDIA, 2022f]. For example, annotation \(\rightarrow\) M: HM[1,2], v: L1[1] below an operator in the de-composition phase indicates to our code generator that the parts (a.k.a tiles) of matrix \(M\) used for this computation step should be read from the HM memory region and that parts of vector \(v\) should be copied to and accessed from fast L1 memory. The annotation also indicates that M should be stored using a row-major memory layout (as we use [1,2] and not [2,1]). The memory regions and layouts are chosen arbitrarily for this example and should be chosen as optimized (auto-tuned) for the particular target architecture and characteristics of the input and output data. Formally, the arrow notation of combination operators is a concise notation to hide MDAs and BUFs for intermediate results (discussed by Rasch [2024], Section C.3, for the interested reader).

Figures 17–22 report the performance of the TVM-generated code, which is in CUDA for GPUs and in OpenCL for CPUs. We observe that we usually achieve with our approach the high performance of TVM and often perform even better. For example, in Figure 21, we achieve a speedup \({\gt}2\times\) over TVM on NVIDIA Ampere GPU for matrix multiplications as used in the inference phase of the ResNet-50 neural network—an actually favorable example for TVM which is designed and optimized toward deep learning computations executed on modern NVIDIA GPUs. Our performance advantage over TVM is because we parallelize and optimize more efficiently reduction-like computations—in the case of MatMul (Figure 14), its \(3\)rd-dimension (a.k.a. \(k\)-dimension). The difficulties of TVM with reduction computations become particularly obvious when computing dot products (Dot) on GPUs (Figure 17): the Dot’s main computation part is a reduction computation (via point-wise addition, see Figure 14), thus requiring reduction-focused optimization, in particular when targeting the highly parallel architecture of GPUs: in the case of Dot (Figure 17), our generated CUDA code exploits parallelization over CUDA blocks, whereas the Ansor-generated TVM code exploits parallelization over threads within in a single block only, because TVM currently cannot use blocks for parallelizing reduction computations [Apache TVM Community, 2022a]. Furthermore, while TVM’s Ansor rigidly parallelizes outer dimensions [Zheng et al., 2020a], our ATF-based tuning process has auto-tuned our tuning parameters D2, S2, R2 in Table 1 to exploit parallelism for inner dimensions, which achieves higher performance for this particular MatMul example used in ResNet-50. Also, for MatMul-like computations, Ansor always caches parts of the input in GPU’s shared memory, and it computes these cached parts always in register memory. In contrast, our caching strategy is auto-tunable (via parameters D3, S3 S5, R3 in Table 1), and ATF has determined to not cache the input matrices into fast memory resources for the MatMul example in ResNet-50. Surprisingly, Ansor does not exploit fast memory resources for Jacobi stencils (Figure 18), as required to achieve high performance for them: our generated and auto-tuned CUDA kernel for Jacobi uses register memory for both inputs (image buffer and filter) when targeting NVIDIA Ampere GPU (small input size), thereby achieving a speedup over TVM+Ansor of \(1.93\times\) for Jacobi. Most likely, Ansor fails to foresee the potential of exploiting fast memory resources for Jacobi stencils, because the Jacobi’s index functions used for memory accesses (Figure 14) are injective. For the MatMul example of ResNet-50’s training phase (Figure 21), we achieve a speedup over TVM on NVIDIA Ampere GPU of \(1.26\times\), because auto-tuning determined to store parts of input matrix \(A\) as transposed into fast memory (via parameter D4 in Table 1). Storing parts of the input/output data as transposed is not considered by Ansor as optimization, perhaps because such optimization must be expressed in TVM’s high-level language, rather than in its scheduling language[Apache TVM Community, 2022c]. For MatVec on NVIDIA Ampere GPU (Figure 17), we achieve a speedup over TVM of \(1.22\times\) for the small input size, by exploiting a so-called swizzle pattern [Phothilimthana et al., 2019]: our ATF tuner has determined to assign threads that are consecutive in CUDA’s x-dimension to the second MDA dimension (via parameters D2, S2, R2 in Table 1), thereby accessing the input matrix in a GPU-efficient manner (a.k.a coalesced global memory accesses [NVIDIA, 2022f]). In contrast, for MatVec computations, Ansor assigns threads with consecutive x-ids always to the first data dimension, in a non-tunable manner, causing lower performance.

Our positive speedups over TVM on CPU are for the same reasons as discussed above for GPU. For example, we achieve a speedup of \({\gt}3\times\) over TVM on Intel Skylake CPU for MCC (Figure 22) as used in the training phase of the MobileNet neural network, because we exploit fast memory resources more efficiently than TVM: our auto-tuning process has determined to use register memory for the MCC’s second input (the filter buffer F, see Fig. 14) and using no fast memory for the first input (image buffer I), whereas Ansor uses shared memory rigidly for both inputs of MCC. Moreover, our auto-tuning process has determined to parallelize the inner dimensions of MCC, while Ansor always parallelizes outer dimensions. We achieve the best speedup over TVM for MCC on an input size taken from TVM’s own tutorials [Apache TVM Documentation, 2022b] (Figure 18), rather than from neural networks (as in Figures 21 and 22). This is because TVM’s MCC size includes large reduction computations, which are not efficiently optimized by TVM (as discussed above).

The TVM compiler achieves higher performance than our approach for some examples in Figures 17–22. However, in most cases, this is for a technical reason only: TVM uses the NVCC compiler for compiling CUDA code, whereas our proof-of-concept code generator currently relies on NVIDIA’s NVRTC library which surprisingly generates less-efficient CUDA assembly than NVCC. In three cases, the higher performance of TVM over our approach is because our ATF was not able to find a better performing tuning configuration than TVM’s Ansor optimization engine during our \(12\)h tuning time; the three cases are: (1) MCC from VGG-16’s inference phase on NVIDIA Ampere GPU (Figure 21), (2) MCC (capsule variant) from VGG-16’s training phase on NVIDIA Ampere GPU (Figure 21), and (3) MCC (capsule variant) from ResNet-50’s training phase on Intel Skylake CPU (Figure 22). However, when we manually set the Ansor-found tuning configurations also for our approach, instead of using the ATF-found configurations, we achieve for these three cases exactly the same high performance as TVM+Ansor, i.e., the well-performing configurations are contained in our search space (Table 1). Most likely, Ansor was able to find this well-performing configuration within the \(12\) h tuning time, because it explores a significantly smaller search space that is particularly designed for deep learning computations. To avoid such tuning issues in our approach, we aim to substantially improve our auto-tuning process in future work: we plan to introduce an analytical cost model that assists (or even replaces) our auto-tuner, as we also outline in Section 8.

Note that the TVM compiler crashes for our data mining example PRL, because TVM has difficulties with computations relying on user-defined combine operators [Apache TVM Community, 2022d].

Figure 23 reports the portability of the TVM compiler. Our portability measurements are based on the Pennycook metric where a value close to 1 indicates high portability and a value close to \(0\) indicates low portability, correspondingly. We observe that except for the example of transposed matrix multiplication \(\texttt{GEMM}^{\texttt{T}}\), we always achieve higher portability than TVM. The higher portability of TVM for \(\texttt{GEMM}^{\texttt{T}}\) is because TVM achieves for this example higher performance than our approach on NVIDIA Volta GPU. However, the higher performance of TVM is only due to the fact that TVM uses NVIDIA’s NVCC compiler for compiling CUDA code, while we currently rely on NVIDIA’s NVRTC library which surprisingly generates less-efficient CUDA assembly, as discussed above.

Listing 1 shows how matrix-vector multiplication (MatVec) is implemented in TVM’s high-level program representation which is embedded into the Python programming language. In line 1, the input size \((I,K)\in\mathbb{N}\times\mathbb{N}\) of matrix \(M\in T^{I\times K}\) (line 2) and vector \(v\in T^{K}\) (line 3) are declared, in the form of function parameters; the matrix and vector are named M and v and both are assumed to contain elements of scalar type \(T=\texttt{float32}\) (floating point numbers). Line 5 defines a so-called reduction axis in TVM, in which all values are combined in line 8 via te.sum (addition). The basic computation part of MatVec—multiplying matrix element M[i,k] with vector element v[k]—is also specified in line 8.

While we consider the MatVec implementations of TVM (Listing 1) and our approach (Figure 6) basically on the same level of abstraction, we consider our approach as more expressive in general. This is because our approach supports multiple reduction dimensions that may rely on different combine operators, e.g., as required for expressing the MBBS example in Figure 14. In contrast, TVM is struggling with different combine operators—adding support for multiple, different reduction dimensions is considered in the TVM community as a non-trivial extension of TVM [Apache TVM Community, 2020, 2022b]. Also, we consider our approach as slightly less error-prone: we automatically compute the expected sizes of matrix \(M\) (as \(I\times K\)) and vector \(v\) (as \(K\)), based on the user-defined input size \((I,K)\) in line 1 and index functions \((i,k)\mapsto(i,k)\) for the matrix and \((i,k)\mapsto(k)\) for the vector in line 8 (the formula for computing the sizes is described by Rasch [2024], Definition 8, for the interested reader). In contrast, TVM redundantly requests these matrix and vector sizes from the user: once in lines 2 and 3 of Listing 1, and again in lines 5 and 7. TVM uses these sizes for generating the function specification of its generated MatVec code, which lets TVM generate incorrect low-level code—without issuing an error message—when the user sets non-matching sizes in lines 2/3 and lines 5/7.

Figures 17–22 report the performance achieved by the PPCG-generated CUDA code for GPUs and of the OpenMP-annotated C code generated by polyhedral compiler Pluto for CPUs. For a fair comparison, we report for both polyhedral compilers their performance achieved for ATF-tuned tile sizes (denoted as PPCG+ATF/Pluto+ATF in the figures), as well as the performance of the two compilers when relying on their internal heuristics instead of auto-tuning (denoted as PPCG and Pluto). In some cases, PPCG’s heuristic crashed with error "too many resources requested for launch," because the heuristic seems to not take into account device-specific constraints, e.g., limited availability of GPUs’ fast memory resources.

We observe from Figures 17–22 that in all cases, our approach achieves better performance than PPCG and Pluto—sometimes by multiple orders of magnitude, in particular for deep learning computations (Figures 21 and 22). This is caused by the rigid optimization goals of PPCG and Pluto, e.g., always parallelizing outer dimensions, which causes severe performance losses. For example, we achieve a speedup over PPCG of \({{\gt}13\times}\) on NVIDIA Ampere GPU and of \({{\gt}60\times}\) over Pluto on Intel Skylake CPU for MCC as used in the inference phase of the real-world ResNet-50 neural network. Compared to PPCG, our better performance for this MCC example is because PPCG has difficulties with efficiently parallelizing computations relying on more than three dimensions. Most likely, this is because CUDA offers per default three dimensions for parallelization (called x, y, z dimension in CUDA). However, MCC relies on seven parallelizable dimensions (as shown in Figure 14), and exploiting the parallelization opportunities of the four further dimensions (as done in our generated CUDA code) is essential to achieve high performance for this MCC example from ResNet-50. Our performance advantage over Pluto for the MCC example is because Pluto parallelizes the outer dimensions of MCC only (whereas our approach has the potential to parallelize all dimensions); however, the dimension has a size of only 1 for this real-world example, resulting in starting only 1 thread in the Pluto-generated OpenMP code.

For dot products Dot (Figure 17), we can observe that PPCG fails to generate parallel CUDA code, because PPCG cannot parallelize and optimize computations which rely solely on combine operators different from concatenation, as we also discuss in Section 6.2. In Section 6.2, we particularly discuss that we do not consider the performance issues of PPCG and Pluto as weaknesses of the polyhedral approach in general, but of the particular polyhedral transformations chosen for PPCG and Pluto.

Note that Pluto crashes for our data mining example (Figure 20), with "Error extracting polyhedra from source file," because the scalar function of this example is too complex for Pluto (it contains if-statements). Moreover, Intel’s icx compiler struggles with compiling the Pluto-generated OpenMP code for quantum chemistry computations (Figure 19): we aborted icx’s compilation process after \(24\) h compilation time. The icx’s issue with the Pluto-generated code is most likely because of too aggressive loop unrolling of Pluto—the Pluto-generated OpenMP code has often a size \({\gt}50\) MB for our real-world quantum chemistry examples.

Since PPCG and Pluto are each designed for particular architectures only, they achieve the lowest portability of \(0\) for all our studies in Figure 23, according to the Pennycook metric. To simplify for PPCG and Pluto the portability comparison with our approach, we compute the Pennycook metric additionally also for two restricted sets of devices: only GPUs to make comparison against our approach easier for PPCG, and only CPUs to make comparison easier for Pluto.

Figures 24–28 report the portability of PPCG when considering only GPUs, as well as the portability of Pluto for only CPUs. We observe that we achieve higher portability for all our studies, as we constantly achieve higher performance than the two polyhedral compilers for the studies.

Note that even when restricting our set of devices to only GPUs for PPCG or only CPUs for Pluto, the two polyhedral compilers still achieve a portability of \(0\) for some examples, because they fail to generate code for them (as discussed above).

Listing 2 shows the input program of polyhedral compilers PPCG and Pluto for MatVec. Both take as input easy-to-implement, straightforward, sequential C code. We consider these two polyhedral compilers as more productive than our approach (as well as scheduling and functional approaches, and also polyhedral compilers that take DSL programs as input, such as TC [Vasilache et al., 2019]), because both compilers fully automatically generate optimized parallel code from unoptimized, sequential program code.

Rasch et al. [2020b,c] show that our approach can achieve the same high user productivity as polyhedral compilers, by using a polyhedral frontend for our approach: we can alternatively take as input the same sequential program code as PPCG and Pluto, instead of programs implemented in our high-level program representation (as in Figure 6). The sequential input program is then transformed via polyhedral tool pet [Verdoolaege and Grosser, 2012] to its polyhedral representation which is then automatically transformed to our high-level program representation, according to the methodology presented by Rasch et al. [2020b,c].

Already experimentally evaluated in previous work [Rasch et al., 2019a] and discussed in general terms in Section 6.3.

Listing 3 shows how MatVec is implemented in Lift. In line 1, type parameters n and m are declared, via the Lift building block nFun. Line 2 declares a function fun that takes as input a matrix of size \(\texttt{m}\times\texttt{n}\) and a vector of size n, both consisting of floating point numbers (float). The computation of MatVec is specified in lines 3 and 4. In line 3, Lift’s map pattern iterates over all rows of the matrix, and the zip pattern in line 4 combines each row pair-wise with the input vector. Afterward, multiplication * is applied to each pair, using Lift’s map pattern again, and the obtained products are finally combined via addition + using Lift’s reduce pattern.

Already for expressing MatVec, we can observe that Lift relies on a vast set of small, functional building blocks (five building blocks for MatVec: nFun, fun, map, zip, and reduce), and the blocks have to be composed and nested in complex ways for expressing computations. Consequently, we consider programming in Lift and Lift-like approaches as complex and their productivity for the user as limited. Moreover, the approaches often need fundamental extension for targeting new kinds of computations, e.g., so-called macro-rules which had to be added to Lift to efficiently target matrix multiplications [Remmelg et al., 2016] and primitives slide and pad together with optimization overlapped tiling for expressing stencil computations [Hagedorn et al., 2018]. This need for extensions limits the expressivity of the Lift language and thus further hinders productivity.

In contrast to Lift, our approach relies on exactly three higher-order functions (Figure 5) to express various kinds of data-parallel computations (Figure 14): (1) inp_view (Definition 7) which prepares the input data; our inp_view function is designed as general enough to subsume, in a structured way, the subset of all Lift patterns intended to change the view on input data, including patterns zip, pad, and slide; (2) md_hom (Definition 3) expresses the actual computation part, and it subsumes the Lift patterns performing actual computations (fun, map, reduce, \(\dotsc\)); (3) out_view (Definition 9) expresses the view on output data and is designed to work similarly as function inp_view (Lemma 2). Our three functions are always composed straightforwardly, in the same, fixed order (Figure 5), and they do not rely on complex function nesting for expressing computations.

Note that even though our language is designed as minimalistic, it should cover the expressivity of the Lift language²² and beyond: for example, we are currently not aware of any Lift program being able to express the prefix-sum examples in Figure 14. For the above reasons, we consider programming in our high-level language as more productive for the user than programming in Lift-like, functional-style languages. Furthermore, as discussed in Section 5.2, our approach can take as input also straightforward, sequential program code, which further contributes to the productivity of our approach.

Figures 17–22 report for completeness and also performance results achieved by domain-specific approaches. Since domain-specific approaches are specifically designed and optimized for particular application domains and often also architectures (e.g., only linear algebra routines on only GPU), we consider comparing to them as most challenging for us: our approach is designed and optimized for data-parallel computations in general, from arbitrary application domains (the same as also polyhedral compilers and many functional approaches), and our approach is also designed as generic in the target parallel architecture.

We observe in Figures 17–22 that the domain-specific libraries NVIDIA cuBLAS/cuDNN (for linear algebra routines and convolutions on GPUs) and Intel oneMKL/oneDNN (for linear algebra routines and convolutions on CPUs) sometimes perform better and sometimes worse than our approach.

The better performance of libraries over our approach is most likely²³ because the libraries internally rely on assembly-level optimizations, while we currently focus on the higher CUDA/OpenCL level of abstraction which offers less optimization opportunities [Goto and Geijn, 2008; Lai and Seznec, 2013]. The cuBLASEx extension of cuBLAS achieves in one case—MatMul on NVIDIA Volta GPU for square \(1024\times 1024\) input matrices—significantly higher performance than our approach. The high performance is achieved by cuBLASEx when using its CUBLAS_GEMM_ALGO1_TENSOR_OP algorithm variant, which casts the float-typed inputs implicitly to the half precision type (a.k.a. half or fp16), allowing cuBLASEx to exploit the GPU’s tensor core extension [NVIDIA, 2017]. Thereby, cuBLASEx achieves significantly higher performance than our approach, because tensor cores compute small matrix multiplication immediately in hardware; however, at the cost of a significant precision loss: the half scalar type achieves only half the accuracy achieved by scalar type float. When using cuBLASEx’s default algorithm CUBLAS_GEMM_DEFAULT (rather than algorithm CUBLAS_GEMM_ALGO1_TENSOR_OP), which retains the float type and thus meets the accuracy expected from the computation, we achieve a speedup of \(1.11\times\) over cuBLASEx.²⁴

The reason for the better performance of our approach over NVIDIA and Intel libraries is most likely because our approach allows generating code that is also optimized (auto-tuned) for data characteristics, which is important for high performance [Tillet and Cox, 2017]. In contrast, the vendor libraries usually rely on pre-implemented code that is optimized toward only average high performance for a range of data characteristics (size, memory layout, etc.). By relying on these fixed, pre-implemented code, the libraries avoid the auto-tuning overhead. However, auto-tuning is often amortized, particularly for deep learning computations—the main target of libraries NVIDIA cuDNN und Intel oneDNN—because the auto-tuned implementations are re-used in many program runs. Moreover, we achieve better performance for convolutions (Figure 18), because the libraries re-use optimizations for these computations originally intended for linear algebra routines [Li et al., 2016], whereas our optimization space (Table 1) is designed for data-parallel computations in general and not as specifically oriented toward linear algebra.

Compared to the EKR library (Figure 20), we achieve higher performance, because the EKR’s Java implementation inefficiently handles memory: the library is implemented using Java’s ArrayList data structure which is convenient to use for the Java programmer, but inefficient in terms of performance, because the structure internally performs costly memory re-allocations.

Similar to polyhedral compilers PPCG and Pluto, the domain-specific approaches work for certain architectures only and thus achieve the lowest portability of \(0\) only in Figure 23 for our studies. The domain-specific approaches are also restricted to a narrow set of studies, e.g., only linear algebra routines as NVIDIA cuBLAS and Intel oneMKL or only data mining example PRL as EKR. Consequently, the approaches achieve for these unsupported studies also a portability of only \(0\) in Figures 24–28 in which our portability evaluation is limited to only GPUs or CPUs, respectively, to make comparison against our approach easier for the vendor libraries.

For their target studies, domain-specific approaches achieve high portability. This is because the approaches are specifically designed and optimized toward these studies, e.g., via assembly-level optimizations which are currently beyond the scope of our work and considered as future work for our approach (see Section 8).

Listing 4 shows the implementation of MatVec in domain-specific approach NVIDIA cuBLAS; the implementation of MatVec in other domain-specific approaches, e.g., Intel oneMKL, is analogous to the implementation in Listing 4.

We consider domain-specific approaches as most productive for their target domain: in the case of MatVec, the user simply calls the high-level function cublasSgemv and passes to it the input matrices (omitted via ellipsis in the listing) together with some meta information (memory layout of matrices, etc); cuBLAS then automatically starts the GPU computation for MatVec.

Besides the fact that domain-specific approaches typically target only particular target architectures, a further fundamental productivity issue of domain-specific approaches is that they can only be used for a narrow class of computations only, e.g., only linear algebra routines as NVIDIA cuBLAS and Intel oneMKL. Moreover, in the case of domain-specific libraries from NVIDIA and Intel, it is often up to the user to manually choose among different, semantically equal but differently performing implementations for high performance. For example, the cuBLAS library offers three different routines for computing matrix multiplications: (1) cublasSgemm (part of standard cuBLAS), (2) cublasGemmEx (part of the cuBLASEx extension of cuBLAS), and (3) routine cublasLtMatmul (part of the cuBLASLt extension). These routines often also offer different, so-called algorithms (e.g., \(42\) algorithm variants in the case cuBLASEx) which impact the internal optimization process. When striving for the highest performance potentials of libraries, the user is in charge of naively testing each possible combination of routine and algorithm variant (as we have done in Figures 17–22 to make experimenting challenging for us). In addition, the user must be aware that different combinations of routines and algorithms can produce results of reduced accuracy (as discussed above), which can be critical for accuracy-sensitive use cases.

Scheduling approaches usually achieve high performance. For this, the approaches incorporate human expert knowledge into their optimization process which is based on two major steps: (1) a human expert implements an optimization program (a.k.a schedule) in a so-called scheduling language—the program specifies the basic optimizations to perform, such as tiling and parallelization; (2) an auto-tuning system (or a human hardware expert) chooses values of performance-critical parameter of the optimizations implemented in the schedule, e.g., particular values of tile sizes and concrete numbers of threads.

Our experiments in Section 5 show that compared to the state-of-the-art scheduling approach TVM (using its recent Ansor optimizer [Zheng et al., 2020a] for schedule generation), our approach achieves competitive and sometimes even better performance, e.g., speedups up to \(2.22\times\) on GPU and \(3.55\times\) on CPU over TVM+Ansor for computations taken from TVM’s favorable application domain (deep learning). Section 5 discusses that our better performance is due to the design and structure of our general optimization space (Table 1) which can be efficiently explored — fully automatically — using state-of-the-art auto tuning techniques [Rasch et al., 2021]. We focus on TVM in our experiments (rather than, e.g. Halide) to make experimenting challenging for us: TVM+Ansor has proved to achieve higher performance on GPUs and CPUs than popular state-of-practice approaches [Zheng et al., 2020a], including Halide, pyTorch [Paszke et al., 2019], and the recent FlexTensor optimizer [Zheng et al., 2020b].

Recent approach TensorIR [Feng et al., 2023] is a compiler for deep learning computations that achieves higher performance than TVM on NVIDIA GPUs. However, this performance gain over TVM is mainly achieved by exploiting the domain-specific tensor core [NVIDIA, 2017] extensions of NVIDIA GPUs, which compute in hardware the multiplications of small, low-precision \(4\times 4\) matrices. For this, TensorIR introduces the concept of blocks which represent sub-computations, e.g., multiplying \(4\times 4\) matrices. These blocks are than mapped by TensorIR to domain-specific hardware extensions, which often leads to high performance.

While domain-specific hardware extensions are not targeted in this article, we can naturally exploit them in our approach, similar to TensorIR, as we plan for our future work: the sub-computations targeted by the current hardware extensions, such as matrix multiplication on \(4\times 4\) matrices, can be straightforwardly expressed in our approach (Figure 14). Thus, we can match these sub-expressions in our low-level representation and map them to hardware extensions in our generated code. For this, instead of relying on a full partitioning in our low-level representation (as in Figure 15) such that we can apply scalar function \(f\) to the fully de-composed data (consisting of a single scalar value only in the case of a full partitioning), we plan to rely on a coarser-grained partitioning schema, e.g., down to only \(4\times 4\) matrices (rather than \(1\times 1\) matrices, as in the case of a full partitioning). This allows us replacing scalar function \(f\) (which in the case of matrix multiplication is a simple scalar multiplication \(*\)) with the operation supported by the hardware extension, such as matrix multiplication on \(4\times 4\) matrices. We expect for our future work to achieve the same advantages over TensorIR as over TVM, because apart from supporting domain-specific hardware extensions, TensorIR is very similar to TVM.

While scheduling approaches achieve high performance, they tend to struggle with achieving portability. This is because even though the approaches often offer different, pre-implemented backends (e.g., a CUDA backend to target NVIDIA GPUs and an OpenCL backend for CPUs), they do not propose any structured methodology about how new backends can be added, e.g., for potentially upcoming architectures, with potentially deeper memory and core hierarchies than GPUs and CPUs. This might be particularly critical (or requiring significant development effort) for the application area of deep learning which is the main target of many scheduling approaches, e.g., TVM and TensorIR, and for which new architectures are arising continuously [Hennessy and Patterson, 2019].

In contrast, we introduce in this article a formally precise recipe for correct-by-construction code generation in different backends (including OpenMP, CUDA, and OpenCL), generically in the target architecture: we introduce an architecture-agnostic low-level representation (Section 3) as target for our high-level programs (Section 2), and we describe formally how our high-level programs are automatically lowered to our low-level representation (Section 4), based on the architecture-agnostic optimization space in Table 1. Rasch [2024] (Section E) outlines how executable, imperative-style program code is straightforwardly generated from low-level expressions, which we plan to discuss and illustrate in detail in our future work.

Scheduling approaches rely on a two-step optimization process, as discussed above: implementing a schedule (first step) and choosing optimized values of performance-critical parameters within that schedule (second step). While the second step often can be easily automatized, e.g., via auto-tuning [Chen et al., 2018b], the first step—implementing a schedule—usually has to be conducted manually by the user to achieve high performance, which requires expert knowledge and thus hinders productivity. The lack of formal foundation of many scheduling approaches further complicates implementing schedules for the user, as implementation becomes error prone and hardly predictable. For example, Fireiron’s schedules can achieve high performance, close to GPUs’ peak, but schedules in Fireiron can easily generate incorrect low-level code: Fireiron cannot guarantee that optimizations expressed in its scheduling language are semantics preserving, e.g., based on a formal foundation as done in this work, making programming Fireiron’s schedules error prone and complex for the user. Similarly, TVM is sometimes unable to detect user errors in both its high-level language (as discussed in Section 5.1) as well as scheduling language [Apache TVM Community, 2022e]. Safety in parallel programming is an ongoing major demand, in particular from industry [Khronos, 2022a].

Auto schedulers, such as Halide’s optimization engine [Mullapudi et al., 2016] and TVM’s recent Ansor [Zheng et al., 2020a], aim to automatically generate well-performing, correct schedules for the user. However, a major flaw of the current auto schedulers is that even though they work well for some computations (e.g., from deep learning, as TVM’s Ansor), they may perform worse for others. For example, our approach achieves a speedup over TVM+Ansor of \({\gt}100\times\) already for straightforward dot products (Figure 17). This is because Ansor does not exploit multiple thread blocks and uses only a small number of threads for reduction computations. While such optimization decisions are often beneficial for reductions as used in deep learning (e.g., within the computations of convolutions and matrix multiplications on deep learning workloads, because parallelization can be better exploited for outer loops of these computations), these rigid optimization decisions of Ansor may perform worse in other contexts (e.g., for computing dot product).

To avoid the productivity issues of scheduling approaches, we have designed our optimization process as fully auto-tunable, thereby freeing the user from the burden and complexity of making complex optimization decisions. Our optimization space (Table 1) is designed as generic in the target application area and hardware architecture, thereby achieving high performance for various combinations of data-parallel computations and architectures (Section 5). Correctness of optimizations is ensured in our approach by introducing a formal foundation that enables mathematical reasoning about correctness. Particularly, our optimization process is designed as correct-by-construction, meaning that any valid optimization decisions (i.e., a particular choice of tuning parameters in Table 1 that satisfy the constraints) leads to a correct expression in our low-level expression (as in Figure 15). In contrast, approaches such as introduced by Clément and Cohen [2022] formally validate optimization decisions of scheduling approaches in already generated low-level code. Thereby, such approaches work potentially for arbitrary scheduling approaches (Halide, TVM, \(\dotsc\)), but the approaches cannot save the user at the high abstraction level from implementing incorrect optimizations (e.g., via easy-to-understand, high-level error messages indicating that an invalid optimization decisions is made) or restricting the optimization space otherwise to valid decisions only, e.g., for an efficient auto-tuning process, because the approaches check already generated program code.

Scheduling approaches often also suffer from expressivity issues. For example, Fireiron is limited to computing only matrix multiplications on only NVIDIA GPUs, and TVM does not support computations that rely on multiple combine operators different from concatenation [Apache TVM Community, 2020, 2022b], e.g., as required for expressing the MBBS example in Figure 14. Also, TVM has difficulties with user-defined combine operators [Apache TVM Community, 2022d] and thus crashes for example PRL in Figure 14. In contrast to TVM, we introduce a formal methodology about how to manage different kinds of arbitrary, user-defined combine operators (Section 3), which is considered challenging [Apache TVM Community, 2020].

Polyhedral compilers tend to struggle with achieving their full performance potential. We argue that this performance issue of polyhedral compilers is mainly caused by the following two major reasons.

While we consider the set of polyhedral transformation (so-called affine transformation) as broad, expressive, and powerful, each polyhedral compiler implements a subset of expert-chosen transformations. This subset of transformations, as well as the application order of transformations, are usually fixed in a particular polyhedral compiler and chosen toward specific optimization goals only, e.g., coarse-grained parallelization and locality-aware data accesses (a.k.a. Pluto algorithm [Bondhugula et al., 2008a]), causing the search spaces of polyhedral compilers to be a proper subset of our space in Table 1. Consequently, computations that require for high performance other subsets of polyhedral transformations and/or application orders of transformations (e.g., transformations toward fine-grained parallelization) might not achieve their full performance potential when compiled with a particular polyhedral compiler [Consolaro et al., 2024].

In contrast to the currently existing polyhedral compilers, we have designed our optimization process as generic in goals: for example, our space is designed such that the degree of parallelization (coarse, fine, \(\dotsc\)) is fully auto-tunable for the particular combination of target architecture and computation to optimize. We consider it as an interesting future work to investigate the strength and weaknesses of the polyhedral model for expressing our generic optimization space.

We see the second reason for potential performance issues in polyhedral compilers in their difficulties with reduction-like computations. This is mainly caused by the fact that the polyhedral model captures less semantic information than the high-level program representation introduced in Section 2 of this article: combine operators which are used to combine the intermediate results of computations (e.g., operator \(+\) from Example 2 for combining the intermediate results of the dot products within matrix multiplication) are not explicitly represented in the polyhedral model; the polyhedral model is rather focused on modeling memory accesses and their relative order only. Most likely, these semantic information are missing in the polyhedral model, because polyhedral approaches were originally intended to fully automatically optimize loop-based, sequential code (such as Pluto and PPCG)—extracting combine operators automatically from sequential code is challenging and often even impossible (Rice’s theorem).

In contrast, our proposed high-level representation explicitly captures combine operators (Figure 14), by requesting these operators explicitly from the user. This is important, because the operators are often required for generating code that fully utilizes the highly parallel hardware of state-of-the-art architectures (GPUs, etc.), as discussed in Section 5. Similarly to our approach, polyhedral compiler TC also requests combine operators explicitly from the user. However, TC is restricted to operators + (addition), * (multiplication), min (minimum), and max (maximum) only, thereby TC is not able to express important examples in Figure 14, e.g., PRL which is popular in data mining. Moreover, TC outsources the computation of its combine operators to the NVIDIA CUB library [NVIDIA, 2022a]; most likely as a workaround, because TC relies on the polyhedral model which is not designed to capture and exploit semantic information about combine operators for optimization. Thereby, TC is dependent on external approaches for computing combine operators, which might not always be available (e.g., for upcoming architectures).

Workarounds have been proposed by the polyhedral community to target reduction-like computations [Doerfert et al., 2015; Reddy et al., 2016]. However, these approaches are limited to a subset of computations, e.g., by not supporting user-defined scalar types [Doerfert et al., 2015] (as required for our PRL example in Figure 14), or by being limited to GPUs only [Reddy et al., 2016]. Comparing the semantic information captured in the polyhedral model vs. our MDH-based representation have been the focus of discussions between polyhedral experts and MDH developers [Google SIG MLIR Open Design Meeting, 2020].

The polyhedral approach, in its general form, is a framework offering transformation rules (affine transformations), and each individual polyhedral compiler implements a set of such transformations which are then instantiated (e.g., with particular tile sizes) and applied when compiling a particular application. However, individual polyhedral compilers (e.g., PPCG and Pluto) apply a fixed set of affine transformations, thereby rigidly optimizing for a particular target architecture only, e.g., only GPU (as PPCG) or only CPU (as Pluto), and it remains open which affine transformations have to be used and how for other architectures, e.g., upcoming accelerators for deep learning computations [Hennessy and Patterson, 2019] with potentially more complex memory and core hierarchies than GPUs and CPUs. Moreover, while we introduce an explicit low-level representation (Section 3), the polyhedral approach does not introduce representations on different abstraction levels: the model relies on one representation that is transformed via affine transformations. Apart from the ability of our low-level representation to handle combine operators (which we consider as complex and important), we see the advantages of our explicit low-level representation in, for example, explicitly representing memory regions, which allows formally defining important correctness constraints, e.g., that GPU architectures allow combining the results of threads in designated memory regions only. Furthermore, our low-level representation also allows straightforwardly generating executable code from it (shown by Rasch [2024], Section E, and planned to be discussed thoroughly in future work). In contrast, code generation from the polyhedral model has proven challenging [Bastoul et al., 2022; Vasilache et al., 2022; Grosser et al., 2015].

Most polyhedral compilers achieve high user productivity, by fully automatically parallelizing and optimizing straightforward sequential code (as Pluto and PPCG). Our approach currently relies on a Domain-Specific Language (DSL) for expressing computations, as discussed in Section 2; thus, our approach can be considered as less productive than many polyhedral compilers. However, Rasch et al. [2020b, c] show that DSL programs in our approach can be automatically generated from sequential code (optionally annotated with simple, OpenMP-like directives for expressing combine operators, enabling advanced optimizations), by using polyhedral tool pet [Verdoolaege and Grosser, 2012] as a frontend for our approach. Thereby, we are able to achieve the same, high user productivity as polyhedral compilers. We consider this direction—combing the polyhedral model with our approach—as promising, as it enables benefitting from the advantages of both directions: optimizing sequential programs and making them parallelizable using polyhedral techniques (like loop skewing, as also outlined in Section 6.5), and mapping the optimized and parallelizable code eventually to parallel architectures based on the concepts and methodologies introduced in this article.

Functional approaches tend to struggle with achieving their full performance potential, often caused by the design of their optimization spaces. For example, analogously to our approach, functional approach Lift relies on an internal low-level representation [Steuwer et al., 2017] that is used as target for Lift’s high-level programs. However, Lift’s transformation process, from high level to low level, turned out to be challenging: Lift’s lowering process relies on an infinitely large optimization space—identifying a well-performing configuration within that space is too complex to be done automatically in general, due to the space’s large and complex structure. As a workaround, Lift currently uses approach Elevate [Hagedorn et al., 2020b] to incorporate user knowledge into the optimization process; however, at the cost of productivity, as manually expressing optimization is challenging, particularly for non-expert users.

In contrast, our optimization process is designed as auto-tunable (Table 1), thereby achieving fully automatically high performance, as confirmed in our experiments (Section 5), without involving the user for optimization decisions. In particular, our previous work already showed that our approach—even in its original, proof-of-concept implementation [Rasch et al., 2019a]—can significantly outperform Lift on GPU and CPU [Rasch et al., 2019a]. Our performance advantage over Lift is mainly caused by the design of our optimization process: relying on formally defined tuning parameters (Table 1), rather than on formal transformation rules that span a too large and complex search space (as Lift), thereby contributing to a simpler, fully auto-tunable optimization process.

The current functional approaches usually are designed and optimized toward code generation in a particular programming model only. For example, Lift inherently relies on the OpenCL programming model, because OpenCL works for multiple kinds of architectures: NVIDIA GPU, Intel CPU, and so on. However, we see two major disadvantages in addressing the portability issue via OpenCL only: (1) GPU-specific optimizations (such as shuffle operations [NVIDIA, 2018]) are available only in the CUDA programming model, but not in OpenCL; (2) the set of OpenCL-compatible devices is broad but still limited; in particular, in the new golden age for computer architectures [Hennessy and Patterson, 2019], upcoming architectures are arising continuously and may not support the OpenCL standard. We consider targeting new programming models as challenging for Lift, as its formal low-level representation is inherently designed for OpenCL [Steuwer et al., 2017]; targeting further programming models with Lift would require the design and implementation of new low-level representations, which we do not consider as straightforward.

To allow easily targeting new programming models with our approach, we have designed our formalism as generic in the target model: our low-level representation (Figure 15) and optimization space (Table 1) are designed and optimized toward an abstract system model (Definition 10) which is capable of representing the device models of important programming approaches, including OpenMP, CUDA, and OpenCL (Example 11). Furthermore, we have designed our high- and low-level representations as minimalistic (Figures 6 and 15), e.g., by relying on three higher-order functions only for expressing programs at the high abstraction level, which simplifies and reduces the development effort for implementing code generators for programming models.

In addition, we believe that compared to our approach, the following basic design decisions of Lift (and similar functional approaches) complicate the process of code generation for them and increase the development effort for implementing code generators: (1) relying on a vast set of small patterns for expressing computations, rather than aiming at a minimalistic design as we do (as also discussed in Section 5.3); (2) relying on complex function nestings and compositions for expressing computations, rather than avoiding nesting and relying on a fixed composition structure of functions, as in our approach (Figure 5); (3) requiring new patterns for targeting new classes of data-parallel computations (such as patterns slide and pad for stencils [Hagedorn et al., 2018]), which have to be non-trivially integrated into Lift’s type and optimization system (often via extensions of the systems [Hagedorn et al., 2018, Remmelg et al., 2016]), instead of relying on a fixed set of expressive patterns (Figure 6) and generalized optimizations (Table 1) that work for various kinds of data-parallel computations (Figure 14); (4) expressing high-level and low-level concepts in the same language, instead of separating high-level and low-level concepts for a more structured and thus simpler code generation process (Figure 4). We consider these four design decisions as disadvantageous for code generation, because they require from a code generator handling various kinds of patterns (decision 1), and the patterns need to be translated to significantly different code variants, depending on their nesting level and composition order (decision 2). Moreover, each extension of patterns (decision 3) might affect code generation also for the already supported patterns, because the existing patterns need to be combined with the new ones via composition and nesting (decision 2). We consider mixing up high-level and low-level concepts in the same language (decision 4) as further complicating the code generation process, because code generators cannot be implemented in clear, distinct stages: high-level language\(\rightarrow\)low-level language\(\rightarrow\)executable programcode.

Functional approaches are expressive frameworks—to the best of our knowledge, the majority of these approaches should also be able to express (possibly after some extension) many of the high-level programs that can also be expressed via our high-level representation (e.g., those presented in Figure 14).

A main difference we see between the high-level representations of existing functional approaches and the representation introduced by our approach is that the existing approaches rely on a vast set of higher-order functions for expressing computations; these functions have to be functionally composed and nested in complex ways for expressing computations. For example, expressing matrix multiplication in Lift requires also involving Lift’s pattern transpose (also when operating on non-transposed input matrices) [Remmelg et al., 2016], as per design in Lift, multi-dimensional data is considered as an array of arrays (rather than a MDA, as in our approach as well as polyhedral approaches). In contrast, we aim to keep our high-level language minimalistic, by expressing data-parallel computations using exactly three higher-order functions and which are always used in the same, fixed order (shown in Figure 5). Rasch et al. [2020b, c] confirm that due to the minimalistic and structured design of our high-level representation, programs in our representation can even be systematically generated from straightforward, sequential program code.

Functional approaches also tend to require extension when targeting new application areas, which hinders the expressivity of the frameworks and thus also their productivity. For example, functional approach Lift [Steuwer et al., 2015] required notable extension for targeting, e.g., matrix multiplications (so-called macro-rules had to be added to Lift [Remmelg et al., 2016]) and stencil computations (primitives slide and pad were added, and Lift’s tiling optimization had to be extended toward overlapped tiling [Hagedorn et al., 2018]). In contrast, we have formally defined our class of targeted computations (as MDH functions, Definition 3), and the generality of our approach allows expressing matrix multiplications and stencils out of the box, without relying on domain-specific building blocks.

Domain-specific approaches, such as cuBLAS and cuDNN, usually achieve high performance. This is because the approaches are hand-optimized by performance experts—on the assembly level—to exploit the full performance potential of their target architecture. In our experiments (Section 5), we show that our approach often achieves competitive and sometimes even better performance than domain-specific approaches provided by NVIDIA and Intel, which is mainly caused by their portability issues across different data characteristics, as we discuss in the next paragraph.

Domain-specific approaches usually struggle with achieving portability across different architectures. This is because the approaches are often implemented in architecture-specific assembly code to achieve high performance, but thereby also being limited to their target architecture. The domain-specific approaches often also struggle with achieving portability of performance across different data characteristics (e.g., their sizes): the approaches usually rely on a set of pre-implemented implementations that are each designed and optimized toward average high performance across a range of data characteristic. In contrast, our approach (as well as many scheduling and polyhedral approaches) allow automatically optimizing (auto-tuning) computations for particular data characteristics, which is important for achieving high performance [Tillet and Cox, 2017]. Thereby, our approach often outperforms domain-specific approaches (as confirmed in Section 5), particularly for advanced data characteristics (small, uneven, irregularly shaped,\(\dots\)), e.g., as used in deep learning. The costly time for auto-tuning is well amortized in many application areas, because the auto-tuned implementations are re-used in many program runs. Furthermore, auto-tuning avoids the time-intensive and costly process of hand-optimization by human experts.

Domain-specific approaches usually achieve highest productivity for their target domain (e.g., linear algebra), by providing easy to use high-level abstractions. However, the approaches suffer from significant expressivity issues, because—per design—they are inherently restricted to their target application domain only. Also, the approaches are often inherently bound to only particular architectures, e.g., only GPU (as NVIDIA cuBLAS and cuDNN) or only CPU (as Intel oneMKL and oneDNN). Domain-specific vendor libraries, such as NVIDIA cuBLAS and Intel oneMKL, also tend to offer the user differently performing variants of computations; the variants have to be naively tested by the user when striving for the full performance potentials of approaches (as discussed in Section 5.4), which is cumbersome for the user.