Mat2Boundary: Treating User-Defined Boundary Condition as SpMV for Distributed PDE Solvers on Block-Structured Grids Yanzheng Cai, Mingzhe Zhang, Shengqi Chen, Haoyuan Song, Wenguang Chen∗
arXiv:2605.14780v1 [cs.PL] 14 May 2026
Department of Computer Science and Technology & BRNist, Tsinghua University, Beijing, China. {cyz22,zmz21,csq20,song-hy22}@mails.tsinghua.edu.cn, [email protected]
Abstract—Boundary-condition (BC) handling is a major source of complexity in PDE solvers on structured and block-structured grids, especially for high-order methods and distributed-memory execution. We present M AT 2B OUNDARY, a DSL and compiler for boundary computations that models a broad class of BCs as affine sparse linear operators. This abstraction unifies halo copying, circular and symmetric mappings, zero padding, block-edge synchronization, and user-defined interpolation, while exposing a modular basic sub-matrix interface for declarative composition. To make this representation efficient, M AT 2B OUNDARY combines multi-stage programming and polyhedral analysis to generate matrix-free kernels for structured cases, support user-defined sparse matrices for irregular cases, eliminate redundant boundary work, and synthesize reusable communication schedules for distributed execution. Evaluated on two shallow-water equation solvers on cubed-sphere grids and HPCG, M AT 2B OUNDARY achieves up to 7.6× BC-kernel speedup, reduces BC code by over 70%, and scales to 1,344 CPU cores with 72%-88% efficiency. Index Terms—domain-specific language, finite difference method, finite volume method, stencil, block-structured grid, multi-stage programming, polyhedral compilation
I. I NTRODUCTION Numerical solutions of partial differential equations (PDEs) are of great importance in a wide range of scientific fields, including seismic imaging [1], weather prediction [2], etc. However, manually developing PDE solvers using generalpurpose programming languages (e.g., FORTRAN, C/C++) imposes a substantial programming burden, requiring the implementation of thousands to tens of thousands of lines of code. Many of these PDE solvers employ structured grids, primarily Cartesian grids, for spatial discretizations. Due to the statically determined neighboring relationships of the grid points, numerous programming frameworks and domainspecific languages (DSLs) [3], [4], [5], [6], [7] offer simplified programming interfaces and dedicated optimization passes for applications on structured grids Block-structured grids, which are the concatenations of several structured blocks, are powerful when processing complex and real geometries by offering uniform or quasi-uniform grid cells. They are widely used in real world applications, such as the dynamical core of weather and climate models [8]. Figure 1 demonstrates the cubed-sphere grid for discretizations of the sphere geometry, and we can see obvious discontinuity over block boundaries. While maintaining the computing efficiency and accuracy as structured grids in the center region of each block, ∗ Corresponding author.
(a) A 3D view of the cubedsphere tiling with resolution of 15 × 15.
(b) A closeup view showing the overlap of grid lines from the upper block ghost cells on the neighboring blocks.
Fig. 1: Geometry of the cubed-sphere grid.
block-structured grids involve additional boundary condition (BC) computations. Simple and naive implementations of the boundary algorithm can introduce degradation on the order of accuracy and grid imprinting [9]. Therefore, researchers have proposed several high-order interpolation schemes for the boundary calculation of cubed-sphere grids [10], [11], [12], [13]. The complexity in the programming of these boundary algorithms is reflected in the length of the code. As shown in Table I, different boundary processing logics are of 1.12-9.87x in code length compared to the corresponding inner calculation kernels, and the total BC lines of code (LOCs) account for over 40% of the whole standalone application. The ratio will increase sharply when we apply program representations with much more simplicity to inner calculations on structured grids [5], [6], [3], which can decrease their code length to less than 20% of origin [3]. Furthermore, process-to-process communications should also be considered when we distribute these standalone applications to large clusters, making the boundary conditions the most lengthy and tedious part of the application programs. Previous frameworks [4], [14] and applications [12] have tried to give solutions to different instances of boundary algorithms, but they fail to balance versatility, programmability and scalability. To address these challenges, we introduce a novel DSL called M AT 2B OUNDARY, designed for processing userdefined BCs on distributed structured and block-structured grids. Our contributions are as follows:
Method
HOPE [12]
MCV [10]
Operator
BC type
BC interp./quadrature dimensions
BC LOCs
Reconstruction
iterative
2d/2d
385
39
Border Flux
direct
-/-
178
159
Calc LOCs
Other
-
-/-
-
542
Flux Derivative
direct
1d/-
203
49
Result Averaging
direct
-/ -
168
23
Other
-
-/-
-
317
TABLE I: Statistics of coding length for implementing two standalone PDE solvers on cubed-sphere grids using FORTRAN
We formalize a wide spectrum of boundary algorithms—from simple padding to complex user-defined interpolations—under a generalized SpMV abstraction, enabling modular construction via programmable basic sub-matrices. • By integrating multi-stage programming with polyhedral analysis, our compiler automatically resolves boundary dependencies at compile-time, eliminates redundant iteration spaces, and generates highly optimized, matrix-free kernels for structured sub-regions. • For large-scale execution, we introduce a domain-specific communication backend for low-overhead remote data fetching. IR-reusing techniques are also introduced to achieve moderate compilation time and output code length. • Integrated with a standalone dsl which describes the computation on structured grids [3], [15], our dsl is capable of expressing a much wider range of boundary conditions which appeared in real applications, including different high order finite-volume method (FVM) solvers for solving shallow-water equation (SWE) on cubed-sphere grids and High Performance Conjugate Gradients (HPCG) on structural grids. For BC computations, we achieve a performance gain of 0.93-7.64× compared to original FORTRAN/C++ implementations with only 30% LOCs for different SWE solvers. By changing several lines of code, we can distribute these application on multiple processes, reaching an average 4.75× speed-up for standalone execution compared with original implementations, and a strong scaling efficiency of 72%-88% on a cluster with 1,344 CPU cores. •
(a) Padding Zero
(c) Pure Function
Figure 2 demonstrates common boundary algorithms implemented in real applications, ranging from convolutions on images to PDE solvers on structured or block-structured grids. The algorithms shown in Figures 2a, 2b, and 2d are widely adopted by popular frameworks and applications even beyond the domain of PDE solver. Meanwhile, the algorithms depicted in Figures 2c, 2e, 2f, and 2g are also commonly utilized in different PDE solvers, but the detailed representation and parameters vary significantly through different governing equations, grids and discretization methods. Notably, the degree of freedom for boundary algorithms comes from not only the increase in basic types, but also the flexible combination of all these types. For example, users may introduce the circular
(f) User-Defined Quadrature
(d) Circular/Symmetric
II. BACKGROUND AND M OTIVATION A. Survey on Existing Boundary Algorithms
(b) Halo Copying
(g) User-Defined Interpolation (e) Edge Averaging
Fig. 2: Common BCs for block-structured grids
type to x-direction ghost cells and the padding zero type to y-direction ghost cells when modeling a cylinder geometry. B. User-Defined Interpolation Example We take the finite-volume SWE solver HOPE [12] as an example to demonstrate complexity of general user-defined BCs on block-structured grids. In order to deal with the discontinuity of different parts of the cubed-sphere grid, for each element in global ghost region, users select one or several reconstruction points in the element’s corresponding area,
and view their numeric quadrature as the value for the whole element. Figure 2f shows the Gaussian quadrature points that the HOPE method uses for the reconstruction of ghost cells on the left block. For calculating the reconstruction points, users view them from other parts of the grid (the right block), translate their position through coordinate transformation between different parts, and carry out another interpolation progress from the data of the right part. Despite the difference in dependencies and weights for the interpolation and quadrature procedures, the main complexity of HOPE lies in the use of the ghost cell values in the right part for calculating reconstruction points derived from the left part (shown in Figure 2g), while other solvers [10], [11], [13] only use data cells. This configuration escalates the difficulty of finding exact solutions to the same level as solving a system of linear equations. In practice, HOPE first initializes all of the ghost cells to zero and carries out a fixed number of iterations (10 times) to obtain a fairly approximate solution, which is equivalent to solving the linear system using Jacobi method during each calculation. Furthermore, HOPE proves that the whole iterating process is still a linear operator on normal data cells to calculate ghost cells, therefore it calculate the weights in advance and stores them down in Compressed-Sparse-Row (CSR) format, therefore the whole iterating process can be expressed as Sparse Matrix-Vector multiplication (SpMV). We reference these two algorithms as HOPE(I) (Iterative) and HOPE(M) (Matrix) in later description.
3) Efficiency: Users may declare a complex boundary region (ring-like ghost region) that share the same computation pattern, therefore it is the compiler’s duty to generate code for traversing these regions efficiently. What is more, there exists two specific optimizing opportunities for BC computation • The boundary conditions may consists of simple types which do not need to iterate over a real sparse matrix (pure function, symmetric and circular). Therefore, generate matrix-free code for these BCs can reduce memory cost while improve execution efficiency. • We notice that BC computations do not appear alone. They are always followed by inner calculations. Additionally, the assignments of BC computations are only used inside its corresponding operator. Therefore, non-equivalent transformations could be taken into account for BC computations as long as they do not affect the inner calculation process. This feature can benefit both the configuration simplicity and the execution efficiency. 4) Scalability: For a framework or DSL, scalability is determined by both its execution efficiency and its ease of programming at scale. A major challenge lies in designing abstractions that allow developers to scale their applications with only minor additions with the underlying data distribution configurations. Achieving scalability also requires a highly capable communication backend — one that can perform rapid halo exchanges while retaining the flexibility to handle irregular data retrieval for unstructured boundary conditions.
C. Opportunities and Challenges
III. S P MV R EPRESENTATION OF B OUNDARY C ONDITIONS
In this section, we demonstrate our intuitive ideas for supporting all kinds of boundary conditions. Ideas and their corresponding opportunities and challenges are summarized in four aspects as below: 1) Universality: Combination of different BCs are necessary when dealing with real-world complicated geometry. Therefore, an unified representation for all types of boundary algorithms will boost programmer’s productivity on processing complicated boundary problems. The sparse matrix is a good abstraction to represent BC and has been used in HOPE-matrix method. However, its versatility among general finite-volume methods remains unexplored. 2) Simplicity: As mentioned in Table I, the workload of writing BC code constitutes a significant or even the main portion of the entire programming workload, therefore the simplicity of boundary configuration is vital for productivity. A declarative configuration style is desired, where users only designate the boundary regions and their corresponding computation patterns. The distributed implementation [16] of Devito DSL [5] meet the requirements, but it only support simple halo-copying. The new Python implementation of FV3 dynamic core [14] uses a cubed-sphere specific communication library for the configuration. But when other grids, discretization methods or distribution schemes are introduced, users have to work on the low-level MPI-based library code, which will block the whole development progress.
In this section, we formalize the representation of different boundary condition algorithms and demonstrate their interconnection with sparse matrix-vector multiplication (SpMV). A. Definitions Shown in [12], [11], [10], and in Figure 2f and Figure 2g, the user-defined boundary condition algorithms compute the values on the extended grid cells of each panel in the cubedsphere grid. Therefore, we structure the full region of a M AT 2B OUNDARY’s grid variable y into 2 distinct regions: data and y ghost. We denote the coordinate sets of data region as RD and y the coordinate sets of ghost region as RG . Then we can use a vector of the size |RFy | to represent the output value of a usery y defined boundary condition algorithm, where RFy = RD ∪RG . When in the context of distributed memory parallelism (DMP), each process owns a sub-part of the whole data region, denoted as Rdyp for process p and variable y. We have the y relationship that ∪p Rdyp = RD . The process p also maintains the data in the local halo region Rhy p for later calculation, and ∪p (Rfyp ) = RFy , where Rfyp = Rdyp ∪ Rhy p . Using all these definitions, we can formulate a distributed user-defined boundary condition algorithms as computing y(|Rfy |) = F (x1 , x2 , ...) p
in each process, where the input parameter xi has a global xi data region of RD .
B. Representation A common and powerful paradigm of computing y using x is matrix-vector multiplication, and we list below how to implement all kinds of common boundary conditions using this paradigm: x |) × x(|Rx |) + b(|Ry |) F (x) = A(|Rfy |,|RD D f p
p
Padding with Zeroes in Figure 2a is the simplest boundary condition for the sub-part of the ghost region that process p locally holds. y x For index i ∈ Rhy p ∩ RG and j ∈ RD , Aij = bi = 0 • Halo Copying in Figure 2b is always introduced by distributed memory parallelism and is fully researched and optimized by manual implementations and dsl implementations. In this case, Ry = Rx . For index i ∈ Rhp ∩ RD and j ∈ RD , Aij = δij , bi = 0. • Circular and Symmetric Boundaries in Figure 2d also demand Ry = Rx . For index i ∈ Rhp ∩ RG and j ∈ RD , •
Aij = Boolean(j = f map(i)), bi = 0 where f map(i) is a simple index mapping. When in the one-dimensional case where RD = {0, 1, ..., |RD | − 1}, f map for the left ghost region has the form ( x−1 symmetric , 0 < x < |RD | f map(−x) = |RD | − x circular The Pure Function Boundary in Figure 2c is common in testing the performance of PDE solvers on equations with analytical solution. For example, building a boundary with Dirichlet Conditions. y x , Aij = 0, bi = udf (i) For index i ∈ Rhy p ∩ RG and j ∈ RD • Block Edge Averaging in Figure 2e occurs when there are several grid cells located on different parts of the blockstructured grid that share the same physical location [10], therefore the condition Ry = Rx holds. We pick all these specific cells on the block edges to form a set RSD , Therefore the relationships between these cells can be viewed as a sub-matrix Sync|RSD |,|RSD | . Sync has a great sparsity. Normally each row of the Sync has only 2 or 3 non-zero elements. For index i ∈ Rdp and j ∈ RD , ( Syncij j ∈ RSD Aij = , bi = 0 0, j∈ / RSD •
quadrature integration scheme can be viewed as another linear mapping Agg|RG |,|P ts| for calculating ghost cell values. Multiplying Agg and Recons together, they obtain the final linear mapping as the matrix Itp|RG |,|RD | . We find out that the simplification process of HOPE is much more general than the method itself. As the set P ts does not affect the final matrix Itp, we claim that the different methods, order of accuracy and reconstruction points P ts only affect the detailed value of Itp. Therefore, by asking users to pre-process the matrix Itp, for index i ∈ Rhp ∩ RD and j ∈ RD , M AT 2B OUNDARY represents the boundary as Aij = Itpij , bi = 0 In summary, we illustrate six types of common boundary conditions in block-structured grids and show their implementation with the abstraction of SpMV. Notably, the construction of these matrices heavily relies on conditional branching and case-by-case analysis, and the final result of matrix A is of great sparsity. This motivates our design of the basic submatrix abstraction which is described in detail in Section IV. IV. M AT 2B OUNDARY D ESIGN Based on the observation mentioned above, we propose M AT 2B OUNDARY, our domain-specific language and its corresponding compiler and runtime library for solving boundary problems on distributed block-structured girds. We demonstrate our design using the 3 blocked 2-D Cartesian grid example shown in Figure 3. In contrast with block 1, block 0 has a similar axis order but its shape is trapezoidal, therefore the left ghost region of block 1 need an user-defined interpolation BC. The geometry of block 2 is the same with block 1, but its axes are rotated, therefore the right ghost region BC of block 1 is in the form of a simple mapping. The top and bottom ghost region BC is set to zero. The data region size is set to 5 × 5 and the full region size is set to 7 × 7.
j i
•
The User Defined Interpolation Boundary in Figure 2f and Figure 2g appears in realistic applications of PDE solvers on block-structured grids [12], [10], [11]. The authors of HOPE [12] prove in their article’s Appendix that their boundary algorithm is equivalent to multiplying a matrix of shape (|RG |, |RD |) on variable x to get ghost points on y, where Ry = Rx . They reach this conclusion by defining the set of reconstruction points as P ts and figuring out that their Tensor product polynomial (TPP) reconstruction scheme is equivalent as a linear mapping Recons|P ts|,|RD | for calculating reconstruction points, and their Gaussian
j
Block 0
i j
i
Block 1
Block 2
Fig. 3: Example Block-Structured Grid A. Modular Boundary Matrix Construction 1) Region Representation: Before demonstrating our interface for configuring matrix, we first illustrate our representation for iteration space and structural grid cells. M AT 2B OUNDARY will first suppose the existence of a multi-dimensional Cartesian grid which covers the range of
[Startd , Stopd ) in each dimension d. Then, M AT 2B OUND ARY configure the full region and the data region of each structural part of a block-structured grid as two slices of the global grid, and the ghost region is the difference between full and data. The slice is capable of not only decreasing the area size, but also applying bigger strides on each dimension to get a coarser grid. We use (R startd , R stopd , R stepd ) for describing the slicing in each dimension for the grid. While the start and stop variables are widely used for grid description [7], [5], the step variable is innovative in PDE solver frameworks and dsls. We list two benefits for the introduction of step: • Naturally support grid staggering technique with step = 2. Taking 2-D problems as examples, in our language, cellcentered, edge-centered, and nodal variables of structural grids now have nothing different from each other except for the detailed shape information. This feature is essential for the development of Finite Volume Method (FVM) solvers of PDE. k • Naturally support multi-grid method with step = 2 . 2) Basic Data Types and Operations: M AT 2B OUNDARY supports a wide range of statements in the procedures of calculating each element of y for y = A × x + b. These statements include creating temporary variables, if statements, while-loop statements without break or continue, and assign statements for temporary variables. We use EI and EF for describing the expression data type of the integer and floating parameters in the above statements. The opportunities for extra optimization in the compiletime arise when we find a subset of the expressions which only rely on the coordinate indices of the current y elements and other integer parameters that remain constant over the whole SpMV procedure (R startd , R stopd , R stepd ). We call these expressions as compile-time resolvable expressions and define two new sub-types CI and CB for Integer and Boolean respectively. The supported operations of the subtypes are listed below using BNF description. Note that the expressions loaded from temporary variables are not treated as resolvable. CI < CI | CI = CI
CB ::= |
CB ∧ CB | CB ∨ CB | ¬CB CI × Literal | CI ÷ Literal | CI mod Literal
CI ::= |
CI + CI | max(CI , CI ) | min(CI , CI )
|
Name | Literal
3) Basic Sub-Matrix/Vector API: As directly using the statements of data types EI , EF , CI , CB to implement a whole SpMV for Figure 3 is still tedious for programmers, basic sub-matrix/vector are introduced to decouple the boundary conditions into several sub-problems which are simple enough to omit the detailed SpMV procedure. Table II shows the core apis of a basic sub-matrix. Comparing with a general meaning row-iterator of sparse matrix, the basic sub-matrix needs an additional Boolean expression
Interface
Parameter Type
Description
is_valid
List[CI ] → CB
Is in coverage area (Rgp ∩ RD , e.g.)
iterator
List[CI ] → Iter
A sparse row of matrix (δij , Itpij , e.g.)
n_group
→ EI
Property for distributed communication
TABLE II: Basic Sub-matrix APIs Interface
Parameter
F_LIST
indices: List[CI ] → elements: List[T ]
T
Description return a static number of entries
Tuple[List[CI ], EF , EI ]
column, data and gid
F_NNZ
indices: List[CI ] → nnz: EI , info: Tuple[EI ]
return the length of the dynamic iterator
F_COL
indices: List[CI ], nnz: EI , info: Tuple[EI ], offset: EI → col: List[EI ]
the column of the offset-th entry of the iterator
F_DATA
List[CI ], EI , Tuple[EI ], EI → data: EF
entry property of data
F_GID
List[CI ], EI , Tuple[EI ], EI → gid: EI
entry property for distributed communication
TABLE III: Parameter and Description of the functions for row iterator building
to describe its responsible area in the whole matrix, therefore the return of iterator interface can be extremely simple. We will describe the usage of n_group in Section V-B. The basic sub-vector api is similar with basic sub-matrix. Instead of iterator and n_group interface, basic subvector only requires user to designate an expression with the indices as input. 4) Row Iterator API: As for row iterator, two formats are supported as the return value of the basic sub-matrix api for describing structured and unstructured part of boundary conditions. 1) Users can directly return a compile-time static list of nonzero elements as result, which can cover simple BCs like Circular and Symmetric Boundaries. This interface is called as F_LIST. 2) For the unstructured and user-defined boundary condition patterns, we demonstrate four interfaces, F_NNZ, F_COL and F_DATA and F_GID to describe an iterator over nonzero elements. Table III shows the meanings of the five function interface. As we can construct the row iterator either in compile time or in real execution using the provided function interface, we have the default implementation of SpMV and its detailed procedure can be omitted by user. 5) Example BC representation: Algorithm 1 shows the pseudo-code implementation of the three boundary conditions in Figure 3. The n_group and F_GID interfaces are omitted for simplicity. Notably, our unstructured api can represent generalized row-based sparse matrices as the user implement interfaces for the detailed storage format, and we provide the
basic implementation for CSR and ELL format in advance, therefore users can configure the boundary using one line like BoundaryCSR(ptr_arr, col_arr, data_arr, calc_addr, extract). Algorithm 1 Pseudo-code for the implementation of basic sub-matrix to describe boundaries in Figure 3 1: RD ← (0, 3, 1) × (0, 5, 1) × (0, 5, 1) ▷ 3D grids are used i ← (i, i + 1, 1) × (0, 5, 1) × (0, 5, 1) 2: RD ▷ for blocks of 2d grids i ← (i, i + 1, 1) × (−1, 6, 1) × (−1, 6, 1) 3: RF 1 4: Rright ← (1, 2, 1) × (0, 5, 1) × (5, 6, 1) 1 5: Rlef t ← (1, 2, 1) × (0, 5, 1) × (−1, 0, 1) 6: 7: function BC R IGHT. IS VALID(indices: List[CI ]) 1 ) 8: return I N A REA(indices, Rright 9: function BC L EFT. IS VALID(indices: List[CI ]) 1 10: return I N A REA(indices, Rlef t) 11: function BC Z ERO . IS VALID(indices: List[CI ]) 1 ) ∧ 12: return I N A REA(indices, RF 1 ) ∧ 13: ¬I N A REA(indices, RD 1 14: ¬I N A REA(indices, Rright )∧ 1 ) 15: ¬I N A REA(indices, Rlef t 16: 17: function BC R IGHT. ITERATOR(indices: List[CI ]) 18: function F LIST(indices) 19: return [([2, 9-indices[2], indices[1]], 1)] ▷ indices mapping 20: return wrap to iterator(F LIST) 21: 22: ptr arr ← Array[n row+1] ▷ use CSR matrix to 23: col arr ← Array[nnz] ▷ represent left boundary 24: data arr ← Array[nnz] 25: function BC L EFT. ITERATOR(indices: List[CI ]) 26: function F NNZ(indices) 27: pt ← calc addr(indices, n row) 28: addr 1 ← row ptr arr[pt] 29: addr 2 ← row ptr arr[pt+1] 30: row nnz ← addr 2 - addr 1 31: return row nnz, addr 1 32: function F COL(indices, row nnz, info, offset) 33: addr 1 ← info 34: raw col ← col arr[addr 1+offset] ▷ type is EI 35: col ← extract(raw col) ▷ type is List[EI ] 36: return col 37: function F DATA(indices, row nnz, info, offset) 38: addr 1 ← info 39: return data arr[addr 1+offset] 40: return wrap to iterator(F NNZ, F COL, F DATA) 41: function BC Z ERO . ITERATOR(indices: List[CI ]) 42: function F LIST(indices) 43: return [] 44: return wrap to iterator(F LIST)
B. SpMV Code Generation As mentioned in Section II-C, optimization opportunities arise when BCs have structured sub-parts and can be hardcoded into the program. We use the multi-stage programming techniques here to generate matrix-free code for these BCs. Multi-stage programming (commonly known as staging) has emerged as a promising technique for developing embedded domain-specific languages. In this paradigm, programs are typically divided into two distinct phases: • Stage 1 (compile-time): A partial evaluation of the user code occurs and produces an optimized intermediate program.
Stage 2 (runtime): The generated intermediate program is executed to compute the final results. The SpMV code generation occurs in Stage 1, when our compiler can easily distinguish between the iterator types and generate different code for them. Furthermore, we can discuss the relationships of these overlapping expressions in originSet with each other and get a new non-overlapping expression set newSet, whose additional features are listed below: •
∨newSet ei ≡ T RUE ∀ei ∈ newSet, ∀ej ∈ originSet either ei → ej ≡ T RUE or ei → ¬ej ≡ T RUE
(1) (2)
By iterating over the power set of the given basic submatrices/sub-vectors, we can easily build the newSet by discussing whether the basic sub-matrix/sub-vector belongs to current subsets of the complete collections. Then nonoverlapping expressions are formed and lifted to the beginning of the loop using a huge if-elsif-else statement, and multiple final assignments points to the vector y are generated corresponding with these expressions. This transformation gathers different BC computations together in certain cell areas and provides more opportunities for later ordinary optimizations and specific non-equivalent optimizations in Section IV-C. In real implementations, we find that the size of power set grows exponentially as the number of basic sub-matrices/subvectors increases, and most of the derived Boolean expressions can be proved to be False during the stage-1 compilation, since in real configurations the valid areas of different sub-matrices are merely intersected. Therefore, we utilize the backtrack and discuss algorithm shown in the Mat2Stencil DSL [3] to reduce the compiling cost and the final number of branches, since this algorithm can early-exit when a sub-part of a Boolean expression can be proved to be false by Z3 SMT solver [17] and skip generating code for their branches. Therefore we can directly obtain a SpMV for-loop with only a handful number of branches inside. C. Optimization for Boundary Condition Computation In this section, we introduce our techniques to optimize the generated code of boundary condition in the runtime stage. Our key idea is to discover optimization opportunities from a global view of the whole program, then we can rely on less instructions from users and gain flexibility in programming. The left part of Figure 4 shows the iteration space of the generated boundary condition code. The initial iteration space has redundant visits to the empty part of the sparse matrix A (central red parts) and parts of the ghost regions that may never be used by later calculations (corner gray parts for 2d-5 points stencil). Expressing boundary and inner computations in a unified IR enables us to extract the program’s complete dependency graph using polyhedral analysis and Presburger arithmetic. In the optimization phase, we capture fine-grained read-after-write
j
j
to distributed memory layouts by injecting partition-aware coordinates. This spatial mapping enables the M AT 2B OUND ARY DSL to isolate standard halo-exchange scenarios from complex, unstructured boundary interpolations at the matrix representation level, significantly reducing the programming burden on users. B. Domain-Specific Communication Backend
i
(a) Origin
i
(b) After
Fig. 4: Example on simplification of the iteration space of boundary condition algorithm in Figure 3
(RAW) dependencies from boundary assignments (S1 ) to inner calculations (S2 ), formulated as: D(S1 , S2 ) = {(i1 , i2 )|S2 at iteration i2 depends on S1 at iteration i1 }
(3)
where ik is the iteration vector of statement Sk in loop Lk , with L2 following L1 . For each boundary write Wi , we union the dependencies of all subsequent dependent reads Rj to construct its effective iteration space: [ U (Wi ) = {i|∃j, (i, j) ∈ D(Wi , Rj )} (4)
Evaluating the SpMV abstraction y = Ax+b in a distributed environment necessitates fetching remote data entries for vector x. Instead of relying on general-purpose sparse matrix communication routines—which incur heavy metadata and packing overheads (e.g., using hash tables to deduplicate fetch indices)—we exploit the domain-specific properties of PDE solvers. Leveraging the polyhedral analysis from Section IV-C, which partitions boundary operations into convex subsets, we guarantee that the remote data dependencies for each subset remain strictly compact in physical space. M AT 2B OUNDARY automatically intercepts these remote dependencies and generates an initialization phase that constructs localized bounding boxes around the required remote indices, with the help of n_group and F_GID interfaces implemented by the users. During the execution phase, these bounding boxes drive highly efficient, contiguous buffer packing and unpacking. This specialized communicator yields high-bandwidth bulk transmission while bypassing the expensive irregular indexing typically associated with unstructured BCs.
Rj
This effective space U (Wi ) allows us to safely eliminate unnecessary boundary assignments (gray regions in Figure 4). We then partition the remaining space U (Wi ) into disjoint convex sets. Because calculations within these separate subsets are embarrassingly parallel, their code can be independently instantiated. As shown in Figure 4 (right), all the user-defined interpolation, user-defined simple mapping and padding-zero boundaries are neatly split into simple Cartesian sub-spaces. We emphasize that M AT 2B OUNDARY compiler automatically removes redundancy of the iteration space of boundary condition and generates high-efficient codes of iterating that are comparable with manual implementation. We use FreeTensor [18] for representing BC code and inner calculation code together and for carrying out polyhedral analysis. V. I MPLEMENTATION To seamlessly deploy M AT 2B OUNDARY-generated solvers onto large-scale clusters, we bridge the gap between high-level matrix representations and low-level distributed infrastructure through several implementation strategies. A. Representation Extension for Real-World Workloads Real-world PDE solvers require abstractions beyond simple scalar grids. We extend M AT 2B OUNDARY’s internal representations to natively support multi-dimensional arrays at each grid point, allowing the compiler to co-optimize tightly coupled physical quantities (e.g., multi-component velocity vectors). Furthermore, we map the global iteration spaces
C. Parameterized IR and JIT Specialization A major drawback of fine-grained polyhedral analysis is its prohibitive compilation time, particularly when dealing with complex inner calculations. To achieve a scalable compilation pipeline, M AT 2B OUNDARY employs an IR-reusing technique. During the heavy Ahead-of-Time (AOT) compilation stage, problem sizes and partition topological parameters are treated as symbolic variables, producing a generalized, dependencyresolved IR. During the initialization phase of execution, a lightweight Just-In-Time (JIT) pass substitutes these symbols with actual runtime dimensions and applies fast, localized optimizations such as constant folding. This hybrid approach preserves all advanced polyhedral transformations while amortizing the heavy compilation overhead across cluster nodes. VI. E VALUATION A. Experimental Setup a) Workloads: We evaluate our language, compiler, and runtime with the applications on two grids: cubed-sphere grid, 3-D Cartesian grid with multi-grid solvers. For cubedsphere grids we choose HOPE [12] and MCV [10] methods of SWE equation solvers. The order of accuracy configuration for HOPE and MCV methods is chosen as 5 and 3 relatively. For multi-grid solvers we choose High Performance Conjugate Gradients (HPCG) [19]. We use the original implementations for the HOPE evaluation [20], which contains FORTRAN implementations for HOPE-iter and HOPE-matrix algorithm and an optimized
We evaluate the programmability of M AT 2B OUNDARY by comparing the LOC against manual FORTRAN baselines (Table IV). For boundary condition implementations on the cubed-sphere grid, M AT 2B OUNDARY utilizes only 18.7% to 32.0% of the FORTRAN LOC across different operators, achieving a geometric average of 23.5%. Coupled with the Mat2Stencil DSL [3] for inner calculations, the end-to-end implementations of the two PDE solvers require merely 36% of the total kernel LOC. These results confirm that integrating M AT 2B OUNDARY with existing structured-grid DSLs substantially minimizes the development effort for complex solvers on block-structured grids. Solver Method HOPE(I) [12]
MCV [10]
Operator
Boundary LOCs
Calculation LOCs
Reconstruction
385 vs. 77
39 vs. 18
Border Flux
178 vs. 57
159 vs. 42
Other
-
542 vs. (290+) 68 49 vs. 67
Flux Derivative
203 vs. 38
Result Averaging
168 vs. 43
23 vs. 10
Other
-
317 vs. (290+) 35
TABLE IV: Comparison of lines of code (LOC) for implementing two PDE solvers on cubed-sphere grids using native FORTRAN and our DSLs (M AT 2B OUNDARY and Mat2Stencil). The notation “(290+)” indicates 290 shared LOC for common infrastructure between the two methods.
102
10 1
10 2 Recons(I)Recons(M) BorderF
1.00 x 1.31 x
100
1.00 x 7.64 x
101
1.00 x 1.92 x
FORTRAN/C++ Pytorch Ours
1.00 x 0.96 x 0.92 x 1.00 x 0.16 x 4.46 x
B. Code Length and Development Effort
case, it achieves substantial and consistent speedups across the remaining four types of BC computations, demonstrating a highly favorable overall performance profile.
Time (s)
Pytorch implementation for HOPE-matrix algorithm. We use a hand-written matrix-free HPCG from [21] for the HPCG evaluation. For the MCV method, we use an open-source FORTRAN implementation as the benchmark [22]. b) Platform: We evaluate M AT 2B OUNDARY on a 24node cluster, each with 2 Intel Xeon 6258R 28-core CPUs, and 192 GB of main memory. The nodes are interconnected with each other by 100 Gbps InfiniBand HDR. We use GCC 13.2.0 to compile the programs, using g++ for C++ and gfortran for FORTRAN, including our generated codes and baseline implementations. We use pip to install PyTorch with version 2.6.0, and carry out optimizations of the PyTorch code using the interface torch.compile.
FDiff EdgeAvg
Fig. 5: BC computation time on a single core. Lower is better. 2) End-to-End Comparison: Figure 6 presents the end-toend performance of all evaluated applications using both a single core and a fully subscribed 56-core node. Crucially, transitioning these applications to a distributed execution model across the entire node required minimal developer effort; by modifying fewer than 20 lines of distribution configurations within the kernel codes, we successfully deployed the fully distributed versions. As illustrated in the benchmark results, M AT 2B OUNDARY achieves an average end-to-end speedup of 1.90 × on a single core , which remarkably scales to an average speedup of 4.75 × across the fully subscribed node. Notably, while some applications like HPCG exhibit a slight single-thread overhead, they achieve massive performance gains at scale (1.79 × on 56 cores). This amplified performance at scale is primarily attributed to two architectural advantages: 1) The underlying designs of M AT 2B OUNDARY and the Mat2Stencil DSL exhibit superior parallel efficiency and reduced scheduling overhead in multi-threading environments. 2) Although multi-process execution introduces inherent data copying overheads, it effectively mitigates NonUniform Memory Access (NUMA) bottlenecks and significantly enhances overall data locality.
C. Benchmarking We use the resolution of 6 × 90 × 90 for SWE solvers and 2243 for HPCG application for small-scale benchmarks on a single node. 1) BC computation Comparison: To evaluate the boundary condition (BC) computation cost, we execute 432 time steps for the SWE solvers on a single core and profile the BC computation times, as depicted in Figure 5. The worst-case scenario for M AT 2B OUNDARY occurs during the reconstruction operator of the HOPE-matrix algorithm. In this specific case, the BC workload is overwhelmingly dominated by a large Sparse Matrix-Vector Multiplication (SpMV) in CSR format. Although M AT 2B OUNDARY incurs a marginal 8% performance overhead in this highly specific SpMV-bound
D. Scalability Figure 7 illustrates the strong scaling efficiency of our applications as we scale from a single node (56 cores) to 1,344 CPU cores (a 24× scale-up). The applications exhibit distinct scaling behaviors based on their numerical properties. The MCV application achieves a remarkable 88% parallel efficiency (a 21.21× speedup). This optimal scaling occurs because MCV features the lowest ratio of boundary condition (BC) complexity to inner calculation workload, making it relatively compute-bound. In contrast, the higher-order HOPE methods (HOPE-Iter and HOPE-Matrix) and the communication-heavy HPCG benchmark maintain parallel efficiencies between 72% and 74%. For the HOPE
HOPE(M)
Parallel Efficiency (%)
9.65 x 10.37 x 11.16 x
17.61 x 17.83 x 21.21 x
5.52 x 5.46 x 5.79 x
2.88 x 2.89 x 2.94 x
1.00 x 1.00 x 1.00 x
672
40
1344
0
120
variants, the higher numerical accuracy inherently requires thicker halo exchanges and more complex BCs, which increases inter-node communication at scale. These scaling results reveal a fundamental trade-off for HPC applications: users must balance the desired numerical accuracy (which dictates BC complexity) against the achievable parallel efficiency. M AT 2B OUNDARY effectively supports this design space, keeping framework-level overheads minimal even for communication-heavy workloads like HPCG. E. Evaluating Compile time MCV 12.77 273.65 6.09 29.66 4299.07
HPCG 7.21 111.36 462.18 33.77 906.10
TABLE V: Comparing the compile and running time (in seconds) of M AT 2B OUNDARY-based programs. The running times listed are at the maximum tested scale. To assess the compilation overhead of M AT 2B OUNDARY, we report the end-to-end compilation times for all workloads at their maximum scale (largest problem sizes and core allocations), as detailed in Table V. Thanks to our IR-reusing technique, the stage-1 and stage-2 optimizations are executed strictly once on a 28-core NUMA of a node. For SWE solvers, this approach confines over 85% of the compilation cost to the
168
17.35 x
10.11 x
2.86 x
1.00 x 56
5.37 x
80
103
102
HPCG
HPCG (Perf) HPCG (Eff)
Parallel Efficiency (%)
100
MCV
HOPE(M) 9.67 287.40 11.18 39.26 4372.42
336
Number of Cores / Scale
60
20
(a) SWE with 6 × 1440 × 1440 resolution
Fig. 6: SWE Solvers and HPCG End-to-End benchmarks
HOPE(I) 8.99 648.04 14.28 43.55 4618.37
168
80
104
(b) single node with 56 cores
Stages Stage 1 Execution Stage 2 Optimize Specialization GCC Run
100
HOPE-Iter (Perf) HOPE-Matrix (Perf) MCV (Perf) HOPE-Iter (Eff) HOPE-Matrix (Eff) MCV (Eff)
56
Perf. (Mpt/s)
1.00 x 18.71 x
1.00 x 1.39 x 4.79 x
10 1 HOPE(I)
120
HPCG
101 100
Perf. (Mpt/s)
1.00 x 1.80 x MCV
1.00 x 1.79 x
102
HOPE(M)
101
(a) single thread with 1 core FORTRAN/C++ Pytorch Ours
1.00 x 4.37 x
Perf. (Mpt/s)
10 2 HOPE(I)
1.00 x 2.14 x 2.81 x
10 1
1.00 x 0.79 x
102
100
1.00 x 4.73 x
Perf. (Mpt/s)
101
103
FORTRAN/C++ Pytorch Ours
60 40 20
336
672
Number of Cores / Scale
1344 0
(b) HPCG with 12803 grid size
Fig. 7: Strong Scaling Performance and Parallel Efficiency
initial run. For HPCG, the compilation time is instead dominated by the later specialization stage due to the large number of loops inherent in multi-grid methods. Ultimately, for longrunning real-world applications, this compilation overhead will finally become negligible, confirming the practicality of M AT 2B OUNDARY. VII. R ELATED W ORKS A. Frameworks and DSLs for Block-Structured Grids While numerous frameworks support block-structured grids, their handling of boundary conditions (BCs) remains a bottleneck. PETSc [23], Hypre [7], and ExaStencils [24], [6] offer semi-structured BC interfaces, but are fundamentally tailored for linear solvers rather than Earth system models. Frameworks like Stella [25], Dawn [26], GT4Py [14], [4], and PSyclone [27], [28] target weather applications on cubedsphere grids. However, their reliance on external BC libraries for distributed partitioning tightly couples them to specific discretizations (e.g., FV3 [9] or GungHo [13]), limiting their extensibility. Other tools prioritize structured inner computations over BC flexibility: Devito [5], [16] supports unaligned sparse data, Open-Earth Compiler [29] targets GPU execution, and Mat2Stencil [3] simplifies explicit and implicit numerical methods. Even recent MLIR-based [30] compiler efforts [31] that successfully unify inner-structured calculations still lack a high-level, unified abstraction for complex boundary conditions.
B. Multi-Stage Programming Lightweight Modular Staging (LMS)[32] is one of such frameworks with Scala as its host language, which separates the two stages through only the type definitions: Rep[T] representing a T-typed value at the second stage. Mat2Stencil [3] provides a more flexible implementation in the Python host language. The backtracking-and-discussing procedure is the basis of Mat2Stencil’s ability to represent boundary calculations on standalone structural grids, and we extend this procedure to generate SpMV code given a list of basic submatrices/vectors. C. Polyhedral Analysis Extensive research exists on fine-grained dependence analysis at the statement-level granularity and its application to program transformation. Early studies are summarized in [33] and [34]. Subsequent advances introduced a mathematical theory named polyhedral analysis for systematic dependence analysis. When memory access patterns are formally defined through Presburger formulas, the dependence relationship can be inferred by solving coupled equations and in-equations. Existing frameworks include PPL [35] and isl [36]. Multiple compilers exploit Polyhedral Analysis and carry out optimizations on general C programs [37], [38], [39] or tensor programs [40], [41], [18]. We adopt FreeTensor [18] as our stage 1 code generation target for its integration of isl and implementation of dependenceaware transformations. VIII. D ISCUSSION While M AT 2B OUNDARY successfully modularizes boundary condition computations while maintaining high performance, several core assumptions and design trade-offs warrant further discussion: • Necessity of Multi-Stage Programming. The multi-stage programming paradigm is fundamental to our framework. The ability to resolve specific variables at compile-time is mandatory; otherwise, the intermediate representation (IR) generated by the code generation algorithm (Section IV-B) would suffer from exponential bloat, rendering the compilation computationally intractable. • Generality of the Global Optimization. The global optimization techniques proposed in Section IV-C are highly agnostic to the underlying IR. As long as the chosen IR and its infrastructure support polyhedral analysis alongside userdefined intrinsics, our optimization passes can be seamlessly applied. Once optimized, the inner-calculation kernels can be subsequently lowered into any specialized IR for more aggressive, architecture-specific optimizations. • Trade-offs in Communication-Computation Overlap. Unlike the distributed algorithms implemented in distributed Devito [16], M AT 2B OUNDARY does not explicitly overlap communications with the center-part of inner calculations. While such overlapping is technically feasible and beneficial for extreme-scale parallel efficiency in explicit methods, altering the iteration execution order breaks the semantic
equivalence for implicit numerical solvers. We prioritized correctness and algorithmic generality across diverse PDE solvers over marginal extreme-scale speedups. • Static Grid Assumptions and GPU Portability. Currently, the underlying communication library relies on a global initialization phase. Consequently, M AT 2B OUNDARY assumes a static grid topology and does not support dynamic meshes or Adaptive Mesh Refinement (AMR), which are occasionally used in specific computational fluid dynamics applications. Furthermore, this global initialization scheme currently presents a bottleneck for migrating M AT 2B OUND ARY to GPU architectures. Overcoming this initialization overhead for GPU portability is the primary focus of our future work. IX. C ONCLUSION We propose M AT 2B OUNDARY, a novel DSL that resolves the programmability and performance bottlenecks of boundary conditions (BCs) in distributed block-structured grids. By unifying diverse BCs into a sparse matrix-vector multiplication (SpMV) abstraction and integrating multi-stage programming with polyhedral compilation, our framework automatically generates highly optimized, matrix-free kernels. Evaluations on diverse workloads demonstrate profound impact: up to 7.6x speedups for BC computations, > 70% reduction in lines of code, and 72% − 88% strong-scaling efficiency across 1,344 CPU cores. ACKNOWLEDGEMENT We use Gemini Pro AI [42] to refine the grammar of the whole article and to generate plotting scripts for figures. R EFERENCES [1] G. Mcmechan, “Migration by extrapolation of time-dependent boundary values,” Geophysical Prospecting, vol. 31, pp. 413–420, 04 2006. [Online]. Available: https://doi.org/10.1111/j.1365-2478.1983.tb01060.x [2] W. C. Skamarock, J. B. Klemp, J. Dudhia, D. O. Gill, Z. Liu, J. Berner, W. Wang, J. G. Powers, M. G. Duda, D. M. Barker et al., “A description of the advanced research wrf version 4,” NCAR tech. note ncar/tn-556+ str, vol. 145, no. 10.5065, 2019. [Online]. Available: https://doi.org/10.5065/1DFH-6P97 [3] H. Cao, S. Tang, Q. Zhu, B. Yu, and W. Chen, “Mat2stencil: A modular matrix-based dsl for explicit and implicit matrix-free pde solvers on structured grid,” Proceedings of the ACM on Programming Languages, vol. 7, no. OOPSLA2, pp. 686–715, 2023. [Online]. Available: https://doi.org/10.1145/3622822 [4] E. G. Paredes, L. Groner, S. Ubbiali, H. Vogt, A. Madonna, K. Mariotti, F. Cruz, L. Benedicic, M. Bianco, J. VandeVondele et al., “Gt4py: High performance stencils for weather and climate applications using python,” arXiv preprint arXiv:2311.08322, 2023. [Online]. Available: https://doi.org/10.48550/arXiv.2311.08322 [5] M. Louboutin, M. Lange, F. Luporini, N. Kukreja, P. A. Witte, F. J. Herrmann, P. Velesko, and G. J. Gorman, “Devito (v3. 1.0): an embedded domain-specific language for finite differences and geophysical exploration,” Geoscientific Model Development, vol. 12, no. 3, pp. 1165–1187, 2019. [Online]. Available: https: //doi.org/10.5194/gmd-12-1165-2019 [6] C. Lengauer, S. Apel, M. Bolten, S. Chiba, U. Rüde, J. Teich, A. Größlinger, F. Hannig, H. Köstler, L. Claus et al., “Exastencils: advanced multigrid solver generation,” in Software for Exascale Computing-SPPEXA 2016-2019. Springer International Publishing Cham, 2020, pp. 405–452. [Online]. Available: https://doi.org/10.1007/ 978-3-030-47956-5 14
[7] R. D. Falgout and U. M. Yang, “hypre: A library of high performance preconditioners,” in International Conference on computational science. Springer, 2002, pp. 632–641. [Online]. Available: https://doi.org/10.100 7/3-540-47789-6 66 [8] M. Ji and F. Toepfer, “Dynamical core evaluation test report for noaa’s next generation global prediction system (nggps),” 2016. [Online]. Available: https://doi.org/10.25923/ztzy-qn82 [9] W. M. Putman and S.-J. Lin, “Finite-volume transport on various cubedsphere grids,” Journal of Computational Physics, vol. 227, no. 1, pp. 55– 78, 2007. [Online]. Available: https://doi.org/10.1016/j.jcp.2007.07.022 [10] C. Chen and F. Xiao, “Shallow water model on cubed-sphere by multi-moment finite volume method,” Journal of Computational Physics, vol. 227, no. 10, pp. 5019–5044, 2008. [Online]. Available: https://doi.org/10.1016/j.jcp.2008.01.033 [11] P. A. Ullrich, C. Jablonowski, and B. Van Leer, “High-order finitevolume methods for the shallow-water equations on the sphere,” Journal of Computational Physics, vol. 229, no. 17, pp. 6104–6134, 2010. [Online]. Available: https://doi.org/10.1016/j.jcp.2010.04.044 [12] L. Zhou, W. Xue, and X. Shen, “Hope: an arbitrary-order nonoscillatory finite-volume shallow water dynamical core with automatic differentiation,” Geoscientific Model Development, vol. 18, no. 21, pp. 8175–8201, 2025, https://doi.org/10.5194/gmd-18-8175-2025. [Online]. Available: https://gmd.copernicus.org/articles/18/8175/2025/ [13] A. Staniforth, T. Melvin, and N. Wood, “Gungho! a new dynamical core for the unified model,” in Proceeding of the ECMWF workshop on recent developments in numerical methods for atmosphere and ocean modelling, 2013, pp. 15–29. [14] T. Ben-Nun, L. Groner, F. Deconinck, T. Wicky, E. Davis, J. Dahm, O. D. Elbert, R. George, J. McGibbon, L. Tr”umper et al., “Productive performance engineering for weather and climate modeling with python,” in SC22: International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE, 2022, pp. 1–14. [Online]. Available: https://dl.acm.org/doi/10.5555/3571885.3571982 [15] H. Cao, “Artifact of mat2stencil: A modular matrix-based dsl for explicit and implicit matrix-free pde solvers on structured grid,” Jul. 2023. [Online]. Available: https://doi.org/10.5281/zenodo.8149701 [16] G. Bisbas, R. Nelson, M. Louboutin, F. Luporini, P. H. Kelly, and G. Gorman, “Automated mpi-x code generation for scalable finitedifference solvers,” in 2025 IEEE International Parallel and Distributed Processing Symposium (IPDPS). IEEE, 2025, pp. 689–701. [Online]. Available: https://doi.org/10.1109/IPDPS64566.2025.00067 [17] L. M. de Moura and N. S. Bjørner, “Z3: an efficient SMT solver,” in Tools and Algorithms for the Construction and Analysis of Systems, 14th International Conference, TACAS 2008, Held as Part of the Joint European Conferences on Theory and Practice of Software, ETAPS 2008, Budapest, Hungary, March 29-April 6, 2008. Proceedings, ser. Lecture Notes in Computer Science, C. R. Ramakrishnan and J. Rehof, Eds., vol. 4963. Springer, 2008, pp. 337–340. [Online]. Available: https://doi.org/10.1007/978-3-540-78800-3 24 [18] S. Tang, J. Zhai, H. Wang, L. Jiang, L. Zheng, Z. Yuan, and C. Zhang, “Freetensor: a free-form dsl with holistic optimizations for irregular tensor programs,” in Proceedings of the 43rd ACM SIGPLAN International Conference on Programming Language Design and Implementation, 2022, pp. 872–887. [Online]. Available: https: //doi.org/10.1145/3519939.3523448 [19] J. Dongarra, M. A. Heroux, and P. Luszczek, “High-performance conjugate-gradient benchmark: A new metric for ranking highperformance computing systems,” The International Journal of High Performance Computing Applications, vol. 30, no. 1, pp. 3–10, 2016. [Online]. Available: https://doi.org/10.1177/1094342015593158 [20] L. Zhou, “High order prediction environment,” Jul. 2025. [Online]. Available: https://doi.org/10.5281/zenodo.16635583 [21] Q. Zhu, H. Luo, C. Yang, M. Ding, W. Yin, and X. Yuan, “Enabling and scaling the HPCG benchmark on the newest generation sunway supercomputer with 42 million heterogeneous cores,” in International Conference for High Performance Computing, Networking, Storage and Analysis, SC 2021, St. Louis, Missouri, USA, November 14-19, 2021, B. R. de Supinski, M. W. Hall, and T. Gamblin, Eds. ACM, 2021, p. 57. [Online]. Available: https://doi.org/10.1145/3458817.3476158 [22] L. Zhou. (2022) MCV SW repository. Accessed: April 9, 2026. [Online]. Available: https://gitee.com/DwyaneChou/MCV SW [23] B. Satish, A. Shrirang, A. Mark, B. Jed, B. Peter, B. Kris, D. Lisandro, D. Alp, and E. Victor, “Gropp w., et al,” PETSc Users Manual, 2019.
[24] C. Lengauer, S. Apel, M. Bolten, A. Größlinger, F. Hannig, H. Köstler, U. Rüde, J. Teich, A. Grebhahn, S. Kronawitter et al., “Exastencils: Advanced stencil-code engineering,” in European Conference on Parallel Processing. Springer, 2014, pp. 553–564. [Online]. Available: https://doi.org/10.1007/978-3-319-14313-2 47 [25] T. Gysi, C. Osuna, O. Fuhrer, M. Bianco, and T. C. Schulthess, “Stella: a domain-specific tool for structured grid methods in weather and climate models,” in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, ser. SC ’15. New York, NY, USA: Association for Computing Machinery, 2015. [Online]. Available: https://doi.org/10.1145/2807591.2807627 [26] C. Osuna, T. Wicky, F. Thuering, T. Hoefler, and O. Fuhrer, “Dawn: a high-level domain-specific language compiler toolchain for weather and climate applications,” Supercomputing Frontiers and Innovations, vol. 7, no. 2, pp. 79–97, 2020. [Online]. Available: https://doi.org/10.14529/jsfi200205 [27] Science and Technology Facilities Council (STFC). (2021) PSycloneBench: small benchmarks used to inform the development of the PSyclone Domain-Specific Compiler. Accessed: Apr. 24, 2023. [Online]. Available: https://github.com/stfc/PSycloneBench [28] S. V. Adams, R. W. Ford, M. Hambley, J. Hobson, I. Kavčič, C. M. Maynard, T. Melvin, E. H. M”uller, S. Mullerworth, A. R. Porter et al., “Lfric: Meeting the challenges of scalability and performance portability in weather and climate models,” Journal of Parallel and Distributed Computing, vol. 132, pp. 383–396, 2019. [Online]. Available: https://doi.org/10.1016/j.jpdc.2019.02.007 [29] T. Gysi, C. M”uller, O. Zinenko, S. Herhut, E. Davis, T. Wicky, O. Fuhrer, T. Hoefler, and T. Grosser, “Domain-specific multi-level ir rewriting for gpu: The open earth compiler for gpu-accelerated climate simulation,” ACM Transactions on Architecture and Code Optimization (TACO), vol. 18, no. 4, pp. 1–23, 2021. [Online]. Available: https://doi.org/10.1145/3469030 [30] C. Lattner, M. Amini, U. Bondhugula, A. Cohen, A. Davis, J. Pienaar, R. Riddle, T. Shpeisman, N. Vasilache, and O. Zinenko, “Mlir: Scaling compiler infrastructure for domain specific computation,” in 2021 IEEE/ACM International Symposium on Code Generation and Optimization (CGO). IEEE, 2021, pp. 2–14. [Online]. Available: https://doi.org/10.1109/CGO51591.2021.9370308 [31] G. Bisbas, A. Lydike, E. Bauer, N. Brown, M. Fehr, L. Mitchell, G. Rodriguez-Canal, M. Jamieson, P. H. Kelly, M. Steuwer et al., “A shared compilation stack for distributed-memory parallelism in stencil dsls,” in Proceedings of the 29th ACM International Conference on Architectural Support for Programming Languages and Operating Systems, Volume 3, 2024, pp. 38–56. [Online]. Available: https://doi.org/10.1145/3620666.3651344 [32] T. Rompf and M. Odersky, “Lightweight modular staging: a pragmatic approach to runtime code generation and compiled dsls,” in Proceedings of the ninth international conference on Generative programming and component engineering, 2010, pp. 127–136. [Online]. Available: https://doi.org/10.1145/1868294.1868314 [33] D. A. Padua and M. J. Wolfe, “Advanced compiler optimizations for supercomputers,” Communications of the ACM, vol. 29, no. 12, pp. 1184–1201, 1986. [Online]. Available: https://doi.org/10.1145/7902.790 4 [34] K. Kennedy and J. R. Allen, Optimizing compilers for modern architectures: a dependence-based approach. Morgan Kaufmann Publishers Inc., 2001. [35] R. Bagnara, E. Ricci, E. Zaffanella, and P. M. Hill, “Possibly not closed convex polyhedra and the parma polyhedra library,” in International Static Analysis Symposium. Springer, 2002, pp. 213–229. [Online]. Available: https://doi.org/10.1007/3-540-45789-5 17 [36] S. Verdoolaege, “isl: An integer set library for the polyhedral model,” in International Congress on Mathematical Software. Springer, 2010, pp. 299–302. [Online]. Available: https://doi.org/10.1007/978-3-642-1 5582-6 49 [37] U. Bondhugula, A. Hartono, J. Ramanujam, and P. Sadayappan, “Pluto: A practical and fully automatic polyhedral program optimization system,” in Proceedings of the ACM SIGPLAN 2008 Conference on Programming Language Design and Implementation (PLDI 08), Tucson, AZ (June 2008). Citeseer, vol. 146, 2008. [38] S. Verdoolaege and G. Janssens, “Scheduling for ppcg,” Report CW, vol. 706, 2017. [39] M. M. Strout, M. Hall, and C. Olschanowsky, “The sparse polyhedral framework: Composing compiler-generated inspector-executor code,”
Proceedings of the IEEE, vol. 106, no. 11, pp. 1921–1934, 2018. [Online]. Available: https://doi.org/10.1109/JPROC.2018.2857721 [40] N. Vasilache, O. Zinenko, T. Theodoridis, P. Goyal, Z. DeVito, W. S. Moses, S. Verdoolaege, A. Adams, and A. Cohen, “Tensor comprehensions: Framework-agnostic high-performance machine learning abstractions,” arXiv preprint arXiv:1802.04730, 2018. [Online]. Available: https://doi.org/10.48550/arXiv.1802.04730 [41] R. Baghdadi, J. Ray, M. B. Romdhane, E. Del Sozzo, A. Akkas, Y. Zhang, P. Suriana, S. Kamil, and S. Amarasinghe, “Tiramisu: A polyhedral compiler for expressing fast and portable code,” in 2019 IEEE/ACM International Symposium on Code Generation and Optimization (CGO). IEEE, 2019, pp. 193–205. [Online]. Available: https://doi.org/10.1109/CGO.2019.8661197 [42] G. Team, R. Anil, S. Borgeaud, J.-B. Alayrac, J. Yu, R. Soricut, J. Schalkwyk, A. M. Dai, A. Hauth, K. Millican et al., “Gemini: a family of highly capable multimodal models,” arXiv preprint arXiv:2312.11805, 2023. [Online]. Available: https://doi.org/10.48550/arXiv.2312.11805