GPU Implementation Research Articles

Path integral radiative transfer via polyline representation allowing GPU implementation

Many research areas require simulating particle transport through scattering media. For instance, radiative transfer is useful for computer graphics and for neutron transport in nuclear physics. These transport simulations tend to be computationally expensive for problems involving large amounts of multiple scattering in generic geometries, requiring significant time to compute. Finding a fast solution for these types of problems remains an open area for research. Previous work shows that radiative transfer can be represented as a Feynman Path Integral over all paths between two points in a space. The path integral assigns a weight to each path based on the local curvature of the path, accumulating a transport kernel by summing the weights of all of the paths. Previous work demonstrated a Monte Carlo method for computing the radiative transfer Feynman path integral via repeatedly perturbing paths to generate new paths, using a discrete Frenet-Serret framework and an expensive root-solve calculation. The approach is highly parallelizable on the CPU, however computations require a supercomputer to complete in reasonable time. While a GPU implementation was considered for this previous approach, the root-solve and path representation mapped poorly to the GPU, hindering a reasonable implementation. In the present work, we propose a representation of paths as polylines, and a new curve perturbation method that guarantees production of a new unique path each execution while maintaining the boundary conditions of each path, but without the time-consuming root-solve. The new perturbation method’s structure maps well to the GPU, allowing implementation in CUDA which outperforms the CPU, and allows the utilization of GPU hardware.

Novel Perlin-based Phantoms Using 3D Models of Compressed Breast Shape and Fractal Noise.

Virtual clinical trials (VCTs) have been used widely to evaluate digital breast tomosynthesis (DBT) systems. VCTs require realistic simulations of the breast anatomy (phantoms) to characterize lesions and to estimate risk of masking cancers. This study introduces the use of Perlin-based phantoms to optimize the acquisition geometry of a novel DBT prototype. These phantoms were developed using a GPU implementation of a novel library called Perlin-CuPy. The breast anatomy is simulated using 3D models under mammography cranio-caudal compression. In total, 240 phantoms were created using compressed breast thickness, chest-wall to nipple distance, and skin thickness that varied in a {[35, 75], [59, 130), [1.0, 2.0]} mm interval, respectively. DBT projections and reconstructions of the phantoms were simulated using two acquisition geometries of our DBT prototype. The performance of both acquisition geometries was compared using breast volume segmentations of the Perlin phantoms. Results show that breast volume estimates are improved with the introduction of posterior-anterior motion of the x-ray source in DBT acquisitions. The breast volume is overestimated in DBT, varying substantially with the acquisition geometry; segmentation errors are more evident for thicker and larger breasts. These results provide additional evidence and suggest that custom acquisition geometries can improve the performance and accuracy in DBT. Perlin phantoms help to identify limitations in acquisition geometries and to optimize the performance of the DBT prototypes.

Scalable communication for high-order stencil computations using CUDA-aware MPI

Modern compute nodes in high-performance computing provide a tremendous level of parallelism and processing power. However, as arithmetic performance has been observed to increase at a faster rate relative to memory and network bandwidths, optimizing data movement has become critical for achieving strong scaling in many communication-heavy applications. This performance gap has been further accentuated with the introduction of graphics processing units, which can provide by multiple factors higher throughput in data-parallel tasks than central processing units. In this work, we explore the computational aspects of iterative stencil loops and implement a generic communication scheme using CUDA-aware MPI, which we use to accelerate magnetohydrodynamics simulations based on high-order finite differences and third-order Runge–Kutta integration. We put particular focus on improving intra-node locality of workloads. Our GPU implementation scales strongly from one to 64 devices at 50%–87% of the expected efficiency based on a theoretical performance model. Compared with a multi-core CPU solver, our implementation exhibits 20–60× speedup and 9–12× improved energy efficiency in compute-bound benchmarks on 16 nodes.

GPU implementation of Explicit and Implicit Eulerian methods with TVD schemes for solving 2D solute transport in heterogeneous flows

In this work we present an efficient implementation of Eulerian TVD methods. We apply parallelization strategies based entirely on GPU for the solution of the 2D transport equation in heterogeneous porous media. Additionally, a parallel strategy is proposed for the generation of exponentially correlated lognormally distributed permeability fields in GPU. The programs are developed using C++/CUDA. The implemented methods are used to solve advective dominant problems, in a context of Monte Carlo type simulations to numerically determine the longitudinal and transversal macrodispersion coefficients averaging over 100 simulations for permeability fields for a large range of variances. The following types of transport are considered for testing: pure advection, advection-diffusion and advection-dispersion. The performance in terms of the computation time of explicit and implicit methods are compared. We show that the implemented algorithms allow to efficiently solve problems in computational domains of up to 134.5 million cells in a single GPU.

JefiGPU: Jefimenko's equations on GPU

We have implemented a GPU version of the Jefimenko's equations — JefiGPU. Given the proper distributions of the source terms ρ (charge density) and J (current density) in the source volume, the algorithm gives the electromagnetic fields in the observational region (not necessarily overlaps the vicinity of the sources). To verify the accuracy of the GPU implementation, we have compared the obtained results with that of the theoretical ones. Our results show that the deviations of the GPU results from the theoretical ones are around 5%. Meanwhile, we have also compared the performance of the GPU implementation with a CPU version. Finally, we have studied the parameter dependence of the execution time and memory consumption on one NVIDIA Tesla V100 card. Our code can be consistently coupled to RBG (Relativistic Boltzmann equations on GPUs) and many other GPU-based algorithms in physics. Program summaryProgram Title: JefiGPUCPC Library link to program files:https://doi.org/10.17632/6vnghzknx4.1Developer's repository link:https://github.com/Juenjie/CRBMGCode Ocean capsule:https://codeocean.com/capsule/9007956Licensing provisions: Apache-2.0Programming language: PythonNature of problem: Jefimenko's equations are numerically stable with given sources ρ (charge density) and J (current density). Providing proper sources, Jefimenko's equations can give the electromagnetic fields without boundary conditions. However, the relevant integrations of the Jefimenko's equations involve a 4D space-time which is both memory and time consuming. With the help of the state-of-art GPU technique, these integrations can be evaluated with acceptable accuracy and execution time. In the present work, we have implemented a GPU version of the Jefimenko's equations and tested the accuracy and the performance of the code.Solution method: The GPU code mainly deals with the implementation of the integrations in Jefimenko's equations. Each CUDA kernel corresponds to a spatial point of E and B. The integrations involving the source region and the retarded time are performed in each CUDA kernel. Ray, Numba and Cupy packages are used to manipulate the scaling of GPU allocation and evaluation.

A Unified Software/Hardware Scalable Architecture for Brain-Inspired Computing Based on Self-Organizing Neural Models.

The field of artificial intelligence has significantly advanced over the past decades, inspired by discoveries from the fields of biology and neuroscience. The idea of this work is inspired by the process of self-organization of cortical areas in the human brain from both afferent and lateral/internal connections. In this work, we develop a brain-inspired neural model associating Self-Organizing Maps (SOM) and Hebbian learning in the Reentrant SOM (ReSOM) model. The framework is applied to multimodal classification problems. Compared to existing methods based on unsupervised learning with post-labeling, the model enhances the state-of-the-art results. This work also demonstrates the distributed and scalable nature of the model through both simulation results and hardware execution on a dedicated FPGA-based platform named SCALP (Self-configurable 3D Cellular Adaptive Platform). SCALP boards can be interconnected in a modular way to support the structure of the neural model. Such a unified software and hardware approach enables the processing to be scaled and allows information from several modalities to be merged dynamically. The deployment on hardware boards provides performance results of parallel execution on several devices, with the communication between each board through dedicated serial links. The proposed unified architecture, composed of the ReSOM model and the SCALP hardware platform, demonstrates a significant increase in accuracy thanks to multimodal association, and a good trade-off between latency and power consumption compared to a centralized GPU implementation.

CuNH: Efficient GPU Implementations of Post-Quantum KEM NewHope

Post-quantum cryptography was proposed in the past years due to the foreseeable emergence of quantum computers that are able to break the conventional public key cryptosystems at acceptable costs. However, post-quantum schemes are usually less efficient than conventional ones, which makes them less practical in scenarios with limited resources or high concurrency. Server-side applications always feature multiple users, therefore requiring efficient execution of batch tasks. GPU is intrinsically well-suited to batch tasks owing to its SIMD/SIMT execution fashion, so it naturally helps to achieve high performance. However, a naive GPU-based implementation cannot make the best use of hardware resources of the GPU regardless of task loads. In this article, we propose SIMD parallelization paradigms for fine-grained GPU implementations and then apply them to a post-quantum key encapsulation algorithm called NewHope, where we carefully design every module, especially NTT and inverse NTT, to fit into the SIMD parallelization paradigms. In addition, we employ multi-streaming to improve performance in user's perspective. Finally, our evaluations are made on two testbeds with GPU accelerators NVIDIA GeForce MX150 and GeForce GTX 1650, respectively. The experimental results show that the fine-grained implementations save up to 98 percent latency at low task loads, and their throughputs increase by up to 86 percent at high task loads, when compared with the naive ones in kernel's perspective, and the multi-streaming implementations greatly reduce the latency overhead percentage at high task loads by up to 86 percent, when compared with the fine-grained implementation in user's perspective. Moreover, our fine-grained implementation and multi-streaming implementation are respectively 51.5 and 45.5 percent faster than Gupta et al.'s implementations when compared with it under reasonable assumptions. Furthermore, as lattice-based post-quantum schemes have similar operations, our proposal also easily applies to other lattice-based post-quantum schemes.

Weakly-supervised learning for catheter segmentation in 3D frustum ultrasound.

Accurate and efficient catheter segmentation in 3D ultrasound (US) is essential for ultrasound-guided cardiac interventions. State-of-the-art segmentation algorithms, based on convolutional neural networks (CNNs), suffer from high computational cost and large 3D data size for GPU implementation, which are far from satisfactory for real-time applications. In this paper, we propose a novel approach for efficient catheter segmentation in 3D US. Instead of using Cartesian US, our approach performs catheter segmentation in Frustum US (i.e., the US data before scan conversion). Compared to Cartesian US, Frustum US has a much smaller volume size, therefore the catheter can be segmented more efficiently in Frustum US. However, annotating the irregular and deformed Frustum images is challenging, and it is laborious to obtain the voxel-level annotation. To address this, we propose a weakly supervised learning framework, which requires only bounding-box annotations. The labels of the voxels are generated by incorporating class activation maps with line filtering, which are iteratively updated during the training cycles. Our experimental results show that, compared to Cartesian US, the catheter can be segmented much more efficiently in Frustum US (i.e., 0.25s per volume) with better accuracy. Extensive experiments also validate the effectiveness of the proposed weakly supervised learning method.

Energy-efficient distributed password hash computation on heterogeneous embedded system

This paper presents the improved version of our cool Cracker cluster (cCc), a heterogeneous distributed system for parallel and energy-efficient bcrypt password hash computation. The cluster consists of up to 8 computational units (nodes) with different performances measured in bcrypt hash computations per second [H/s]. In the cluster, nodes are low-power heterogeneous embedded systems with programmable logic containing specialized hash computation accelerators. In the experiments, we used a combination of Xilinx Zynq-series SoC boards and ZTEX 1.15y board which was initially used as a bitcoin miner. Zynq based nodes use the improved version of our custom bcrypt accelerator, which executes the most costly parts of the bcrypt hash computation in programmable logic. The cluster was formed around the famous open-source password cracking software package John the Ripper (abbr. JtR). On the communication layer, we used Message Passing Interface (MPI)library with a standard Ethernet network connecting the nodes. To mitigate the different performances among the cluster nodes and to balance the load, we developed and implemented password candidate distribution scheme based on the passwords' probability distribution, i.e. the order of appearance in the dictionary. We tested individual nodes and the cluster as a whole, trying different combinations of nodes and evaluating our distribution scheme for password candidates. We also compared our cluster with various GPU implementations in terms of performance, energy-efficiency, and price-efficiency. We show that our solution outperforms other platforms such as high-end GPUs, by a factor of at least 3 in terms of energy-efficiency and thus producing less overall cost of password attack than other platforms. In terms of the total operational costs, our cluster pays off after 4500 cracked passwords for a bcrypt hash with cost parameter 12, which makes it more appealing for real-world password-based system attacks. We also demonstrate the scalability of our cCc cluster.

Highly Efficient and Accurate Deep Learning-Based Classification of MRI Contrast on a CPU and GPU.

Classifying MR images based on their contrast mechanism can be useful in image segmentation where additional information from different contrast mechanisms can improve intensity-based segmentation and help separate the class distributions. In addition, automated processing of image type can be beneficial in archive management, image retrieval, and staff training. Different clinics and scanners have their own image labeling scheme, resulting in ambiguity when sorting images. Manual sorting of thousands of images would be a laborious task and prone to error. In this work, we used the power of transfer learning to modify pretrained residual convolution neural networks to classify MRI images based on their contrast mechanisms. Training and validation were performed on a total of 5169 images belonging to 10 different classes and from different MRI vendors and field strengths. Time for training and validation was 36min. Testing was performed on a different data set with 2474 images. Percentage of correctly classified images (accuracy) was 99.76%. (A deeper version of the residual network was trained for 103min and showed slightly lower accuracy of 99.68%.) In consideration of model deployment in the real world, performance on a single CPU computer was compared with GPU implementation. Highly accurate classification, training, and testing can be achieved without use of a GPU in a relatively short training time, through proper choice of a convolutional neural network and hyperparameters, making it feasible to improve accuracy by repeated training with cumulative training sets. Techniques to improve accuracy further are discussed and demonstrated. Derived heatmaps indicate areas of image used in decision making and correspond well with expert human perception. The methods used can be easily extended to other classification tasks with minimal changes.

High performance sparse multifrontal solvers on modern GPUs

We have ported the numerical factorization and triangular solve phases of the sparse direct solver STRUMPACK to GPU. STRUMPACK implements sparse LU factorization using the multifrontal algorithm, which performs most of its operations in dense linear algebra operations on so-called frontal matrices of various sizes. Our GPU implementation off-loads these dense linear algebra operations, as well as the sparse scatter–gather operations between frontal matrices. For the larger frontal matrices, our GPU implementation relies on vendor libraries such as cuBLAS and cuSOLVER for NVIDIA GPUs and rocBLAS and rocSOLVER for AMD GPUs. For the smaller frontal matrices we developed custom CUDA and HIP kernels to reduce kernel launch overhead. Overall, high performance is achieved by identifying submatrix factorizations corresponding to sub-trees of the multifrontal assembly tree which fit entirely in GPU memory. The multi-GPU setting uses SLATE (Software for Linear Algebra Targeting Exascale) as a modern GPU-aware replacement for ScaLAPACK. On 4 nodes of SUMMIT the code runs ∼10× faster when using all 24 V100 GPUs compared to when it only uses the 168 POWER9 cores. On 8 SUMMIT nodes, using 48 V100 GPUs, the sparse solver reaches over 50TFlop/s. Compared to SuperLU, on a single V100, for a set of 17 matrices our implementation is faster for all but one matrix, and is on average 5× (median 4×) faster

Analysis of the Leaky Integrate-and-Fire neuron model for GPU implementation

Understanding how neurons perform, when they are organized in interacting networks, is a key to understanding how the brain performs complex functions. Different models that approximate the behavior of interconnected neurons have been proposed in the literature. Implementing these models to simulate neuron behavior at an appropriately detailed level to observe collective phenomena is computationally intensive. In this study we analyze the coupled Leaky Integrate-and-Fire model and report on the issues that affect performance when the model is implemented on a GPU. We conclude that the problem is heavily memory-bound. Advances in memory technology at the hardware level seem to be the deciding factor to achieve better performance on the GPU. Our results show that using an NVidia K40 GPU a modest 2x speedup can be achieved compared to a parallel implementation running on a modern multi-core CPU. However, a substantial speedup of 11.1x can be achieved using an NVidia V100 GPU, mainly due to the improvements in its memory subsystem.

How separable median filters can get better results than full 2D versions

The well-known method of median filtering is used both in a wide range of application frameworks and as a standalone filter. Small-window median filters can highly reduce the power of salt and pepper or additive Gaussian noise and minimize edge blurring. Large-window filters are also used, for example to estimate the background of images. Currently, the window size should not be an issue as a constant time algorithm and several implementations, including GPU-based codes, have been proposed in recent years. Unfortunately, none of these constant time implementations manage to fully exploit the capabilities of modern GPUs, and thus the throughputs of large-window median filters remain far below the peak throughputs allowed by recent GPU models. This paper aims at showing that a separable approximation of a 2D median filtering is often stronger than its full 2D implementation. Statistical and theoretical analysis are conducted to show and explain this fact which had so far remained unobserved. It is confirmed by experimentations on a dataset composed of 10,000 images, corrupted by different levels of salt and pepper noise. Separable and full 2D median filter algorithms are compared with several metrics, notably PSNR and MSSIM. In addition, a GPU implementation of 2D separable median filters is also proposed. This implementation is able to output up to 125 billion pixels per second on a recent Volta V100, which significantly outperforms existing implementations like Nvidia’s NPP library or Green’s code, resulting in the fastest median filtering solution to date.

An evaluation of fast segmented sorting implementations on GPUs

Many real-world problems require a sorting operation as part of their efficient solution. Some examples of this are real-time plasma diagnostic, image re-ranking, and suffix array construction. These problems usually involve a large amount of data, so their solutions need a particular application of the sorting procedure, consisting of sorting several arrays in matrix rows or array segments, an operation called segmented sorting. Previous studies showed that a merge sort-based strategy and a strategy called fix sort executed this operation on GPUs with good performance for different array sizes. In this work, we compare the fastest segmented sorting GPU implementations on seven different GPU models with various input data scenarios, including scenarios with varying numbers of segments, segment sizes, and considering segments of the same and different sizes. We first performed algorithm analysis to explain how the number of segments affects each implementation’s performance. Then, we perform an S-curve analysis and observe that, even though each strategy might be the fastest option for a subset of the sorting scenarios, some approaches may cause very high slowdowns on specific scenarios. We also compare the strategies using heat maps, show that their performance depends on the array size and number of segments, and propose a recommendation map to support selecting the best overall implementation based on the size and number of segments. Our experimental results show that choosing a strategy based on our recommendation map leads to the best strategy on 47.57% of the cases and a maximum slowdown of less than 1.5 times in 93.58% of the cases. Moreover, on average, the recommended strategy is only 1.11× worse than the optimum one.Finally, we evaluated how each strategy behaves when sorting arrays with different and equal segment sizes and showed that the fix sort-based approaches take roughly the same time to sort arrays with equal or different segment sizes, while the approach proposed by Hou et al. usually takes longer to sort arrays with different segment sizes than arrays with equal segment sizes.

Universal Background Subtraction Based on Arithmetic Distribution Neural Network.

We propose a universal background subtraction framework based on the Arithmetic Distribution Neural Network (ADNN) for learning the distributions of temporal pixels. In our ADNN model, the arithmetic distribution operations are utilized to introduce the arithmetic distribution layers, including the product distribution layer and the sum distribution layer. Furthermore, in order to improve the accuracy of the proposed approach, an improved Bayesian refinement model based on neighboring information, with a GPU implementation, is incorporated. In the forward pass and backpropagation of the proposed arithmetic distribution layers, histograms are considered as probability density functions rather than matrices. Thus, the proposed approach is able to utilize the probability information of the histogram and achieve promising results with a very simple architecture compared to traditional convolutional neural networks. Evaluations using standard benchmarks demonstrate the superiority of the proposed approach compared to state-of-the-art traditional and deep learning methods. To the best of our knowledge, this is the first method to propose network layers based on arithmetic distribution operations for learning distributions during background subtraction.

NL-CALIC Soft Decoding Using Strict Constrained Wide-Activated Recurrent Residual Network.

In this work, we propose a normalized Tanh activate strategy and a lightweight wide-activate recurrent structure to solve three key challenges of the soft-decoding of near-lossless codes: 1. How to add an effective strict constrained peak absolute error (PAE) boundary to the network; 2. An end-to-end solution that is suitable for different quantization steps (compression ratios). 3. Simple structure that favors the GPU and FPGA implementation. To this end, we propose a Wide-activated Recurrent structure with a normalized Tanh activate strategy for Soft-Decoding (WRSD). Experiments demonstrate the effectiveness of the proposed WRSD technique that WRSD outperforms better than the state-of-the-art soft decoders with less than 5% number of parameters, and every computation node of WRSD requires less than 64KB storage for the parameters which can be easily cached by most of the current consumer-level GPUs. Source code is available at https://github.com/dota-109/WRSD.

The Performance Analysis of the Thermal Discrete Element Method Computations on the GPU

The paper presents a GPU implementation of the thermal discrete element method (TDEM) and the comparative analysis of its performance. Several discrete element models for granular flows, the bonded particle model and the TDEM are considered for quantitative comparison of computational performance. The performance measured on NVIDIA(R) Tesla™ P100 GPU is compared with that attained by running the same OpenCL code on Intel(R) Xeon™ E5-2630 CPU with 20 cores. The presented GPU implementation of the TDEM increases the computing time of the bonded particle model only up to 30.6 % of the computing time of the simplest DEM model, which is an acceptable decrease in the performance required for solving coupled thermomechanical problems.

Collective proposal distributions for nonlinear MCMC samplers: Mean-field theory and fast implementation

Over the last decades, various “non-linear” MCMC methods have arisen. While appealing for their convergence speed and efficiency, their practical implementation and theoretical study remain challenging. In this paper, we introduce a non-linear generalization of the Metropolis-Hastings algorithm to a proposal that depends not only on the current state, but also on its law. We propose to simulate this dynamics as the mean field limit of a system of interacting particles, that can in turn itself be understood as a generalisation of the Metropolis-Hastings algorithm to a population of particles. Under the double limit in number of iterations and number of particles we prove that this algorithm converges. Then, we propose an efficient GPU implementation and illustrate its performance on various examples. The method is particularly stable on multimodal examples and converges faster than the classical methods.

Comparative Study of CUDA GPU Implementations in Python With the Fast Iterative Shrinkage-Thresholding Algorithm for LASSO

Comparative Study of CUDA GPU Implementations in Python With the Fast Iterative Shrinkage-Thresholding Algorithm for LASSO

Dose kernel decomposition for spot-based radiotherapy treatment planning.

Pre-calculation of accurate dose deposition kernels for treatment planning of spot-based radiotherapies, such as Gamma Knife (GK) and Gamma Pod (GP), can be very time-consuming and may require large data storage with an enormous number of possible spots. We proposed a novel kernel decomposition (KD) model to address accurate and fast (real-time) dose calculation with reduced data storage requirements for spot-based treatment planning. The application of the KD model was demonstrated for clinical GK and GP radiotherapy platforms. The dose deposition kernel at each spot (shot position) is modeled as the product of a shift-invariant kernel based on a reference kernel and spatially variant scale factor. The reference kernel, one for each collimator, is defined at the center of the commissioning phantom for GK and at the center of the treatment target for GP and calculated using the Monte Carlo (MC) method. The spatially variant scale factor is defined as the ratio of the mean tissue maximum ratio (TMR) at the candidate shot position to that at the reference kernel position, and the mean TMR map is calculated within the entire volume through parallel beam ray tracing on the density image followed by averaging over all source directions. The proposed KD dose calculations were compared with the MC method and with the GK and GP treatment planning system (TPS) computations for various shot positions and collimator sizes utilizing a phantom and 14 and 12 clinical plans for GK and GP, respectively. For the phantom study, the KD Gamma index (3%/1mm) passing rates were greater than 99% (median 100%) relative to the MC doses, except for the shots close to the boundary. The passing rates dropped below 90% for 8mm (16mm) shots positioned within ∼1cm (∼2cm) of the boundary. For the clinical GK plans, the KD Gamma passing rates were greater than 99% (median 100%) compared to the MC and greater than 92% (median 99%) compared to the TPS. For the clinical GP plans, the KD Gamma passing rates were greater than 95% (median 98%) compared to the MC and greater than 91% (median 97%) compared to the TPS. The scale factors were calculated in sub-seconds with GPU implementation and only need to be calculated once before treatment plan optimization. The calculation of the dose kernel was also within sub-seconds without requiring beam-by-beam calculation commonly done in the TPS. The proposed model can provide an accurate dose and enables real-time dose and derivative calculations by kernel shifting and scaling without pre-calculating or requiring large data storage for GK and GP dose deposition kernels during treatment planning. This model could be useful for spot-based radiotherapy treatment planning by allowing an efficient global fine search for optimal spots.