Computational Kernels Research Articles

Programming in OpenCL and its advantages in a GPU Framework

Abstract: OpenCL is a framework used for building applications that mostly run on heterogenous platforms containing CPUs, GPUs, and DSPs. OpenCL provides an interface for parallel programming that can be used to take advantage of the GPUs high parallel computing power. Programmers who need complete control over the parallelization process and who are required to write portable heterogeneous code mostly use OpenCL. OpenCL views a processing unit as a collection of compute units which in turn are made up work items. According to OpenCL’s workflow and memory hierarchy, each work item is a thread as far in terms of control and memory model. A collection of work items is called a work group which is mapped to a compute unit. The language used to write “compute kernels” is called kernel language. OpenCL uses C/C++ to carry over the kernel computations done on the device. The host code specifies the kernel specifications that is needed for the computation of the device which includes creating buffers, calling kernels, mapping the memory back to CPU from device, etc. OpenCL also has specific optimization techniques that helps improve parallelization while computing on a GPU which results in better performance numbers

WeaveNet: Solution for Variable Input Sparsity Depth Completion

LIDARs produce depth measurements, which are relatively sparse when compared with cameras. Current state-of-the-art solutions for increasing the density of LIDAR-derived depth maps rely on training the models for specific input measurement density. This assumption can easily be violated. The goal of this work was to develop a solution capable of producing reasonably accurate depth predictions while using input with a very wide range of depth information densities. To that end, we defined a WeaveBlock capable of efficiently propagating depth information. To achieve this goal, WeaveBlocks utilize long and narrow horizontal and vertical convolution kernels together with MobileNet-inspired pointwise convolutions serving as computational kernels. In this paper, we present the WeaveNet architecture for guided (LIDAR and camera) and unguided (LIDAR only) depth completion as well as a non-standard network training procedure. We present the results of the network on the KITTI test and validation sets. We analyze the network performance at various levels of input sparsity by randomly removing between 0% and 99% of the LIDAR points from the network inputs, and in each case, we obtain reasonable quality output. Additionally, we show that our trained network weights can easily be reused with a different LIDAR sensor.

An efficient full-wavefield computational approach for seismic testing in a layered half-space

Inversion of subsurface shear-wave velocity profile based on effective dispersion curve or full wavefield is still rare due to the lack of efficient forward computation kernel for the full wavefield in a layered medium. To facilitate exploiting the complete signal content in seismic testing, an efficient new computational framework for modeling synthetic full wavefields in horizontally layered systems is proposed. Dynamic response of an isotropic elastic layered system subjected to vertical loading is implemented and thoroughly validated. The new procedure is based the Fourier-Bessel series expansion (not the traditional Fourier-Bessel transform), and it is at least two-orders faster than the current full wavefield computation method based on integral transformation. It considers survey parameters and models all wave phenomena, including near-field effect and leaky waves, which cannot be captured by the modal summation method. Examples of full wavefield simulation are presented to demonstrate different levels of information provided by the experimentally available fundamental mode, broadband effective dispersion curve, and the full velocity spectrum. Different earth models may have similar fundamental modes or even effective dispersion curves; however, their difference manifests in different higher-mode contributions and can only be exploited by utilizing the entire velocity spectra including the vertical and radial components. The coupling between the vertical and radial components in terms of their spectral ratio and phase difference can provide additional constraints (data) in seismic testing.

Implementation of symmetry-adapted perturbation theory based on density functional theory and using hybrid exchange-correlation kernels for dispersion terms.

We report the implementation of a symmetry-adapted perturbation theory algorithm based on a density functional theory [SAPT(DFT)] description of monomers. The implementation adopts a density-fitting treatment of hybrid exchange-correlation kernels to enable the description of monomers with hybrid functionals, as in the algorithm by Bukowski, Podeszwa, and Szalewicz [Chem. Phys. Lett. 414, 111 (2005)]. We have improved the algorithm by increasing numerical stability with QR factorization and optimized the computation of the exchange-correlation kernel with its 2-index density-fitted representation. The algorithm scales as O(N5) formally and is usable for systems with up to ∼3000 basis functions, as demonstrated for the C60-buckycatcher complex with the aug-cc-pVDZ basis set. The hybrid-kernel-based SAPT(DFT) algorithm is shown to be as accurate as SAPT(DFT) implementations based on local effective exact exchange potentials obtained from the local Hartree-Fock (LHF) method while avoiding the lower-scaling [O(N4)] but iterative and sometimes hard-to-converge LHF process. The hybrid-kernel algorithm outperforms Hartree-Fock-based SAPT (SAPT0) for the S66 test set, and its accuracy is comparable to the many-body perturbation theory based SAPT2+ approach, which scales as O(N7), although SAPT2+ exhibits a more narrow distribution of errors.

Deep Learning Approaches to Source Code Analysis for Optimization of Heterogeneous Systems: Recent Results, Challenges and Opportunities

To cope with the increasing complexity of digital systems programming, deep learning techniques have recently been proposed to enhance software deployment by analysing source code for different purposes, ranging from performance and energy improvement to debugging and security assessment. As embedded platforms for cyber-physical systems are characterised by increasing heterogeneity and parallelism, one of the most challenging and specific problems is efficiently allocating computational kernels to available hardware resources. In this field, deep learning applied to source code can be a key enabler to face this complexity. However, due to the rapid development of such techniques, it is not easy to understand which of those are suitable and most promising for this class of systems. For this purpose, we discuss recent developments in deep learning for source code analysis, and focus on techniques for kernel mapping on heterogeneous platforms, highlighting recent results, challenges and opportunities for their applications to cyber-physical systems.

Towards Efficient Sparse Matrix Vector Multiplication on Real Processing-In-Memory Architectures

Several manufacturers have already started to commercialize near-bank Processing-In-Memory (PIM) architectures, after decades of research efforts. Near-bank PIM architectures place simple cores close to DRAM banks. Recent research demonstrates that they can yield significant performance and energy improvements in parallel applications by alleviating data access costs. Real PIM systems can provide high levels of parallelism, large aggregate memory bandwidth and low memory access latency, thereby being a good fit to accelerate the Sparse Matrix Vector Multiplication (SpMV) kernel. SpMV has been characterized as one of the most significant and thoroughly studied scientific computation kernels. It is primarily a memory-bound kernel with intensive memory accesses due its algorithmic nature, the compressed matrix format used, and the sparsity patterns of the input matrices given. This paper provides the first comprehensive analysis of SpMV on a real-world PIM architecture, and presents SparseP, the first SpMV library for real PIM architectures. We make two key contributions. First, we design efficient SpMV algorithms to accelerate the SpMV kernel in current and future PIM systems, while covering a wide variety of sparse matrices with diverse sparsity patterns. Second, we provide the first comprehensive analysis of SpMV on a real PIM architecture. Specifically, we conduct our rigorous experimental analysis of SpMV kernels in the UPMEM PIM system, the first publicly-available real-world PIM architecture. Our extensive evaluation provides new insights and recommendations for software designers and hardware architects to efficiently accelerate the SpMV kernel on real PIM systems. For more information about our thorough characterization on the SpMV PIM execution, results, insights and the open-source SparseP software package [21], we refer the reader to the full version of the paper [3, 4]. The SparseP software package is publicly and freely available at https://github.com/CMU-SAFARI/SparseP.

Design of Processing-in-Memory With Triple Computational Path and Sparsity Handling for Energy-Efficient DNN Training

As machine learning (ML) and artificial intelligence (AI) have become mainstream technologies, many accelerators have been proposed to cope with their computation kernels. However, they access the external memory frequently due to the large size of deep neural network model, suffering from the von Neumann bottleneck. Moreover, as privacy issue is becoming more critical, on-device training is emerging as its solution. However, on-device training is challenging because it should perform the training under a limited power budget, which requires a lot more computations and memory accesses than the inference. In this paper, we present an energy-efficient processing-in-memory (PIM) architecture supporting end-to-end on-device training named T-PIM. Its macro design includes an 8T-SRAM cell-based PIM block to compute in-memory <italic xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">AND</i> operation and three computational datapaths for end-to-end training. Each of three computational paths integrates arithmetic units for forward propagation, backward propagation, and gradient calculation and weight update, respectively, allowing the weight data stored in the memory stationary. T-PIM also supports variable bit precision to cover various ML scenarios. It can use fully variable input bit precision and 2-bit, 4-bit, 8-bit, and 16-bit weight bit precision for the forward propagation and the same input bit precision and 16-bit weight bit precision for the backward propagation. In addition, T-PIM implements sparsity handling schemes that skip the computation for input data and turn off the arithmetic units for weight data to reduce both unnecessary computations and leakage power. Finally, we fabricate the T-PIM chip on a 5.04mm <sup xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">2</sup> die in a 28-nm CMOS logic process. It operates at 50–280MHz with the supply voltage of 0.75–1.05V, dissipating 5.25–51.23mW power in inference and 6.10–37.75mW in training. As a result, it achieves 17.90–161.08TOPS/W energy efficiency for the inference of 1-bit activation and 2-bit weight data, and 0.84–7.59TOPS/W for the training of 8-bit activation/error and 16-bit weight data. In conclusion, T-PIM is the first PIM chip that supports end-to-end training, demonstrating 2.02 times performance improvement over the latest PIM that partially supports training.

Foreword to the Special Section on Smart Tools and Applications in Graphics (STAG 2021)

The Special Section contains extended and revised versions of the best papers presented at the 8th Conference on Smart Tools and Applications in Graphics (STAG 2021), held virtually on October 28–29, 2021. Five papers were selected by appointed members from the Program Committee; extended versions were submitted and further reviewed by external experts. The result is a collection of paper spanning different Visual Computing domains: from vector graphics on surfaces to image vectorization, from 3D geometric modelling of geological data to computation of the geometric kernel of polyhedra, up to lens-based scene exploration of annotated visual data.

Auto-Selection of an Optimal Sparse Matrix Format in the Neuro-Simulator ANNarchy.

Modern neuro-simulators provide efficient implementations of simulation kernels on various parallel hardware (multi-core CPUs, distributed CPUs, GPUs), thereby supporting the simulation of increasingly large and complex biologically realistic networks. However, the optimal configuration of the parallel hardware and computational kernels depends on the exact structure of the network to be simulated. For example, the computation time of rate-coded neural networks is generally limited by the available memory bandwidth, and consequently, the organization of the data in memory will strongly influence the performance for different connectivity matrices. We pinpoint the role of sparse matrix formats implemented in the neuro-simulator ANNarchy with respect to computation time. Rather than asking the user to identify the best data structures required for a given network and platform, such a decision could also be carried out by the neuro-simulator. However, it requires heuristics that need to be adapted over time for the available hardware. The present study investigates how machine learning methods can be used to identify appropriate implementations for a specific network. We employ an artificial neural network to develop a predictive model to help the developer select the optimal sparse matrix format. The model is first trained offline using a set of training examples on a particular hardware platform. The learned model can then predict the execution time of different matrix formats and decide on the best option for a specific network. Our experimental results show that using up to 3,000 examples of random network configurations (i.e., different population sizes as well as variable connectivity), our approach effectively selects the appropriate configuration, providing over 93% accuracy in predicting the suitable format on three different NVIDIA devices.

MemHC: An Optimized GPU Memory Management Framework for Accelerating Many-body Correlation

The many-body correlation function is a fundamental computation kernel in modern physics computing applications, e.g., Hadron Contractions in Lattice quantum chromodynamics (QCD). This kernel is both computation and memory intensive, involving a series of tensor contractions, and thus usually runs on accelerators like GPUs. Existing optimizations on many-body correlation mainly focus on individual tensor contractions (e.g., cuBLAS libraries and others). In contrast, this work discovers a new optimization dimension for many-body correlation by exploring the optimization opportunities among tensor contractions. More specifically, it targets general GPU architectures (both NVIDIA and AMD) and optimizes many-body correlation’s memory management by exploiting a set of memory allocation and communication redundancy elimination opportunities: first, GPU memory allocation redundancy : the intermediate output frequently occurs as input in the subsequent calculations; second, CPU-GPU communication redundancy : although all tensors are allocated on both CPU and GPU, many of them are used (and reused) on the GPU side only, and thus, many CPU/GPU communications (like that in existing Unified Memory designs) are unnecessary; third, GPU oversubscription: limited GPU memory size causes oversubscription issues, and existing memory management usually results in near-reuse data eviction, thus incurring extra CPU/GPU memory communications. Targeting these memory optimization opportunities, this article proposes MemHC, an optimized systematic GPU memory management framework that aims to accelerate the calculation of many-body correlation functions utilizing a series of new memory reduction designs. These designs involve optimizations for GPU memory allocation, CPU/GPU memory movement, and GPU memory oversubscription, respectively. More specifically, first, MemHC employs duplication-aware management and lazy release of GPU memories to corresponding host managing for better data reusability. Second, it implements data reorganization and on-demand synchronization to eliminate redundant (or unnecessary) data transfer. Third, MemHC exploits an optimized Least Recently Used (LRU) eviction policy called Pre-Protected LRU to reduce evictions and leverage memory hits. Additionally, MemHC is portable for various platforms including NVIDIA GPUs and AMD GPUs. The evaluation demonstrates that MemHC outperforms unified memory management by \( 2.18\times \) to \( 10.73\times \) . The proposed Pre-Protected LRU policy outperforms the original LRU policy by up to \( 1.36\times \) improvement. 1

Quantitative uncertainty estimation in biophysical models of fish larval connectivity in the Florida Keys

Abstract The impacts of seven uncertain biological parameters on simulated larval connectivity in the Florida Keys are investigated using Polynomial chaos surrogates. These parameters describe biological traits and behaviours—such as mortality, swimming abilities, and orientation—and modulate larval settlement as well as dispersal forecasts. However, these parameters are poorly constrained by observations and vary naturally between individual larvae. The present investigation characterizes these input uncertainties with probability density functions informed by previous studies of Abudefduf saxatilis. The parametric domain is sampled via ensemble calculations, then a polynomial-based surrogate is built to explicitly approximate the dependence of the model outputs on the uncertain model inputs, which enables a robust statistical analysis of uncertainties. This approach allows the computation of probabilistic dispersal kernels that are further analyzed to understand the impact of the parameter uncertainties. We find that the biological input parameters influence the connectivity differently depending on dispersal distance and release location. The global sensitivity analysis shows that the interactions between detection distance threshold, orientation ontogeny, and orientation accuracy, are the dominant contributors to the uncertainty in settlement abundance in the Florida Keys. Uncertainties in swimming speed and mortality, on the other hand, seem to contribute little to dispersal uncertainty.

Numerical evaluation of time-distance helioseismic sensitivity kernels in spherical geometry

Context. Helioseismic analysis of large-scale flows and structural inhomogeneities in the Sun requires the computation of sensitivity kernels that account for the spherical geometry of the Sun, as well as systematic effects such as line-of-sight projection. Aims. I aim to develop a code to evaluate helioseismic sensitivity kernels for flows using line-of-sight projected measurements. Methods. I decomposed the velocity field in a basis of vector spherical harmonics and computed the kernel components corresponding to the coefficients of velocity in this basis. The kernels thus computed are radial functions that set up a 1.5D inverse problem to infer the flow from surface measurements. I demonstrate that using the angular momentum addition formalism lets us express the angular dependence of the kernels as bipolar spherical harmonics, which may be evaluated accurately and efficiently. Results. Kernels for line-of-sight projected measurements may differ significantly from those that don’t account for projection. Including projection in our analysis does not increase the computational time significantly. We demonstrate that it is possible to evaluate kernels for pairs of points that are related through a rotation by linearly transforming the terms that enter the expression of the kernel, and that this result holds even for line-of-sight projected kernels. Conclusions. I developed a Julia code that may be used to evaluate sensitivity kernels for seismic wave travel times computed using line-of-sight projected measurements, which is made freely available under the MIT license.

Detection and computation of conservative kernels of models consisting of freeform curves and surfaces, using inequality constraints

We present an algorithm to compute a tight-as-needed conservative approximation of the kernel domain of freeform curves in R2 and freeform surfaces in R3. Inequality constraints to detect the interior of the kernel domain are formulated as multivariates, and solved with a subdivision-based approach to find the domains in R2 or R3 that satisfy the inequalities and are in the kernel. The convex hull of the computed domains is also included in the kernel, and adopted as the approximated kernel domain. We can apply the presented algorithm to detect the kernel domain of not only C1 continuous closed regular curves and surfaces, but also the kernel domains of multiple piecewise C1 continuous regular freeform curves and surfaces. Further, the presented algorithm can be applied to find the gamma-kernel as well as the kernel domain of open curves and surfaces, under some assumptions. We demonstrate our experimental result using various freeform curves and surfaces, and compare it with the kernel computation algorithm presented in Elber et al. (2006).

A Multi-Bit Non-Volatile Compute-in-Memory Architecture with Quantum-Dot Transistor Based Unit

The recent advance of artificial intelligence (AI) has shown remarkable success for numerous tasks, such as cloud computing, deep-learning, neural network and so on. Most of those applications rely on fast computation and large storage, which brings various challenges to the hardware platform. The hardware performance is the bottle neck to break through and therefore, there is a lot of interest in exploring new solutions for computation architecture in recent years. Compute-in-memory (CIM) has drawn attention to the researchers and it is considered as one of the most promising candidates to solve the above challenges. Computing-In-memory is an emerging technique to fulfill the fast-growing demand for high-performance data processing. This technique offers fast processing, low power and high performance by blurring the boundary between processing cores and memory units. One key aspect of CIM is performing matrix-vector multiplication (MVM) or dot product operation through intertwining of processing and memory elements. As the primary computational kernel in neural networks, dot product operation is targeted to be improved in terms of its performance. In this paper, we present the design, implementation and analysis of quantum-dot transistor (QDT) based CIM, from the multi-bit multiplier to the dot product unit, and then the in-memory computing array.

SparseP

Several manufacturers have already started to commercialize near-bank Processing-In-Memory (PIM) architectures, after decades of research efforts. Near-bank PIM architectures place simple cores close to DRAM banks. Recent research demonstrates that they can yield significant performance and energy improvements in parallel applications by alleviating data access costs. Real PIM systems can provide high levels of parallelism, large aggregate memory bandwidth and low memory access latency, thereby being a good fit to accelerate the Sparse Matrix Vector Multiplication (SpMV) kernel. SpMV has been characterized as one of the most significant and thoroughly studied scientific computation kernels. It is primarily a memory-bound kernel with intensive memory accesses due its algorithmic nature, the compressed matrix format used, and the sparsity patterns of the input matrices given. This paper provides the first comprehensive analysis of SpMV on a real-world PIM architecture, and presents SparseP, the first SpMV library for real PIM architectures. We make three key contributions. First, we implement a wide variety of software strategies on SpMV for a multithreaded PIM core, including (1) various compressed matrix formats, (2) load balancing schemes across parallel threads and (3) synchronization approaches, and characterize the computational limits of a single multithreaded PIM core. Second, we design various load balancing schemes across multiple PIM cores, and two types of data partitioning techniques to execute SpMV on thousands of PIM cores: (1) 1D-partitioned kernels to perform the complete SpMV computation only using PIM cores, and (2) 2D-partitioned kernels to strive a balance between computation and data transfer costs to PIM-enabled memory. Third, we compare SpMV execution on a real-world PIM system with 2528 PIM cores to an Intel Xeon CPU and an NVIDIA Tesla V100 GPU to study the performance and energy efficiency of various devices, i.e., both memory-centric PIM systems and conventional processor-centric CPU/GPU systems, for the SpMV kernel. SparseP software package provides 25 SpMV kernels for real PIM systems supporting the four most widely used compressed matrix formats, i.e., CSR, COO, BCSR and BCOO, and a wide range of data types. SparseP is publicly and freely available at https://github.com/CMU-SAFARI/SparseP. Our extensive evaluation using 26 matrices with various sparsity patterns provides new insights and recommendations for software designers and hardware architects to efficiently accelerate the SpMV kernel on real PIM systems.

Sortledton

Despite the wide adoption of graph processing across many different application domains, there is no underlying data structure that can serve a variety of graph workloads (analytics, traversals, and pattern matching) on dynamic graphs with transactional updates. In this paper, we present Sortledton, a universal graph data structure that addresses the open problem by being carefully optimizing for the most relevant data access patterns used by graph computation kernels. It can support millions of transactional updates per second, while providing competitive performance (1.22x on average) for the most common graph workloads to the best-known baseline for static graphs - CSR. With this, we improve the ingestion throughput over state-of-the-art dynamic graph data structures, while supporting a wider range of graph computations under transactional guarantees, with a much simpler design and significantly smaller memory footprint (2.1x that of CSR).

PFOB sonosensitive microdroplets: determining their interaction radii with focused ultrasound using MR thermometry and a Gaussian convolution kernel computation

Purpose: Micron-sized perfluorocarbon droplet adjuvants to focused ultrasound therapies allow lower applied power, circumvent unwanted prefocal heating, and enhance thermal dose in highly perfused tissues. The heat enhancement has been shown to saturate at increasing concentrations. Experiments were performed to empirically model the saturating heating effects during focused ultrasound. Materials and methods: The measurements were made at varying concentrations using magnetic resonance thermometry and focused ultrasound by circulating droplets of mean diameter 1.9 to 2.3 µm through a perfused phantom. A simulation was performed to estimate the interaction radius size, empirically. Results: The interaction radius, representing the radius of a sphere encompassing 90% of the probability for the transformation of acoustic energy into heat deposition around a single droplet, was determined experimentally from ultrasonic absorption coefficient measurements The simulations suggest the interaction radius was approximately 12.5-fold larger than the geometrical radius of droplets, corresponding to an interaction volume on the order of 2000 larger than the geometrical volume. Conclusions: The results provide information regarding the dose–response relationship from the droplets, a measure with 15% precision of their interaction radii with focused ultrasound, and subsequent insights into the underlying physical heating mechanism.

Сравнение инструментов высокоуровневого синтеза и конструирования цифровой аппаратуры

Application-specific systems with FPGA accelerators are often designed using high-level synthesis or hardware construction tools. Nowadays, there are many frameworks available, both open-source and commercial. In this work, we attempt to fairly compare several existing solutions (languages and tools), including Verilog (our baseline), Chisel, Bluespec SystemVerilog (Bluespec Compiler), DSLX (XLS), MaxJ (MaxCompiler), and C (Bambu and Vivado HLS). Our analysis has been carried out using a representative example of 8×8 inverse discrete cosine transform (IDCT), a widely used algorithm engaged in, among others, JPEG and MPEG decoders. The metrics under consideration include: (a) the degree of automation (how much less code is required compared to Verilog), (b) the controllability (possibility to achieve given design characteristics, namely a given ratio of the performance and area), and (c) the flexibility (ease of design modification to achieve certain characteristics). Rather than focusing on computational kernels only, we have developed AXI-Stream wrappers for the synthesized implementations, which allows adequately evaluating characteristics of the designs when they are used as parts of real computer systems. Our study shows clear examples of what impact specific optimizations (tool settings and source code modifications) have on the overall system performance and area. It emphasizes how important is to be able to control the balance between the communication interface utilization and the computational kernel performance and delivers clear guidelines for the next generation tools for designing FPGA accelerator based systems.

On exponentially convergent parameter estimation with lack of persistency of excitation

Given a linear regression model, the problem of exponentially convergent identification of a vector of p constant parameters is addressed, without the assumption that the regressor is persistently exciting. A new definition of lack of persistency of excitation of order q , with 0 ≤ q ≤ p , is introduced, which reduces to the well-known persistency of excitation condition in the special case q = 0 . This definition implies that the kernel of the Gram matrix of the regressor has constant dimension q . Gradient and least squares identification algorithms may be modified on the basis of such a kernel computation, so that exponential convergence and robustness of parameter estimation errors is achieved both for gradient and least squares algorithms, even without persistency of excitation.

Fast Kernel Smoothing in R with Applications to Projection Pursuit

This paper introduces the R package FKSUM, which offers fast and exact evaluation of univariate kernel smoothers. The main kernel computations are implemented in C++, and are wrapped in simple, intuitive and versatile R functions. The fast kernel computations are based on recursive expressions involving the order statistics, which allows for exact evaluation of kernel smoothers at all sample points in log-linear time. In addition to general purpose kernel smoothing functions, the package offers purpose built and readyto-use implementations of popular kernel-type estimators. On top of these basic smoothing problems, this paper focuses on projection pursuit problems in which the projection index is based on kernel-type estimators of functionals of the projected density.