Source-linked AI summary
Architecture and performance of Devito, a system for automated stencil computation
Fabio Luporini, Michael Lange, Mathias Louboutin, Navjot Kukreja, Jan Hückelheim, Charles Yount, Philipp Witte, Paul H. J. Kelly, Felix J. Herrmann, Gerard J. Gorman
TL;DR
High-performance stencil applications need abstractions that preserve productivity while producing efficient, portable code. This paper presents Devito, a Python-based symbolic framework and compiler that lowers equations through intermediate representations into optimized C++, and demonstrates production-level performance on state-of-the-art architectures. The supported scope nevertheless includes practical boundaries such as arithmetic-intensity requirements for some memory-saving transformations and incomplete distributed-memory and convenient staggered-grid support.
Problem
High-performance software interweaves scientific methods with hardware-specific optimization, making code difficult to maintain, debug, and port across architectures.
Method
Devito uses a Python-based symbolic DSL, compiler intermediate representations, and backend optimizations to generate executable C++ for stencil computations.
Results
Devito generates production-level code with compelling performance on state-of-the-art architectures, while its compilation and code-generation overheads are generally quick on evaluated systems.
Takeaways & Limitations
Devito supports general loop nests and heterogeneous expression sequences while composing with Python workflows and existing components such as SymPy, NumPy, and YASK.
Takeaways & Limitations
Some memory-saving transformations are useful mainly for arithmetic-intensive kernels, while distributed-memory parallelism and convenient staggered-grid support remain planned or limited.
Abstract
from arXiv · showhide
Stencil computations are a key part of many high-performance computing applications, such as image processing, convolutional neural networks, and finite-difference solvers for partial differential equations. Devito is a framework capable of generating highly-optimized code given symbolic equations expressed in Python, specialized in, but not limited to, affine (stencil) codes. The lowering process---from mathematical equations down to C++ code---is performed by the Devito compiler through a series of intermediate representations. Several performance optimizations are introduced, including advanced common sub-expressions elimination, tiling and parallelization. Some of these are obtained through well-established stencil optimizers, integrated in the back-end of the Devito compiler. The architecture of the Devito compiler, as well as the performance optimizations that are applied when generating code, are presented. The effectiveness of such performance optimizations is demonstrated using operators drawn from seismic imaging applications.
1 INTRODUCTION
Devito addresses the complexity of high-performance scientific software by separating mathematical specification from implementation and optimization. It provides a Python-based DSL and compiler that generate optimized C++ for stencil computations across varied applications and architectures.
- DSLs separate domain knowledge, numerical methods, and hardware-specific optimization, addressing monolithic HPC code that is difficult to maintain, debug, and port.
- Devito expresses finite-difference methods through high-level mathematical syntax while supporting general affine stencil computations and arbitrarily nested, possibly irregular, loops.
- Python provides access to scientific libraries, while Devito transforms symbolic input into optimized C++ that can also support incremental modernization of legacy software.
- Multiple intermediate representations and optimization passes enable complex transformations and programmatic performance portability across target architectures.
- Devito combines symbolic FLOP reduction with loop-level optimizations implemented by its own optimizer or integrated stencil compilers such as YASK.
- The evaluation uses test cases inspired by real-world seismic-imaging problems.
2 RELATED WORK
Devito builds on DSL research for numerical methods and distinguishes itself through compiler-based optimization, symbolic mathematics, Python composability, and graceful degradation. Its broader goal is to combine high-level usability with efficient, portable stencil code.
- Earlier DSL systems target broad mathematics or specialized PDE methods, including finite-element and finite-volume workflows.
- Devito distinguishes itself through graceful degradation, symbolic mathematics, compiler-based rather than template-based generation, and a native Python interface.
- Stencil tools include DSLs, polyhedral compilers, and user-driven generators, which Devito aims to integrate for performance portability.
- Unlike software specialized to a small set of seismic wave equations, Devito supports arbitrary wave-equation forms with greater portability and composability.
- Devito’s compiler targets FLOP reduction, data locality, and parallelism, including common-subexpression elimination, factorization, code motion, and locality optimizations.
3 SPECIFICATION OF A FINITE-DIFFERENCE METHOD WITH DEVITO
Devito lets users express finite-difference and general stencil computations symbolically, then constructs executable solvers from grids, symbolic functions, discretized equations, boundary conditions, and sparse operations. Its API supports concise stencil definitions, automatic equation rearrangement, source injection, and time-history control.
- Devito’s DSL uses mathematical notation and SymPy to express finite-difference, convolution, and basic linear-algebra operations.
- The workflow begins by defining a computational Grid that contains model parameters, wavefields, and sources over specified spatial dimensions.
- Function, TimeFunction, and sparse-function objects represent spatial fields, time-dependent wavefields, and quantities defined only on selected grid points.
- Devito provides shorthand derivatives and Laplacians, allowing common finite-difference stencils and complete wave equations to be written concisely.
- The solve function automatically rearranges a PDE expression into an explicit update rule for the next wavefield time step.
- SparseFunctions inject or sample values at selected coordinates, using interpolation and its adjoint to maintain a consistent discrete adjoint.
- Boundary conditions can be imposed through equation keywords for simple cases or written explicitly for more exotic schemes.
- TimeFunction storage defaults to a time-order-sized rolling buffer, while ConditionalDimension can save periodic time slices or longer histories.
4 THE DEVITO COMPILER
Devito lowers symbolic equations through multiple compiler passes that analyze domains and dependencies, group equations into Clusters, and prepare optimized executable code. Its Operator coordinates code generation, JIT compilation, and execution.
- Operator and equation lowering: The Operator generates low-level code, JIT-compiles it, executes it, and can be reused with different input data after one compilation.It accepts symbolic equations, substitutions, and optimization levels; compilation is triggered by the first execution.
- Equation lowering: Equation lowering converts indexed Functions into arrays, applies substitutions, and aligns accesses with halo and padding regions.Substitutions can improve constant folding but reduce Operator genericity, while unspecified symbol values remain runtime inputs.
- Equation analysis: Local analysis extracts input and output Functions, orders Dimensions, and constructs iteration and data spaces for each equation.These spaces describe iteration domains and accessed data domains, including dimensions that appear in only one of them.
- Clustering: Clusters group equations with matching iteration spaces and control flow when no dimension-carried true anti-dependencies prevent grouping.Dependence analysis resolves iteration directions and handles reductions specially, while conservative testing assumes a dependency when inconclusive.
- Clustering: Clustering algorithms output an ordered sequence of Clusters while enforcing directions, checking anti-dependencies, and handling control flow.The algorithms may preserve original directions when conflicting requirements would otherwise arise.
4.4 Symbolic optimization
The Devito Symbolic Engine reduces the arithmetic strength of Clusters through a sequence of symbolic optimization passes.
- Symbolic optimization: The DSE applies standard common sub-expression elimination and more advanced rewrites individually to each Cluster.Its output may change the number of Clusters, equations, and arithmetic-operation sequence.
4.5 IET construction
IET construction lowers Clusters into an Iteration/Expression Tree whose nested Iterations represent loops and whose Expressions contain equations.
- IET construction: The IET represents each Cluster in a loop nest, while suitable nesting of Iterations constructs the generated loop structure.Some outer loops may be shared by multiple Clusters.
- IET construction: The scheduling algorithm compares each Cluster’s iteration space with the current schedule to build or reuse nested Iterations.In the running example, the schedule handles stencil, save, and injection Clusters in sequence.
4.6 IET analysis
IET analysis attaches sequential, parallel, and vectorizable properties to Iterations so later loop-optimization passes can use them.
- IET analysis: The analysis determines whether Iterations are sequential or parallel using data-dependence information represented by distance vectors.These properties are attached to IET nodes before later loop optimization.
- IET analysis: The cluster-scheduling algorithm converts a sequence of Clusters into an Iteration/Expression Tree and can share outer Iterations across Clusters.Its input is a sequence of Clusters and its output is an IET.
- IET analysis: The resulting IET contains nested loops over dimensions such as t, x, and p_q, with expressions and conditional structure embedded in the tree.The graphical representation shows separate loop and expression nodes for the running example.
4.7 IET optimization
The IET optimization macro-pass transforms Devito’s intermediate expression tree for runtime performance and faster JIT compilation, while supporting modular backend specialization. It introduces loop optimizations directly or delegates them to backend engines such as YASK.
- The macro-pass optimizes the IET for runtime performance and rapid JIT compilation by the underlying C compiler.It introduces loop blocking, remainder-loop minimization, SIMD vectorization, OpenMP shared-memory parallelism, and software prefetching.
- The IET is completed with variable declarations, headers, and profiling instrumentation before C code generation and JIT compilation.Declarations respect scope and OpenMP private/shared semantics; generated operators are cached and compiled into shared objects.
- Backends specialize data types and Operator passes while preserving software modularity and composability.This infrastructure is intended to integrate specialized stencil tools and support future architectures such as GPUs.
- Devito provides a core backend using the DLE and a yask backend using YASK to generate optimized C++ for Intel Xeon and Xeon Phi.The yask backend uses YASK for loop optimization, data management, JIT compilation, and execution; both backends share the pipeline until loop optimization.
5 AUTOMATED PERFORMANCE OPTIMIZATIONS
Devito combines symbolic transformations with loop-level optimization to reduce arithmetic, improve locality, and expose parallelism. Its DLE and YASK paths provide complementary stencil optimizations, including SIMD layouts, blocking, and vector-folding.
- Devito’s optimization pipeline combines symbolic arithmetic reduction with loop transformations for data locality and parallelism.Optimization modes compose predefined sequences of DSE, DLE, and YLE passes.
- DSE applies common sub-expression elimination, factorization, extraction, and shift-invariant detection to reduce expression strength and repeated work.The available CSE implementations trade recognition power against turnaround time and preservation of factorization opportunities.
- The DLE applies SIMD vectorization, loop blocking, and parallelism, but its SIMD strategy can produce many unaligned vector loads and stores.Blocking currently targets fully parallel iterations, and block shapes are selected through empirical auto-tuning.
- Vector-folding packs neighboring values into multidimensional tiles, increasing overlap and reuse between stencil iterations and reducing memory-bandwidth demand for bandwidth-bound stencils.The described eight-element layouts include 1D, 2D, and 3D folds; the 1 × 1 × 8 layout is traditional inline vectorization.
- Devito integrates YASK so symbolic processing can feed YASK’s low-level stencil optimizations through a dedicated backend.In Devito v3.1, the yask backend supported roughly 70% of the Devito API.
6 THE SHIFT-INVARIANTS ELIMINATION ALGORITHM
Shift-invariants are cross-iteration redundancies identified by translated operand structure and replaced with shared pivots to reduce operation counts. Devito combines extraction, detection, and blocked temporary storage, trading memory for arithmetic savings.
- Shift-invariants are cross-iteration redundancies whose presence depends on the PDE’s differential operators and discretization scheme.They matter particularly in kernels with high arithmetic intensity, while being less consequential in memory-bound operators.
- Candidate extraction selects maximal subexpressions using empirically chosen operation-count thresholds, including a default time-invariant threshold of Thr0 = 10.The extraction pass is reused by shift-invariant elimination.
- Two expressions are shifted when corresponding operands share compatible coefficients, labels, dimensions, and operators under a common translation of displacements.The relation excludes mismatched labels, dimensions, or non-translatable displacement patterns.
- Less than 2 seconds were required to detect shift-invariants for the challenging tti test case with so=16 despite O(n^2) complexity.The reported timing was measured on an Intel Xeon E5-2620 v4.
- The detection algorithm compares candidate expressions, prunes the search, and constructs pivots that replace members of each shift-invariant set.For a set of k expressions, pivot reuse yields an operation-count reduction proportional to k.
- SIE trades arithmetic operations for temporary-array memory, so its runtime benefit is mainly expected in arithmetic-intensive kernels.Whole-grid temporary arrays rarely help unless operation-count reductions are exceptionally high.
- The SIE algorithm uses cross-loop-nest blocking so array temporaries scale with block shape rather than spanning the entire grid.One loop nest produces temporaries and a subsequent loop nest consumes them; the DLE supplies blocking across loop sequences.
7 PERFORMANCE EVALUATION
The evaluation measures generated-code performance across architectures, optimization modes, and seismic wave-equation test cases. Results show backend- and architecture-dependent behavior, with YASK improving isotropic performance and DSE choices strongly affecting TTI runtime and operation counts.
- Experimental methodology: The experiments use enriched roofline plots, joint GFlops/s and runtime analysis, and attainable-bandwidth and floating-point ceilings to assess code quality.Operational intensity uses a realistic model that reloads time-invariant functions at each time iteration.
- Optimization modes: Devito optimization modes vary compiler time, operation-count reduction, and temporary-memory allocation while sharing OpenMP parallelism, SIMD vectorization, and loop blocking.Basic, advanced, and aggressive modes differ primarily in their symbolic transformations.
- Isotropic performance: At space order 16, isotropic core performance on skl8180 falls from 59% to 36% of attainable machine peak as unaligned memory streams and SIMD pack/unpack overhead increase.Profiling identifies memory pressure, split-load handling, and auto-vectorized SIMD quality as major issues; disabling DSE barely changes runtime.
- Isotropic performance: On knl7250, YASK is roughly 3× faster than core at space order 4 and more than 4× faster at space order 12.Vector folding, software prefetching, and hierarchical OpenMP help YASK address core’s memory and parallelism limitations.
- TTI performance: TTI benefits substantially from DSE because symbolic transformations reduce operation count and register pressure, although its higher arithmetic intensity and load count limit performance relative to isotropic.The evaluation covers multiple space orders on skl8180 and knl7250 using two cubic grids.
- TTI performance: At space order 8, advanced becomes slower than aggressive on both architectures and grids; beyond that point, aggressive’s operation-count reduction outweighs its additional memory and cache-communication overhead.The operation-count reduction grows almost quadratically with space order.
- TTI performance: On knl7250, advanced reaches approximately 17% of attainable peak at space order 4 with a 512^3 grid, while exceeding MCDRAM capacity at the larger grid contributes to a runtime increase from 34s to 173s.The larger grid has a 25.5GB working set versus roughly 7.5GB for the smaller grid.
- Overhead summary: On skl8180, generating C code for TTI at space order 16 takes around 3 seconds and compiling it takes less than 7 seconds, while knl7250 compilation takes slightly more than one minute.Auto-tuning takes about 3 minutes on skl8180 and 15 minutes on knl7250, but these costs are amortized in production runs.
8 FURTHER WORK
Devito’s further work targets broader execution support and additional optimization capabilities. Priorities include distributed-memory parallelism, easier staggered-grid use, completing the yask backend, and adding advanced optimizations and GPU support.
- Distributed-memory parallelism via MPI is planned for applications requiring solvers to run across multiple compute nodes.The immediate plan is to leverage yask’s MPI support, with possible later integration into Devito core.
- Staggered grids are needed for many finite-difference discretizations, including elastic wave propagation, but Devito currently exposes them only through a low-level API.Basic support exists in Devito v3.1; making it more convenient is planned.
- The yask backend is not feature-complete because it cannot run tti equations with array temporaries.The paper identifies extending this capability as a high-priority task because tti is an advanced industrial wave-propagation model.
- Planned optimizations include time tiling, on-the-fly data compression, and mixed-precision arithmetic using application knowledge.An ops backend is also being developed for GPU code generation and distributed-memory parallelism via MPI.
9 CONCLUSIONS
Devito automates high-performance stencil computation through a Python interface that expresses finite-difference approximations and more general loop nests. Its experiments show production-level generated code with compelling performance on modern architectures.
- Devito automates high-performance stencil computations while supporting Python-based finite-difference expressions and arbitrary loop nests.Operators can evaluate heterogeneous expression sequences arising in finite-difference solvers, linear algebra, or interpolation.
- The compiler generates production-level code with compelling performance on state-of-the-art architectures.