GPU Algorithm Research Articles

Ultra-fast computation of fractal dimension for RGB images

The fractal dimension (FD) is a quantitative parameter widely used to analyze digital images in many application fields such as image segmentation, feature extraction, object recognition, texture analysis, and image compression and denoising, among many others. A variety of algorithms have been previously proposed for estimating the FD, however most of them are limited to binary or gray-scale images only. In recent years, several authors have proposed algorithms for computing the FD of color images. Nevertheless, almost all these methods are computationally inefficient when analyzing large images. Nowadays, color images can be very large in size, and there is a growing trend toward even larger datasets. This implies that the time required to calculate the FD of such datasets can become extremely long. In this paper we present a very efficient GPU algorithm, implemented in CUDA, for computing the FD of RGB color images. Our solution is an extension to RGB of the differential box-counting (DBC) algorithm for gray-scale images. Our implementation simplifies the box-counting computation to very simple operations which are easily combined across iterations. We evaluated our algorithm on two distinct hardware/software platforms using a set of images of increasing size. The performance of our method was compared against two recent FD algorithms for RGB images: a fast box-merging GPU algorithm, and the most advanced approach based on extending the DBC method. The results showed that our GPU algorithm performed very well and achieved speedups of up to 7.9× and 6172.6× regarding these algorithms, respectively. In addition, our algorithm achieved average error rates similar to those obtained by the two reference algorithms when estimating the FD for synthetic images with known FD values, and even outperformed them when processing large images. These results suggest that our GPU algorithm offers a highly reliable and ultra-fast solution for estimating the FD of color images.

Implementation and analysis of GPU algorithms for Vecchia Approximation

Implementation and analysis of GPU algorithms for Vecchia Approximation

A GPU Algorithm for Solving the Positions of New Pulsars

Timing newly discovered pulsars requires gradually building up a timing model that connects observations taken days to months apart. This sometimes can be challenging when our initial knowledge of the pulsar’s position is arcminutes off from its true position. Such a position error leads to significant arrival time shifts as a result of the Earth’s orbital motion. Traditional down-hill fitting timing algorithms become ineffective when our model predicts the wrong pulse rotations for our next observation. For some pulsars whose model prediction is not too far off, the correct rotation number could be found by trial-and-error methods. For the remaining challenging pulsars, a more generalized method is called for. This paper proposes a GPU-based algorithm that could exhaustively search a large area of trail positions for probable timing solutions. This could help find phase-connected timing solutions for new pulsars using brute force.

GPU-accelerated parallel all-pair shortest path routing within stochastic road networks

All-pair shortest path routing within stochastic road networks is often more complicated and computationally challenging than routing in deterministic networks because uncertainties in travel time are introduced. We develop and test a GPU-enabled approach to resolve this computational challenge. This approach is based on a simulation framework that includes four steps: simulation of road networks impedances, computation of candidate paths using an all-pair shortest path algorithm, travel time estimation across candidate paths, and all-pair optimal path calculations. We developed GPU algorithms to accelerate the second and third steps, which are computationally demanding for large road networks. We used the road network of Wuhan, China as a case study. Our experiments indicate that the GPU approach leads to thousands of times in improvement of acceleration, which makes all-pair shortest path routing within stochastic road networks computationally feasible.

Efficient Debris-flow Simulation for Steep Terrain Erosion

Erosion simulation is a common approach used for generating and authoring mountainous terrains. While water is considered the primary erosion factor, its simulation fails to capture steep slopes near the ridges. In these low-drainage areas, erosion is often approximated with slope-reducing erosion, which yields unrealistically uniform slopes. However, geomorphology observed that another process dominates the low-drainage areas: erosion by debris flow, which is a mixture of mud and rocks triggered by strong climatic events. We propose a new method to capture the interactions between debris flow and fluvial erosion thanks to a new mathematical formulation for debris flow erosion derived from geomorphology and a unified GPU algorithm for erosion and deposition. In particular, we observe that sediment and debris deposition tend to intersect river paths, which motivates the design of a new, approximate flow routing algorithm on the GPU to estimate the water path out of these newly formed depressions. We demonstrate that debris flow carves distinct patterns in the form of erosive scars on steep slopes and cones of deposited debris competing with fluvial erosion downstream.

Meta-Meshing and Triangulating Lattice Structures at a Large Scale

Lattice structures have been widely used in applications due to their superior mechanical properties. To fabricate such structures, a geometric processing step called triangulation is often employed to transform them into the STL format before sending them to 3D printers. Because lattice structures tend to have high geometric complexity, this step usually generates a large amount of triangles, a memory and compute-intensive task. This problem manifests itself clearly through large-scale lattice structures that have millions or billions of struts. To address this problem, this paper proposes to transform a lattice structure into an intermediate model called meta-mesh before undergoing real triangulation. Compared to triangular meshes, meta-meshes are very lightweight and much less compute-demanding. The meta-mesh can also work as a base mesh reusable for conveniently and efficiently triangulating lattice structures with arbitrary resolutions. A CPU+GPU asynchronous meta-meshing pipeline has been developed to efficiently generate meta-meshes from lattice structures. It shifts from the thread-centric GPU algorithm design paradigm commonly used in CAD to the recent warp-centric design paradigm to achieve high performance. This is achieved by a new data compression method, a GPU cache-aware data structure, and a workload-balanced scheduling method that can significantly reduce memory divergence and branch divergence. Experimenting with various billion-scale lattice structures, the proposed method is seen to be two orders of magnitude faster than previously achievable.

Lattice Boltzmann simulation of dissolution patterns in porous media: Single porosity versus dual porosity media

Understanding the influence of porous media structure, particularly dual porosity, on solvent transport and pore geometry evolution during chemical reactions is a complex and critical area of study. This research leverages the lattice Boltzmann method to investigate how the presence of aggregates in a medium affects solvent transport and pore space development, focusing on distinct dissolution regimes: face and wormhole dissolution. The study addresses the challenge of managing variable pore sizes in dual porosity media by developing specialized GPU algorithms, which efficiently handle fine grids and complex pore spaces. The findings reveal that dual porosity significantly enhances dissolution rates in both the face and wormhole dissolution regimes. Intriguingly, while the pattern of face dissolution remains largely unchanged, dual porosity markedly alters the pattern of wormhole dissolution. In dual-porosity media, the wormholes tend to be narrower and more elongated compared to the wider wormholes observed in single-porosity media. This variation is attributed to the reaction area dynamics, where the reduced reactive surface area along the main wormhole path in dual-porosity media results in less solvent engagement in the reaction processes. Moreover, the research provides insights into the microscale interactions in porous media, emphasizing how variations in microscale porosity can have substantial impacts on the overall dissolution dynamics. The study results are not only significant for understanding the fundamental aspects of chemical dissolution in porous media but also have practical implications in fields such as geo-energy and groundwater remediation. These findings help optimizing reaction processes in complex and heterogeneous porous systems, highlighting the need for detailed consideration of microstructural characteristics in modeling and industrial applications.

An efficiency-improved GPU algorithm for the 2 + 2 + 1 method in nonlinear beamforming

Abstract Nonlinear beamforming (NLBF) has emerged as a highly effective technology for enhancing seismic data quality. The crux of NLBF's success lies in its ability to robustly estimate local traveltime operators directly from input data, a process that entails solving millions or even billions of nonlinear optimization problems per input gather. Among the solvers utilized for estimating these operators is the 2 + 2 + 1 method, for which we have previously introduced algorithmic implementations on both the CPU and GPU platforms. In this paper, we present an efficiency-improved GPU algorithm for the 2 + 2 + 1 method, particularly beneficial when dealing with small data apertures in NLBF. Our enhanced GPU algorithm brings significant improvements in computation efficiency through several strategic measures, which include leveraging Horner's method to minimize the mathematical overhead of traveltime calculation, implementing a GPU-friendly data reduction algorithm to exploit GPU computational power, and optimizing shared GPU memory usage as the primary workspace whenever feasible. To demonstrate the tangible efficiency enhancement achieved by our new GPU algorithm, via two illustrative examples, we compare its performance with that of our previous implementation.

Stochastic seismic acoustic impedance inversion via a Markov-chain Monte Carlo method using a single GPU card

Seismic acoustic impedance inversion plays an important role in understanding subsurface structures and obtaining subsurface properties. The stochastic approach is one of the methods used for impedance inversion, and it aims to produce more reliable results by accounting for modeling uncertainty. Stochastic inversion represents the uncertainty of a subsurface model as a probability distribution and uses this distribution to estimate model parameters. In this study, seismic acoustic impedance inversion was performed using a Markov-chain Monte Carlo approach, which is a sampling method used in the stochastic process. Utilizing Bayesian inference based on prior information and observed data, we implemented acceptance probabilities and performed acoustic impedance inversion for field data through iterative calculations. We employed a single GPU card to execute the inversion algorithm and conducted a comparative analysis with the results obtained using a cluster composed of multiple CPU cores. Through a computational speed analysis, the efficiency of the algorithm using a GPU was verified, while uncertainty analysis was employed for algorithm validation. Through these analyses, we confirmed the feasibility of applying our developed GPU algorithm to tasks that require the inversion of extensive data, such as high-resolution seismic exploration and CO2 storage monitoring.

A parallel acceleration GPU algorithm for large deformation of thin shell structures based on peridynamics

A parallel acceleration GPU algorithm for large deformation of thin shell structures based on peridynamics

GPU Algorithms for Structured Sparse Matrix Multiplication with Diagonal Storage Schemes

Matrix–matrix multiplication is of singular importance in linear algebra operations with a multitude of applications in scientific and engineering computing. Data structures for storing matrix elements are designed to minimize overhead information as well as to optimize the operation count. In this study, we utilize the notion of the compact diagonal storage method (CDM), which builds upon the previously developed diagonal storage—an orientation-independent uniform scheme to store the nonzero elements of a range of matrices. This study exploits both these storage schemes and presents efficient GPU-accelerated parallel implementations of matrix multiplication when the input matrices are banded and/or structured sparse. We exploit the data layouts in the diagonal storage schemes to expose a substantial amount of fine-grained parallelism and effectively utilize the GPU shared memory to improve the locality of data access for numerical calculations. Results from an extensive set of numerical experiments with the aforementioned types of matrices demonstrate orders-of-magnitude speedups compared with the sequential performance.

GPU implementation of the discrete unified gas kinetic scheme for low-speed isothermal flows

In this paper, two GPU parallel algorithms are proposed for the discrete unified gas kinetic scheme (DUGKS) for simulating low-speed isothermal flows. Algorithm-I uses a two-level fine-grain technique for the parallelization of physical spatial space, while Algorithm-II adopts this technique for both physical spatial and particle velocity spaces. To evaluate the performance of the proposed algorithms, several typical benchmark problems are simulated, including the two-dimensional (2D) and three-dimensional (3D) lid-driven cavity flows, the micro channel and cavity flows. Numerical results show that our GPU algorithms can achieve satisfactory computational efficiency. For Algorithm-I, the speedup can reach 250 and 338 on a Tesla V100 GPU card for the 2D and 3D continuum cavity flows, respectively, and a hundredfold acceleration can be obtained for the rarefied cases. While for Algorithm-II, a speedup of about 70 can be attained for rarefied cases. However, it is not applied to continuum problems that only require a small number of velocity points. Moreover, comparisons between the two GPU algorithms are also conducted for the rarefied flows with various grid meshes and velocity directions. The results show that Algorithm-I performs better when physical mesh size is large, while Algorithm-II can provide higher efficiency for a coarser mesh with medium number of discrete velocities. Special attention is also paid to comparisons between Algorithm-I and MPI parallelization with 128 CPU cores based on physical space discretization approach, and it is found that Algorithm-I has a clear advantage on V100 GPU when dealing with sparse physical grids in both continuum and rarefied cases.

Second-order Stencil Descent for Interior-point Hyperelasticity

In this paper, we present a GPU algorithm for finite element hyperelastic simulation. We show that the interior-point method, known to be effective for robust collision resolution, can be coupled with non-Newton procedures and be massively sped up on the GPU. Newton's method has been widely chosen for the interior-point family, which fully solves a linear system at each step. After that, the active set associated with collision/contact constraints is updated. Mimicking this routine using a non-Newton optimization (like gradient descent or ADMM) unfortunately does not deliver expected accelerations. This is because the barrier functions employed in an interior-point method need to be updated at every iteration to strictly confine the search to the feasible region. The associated cost (e.g., per-iteration CCD) quickly overweights the benefit brought by the GPU, and a new parallelism modality is needed. Our algorithm is inspired by the domain decomposition method and designed to move interior-point-related computations to local domains as much as possible. We minimize the size of each domain (i.e., a stencil) by restricting it to a single element, so as to fully exploit the capacity of modern GPUs. The stencil-level results are integrated into a global update using a novel hybrid sweep scheme. Our algorithm is locally second-order offering better convergence. It enables simulation acceleration of up to two orders over its CPU counterpart. We demonstrate the scalability, robustness, efficiency, and quality of our algorithm in a variety of simulation scenarios with complex and detailed collision geometries.

Research on the Construction of Financial Computing Model Based on BSDE Algorithm

Economic and social development has made financial engineering an increasingly important research area, and more and more financial problems cannot be solved directly by analytical formulas. In view of this, algorithms that apply computer technology to financial engineering have emerged. In this study, the Backward Stochastic Differential Equation (BSDE) algorithm is used to investigate and analyse the problem of option pricing calculation in finance. In the research process, GBSDE-Theta parallel algorithm composed of BSDE-Theta algorithm and GPU algorithm uses the new algorithm to establish a computing model in the financial engineering field, which applies to the calculation of enterprise option pricing. The research results show that compared with the basic algorithm, the actual option values of the option pricing data obtained by using the GBSDE-Theta parallel algorithm are more closely matched. The computational model can achieve a speedup ratio of about 230 times of the serial version with the number of time steps [Formula: see text] and the number of simulated paths 80,000. About the relative error of the GBSDE-Theta algorithm, there are 80 points within 3% and only 16 points over 3.00%, which is a relatively small error. The above results show that the financial computing system obtained in this study is highly feasible and effective, and can provide a new research idea for the progress and development of other computations in the financial field.

End-to-end GPU acceleration of low-order-refined preconditioning for high-order finite element discretizations

In this article, we present algorithms and implementations for the end-to-end GPU acceleration of matrix-free low-order-refined preconditioning of high-order finite element problems. The methods described here allow for the construction of effective preconditioners for high-order problems with optimal memory usage and computational complexity. The preconditioners are based on the construction of a spectrally equivalent low-order discretization on a refined mesh, which is then amenable to, for example, algebraic multigrid preconditioning. The constants of equivalence are independent of mesh size and polynomial degree. For vector finite element problems in H(curl) and H(div) (e.g., for electromagnetic or radiation diffusion problems), a specially constructed interpolation–histopolation basis is used to ensure fast convergence. Detailed performance studies are carried out to analyze the efficiency of the GPU algorithms. The kernel throughput of each of the main algorithmic components is measured, and the strong and weak parallel scalability of the methods is demonstrated. The different relative weighting and significance of the algorithmic components on GPUs and CPUs is discussed. Results on problems involving adaptively refined nonconforming meshes are shown, and the use of the preconditioners on a large-scale magnetic diffusion problem using all spaces of the finite element de Rham complex is illustrated.

Multiparallel MMT: Faster ISD Algorithm Solving High-Dimensional Syndrome Decoding Problem

The hardness of the syndrome decoding problem (SDP) is the primary evidence for the security of code-based cryptosystems, which are one of the finalists in a project to standardize post-quantum cryptography conducted by the U.S. National Institute of Standards and Technology (NIST-PQC). Information set decoding (ISD) is a general term for algorithms that solve SDP efficiently. In this paper, we conducted a concrete analysis of the time complexity of the latest ISD algorithms under the limitation of memory using the syndrome decoding estimator proposed by Esser et al. As a result, we present that theoretically nonoptimal ISDs, such as May-Meurer-Thomae (MMT) and May-Ozerov, have lower time complexity than other ISDs in some actual SDP instances. Based on these facts, we further studied the possibility of multiple parallelization for these ISDs and proposed the first GPU algorithm for MMT, the multiparallel MMT algorithm. In the experiments, we show that the multiparallel MMT algorithm is faster than existing ISD algorithms. In addition, we report the first successful attempts to solve the 510-, 530-, 540- and 550-dimensional SDP instances in the Decoding Challenge contest using the multiparallel MMT.

Parallel border tracking in binary images for multicore computers

Border tracking in binary images is an important operation in many computer vision applications. The problem consists in finding borders in a 2D binary image (where all of the pixels are either 0 or 1). There are several algorithms available for this problem, but most of them are sequential. In a former paper, a parallel border tracking algorithm was proposed. This algorithm was designed to run in Graphics Processing units, and it was based on the sequential algorithm known as the Suzuki algorithm. In this paper, we adapt the previously proposed GPU algorithm so that it can be executed in multicore computers. The resulting algorithm is evaluated against its GPU counterpart. The results show that the performance of the GPU algorithm worsens (or even fails) for very large images or images with many borders. On the other hand, the proposed multicore algorithm can efficiently cope with large images.

Improving the Scalability of GPU Synchronization Primitives

General-purpose GPU applications increasingly use synchronization to enforce ordering between many threads accessing shared data. Accordingly, recently there has been a push to establish a common set of GPU synchronization primitives. However, the expressiveness of existing GPU synchronization primitives is limited. In particular the expensive GPU atomics often used to implement fine-grained synchronization make it challenging to implement efficient algorithms. Consequently, as GPU algorithms scale to millions or billions of threads, existing GPU synchronization primitives either scale poorly or suffer from livelock or deadlock issues because of heavy contention between threads accessing shared synchronization objects. We seek to overcome these inefficiencies by designing more efficient, scalable GPU barriers and semaphores. In particular, we show how multi-level sense reversing barriers and priority mechanisms for semaphores can be designed with the GPUs unique processing model in mind to improve performance and scalability of GPU synchronization primitives. Our results show that the proposed designs significantly improve performance compared to state-of-the-art solutions like CUDA Cooperative Groups and optimized CPU-style synchronization algorithms at medium and high contention levels, scale to an order of magnitude more threads, and avoid livelock in these situations unlike prior open source algorithms. Overall, across three modern GPUs the proposed barrier algorithm improves performance by an average of 33% over a GPU tree barrier algorithm and improves performance by an average of 34% over CUDA Cooperative Groups for five full-sized benchmarks at high contention levels; the new semaphore algorithm improves performance by an average of 83% compared to prior GPU semaphores.

Max-Tree Computation on GPUs

In Mathematical Morphology, the max-tree is a region-based representation that encodes the inclusion relationship of the threshold sets of an image. This tree has proved useful in numerous image processing applications. For the last decade, work has led to improving the construction time of this structure; mixing algorithmic optimizations, parallel and distributed computing. Nevertheless, there is still no algorithm that benefits from the computing power of the massively parallel architectures. In this work, we propose the first GPU algorithm to compute the max-tree. The proposed approach leads to significant speed-ups, and is up to one order of magnitude faster than the current State-of-the-Art parallel CPU algorithms. This work paves the way for a max-tree integration in image processing GPU pipelines and real-time image processing based on Mathematical Morphology. It is also a foundation for porting other image representations from Mathematical Morphology on GPUs.

DAMO: Deep Agile Mask Optimization for Full-Chip Scale

Continuous scaling of the very-large-scale integration system leaves a significant challenge on manufacturing; thus optical proximity correction (OPC) is widely applied in conventional design flow for manufacturability optimization. Traditional techniques conduct OPC by leveraging a lithography model but may suffer from prohibitive computational overhead. In addition, most of them focus on optimizing a single and local clip instead of addressing how to tackle the full-chip scale. In this article, we present DAMO, a high-performance and scalable deep-learning-enabled OPC system for full-chip scale. It is an end-to-end mask optimization paradigm that contains a deep lithography simulator (DLS) for lithography modeling and a deep mask generator (DMG) for mask pattern generation. Moreover, a novel layout splitting algorithm customized for DAMO is proposed, composed of DBSCAN clustering and KMeans++ clustering, to handle the full-chip OPC problem. Further, graph-based computation and parallelism techniques are proposed to deploy our GPU algorithms to accelerate computations. Extensive experiments show that DAMO outperforms state-of-the-art OPC solutions in both academia and industrial commercial toolkit.