Stencil Codes Research Articles

Analytical performance estimation during code generation on modern GPUs

Automatic code generation is frequently used to create implementations of algorithms specifically tuned to particular hardware and application parameters. The code generation process involves the selection of adequate code transformations, tuning parameters, and parallelization strategies. We propose an alternative to time-intensive autotuning, scenario-specific performance models, or black-box machine learning to select the best-performing configuration.This paper identifies the relevant performance-defining mechanisms for memory-intensive GPU applications through a performance model coupled with an analytic hardware metric estimator. This enables a quick exploration of large configuration spaces to identify highly efficient code candidates with high accuracy.We examine the changes of the A100 GPU architecture compared to the predecessor V100 and address the challenges of how to model the data transfer volumes through the new memory hierarchy.We show how our method can be coupled to the “pystencils” stencil code generator, which is used to generate kernels for a range-four 3D-25pt stencil and a complex two-phase fluid solver based on the Lattice Boltzmann Method. For both, it delivers a ranking that can be used to select the best-performing candidate.The method is not limited to stencil kernels but can be integrated into any code generator that can generate the required address expressions.

Accelerating high‐order stencils on GPUs

SummaryFinite‐difference methods based on high‐order stencils are commonly used for modeling of seismic wave propagation, weather forecasting, computational fluid dynamics, convolutional neural networks, and others. Nowadays, the community commonly employs graphics processing units (GPUs) to accelerate such stencil computations. As a result, knowing how to write efficient stencil computations for GPUs is of significant interest. While high‐performance, low‐order stencils on GPUs have been studied extensively in the literature, not all proposed approaches work well for high‐order stencils. Furthermore, coping with boundary conditions used with stencils for seismic modeling makes it challenging to efficiently exploit thread‐level parallelism on GPUs. In this article, we describe several implementations of a 25‐point stencil. We evaluate our stencil code shapes, memory hierarchy usage, data access patterns, and other performance attributes on several modern GPUs and compare them with machine rooflines. On average, our top‐performing kernels achieve six times the performance of a 25‐point stencil code developed in C and mapped to GPUs using OpenACC. Several of our implementations have excellent performance portability across multiple generations of both NVIDIA and AMD GPUs.

Efficient Nearest-Neighbor Data Sharing in GPUs

Stencil codes (a.k.a. nearest-neighbor computations) are widely used in image processing, machine learning, and scientific applications. Stencil codes incur nearest-neighbor data exchange because the value of each point in the structured grid is calculated as a function of its value and the values of a subset of its nearest-neighbor points. When running on Graphics Processing Unit (GPUs), stencil codes exhibit a high degree of data sharing between nearest-neighbor threads. Sharing is typically implemented through shared memories, shuffle instructions, and on-chip caches and often incurs performance overheads due to the redundancy in memory accesses. In this article, we propose Neighbor Data (NeDa), a direct nearest-neighbor data sharing mechanism that uses two registers embedded in each streaming processor (SP) that can be accessed by nearest-neighbor SP cores. The registers are compiler-allocated and serve as a data exchange mechanism to eliminate nearest-neighbor shared accesses. NeDa is embedded carefully with local wires between SP cores so as to minimize the impact on density. We place and route NeDa in an open-source GPU and show a small area overhead of 1.3%. The cycle-accurate simulation indicates an average performance improvement of 21.8% and power reduction of up to 18.3% for stencil codes in General-Purpose Graphics Processing Unit (GPGPU) standard benchmark suites. We show that NeDa’s performance is within 13.2% of an ideal GPU with no overhead for nearest-neighbor data exchange.

A Data-Centric Directive-Based Framework to Accelerate Out-of-Core Stencil Computation on a GPU

Graphics processing units (GPUs) are highly efficient architectures for parallel stencil code; however, the small device (i.e., GPU) memory capacity (several tens of GBs) necessitates the use of out-of-core computation to process excess data. Great programming effort is needed to manually implement efficient out-of-core stencil code. To relieve such programming burdens, directive-based frameworks emerged, such as the pipelined accelerator (PACC); however, they usually lack specific optimizations to reduce data transfer. In this paper, we extend PACC with two data-centric optimizations to address data transfer problems. The first is a direct-mapping scheme that eliminates host (i.e., CPU) buffers, which intermediate between the original data and device buffers. The second is a region-sharing scheme that significantly reduces host-to-device data transfer. The extended PACC was applied to an acoustic wave propagator, automatically extending the length of original serial code 2.3-fold to obtain the out-of-core code. Experimental results revealed that on a Tesla V100 GPU, the generated code ran 41.0, 22.1, and 3.6 times as fast as implementations based on Open Multi-Processing (OpenMP), Unified Memory, and the previous PACC, respectively. The generated code also demonstrated usefulness with small datasets that fit in the device capacity, running 1.3 times as fast as an in-core implementation.

A domain-specific language and matrix-free stencil code for investigating electronic properties of Dirac and topological materials

We introduce PVSC-DTM (Parallel Vectorized Stencil Code for Dirac and Topological Materials), a library and code generator based on a domain-specific language tailored to implement the specific stencil-like algorithms that can describe Dirac and topological materials such as graphene and topological insulators in a matrix-free way. The generated hybrid-parallel (MPI+OpenMP) code is fully vectorized using Single Instruction Multiple Data (SIMD) extensions. It is significantly faster than matrix-based approaches on the node level and performs in accordance with the roofline model. We demonstrate the chip-level performance and distributed-memory scalability of basic building blocks such as sparse matrix-(multiple-) vector multiplication on modern multicore CPUs. As an application example, we use the PVSC-DTM scheme to (i) explore the scattering of a Dirac wave on an array of gate-defined quantum dots, to (ii) calculate a bunch of interior eigenvalues for strong topological insulators, and to (iii) discuss the photoemission spectra of a disordered Weyl semimetal.

Efficient Acceleration of Stencil Applications through In-Memory Computing.

The traditional computer architectures severely suffer from the bottleneck between the processing elements and memory that is the biggest barrier in front of their scalability. Nevertheless, the amount of data that applications need to process is increasing rapidly, especially after the era of big data and artificial intelligence. This fact forces new constraints in computer architecture design towards more data-centric principles. Therefore, new paradigms such as in-memory and near-memory processors have begun to emerge to counteract the memory bottleneck by bringing memory closer to computation or integrating them. Associative processors are a promising candidate for in-memory computation, which combines the processor and memory in the same location to alleviate the memory bottleneck. One of the applications that need iterative processing of a huge amount of data is stencil codes. Considering this feature, associative processors can provide a paramount advantage for stencil codes. For demonstration, two in-memory associative processor architectures for 2D stencil codes are proposed, implemented by both emerging memristor and traditional SRAM technologies. The proposed architecture achieves a promising efficiency for a variety of stencil applications and thus proves its applicability for scientific stencil computing.

Tiling Optimizations for Stencil Computations Using Rewrite Rules in L ift

Stencil computations are a widely used type of algorithm, found in applications from physical simulations to machine learning. Stencils are embarrassingly parallel, therefore fit on modern hardware such as Graphic Processing Units perfectly. Although stencil computations have been extensively studied, optimizing them for increasingly diverse hardware remains challenging. Domain-specific Languages (DSLs) have raised the programming abstraction and offer good performance; however, this method places the burden on DSL implementers to write almost full-fledged parallelizing compilers and optimizers. Lift has recently emerged as a promising approach to achieve performance portability by using a small set of reusable parallel primitives that DSL or library writers utilize. L ift ’s key novelty is in its encoding of optimizations as a system of extensible rewrite rules which are used to explore the optimization space. This article demonstrates how complex multi-dimensional stencil code and optimizations are expressed using compositions of simple 1D L ift primitives and rewrite rules. We introduce two optimizations that provide high performance for stencils in particular: classical overlapped tiling for multi-dimensional stencils and 2.5D tiling specifically for 3D stencils. We provide an in-depth analysis on how the tiling optimizations affects stencils of different shapes and sizes across different applications. Our experimental results show that our approach outperforms existing compiler approaches and hand-tuned codes.

Impact study of data locality on task-based applications through the Heteroprio scheduler.

The task-based approach has emerged as a viable way to effectively use modern heterogeneous computing nodes. It allows the development of parallel applications with an abstraction of the hardware by delegating task distribution and load balancing to a dynamic scheduler. In this organization, the scheduler is the most critical component that solves the DAG scheduling problem in order to select the right processing unit for the computation of each task. In this work, we extend our Heteroprio scheduler that was originally created to execute the fast multipole method on multi-GPUs nodes. We improve Heteroprio by taking into account data locality during task distribution. The main principle is to use different task-lists for the different memory nodes and to investigate how locality affinity between the tasks and the different memory nodes can be evaluated without looking at the tasks’ dependencies. We evaluate the benefit of our method on two linear algebra applications and a stencil code. We show that simple heuristics can provide significant performance improvement and cut by more than half the total memory transfer of an execution.

A large-scale framework for symbolic implementations of seismic inversion algorithms in Julia

Writing software packages for seismic inversion is a very challenging task because problems such as full-waveform inversion or least-squares imaging are algorithmically and computationally demanding due to the large number of unknown parameters and the fact that waves are propagated over many wavelengths. Therefore, software frameworks need to combine versatility and performance to provide geophysicists with the means and flexibility to implement complex algorithms that scale to exceedingly large 3D problems. Following these principles, we have developed the Julia Devito Inversion framework, an open-source software package in Julia for large-scale seismic modeling and inversion based on Devito, a domain-specific language compiler for automatic code generation. The framework consists of matrix-free linear operators for implementing seismic inversion algorithms that closely resemble the mathematical notation, a flexible resilient parallelization, and an interface to Devito for generating optimized stencil code to solve the underlying wave equations. In comparison with many manually optimized industry codes written in low-level languages, our software is built on the idea of independent layers of abstractions and user interfaces with symbolic operators. Through a series of numerical examples, we determined that this allows users to implement a series of increasingly complex algorithms for waveform inversion and imaging as simple Julia scripts that scale to large-scale 3D problems. This illustrates that software based on the paradigms of abstract user interfaces and automatic code generation and makes it possible to manage the complexity of the algorithms and performance optimizations, thus providing a high-performance research and production framework.

OpenFPM: A scalable open framework for particle and particle-mesh codes on parallel computers

Scalable and efficient numerical simulations continue to gain importance, as computation is firmly established as the third pillar of discovery, alongside theory and experiment. Meanwhile, the performance of computing hardware grows through increasingly heterogeneous parallelism, enabling simulations of ever more complex models. However, efficiently implementing scalable codes on heterogeneous, distributed hardware systems becomes the bottleneck. This bottleneck can be alleviated by intermediate software layers that provide higher-level abstractions closer to the problem domain, reducing development times and allowing computational scientists to focus. Here, we present OpenFPM, an open and scalable framework that provides an abstraction layer for numerical simulations using particles and/or meshes. OpenFPM provides transparent and scalable infrastructure for shared-memory and distributed-memory implementations of particles-only and hybrid particle-mesh simulations of both discrete and continuous models, as well as non-simulation codes. This infrastructure is complemented with frequently used numerical routines, as well as interfaces to third-party libraries. We present the architecture and design of OpenFPM, detail the underlying abstractions, and benchmark the framework in applications ranging from Smoothed-Particle Hydrodynamics (SPH) to Molecular Dynamics (MD), Discrete Element Methods (DEM), Vortex Methods, stencil codes (finite differences), and high-dimensional Monte Carlo sampling (CMA-ES), comparing it to the current state of the art and to existing software frameworks. Program summaryProgram Title: OpenFPMProgram Files doi:http://dx.doi.org/10.17632/4yrp8nbm7c.1Licensing provisions: GPLv3Programming language: C++Nature of problem: Writing numerical simulation programs that use meshes, particles, or any combination of the two typically requires long development times, in particular if the code is to scale efficiently on parallel distributed-memory computers. The long development times incur high financial and project-time costs and often lead to sub-optimal program performance as shortcuts are taken. Yet, a large portion of the functionality is common across programs and could be automated or provided as reusable software components, leading to large savings in project costs and potentially improved software performance.Solution method: OpenFPM provides a scalable, highly efficient software platform for numerical simulations using meshes, particles, or any combination of the two on parallel computers. It is based on a well-known set of abstract data types and operators that suffice to express any such simulation, regardless of the application domain. OpenFPM provides reusable, tested, and internally parallelized software components that reduce development times and make parallel computing accessible to computational scientists without extensive knowledge in parallel programming.Additional comments including restrictions and unusual features: OpenFPM is a software library based on which users can implement their simulation codes at a fraction of the development cost. All parallelization and memory handling is transparently done by the library. As its main innovation, OpenFPM makes use of C++ Template Meta Programming in order to enable simulations in arbitrary-dimensional spaces, distribution of arbitrary user-defined C++ objects, and compile-time code optimization and targeting for specific hardware platforms. OpenFPM-based simulations can directly output VTK files for visualization of results and HDF5 files for data archiving.

Polyhedral Search Space Exploration in the ExaStencils Code Generator

Performance optimization of stencil codes requires data locality improvements. The polyhedron model for loop transformation is well suited for such optimizations with established techniques, such as the PLuTo algorithm and diamond tiling. However, in the domain of our project ExaStencils, stencil codes, it fails to yield optimal results. As an alternative, we propose a new, optimized, multi-dimensional polyhedral search space exploration and demonstrate its effectiveness: we obtain better results than existing approaches in several cases. We also propose how to specialize the search for the domain of stencil codes, which dramatically reduces the exploration effort without significantly impairing performance.

Massively Parallel Stencil Code Solver with Autonomous Adaptive Block Distribution

In the last decades, simulations have been established in several fields of science and industry to study various phenomena by solving, inter alia, partial differential equations. For an efficient use of current and future high performance computing systems, with many thousands of computation ranks, high node-level performance, scalable communication, and the omission of unnecessary calculations are of high priority in the development of new solvers. The challenge of contemporary simulation applications is to bridge the gap between the scales of the various physical processes. We introduce the NAStJA framework, a block-based MPI parallel solver for arbitrary algorithms, based on stencil code or other regular grid methods. NAStJA decomposes the domain of spatially complex structures into small cuboid blocks. A special feature of NAStJA is the dynamic block adaption which modifies the calculation domain around the region where the computation currently takes place, and hence avoids unnecessary calculations. This often occurs, inter alia, in phase-field simulations. Block creation and deletion is managed autonomously within local neighborhoods. A basic load balancing mechanism allows a re-distribution of newly created blocks to the involved computing ranks. The use of a multi-hop network, to distribute information to the entire domain, avoids collective all-gather communications. Thus, we can demonstrate excellent scaling. The present scaling tests substantiate the enormous advantage of this adaptive method. For certain simulation scenarios, we can show that the calculation effort and memory consumption can be reduced to only 3.5 percent, compared to the classical full-domain reference simulation. The overhead of 70–100 percent for the dynamic adapting block creation is significantly lower than the gain. The approach is not restricted to phase-field simulations, and can be employed in other domains of computational science to exploit sparsity of computing regions.

Automatic Data Layout Transformations in the ExaStencils Code Generator

Performance optimizations should focus not only on the computations of an application, but also on the internal data layout. A well-known problem is whether a struct of arrays or an array of structs results in a higher performance for a particular application. Even though the switch from the one to the other is fairly simple to implement, testing both transformations can become laborious and error-prone. Additionally, there are more complex data layout transformations, such as a color splitting for multi-color kernels in the domain of stencil codes, that are manually difficult. As a remedy, we propose new flexible layout transformation statements for our domain-specific language ExaSlang that support arbitrary affine transformations. Since our code generator applies them automatically to the generated code, these statements enable the simple adaptation of the data layout without the need for any other modifications of the application code. This constitutes a big advance in the ease of testing and evaluating different memory layout schemes in order to identify the best.

Thoroughly Exploring GPU Buffering Options for Stencil Code by Using an Efficiency Measure and a Performance Model

Stencil computations form the basis for computer simulations across almost every field of science, such as computational fluid dynamics, data mining, and image processing. Their mostly regular data access patterns potentially enable them to take advantage of the high computation and data bandwidth of GPUs, but only if data buffering and other issues are handled properly. Finding a good code generation strategy presents a number of challenges, one of which is the best way to make use of memory. GPUs have several types of on-chip storage including registers, shared memory, and a read-only cache. The choice of type of storage and how it’s used, a buffering strategy , for each stencil array ( grid function , [GF]) not only requires a good understanding of its stencil pattern, but also the efficiency of each type of storage for the GF, to avoid squandering storage that would be more beneficial to another GF. For a stencil computation with $N$ GFs, the total number of possible assignments is $\beta ^{N}$ where $\beta$ is the number of buffering strategies. Our code-generation framework supports five buffering strategies ( $\beta =5$ ). Large, complex stencil kernels may consist of dozens of GFs, resulting in significant search overhead. In this work, we present an analytic performance model for stencil computations on GPUs and study the behavior of read-only cache and L2 cache. Next, we propose an efficiency-based assignment algorithm which operates by scoring a change in buffering strategy for a GF using a combination of (a) the predicted execution time and (b) on-chip storage usage. By using this scoring, an assignment for $N$ GFs can be determined in $(\beta -1)N(N+1)/2$ steps. Results show that the performance model has good accuracy and that the assignment strategy is highly efficient.

Automatic Matching of Legacy Code to Heterogeneous APIs

Heterogeneous accelerators often disappoint. They provide the prospect of great performance, but only deliver it when using vendor specific optimized libraries or domain specific languages. This requires considerable legacy code modifications, hindering the adoption of heterogeneous computing. This paper develops a novel approach to automatically detect opportunities for accelerator exploitation. We focus on calculations that are well supported by established APIs: sparse and dense linear algebra, stencil codes and generalized reductions and histograms. We call them idioms and use a custom constraint-based Idiom Description Language (IDL) to discover them within user code. Detected idioms are then mapped to BLAS libraries, cuSPARSE and clSPARSE and two DSLs: Halide and Lift. We implemented the approach in LLVM and evaluated it on the NAS and Parboil sequential C/C++ benchmarks, where we detect 60 idiom instances. In those cases where idioms are a significant part of the sequential execution time, we generate code that achieves 1.26x to over 20x speedup on integrated and external GPUs.

Special issue on high performance computing in modeling and simulation

Special issue on high performance computing in modeling and simulation

A high‐level library for multidimensional arrays programming in computational science

SummaryThis paper describes ExaShark, a hybrid n‐dimensional array toolkit offered as a high‐level library for scientists to compute large‐scale simulations. It offers a global‐array–like interface while its runtime can be configured to use shared memory threading techniques, inter‐node distribution techniques, or combinations of both. ExaShark takes advantage of the latest HPC technologies, helping to scale to future generation systems. It has been used to develop several scientific applications including stencil codes, solvers, and matrix factorization algorithms. These applications are used to demonstrate that it improves on the state of the art by providing a user‐friendly, generic API without sacrificing performance.

Strategy for data-flow synchronizations in stencil parallel computations on multi-/manycore systems

In this paper, an innovative strategy for the data-flow synchronization in shared-memory systems is proposed. This strategy assumes to synchronize only interdependent threads instead of using the barrier approach that—in contrast to our approach—synchronize all threads. We demonstrate the adaptation of the data-flow synchronization strategy to two complex scientific applications based on stencil codes. An algorithm for the data-flow synchronization is developed and successfully used for both applications. The proposed approach is evaluated for various Intel microarchitectures released in the last 5 years, including the newest processors: Skylake and Knights Landing. The important part of this assessment is the performance comparison of the proposed data-flow synchronization with the OpenMP barrier. The experimental results show that the performance of the studied applications can be accelerated up to 1.3 times using the proposed data-flow synchronizations strategy.

Diamond Tiling: Tiling Techniques to Maximize Parallelism for Stencil Computations

Most stencil computations allow tile-wise concurrent start, i.e., there always exists a face of the iteration space and a set of tiling directions such that all tiles along that face can be started concurrently. This provides load balance and maximizes parallelism. However, existing automatic tiling frameworks often choose hyperplanes that lead to pipelined start-up and load imbalance. We address this issue with a new tiling technique, called diamond tiling, that ensures concurrent start-up as well as perfect load-balance whenever possible. We first provide necessary and sufficient conditions for a set of tiling hyperplanes to allow concurrent start for programs with affine data accesses. We then provide an approach to automatically find such hyperplanes. Experimental evaluation on a 12-core Intel Westmere shows that diamond tiled code is able to outperform a tuned domain-specific stencil code generator by 10 to 40 percent, and previous compiler techniques by a factor of 1.3 $\times$ to 10.1 $\times$ .

Loop Tiling in Large-Scale Stencil Codes at Run-Time with OPS

The key common bottleneck in most stencil codes is data movement, and prior research has shown that improving data locality through optimisations that schedule across loops do particularly well. However, in many large PDE applications it is not possible to apply such optimisations through compilers because there are many options, execution paths and data per grid point, many dependent on run-time parameters, and the code is distributed across different compilation units. In this paper, we adapt the data locality improving optimisation called iteration space slicing for use in large OPS applications both in shared-memory and distributed-memory systems, relying on run-time analysis and delayed execution. We evaluate our approach on a number of applications, observing speedups of 2$\times$ on the Cloverleaf 2D/3D proxy application, which contain 83/141 loops respectively, $3.5\times$ on the linear solver TeaLeaf, and $1.7\times$ on the compressible Navier-Stokes solver OpenSBLI. We demonstrate strong and weak scalability up to 4608 cores of CINECA's Marconi supercomputer. We also evaluate our algorithms on Intel's Knights Landing, demonstrating maintained throughput as the problem size grows beyond 16GB, and we do scaling studies up to 8704 cores. The approach is generally applicable to any stencil DSL that provides per loop data access information.