GPU-oriented numerical algorithm to estimate formation factor of porous materials

Computer architectures are becoming more and more complex due to non-standardized memory accesses and hierarchical caches. It is very difficult for scientists and engineers to optimize their code to extract potential performance improvements on these architectures. Due to this, automatic performance tuning (AT) technology, hence, is a key technology to reduce cost of development for high performance numerical software.In this talk, the following two aims are folded. First, we introduce current AT studies. We focus on AT technology for numerical computations in viewpoint of numerical libraries, languages, code generators, and OS run-time software.Second, we explain ABCLibScript [1], which is an auto-tuning description language for C and Fortran90 for numerical computations to numerical software developers. ABCLibScript provides automatic code generation functions for dedicated code optimization, such as loop unrolling, algorithm selection, and varying of specified variables described by the user. We also explain HxABCLibScript[2], which is an AT language with extended function from original ABCLibScript to heterogeneous computer environment, which includes CPU and GPU (Graphics Processing Unit). The description of HxABCLibScript can free from selection of CPU and GPU switching to the arbitrary parts of program from users.The preliminary results show that the function of HxABCLibScript was highly efficient for simple kernels of typical numerical computations, such as a matrix-matrix multiplication, or a stencil computation from the Poisson's equation solver. The automatically generated codes from the description of HxABCLibScript can select the best computer resources between CPU and GPU according to problem size or the number of iterations on the program.

In this work we present a programming philosophy and a template code for achieving computational scalability when using multiple graphics processing units (GPUs) in the numerical solution of any mathematical system of equations found in continuum dynamic simulations. The programming philosophy exploits the principal characteristics of the GPU hardware, with emphasis in the delivering of threads with massive memory fetches, intense calculations using local registers and limited writes to global memory. The philosophy requires explicit formulas for calculations for which domain decomposition is trivial. The domains are decomposed in regions that use the local central processing unit (CPU) to communicate common interfaces using the message passing interface (MPI). A template code for the heat equation is established and tested for scalability. The novelty is that we show a series of codes, constructed from the basic template, that solve all the basic model equations found in continuum dynamics, and present illustrative results. The model equations are the heat equation, the Poisson equation, the shallow-water equations, the flow in porous media equations and the vorticity equations.

The Legendre spectral Galerkin method of self-adjoint second order elliptic equations usually results in a linear system with a dense and ill-conditioned coefficient matrix. In this paper, the linear system is solved by a preconditioned conjugate gradient (PCG) method where the preconditioner M is constructed by approximating the variable coefficients with a (T+1)-term Legendre series in each direction to a desired accuracy. A feature of the proposed PCG method is that the iteration step increases slightly with the size of the resulting matrix when reaching a certain approximation accuracy. The efficiency of the method lies in that the system with the preconditioner M is approximately solved by a one-step method based on the incomplete LU factorization technique with no fill-in, denoted by ILU(0). The ILU(0) factorization of \(M\in {\mathbb {R}}^{(N-1)^d\times (N-1)^d}\) can be computed using \({\mathcal {O}}(T^{2d} N^d)\) operations, and the number of nonzeros in the factorization factors is of \({\mathcal {O}}(T^{d} N^d)\), \(d=1,2,3\). A conclusion of the algorithm is to fast solve the resulting system from the Legendre Galerkin spectral method for Poisson equations with Dirichlet boundary conditions, which has a complexity of \({\mathcal {O}}(N^d)\). To further speed up the PCG method, an algorithm is developed for fast matrix-vector multiplications by the resulting matrix of Legendre-Galerkin spectral discretization, without the need to explicitly form it. The complexity of the fast matrix-vector multiplications is of \({\mathcal {O}}(N^d (\log _2 N)^{2})\). In view that T is independent of N in one dimension and is set to be of order \({\mathcal {O}}(\log _2 N)\) in two and three dimensions, the PCG method has a \({\mathcal {O}}(N^d (\log _2 N)^{2d})\) quasi-optimal complexity for a d dimensional domain with \((N-1)^d\) unknows, \(d=1,2,3\). In addition, a fast direct solver for the three-dimensional Poisson equation is developed, which is of \({\mathcal {O}}(N^{3} (\log _2 N)^{2})\) and improves the existing results on the computational complexity. Numerical examples are given to demonstrate the efficiency of proposed preconditioners and the algorithm for fast matrix-vector multiplications.

In this paper, through solving Poisson equations, a precise model for double-gates and dual-material carbon nanotube transistors is presented that is more comprehensive than the previous models. Thus, considering a double-gate transistor with a two-material gate, modeling and solving the Poisson equation are taken into consideration in two dimensions for a more accurate analysis of the inner potential of the channel. In this method, the inside electrostatic potential of the channel is obtained by summing two parts of the long channel potential. Also, the dependence of the inside load of the channel on its inside potential is considered in a way that gives a more accurate answer to the potential of the transistor. The potential is generally considered based on the sum of the one-dimensional potential in the channel (without any dependence on the surface potential) and the lateral potential changes in two dimensions. The charge in the Poisson relation is also considered to evaluate the dependence on the potential within the channel which provides a more complete analysis of the Poisson equation.

A parallel numerical acoustic simulation on a GPU using an edge-based smoothed finite element method

The Markov-modulated Poisson process (MMPP) is a generalization of the Poisson process and is commonly used in modeling the input process of communication systems such as data traffic systems and ATM networks. In this paper, we give fast algorithms for solving queueing systems and manufacturing systems with MMPP inputs. We consider queueing systems where the input of the queues is a superposition of the MMPP which is still an MMPP. The generator matrices of these processes are tridiagonal block matrices with each diagonal block being a sum of tensor products of matrices. We are interested in finding the steady state probability distributions of these processes which are the normalized null vectors of their generator matrices. Classical iterative methods, such as the block Gauss-Seidel method, are usually employed to solve for the steady state probability distributions. They are easy to implement, but their convergence rates are slow in general. The number of iterations required for convergence increases like O(m). where m is the size of the waiting spaces in the queues. Here, we propose to use the preconditioned conjugate gradient method. We construct our preconditioners by taking circulant approximations of the tensor blocks of the generator matrices. We show that the number of iterations required for convergence increases at most like O(log2 m) for large m. Numerical results are given to illustrate the fast convergence. As an application, we apply the MMPP to model unreliable manufacturing systems. The production process consists of multiple parallel machines which produce one type of product. Each machine has exponentially distributed up time, down time, and processing time for one unit of product. The interarrival of a demand is exponentially distributed and finite backlog is allowed. We consider hedging point policy as the production control. The average running cost of the system can be written in terms of the steady state probability distribution. Our numerical algorithm developed for the queueing systems can be applied to obtain the steady state distribution for the system and hence the optimal hedging point. Furthermore, our method can be generalized to handle the case when the machines have a more general type of repairing process distribution such as the Erlangian distribution.

Anonoverlapping domain decomposition approach is used on uniform and matching grids to first define and then to compute the orthogonal spline collocation solution of the Dirichlet boundary value problem for Poisson's equation on an L-shaped region. We prove existence and uniqueness of the collocation solution and derive optimal order Hs-norm error bounds for s=0,1,2$. The collocation solution on two interfaces is computed using the preconditioned conjugate gradient method, and the collocation solution on three squares is computed by a matrix decomposition method that uses fast Fourier transforms. The total cost of the algorithm is O(N2 \log N), where the number of unknowns in the collocation solution is O(N2 )$.

A family of factored complex discretizations of the Laplacian is proposed which serve as a basis for fast, second-order accurate Poisson solvers on general two-dimensional regions. One direct implementation of these discretizations is given: a variant of the marching method which is much more stable than the latter, and thus is applicable to grids with relatively large numbers of discretization steps in each direction, without resorting to domain decomposition and multiple shooting. The gain in stability is due to solving initial value problems for twofirst-order difference equations, rather than onesecond-order equation as in the conventional marching method. These initial value problems can be interpreted as backsolves in a sparse Choleski decomposition of the coefficient matrix, induced by the factored discretization operator. The stability and accuracy of the method can be further controlled by the choice of the parameter on which the discretizations depend. The marching algorithms were tested successfully for various geometries. The factored discretization also lend themselves well to an iterative implementation by the preconditioned conjugate gradient method.

Parallel direct Poisson solver for discretisations with one Fourier diagonalisable direction

This paper presents a computational acceleration of image inpainting using parallel processing based on Graphics Processing Unit (GPU) Compute Unified Device Architecture (CUDA). We use parabolic partial differential equation (PDE) called heat equation as the model equation. The heat equation is discretized numerically using Finite Difference method. Semi-algebraic equation that formed then solved by using Alternating-Direction Implicit (ADI) scheme. The numerical algorithm is implemented in GPU CUDA parallel computing to speed up the computational time. The computational process of the inpainting can be done using larger time-step. The computational time can be accelerated to 5.86 times faster using an image with 2736×1824 resolution.

Three-dimensional simulation studies on electrostatic predictions for carbon nanotube field effect transistors

Cloud environments today are increasingly featuring hybrid nodes containing multicore CPU processors and a diverse mix of accelerators such as Graphics Processing Units (GPUs), Intel Xeon Phi co-processors, and Field-Programmable Gate Arrays (FPGAs) to facilitate easier migration to them of HPC workloads. While virtualization of accelerators in clouds is a leading research challenge, we address the programming challenges that assail execution of large instances of data-parallel applications using these accelerators in this paper. In a typical hybrid node in a cloud, the tight integration of accelerators with multicore CPUs via PCI-E communication links contains inherent limitations such as limited main memory of accelerators and limited bandwidth of the PCI-E communication links. These limitations poses formidable programming challenges to execution of large problem sizes on these accelerators. In this paper, we describe a library containing interfaces (HCLOOC) that addresses these challenges. It employs optimal software pipelines to overlap data transfers between host CPU and the accelerator and computations on the accelerator. It is designed using the fundamental building blocks, which are OpenCL command queues for FPGAs, Intel offload streams for Intel Xeon Phis, and CUDA streams and events that allow concurrent utilization of the copy and execution engines provided in NVidia GPUs. We elucidate the key features of our library using an out-of-core implementation of matrix multiplication of large dense matrices on a hybrid node, an Intel Haswell multicore CPU server hosting three accelerators that includes NVidia K40c GPU, Intel Xeon Phi 3120P, and a Xilinx FPGA. Based on experiments with the GPU, we show that our out-of-core implementation achieves 82% of peak double-precision floating performance of the GPU and a speedup of 2.7 times over the NVidia’s out-of-core matrix multiplication implementation (CUBLAS-XT). We also demonstrate that our implementation exhibits 0% drop in performance when the problem size exceeds the main memory of the GPU. We observe this 0% drop also for our implementation for Intel Xeon Phi and Xilinx FPGA.

This paper presents a progressive approach for optimizing the batched LU factorization on graphics processing units (GPUs). The paper shows that the reliance on level-3 BLAS routines for performance does not really pay off, and that it is indeed important to pay attention to the memory-bound part of the algorithm, especially when the problem size is very small. In this context, we develop a size-aware multi-level blocking technique that utilizes different granularities for kernel fusion according to the problem size. Our experiments, which are conducted on a Tesla V100 GPU, show that the multi-level blocking technique achieves speedups for single/double precisions that are up to 3.28×/2.69× against the generic LAPACK-style implementation. It is also up to 8.72×/7.2× faster than the cuBLAS library for single and double precisions, respectively. The developed solution is integrated into the open-source MAGMA library.

Most Computational Fluid Dynamics (CFD) applications solve the Naiver-Stokes equations using various discretization methods like Finite Difference Method (FDM), Finite Element Method (FEM), Finite Volume Method (FVM), etc. The most compute-intensive process in CFD algorithms is the pressure Poisson solver. Solving the pressure Poisson equation requires the solution of a set of linear simultaneous equations. As the problem size/complexity increases, so does the effort and time invested in solving the equations. This study presents fast and robust iterative solvers and the acceleration of the same using the Graphics Processing Unit (GPU) architecture. Serial versions of the codes are implemented in the C language and the parallelization on GPU is achieved using Compute Unified Device Architecture (CUDA). For this study, a two-dimensional heat conduction problem is considered. FDM is used for discretization. Iterative solvers Conjugate Gradient and Multi Grid methods have been implemented in serial and parallel. The performance enhancement with single GPU and CPU system over a single CPU in terms of computing time is reported. It is seen that multigrid algorithms are superior in convergence and nearly 13 times speed-up has been obtained through GPU acceleration using CUDA, for a grid size of 2048 X 2048. The tests have been done on a cluster with Intel(R) Xeon(R) CPU E5-2670 CPU and NVIDIA K10 GPU.

We study a semi-implicit time-difference scheme for magnetohydrodynamics of a viscous and resistive incompressible fluid in a bounded smooth domain with a perfectly conducting boundary.In the scheme, the velocity and magnetic fields are updated by solving simple Helmholtz equations.Pressure is treated explicitly in time, by solving Poisson equations corresponding to a recently developed formula for the Navier-Stokes pressure involving the commutator of Laplacian and Leray projection operators.We prove stability of the time-difference scheme, and deduce a local-time wellposedness theorem for MHD dynamics extended to ignore the divergence-free constraint on velocity and magnetic fields.These fields are divergence-free for all later time if they are initially so.