Graphics Processing Units Algorithm Research Articles

Many-core algorithms for high-dimensional gradients on phylogenetic trees.

Advancements in high-throughput genomic sequencing are delivering genomic pathogen data at an unprecedented rate, positioning statistical phylogenetics as a critical tool to monitor infectious diseases globally. This rapid growth spurs the need for efficient inference techniques, such as Hamiltonian Monte Carlo (HMC) in a Bayesian framework, to estimate parameters of these phylogenetic models where the dimensions of the parameters increase with the number of sequences N. HMC requires repeated calculation of the gradient of the data log-likelihood with respect to (wrt) all branch-length-specific (BLS) parameters that traditionally takes O(N2) operations using the standard pruning algorithm. A recent study proposes an approach to calculate this gradient in O(N), enabling researchers to take advantage of gradient-based samplers such as HMC. The CPU implementation of this approach makes the calculation of the gradient computationally tractable for nucleotide-based models but falls short in performance for larger state-space size models, such as Markov-modulated and codon models. Here, we describe novel massively parallel algorithms to calculate the gradient of the log-likelihood wrt all BLS parameters that take advantage of graphics processing units (GPUs) and result in many fold higher speedups over previous CPU implementations. We benchmark these GPU algorithms on three computing systems using three evolutionary inference examples exploring complete genomes from 997 dengue viruses, 62 carnivore mitochondria and 49 yeasts, and observe a >128-fold speedup over the CPU implementation for codon-based models and >8-fold speedup for nucleotide-based models. As a practical demonstration, we also estimate the timing of the first introduction of West Nile virus into the continental Unites States under a codon model with a relaxed molecular clock from 104 full viral genomes, an inference task previously intractable. We provide an implementation of our GPU algorithms in BEAGLE v4.0.0 (https://github.com/beagle-dev/beagle-lib), an open-source library for statistical phylogenetics that enables parallel calculations on multi-core CPUs and GPUs. We employ a BEAGLE-implementation using the Bayesian phylogenetics framework BEAST (https://github.com/beast-dev/beast-mcmc).

GPU accelerated volumetric lattice Boltzmann model for image-based hemodynamics in portal hypertension

We are developing a new technique for monitoring portal hypertension by the pressure gradient between the portal vein and the inferior vena cava (PPG) based on non-invasive measurements (MRI images). Massive parametrization and classification are required to investigate the underlying relationship between the porosity and the stages of liver cirrhosis numerically, and the hepatic-portal venous system is a multi-scale system. Both of them need high computational costs. The suitability of the lattice Boltzmann method for GPU (Graphics Processing Unit) parallel computation provides an opportunity to overcome it. In this paper, we perform GPU parallelization and optimization for the volumetric lattice Boltzmann method with arbitrary geometry based on images. Three application cases, including pipe flow, hemodynamics in the portal venous system, and hemodynamics in a simple hepatic-portal venous system, are employed to prove the method can be applied in the hepatic-portal venous system based on accuracy and efficiency. The reliability of the model is qualitatively validated by the analytical solution of velocity and pressure difference distribution of pipe flow and quantitatively confirmed by the pulsatility of velocity and pressure difference that can be neglected in the portal venous system. The performance of the application cases is examined with Intel Broadwell E5-2683 v3@ 2.30 GHz (CPU) and NVIDIA Tesla V100 16GB (GPU). It shows the GPU algorithm for sparse geometry (SPARSE) has a similar speed to the regular GPU algorithm for dense geometry (DENSE) when the fluid volume fraction (q) is close to 1. And SPARSE speeds up to 2.2 times compared with DENSE when q is in the range of 0.19∼0.27. Meanwhile, the saving ratio of memory cost depends on q. For a numerical case in the hepatic-portal venous system, i.e., a large-scale system, parallel execution can be converged around half an hour with SPARSE, while the memory spills the limitation with DENSE with a single GPU. Hence, multi-GPU implementation is applied to release the limitation, and it can improve performance by increasing the number of GPU cards. In summary, the method presented in the paper is feasibility applied in the hepatic-portal venous system, laying the foundation of the new technique development.

Large‐scale 3D fast Fourier transform computation on a GPU

AbstractWe propose a novel graphics processing unit (GPU) algorithm that can handle a large‐scale 3D fast Fourier transform (i.e., 3D‐FFT) problem whose data size is larger than the GPU's memory. A 1D FFT‐based 3D‐FFT computational approach is used to solve the limited device memory issue. Moreover, to reduce the communication overhead between the CPU and GPU, we propose a 3D data‐transposition method that converts the target 1D vector into a contiguous memory layout and improves data transfer efficiency. The transposed data are communicated between the host and device memories efficiently through the pinned buffer and multiple streams. We apply our method to various large‐scale benchmarks and compare its performance with the state‐of‐the‐art multicore CPU FFT library (i.e., fastest Fourier transform in the West [FFTW]) and a prior GPU‐based 3D‐FFT algorithm. Our method achieves a higher performance (up to 2.89 times) than FFTW; it yields more performance gaps as the data size increases. The performance of the prior GPU algorithm decreases considerably in massive‐scale problems, whereas our method's performance is stable.

Fast computation of fractal dimension for 2D, 3D and 4D data

The box-counting (BC) algorithm is one of the most popular methods for calculating the fractal dimension (FD) of binary data. FD analysis has many important applications in the biomedical field, such as cancer detection from 2D computed axial tomography images, Alzheimer’s disease diagnosis from magnetic resonance 3D volumetric data, and consciousness states characterization based on 4D data extracted from electroencephalography (EEG) signals, among many others. Currently, these kinds of applications use data whose size and amount can be very large, with high computation times needed to calculate the BC of the whole datasets. In this study we present a very efficient parallel implementation of the BC algorithm for its execution on Graphics Processing Units (GPU). Our algorithm can process 2D, 3D and 4D data and we tested it on two platforms with different hardware configurations. The results showed speedups of up to 92.38 × (2D), 57.27 × (3D) and 75.73 × (4D) with respect to the corresponding CPU single-thread implementations of the same algorithm. Against an OpenMP multi-thread CPU implementation, our GPU algorithm achieved speedups of up to 16.12 × (2D), 6.86 × (3D) and 7.49 × (4D). We have also compared our algorithm to a previous GPU implementation of the BC algorithm in 3D, achieving a speedup of up to 4.79 × . Finally, as a practical application of our GPU BC algorithm a study comparing the FD of 4D data extracted from the EEGs of a schizophrenia patient and a healthy subject was performed. The computation time for processing 40 4D matrices was reduced from three hours (sequential CPU) to less than three minutes with our GPU algorithm.

Focused wave interaction with a partially-immersed rectangular box using 2-D incompressible SPH on a GPU comparing with experiment and linear theory

Wave interaction with bodies spanning the water surface is a common problem. A basic rectangular box of small draft is considered in two dimensions and focused waves are used to minimise reflections. Overtopping occurs in some cases making the inherently mass conservative smoothed particle hydrodynamics (SPH) method appropriate. This is applied in incompressible (divergence-free) form with a 2-D graphics processing unit (GPU) version of the Chow et al. (2019) code with some improvements. A newly parallelised GPU algorithm for the population of a compressed sparse row matrix storage array for the pressure Poisson equation is presented and provides a 327–446 times function speedup and overall simulation time speedups of 1.45–1.78 times with over 5 million particles. The inclusion of the parallel function leads to a completely parallel ISPH algorithm. Peak pressures on the base and front face are compared with experiment and linear (potential-flow) theory. The experiments used periodic focused waves which showed some variation in form about a peak crest although crest elevations were repeatable, and these are reproduced in the model. Converged incompressible SPH (ISPH) values are in approximate agreement with both. Overtopping of the box shows qualitative agreement with experiment. While linear theory cannot account for overtopping or viscous (eddy-shedding) effects, submerged pressure prediction provides a useful approximation. The ISPH model is single phase and limited video evidence suggests that air entrapment can occur in the initial stages of overtopping near the front face. Quite complex vorticity generation and eddy shedding is predicted with free-surface interaction. The study shows the capability and limitations of single phase ISPH for complex wave-body interaction.

An investigation of hybrid CPU-GPU solvers for supersonic reacting flow simulation with detailed chemical kinetics

Accurate simulations of supersonic reacting flow require the use of detailed chemical kinetics to capture the complex phenomena such as auto-ignition and broken, thickened flamelets. However, the detailed chemical mechanisms for hydrocarbon fuels used in supersonic combustion typically have a great number of species and reactions. This paper reports on the implementation of a hybrid CPU-GPU solver by using the graphics processing unit (GPU) as a coprocessor to accelerate the computation of chemical source terms in modeling supersonic reacting flows with detailed chemical kinetics. The chemical kinetics are solved explicitly via a fully-coupled method with an implicit scheme. Three parallel strategies with different storage schemes were explored in implementing the highly parallel GPU algorithms for accelerating combustion modeling. Significant acceleration was obtained by using these parallel algorithms for chemical source terms evaluation. In order to verify the viability of the solver, two supersonic combustion flow test cases were conducted. The performance of these algorithms is evaluated by modeling four chemical reaction mechanisms on different grids compared with the CPU-only solver. The computation cost of these algorithms is basically proportional to the size of the grids. The performance gain of these parallel algorithms depends on the size of the reaction mechanism. A maximum speedup factor of over 80 was obtained by modeling the mechanism consisting of 348 species and 2163 reactions. The significant performance improvement provided by these parallel algorithms can provide a significant perspective for designing GPU-accelerated algorithms for applications in simulating supersonic reacting flows.

The transformational role of GPU computing and deep learning in drug discovery

Deep learning has disrupted nearly every field of research, including those of direct importance to drug discovery, such as medicinal chemistry and pharmacology. This revolution has largely been attributed to the unprecedented advances in highly parallelizable graphics processing units (GPUs) and the development of GPU-enabled algorithms. In this Review, we present a comprehensive overview of historical trends and recent advances in GPU algorithms and discuss their immediate impact on the discovery of new drugs and drug targets. We also cover the state-of-the-art of deep learning architectures that have found practical applications in both early drug discovery and consequent hit-to-lead optimization stages, including the acceleration of molecular docking, the evaluation of off-target effects and the prediction of pharmacological properties. We conclude by discussing the impacts of GPU acceleration and deep learning models on the global democratization of the field of drug discovery that may lead to efficient exploration of the ever-expanding chemical universe to accelerate the discovery of novel medicines.

Perfect Hashing for State Space Exploration on the GPU

This paper exploits parallel computing power of graphics cards to accelerate state space search. We illustrate that modern graphics processing units (GPUs) have the potential to speed up breadth-first search significantly. For a bitvector representation of the search frontier, GPU algorithms with one and two bits per state are presented. Efficient perfect hash functions and their inverse are explored in order to achieve enhanced compression. We report maximal speed-ups of up to a factor of 27 wrt. single core CPU computation.

Exploring Means to Enhance the Efficiency of GPU Bitmap Index Query Processing

Once exotic, computational accelerators are now commonly available in many computing systems. Graphics processing units (GPUs) are perhaps the most frequently encountered computational accelerators. Recent work has shown that GPUs are beneficial when analyzing massive data sets. Specifically related to this study, it has been demonstrated that GPUs can significantly reduce the query processing time of database bitmap index queries. Bitmap indices are typically used for large, read-only data sets and are often compressed using some form of hybrid run-length compression. In this paper, we present three GPU algorithm enhancement strategies for executing queries of bitmap indices compressed using word aligned hybrid compression: (1) data structure reuse (2) metadata creation with various type alignment and (3) a preallocated memory pool. The data structure reuse greatly reduces the number of costly memory system calls. The use of metadata exploits the immutable nature of bitmaps to pre-calculate and store necessary intermediate processing results. This metadata reduces the number of required query-time processing steps. Preallocating a memory pool can reduce or entirely remove the overhead of memory operations during query processing. Our empirical study showed that performing a combination of these strategies can achieve 32.4times to 98.7times speedup over the current state-of-the-art implementation. Our study also showed that by using our enhancements, a common gaming GPU can achieve a 15.0times speedup over a more expensive high-end CPU.

GPU algorithms for density matrix methods on MOPAC: linear scaling electronic structure calculations for large molecular systems.

Purification of the density matrix methods should be employed when dealing with complex chemical systems containing many atoms. The running times for these methods scale linearly with the number of atoms if we consider the sparsity from the density matrix. Since the efficiency expected from those methods is closely tied to the underlying parallel implementations of the linear algebra operations (e.g., P2 = P × P), we proposed a central processing unit (CPU) and graphics processing unit (GPU) parallel matrix-matrix multiplication in SVBR (symmetrical variable block row) format for energy calculations through the SP2 algorithm. This algorithm was inserted in MOPAC's MOZYME method, using the original LMO Fock matrix assembly, and the atomic integral calculation implemented on it. Correctness and performance tests show that the implemented SP2 is accurate and fast, as the GPU is able to achieve speedups up to 40 times for a water cluster system with 42,312 orbitals running in one NVIDIA K40 GPU card compared to the single-threaded version. The GPU-accelerated SP2 algorithm using the MOZYME LMO framework enables the calculations of semiempirical wavefunction with stricter SCF criteria for localized charged molecular systems, as well as the single-point energies of molecules with more than 100.000 LMO orbitals in less than 1h. Graphical abstract Parallel CPU and GPU purification algorithms for electronic structure calculations were implemented in MOPAC's MOZYME method. Some matrices in these calculations, e.g., electron density P, are compressed, and the developed linear algebra operations deal with non-zero entries only. We employed the NVIDIA/CUDA platform to develop GPU algorithms, and accelerations up to 40 times for larger systems were achieved.

Parallel acceleration of CPU and GPU range queries over large data sets

Data management systems commonly use bitmap indices to increase the efficiency of querying scientific data. Bitmaps are usually highly compressible and can be queried directly using fast hardware-supported bitwise logical operations. The processing of bitmap queries is inherently parallel in structure, which suggests they could benefit from concurrent computer systems. In particular, bitmap-range queries offer a highly parallel computational problem, and the hardware features of graphics processing units (GPUs) offer an alluring platform for accelerating their execution.In this paper, we present four GPU algorithms and two CPU-based algorithms for the parallel execution of bitmap-range queries. We show that in 98.8% of our tests, using real and synthetic data, the GPU algorithms greatly outperform the parallel CPU algorithms. For these tests, the GPU algorithms provide up to 54.1 × speedup and an average speedup of 11.5× over the parallel CPU algorithms. In addition to enhancing performance, augmenting traditional bitmap query systems with GPUs to offload bitmap query processing allows the CPU to process other requests.

Optimized Block-Based Algorithms to Label Connected Components on GPUs

Connected Components Labeling (CCL) is a crucial step of several image processing and computer vision pipelines. Many efficient sequential strategies exist, among which one of the most effective is the use of a block-based mask to drastically cut the number of memory accesses. In the last decade, aided by the fast development of Graphics Processing Units (GPUs), a lot of data parallel CCL algorithms have been proposed along with sequential ones. Applications that entirely run in GPU can benefit from parallel implementations of CCL that allow to avoid expensive memory transfers between host and device. In this paper, two new eight-connectivity CCL algorithms are proposed, namely Block-based Union Find (BUF) and Block-based Komura Equivalence (BKE). These algorithms optimize existing GPU solutions introducing a block-based approach. Extensions for three-dimensional datasets are also discussed. In order to produce a fair comparison with previously proposed alternatives, YACCLAB, a public CCL benchmarking framework, has been extended and made suitable for evaluating also GPU algorithms. Moreover, three-dimensional datasets have been added to its collection. Experimental results on real cases and synthetically generated datasets demonstrate the superiority of the new proposals with respect to state-of-the-art, both on 2D and 3D scenarios.

Heterogeneous CPU + GPU Algorithm for Variational Two-Electron Reduced-Density Matrix-Driven Complete Active-Space Self-Consistent Field Theory.

We present a heterogeneous central processing unit (CPU) + graphical processing unit (GPU) algorithm for the direct variational optimization of the two-electron reduced-density matrix (2RDM) under two-particle N-representability conditions. This variational 2RDM (v2RDM) approach is the driver for a polynomially scaling approximation to configuration-interaction-driven complete active-space self-consistent field (CASSCF) theory. For v2RDM-based CASSCF computations involving an active space consisting of 50 electrons in 50 orbitals, we observe a speedup of a factor of 3.7 when the code is executed on a combination of an NVIDIA TITAN V GPU and an Intel Core i7-6850k CPU, relative to the case when the code is executed on the CPU alone. We use this GPU-accelerated v2RDM-CASSCF algorithm to explore the electronic structure of the 3,k-circumacene and 3,k-periacene series (k = 2-7) and compare indicators of polyradical character in the lowest-energy singlet states to those observed for oligoacene molecules. The singlet states in larger circumacene and periacene molecules display the same polyradical characteristics observed in oligoacenes, with the onset of this behavior occurring at smallest k for periacenes, followed by the circumacenes and then the oligoacenes. However, the unpaired electron density that accumulates along the zigzag edge of the circumacenes is slightly less than that which accumulates in the oligoacenes, while periacenes clearly exhibit the greatest buildup of unpaired electron density in this region.

A GPU Algorithm for Agent-Based Models to Simulate the Integration of Cell Membrane Signals

Simulation of complex biological systems with agent-based models is becoming more relevant with the increase in Graphics Processing Unit (GPU) power. In those simulations, up to millions of virtual cells are individually computed, involving daunting processing times. An important part of computational models is the algorithm that manages how agents perceive their surroundings. This can be particularly problematic in three-dimensional environments where agents have deformable virtual membranes. This article presents a GPU algorithm that gives the possibility for agents to integrate the signals scattered on their virtual membrane. It is detailed to be coded in languages like OpenCL or Cuda. Its performances are tested to show its speed with current GPU devices. Finally, it was implemented inside an existing software to test and illustrate the possibilities it offers.

GPU acceleration of an iterative scheme for gas-kinetic model equations with memory reduction techniques

This paper presents a Graphics Processing Unit (GPU) acceleration of an iteration-based discrete velocity method (DVM) for gas-kinetic model equations. Unlike the previous GPU parallelization of explicit kinetic schemes, this work is based on a fast converging iterative scheme. The memory reduction techniques previously proposed for DVM are applied for GPU computing, enabling full three-dimensional (3D) solutions of kinetic model equations in the contemporary GPUs usually with a limited memory capacity that otherwise would need terabytes of memory. The GPU algorithm is validated against the direct simulation Monte Carlo (DSMC) simulation of the 3D lid-driven cavity flow and the supersonic rarefied gas flow past a cube with the phase-space grid points up to 0.7 trillion. The computing performance profiling on three models of GPUs shows that the two main kernel functions can utilize 56%∼79% of the GPU computing and memory resources. The performance of the GPU algorithm is compared with a typical parallel CPU implementation of the same algorithm using the Message Passing Interface (MPI). The comparison shows that the GPU program on K40 and K80 achieves 1.2∼2.8 and 1.2∼2.4 speedups for the 3D lid-driven cavity flow, respectively, compared with the MPI parallelized CPU program running on 96 CPU cores.

GPU parameter tuning for tall and skinny dense linear least squares problems

ABSTRACTLinear least squares problems (LLSPs) routinely arise in many scientific and engineering problems. One of the fastest ways to solve LLSPs involves performing calculations in parallel on graphics processing units (GPUs). However, GPU algorithms are typically designed for one GPU architecture and may be suboptimal or unusable on another GPU. To design optimal algorithms for any GPU with little need for modifying code, tuneable parameters can simplify the transition of GPU algorithms to different GPU architectures. In this paper, we investigate the benefits of using derivative-free optimization (DFO) and simulation optimization (SO) to systematically optimize tuneable parameters for a GPU or hybrid CPU/GPU LLSP solvers. Computational experiments show that both DFO and SO can be effective tools for determining optimal tuning parameters that can speed up the performance of the popular LLSP solver MAGMA by about 1.8x, compared to MAGMA's default parameters for large tall and skinny matrices. Using DFO solvers, we were able to identify optimal parameters after enumerating an order of magnitude fewer parameter combinations than with direct enumeration. Additionally, the proposed approach is faster than a state-of-the-art autotuner and provides better tuning parameters.

GPU-DAEMON: GPU algorithm design, data management & optimization template for array based big omics data

In the age of ever increasing data, faster and more efficient data processing algorithms are needed. Graphics Processing Units (GPU) are emerging as a cost-effective alternative architecture for high-end computing. The optimal design of GPU algorithms is a challenging task which requires thorough understanding of the high performance computing architecture as well as the algorithmic design. The steep learning curve needed for effective GPU-centric algorithm design and implementation requires considerable expertise, time, and resources. In this paper, we present GPU-DAEMON, a GPU Data Management, Algorithm Design and Optimization technique suitable for processing array based big omics data. Our proposed GPU algorithm design template outlines and provides generic methods to tackle critical bottlenecks which can be followed to implement high performance, scalable GPU algorithms for given big data problem. We study the capability of GPU-DAEMON by reviewing the implementation of GPU-DAEMON based algorithms for three different big data problems. Speed up of as large as 386x (over the sequential version) and 50x (over naive GPU design methods) are observed using the proposed GPU-DAEMON. GPU-DAEMON template is available at https://github.com/pcdslab/GPU-DAEMON and the source codes for GPU-ArraySort, G-MSR and GPU-PCC are available at https://github.com/pcdslab.

Accelerating GPU betweenness centrality

Graphs that model social networks, numerical simulations, and the structure of the Internet are enormous and cannot be manually inspected. A popular metric used to analyze these networks is Betweenness Centrality (BC), which has applications in community detection, power grid contingency analysis, and the study of the human brain. However, these analyses come with a high computational cost that prevents the examination of large graphs of interest. Recently, the use of Graphics Processing Units (GPUs) has been promising for efficient processing of unstructured data sets. Prior GPU implementations of BC suffer from large local data structures and inefficient graph traversals that limit scalability and performance. Here we present a hybrid GPU implementation that provides good performance on graphs of arbitrary structure rather than just scale-free graphs as was done previously. Our methods achieve up to 13× speedup on high-diameter graphs and an average of 2.71× speedup overall compared to the best existing GPU algorithm. We also observe near linear speedup when running BC on 192 GPUs.

An efficient GPU algorithm for tetrahedron-based Brillouin-zone integration

We report an efficient algorithm for calculating momentum-space integrals in solid state systems on modern graphics processing units (GPUs). Our algorithm is based on the tetrahedron method, which we demonstrate to be ideally suited for execution in a GPU framework. In order to achieve maximum performance, all floating point operations are executed in single precision. For benchmarking our implementation within the CUDA programming framework we calculate the orbital-resolved density of states in an iron-based superconductor. However, our algorithm is general enough for the achieved improvements to carry over to the calculation of other momentum integrals such as, e.g. susceptibilities. If our program code is integrated into an existing program for the central processing unit (CPU), i.e. when data transfer overheads exist, speedups of up to a factor ∼130 compared to a pure CPU implementation can be achieved, largely depending on the problem size. In case our program code is integrated into an existing GPU program, speedups over a CPU implementation of up to a factor ∼165 are possible, even for moderately sized workloads.

Simulating thermal effects of MR-guided focused ultrasound in cortical bone and its surrounding tissue.

Magnetic resonance-guided focused ultrasound (MRgFUS) is emerging as a treatment alternative for osteoid osteoma and painful bone metastases. This study describes a new simulation platform that predicts the distribution of heat generated by MRgFUS when applied to bone tissue. Calculation of the temperature distribution was performed using two mathematical models. The first determined the propagation and absorption of acoustic energy through each medium, and this was performed using a multilayered approximation of the Rayleigh integral method. The ultrasound energy distribution derived from these equations could then be converted to heat energy, and the second mathematical model would then use the heat generated to determine the final temperature distribution using a finite-difference time-domain application of Pennes' bio-heat transfer equation. Anatomical surface geometry was generated using a modified version of a mesh-based semiautomatic segmentation algorithm, and both the acoustic and thermodynamic models were calculated using a parallelized algorithm running on a graphics processing unit (GPU) to greatly accelerate computation time. A series of seven porcine experiments were performed to validate the model, comparing simulated temperatures to MR thermometry and assessing spatial, temporal, and maximum temperature accuracy in the soft tissue. The parallelized algorithm performed acoustic and thermodynamic calculations on grids of over 108 voxels in under 30s for a simulated 20s of heating and 40s of cooling, with a maximum time per calculated voxel of less than 0.3μs. Accuracy was assessed by comparing the soft tissue thermometry to the simulation in the soft tissue adjacent to bone using four metrics. The maximum temperature difference between the simulation and thermometry in a region of interest around the bone was measured to be 5.43±3.51°C average absolute difference and a percentage difference of 16.7%. The difference in heating location resulted in a total root-mean-square error of 4.21±1.43mm. The total size of the ablated tissue calculated from the thermal dose approximation in the simulation was, on average, 67.6% smaller than measured from the thermometry. The cooldown was much faster in the simulation, where it decreased by 14.22±4.10°C more than the thermometry in 40s after sonication ended. The use of a Rayleigh-based acoustic model combined with a discretized bio-heat transfer model provided a rapid three-dimensional calculation of the temperature distribution through bone and soft tissue during MRgFUS application, and the parallelized GPU algorithm provided the computational speed that would be necessary for an intraoperative treatment planning software platform.