Cyclic Reduction Algorithm Research Articles

Parallel cyclic reduction of padded bordered almost block diagonal matrices

The solution of linear systems is a crucial and indispensable technique in the field of numerical analysis. Among linear system solvers, the cyclic reduction algorithm stands out for its natural inclination to parallelization. So far, the cyclic reduction has been applied primarily to linear systems with almost block diagonal matrices. Some of its variants widen the usage to almost block diagonal matrices with a last block of rows introducing a set of non-zero elements in the first group of columns. In this work, we extend cyclic reduction to matrices with additional non-zero elements below and to the right without any limitations. These matrices, called padded bordered almost block diagonal, arise from the discretization of optimal control problems featuring arbitrary boundary conditions and free parameters. Nonetheless, they also appear in two-point boundary value problems with free parameters. The proposed algorithm is based on the LU factorizations, and it is designed to be executed in parallel on multi-thread architectures. The algorithm performance is assessed through numerical experiments with different matrix sizes and threads. The computation times and speedups obtained with the parallel implementation indicate that the suggested algorithm is a robust solution for solving padded bordered almost block diagonal linear systems. Furthermore, its structure makes it suitable for the use of different matrix factorization techniques, such as QR or SVD. This flexibility enables tailored customization of the algorithm on the basis of the specific application requirements.

2차원 압축성 유동 해석자의 OpenACC를 이용한 GPU 활용

A two-dimensional compressible flow solver was ported to GPU using OpenACC. The performance of this program was tested in three cases: viscous flow on a flat plate, inviscid transonic flow in a channel with a bump, and inviscid transonic flow around n0012 airfoil. The GPU program running on RTX 3090 showed up to 20 times of speedup compared to the original program using a single core of Ryzen 5950X, and up to 10 times of speedup compared to the program that uses all resources on the CPU. Among the computation steps of the program, implicit residual smoothing was the worst performing step on the GPU. To address this problem, skipping residual smoothing and using the Cyclic Reduction algorithm were tried as alternatives. These approaches improved speedup in almost every case, increasing the maximum speedups to 24 times compared to the single core program and 11 times compared to the 16-core program. The results show that GPU is a suitable tool for accelerating CFD solvers due to its memory bandwidth and parallel computing capability and using OpenACC is a viable method for porting programs to GPU.

An Efficient ADI Method for Transient Thermal Simulation of Liquid-Cooled 3-D ICs

In this article, an efficient approach based on the parallel alternating direction implicit (ADI) method is presented for the transient thermal simulation of liquid-cooled 3-D integrated circuits (ICs). To accurately characterize the thermal behaviors, a compact thermal model based on finite-difference approximation is established. The ADI method with linear time complexity is utilized to reduce computational cost, which can easily deal with large-scale problems. To suppress the oscillations of transient results in the case of large time steps, an effective technique using the exponentially expanding subtime steps is introduced, which improves the accuracy of results significantly. In addition, taking the advantage of parallel computing, the transient simulation is further accelerated by the parallel cyclic reduction (PCR) algorithm implemented on graphics processing units (GPUs). The accuracy and efficiency of the proposed method are demonstrated by several numerical examples.

High-performance computation of pricing two-asset American options under the Merton jump-diffusion model on a GPU

This paper is concerned with fast, parallel and numerically accurate pricing of two-asset American options under the Merton jump-diffusion model, which gives rise to a two-dimensional partial integro-differential complementarity problem (PIDCP) with a nonlocal two-dimensional integral term. Following method-of-lines approach, the solution to the PIDCP can be computed quite accurately by a robust numerical technique that combines Ikonen–Toivanen splitting with an alternating direction implicit scheme. However, we observed that computing the numerical solution with this technique becomes extremely time consuming, mainly due to the handling of the integral term. In this paper we parallelize this technique by applying a parallel fast Fourier transformation algorithm to all matrix-vector multiplications involving the huge and dense integral approximation matrix by exploiting its block Toeplitz with Toeplitz block structure. We also parallelize other computationally intensive steps of this technique by applying a recently developed parallel cyclic reduction algorithm for pentadiagonal systems. Our solutions computed on a graphics processing unit (GPU) using CUDA® platform are compared for accuracy with those available in the literature. It is observed that by solving the PIDCP parallelly we could bring down the computational times from several hours to a few seconds in certain cases in our experiments.

Highly efficient parallel algorithms for solving the Bates PIDE for pricing options on a GPU

In this paper we investigate faster and memory efficient parallel techniques to numerically solve the Bates model for European options. We have followed method-of-linesapproach and implemented the numerical algorithms on a graphics processing unit (GPU). Two second order finite difference (FD) schemes are taken into account that yield discretization matrices with tridiagonal and pentadiagonal block structures. Three recent adaptations of an alternating direction implicit scheme are employed for time-stepping. Spatial and temporal errors corresponding to our chosen FD and time-stepping schemes are numerically studied. For parallel computation of solutions we have applied the well-know parallel cyclic reduction (PCR) algorithm for tridiagonal systems and our novel PCR algorithm for pentadiagonal systems. Ample numerical experiments are performed to study speed and accuracy on three platforms: single GPU using CUDA, multi-core CPU using OpenMP and an efficient sequential algorithm on a single core using MATLAB, where substantial speed-up is observed on the GPU. Sensitivities of computational times of the sequential algorithm in MATLAB with respect to certain parameters in the Bates model are also analysed.

A Parallel Cyclic Reduction Algorithm for Pentadiagonal Systems with Application to a Convection-Dominated Heston PDE

A Parallel Cyclic Reduction Algorithm for Pentadiagonal Systems with Application to a Convection-Dominated Heston PDE

ПАРАЛЛЕЛЬНЫЙ АЛГОРИТМ НА CUDA ДЛЯ РЕШЕНИЯ ЗАДАЧ МНОГОФАЗНОЙ ФИЛЬТРАЦИИ МНОГОКОМПОНЕНТНОЙ ЖИДКОСТИ В ПОРИСТЫХ СРЕДАХ

The paper discusses the implementation of alternating direction implicit parallel algorithm with CUDA technology to solve the problems of multiphase filtration of a multicomponent fluid in porous media. To solve the tridiagonal equations systems of alternating direction implicit method, cyclic and parallel cyclic reduction methods were used. The results of the study showed that the implementation of cyclic and parallel cyclic reduction algorithms on modern GPUs is much more efficient than on CPUs.

A parallel hybrid implementation of the 2D acoustic wave equation

Abstract In this paper, we propose a hybrid parallel programming approach for a numerical solution of a two-dimensional acoustic wave equation using an implicit difference scheme for a single computer. The calculations are carried out in an implicit finite difference scheme. First, we transform the differential equation into an implicit finite-difference equation and then using the alternating direction implicit (ADI) method, we split the equation into two sub-equations. Using the cyclic reduction algorithm, we calculate an approximate solution. Finally, we change this algorithm to parallelize on graphics processing unit (GPU), GPU + Open Multi-Processing (OpenMP), and Hybrid (GPU + OpenMP + message passing interface (MPI)) computing platforms. The special focus is on improving the performance of the parallel algorithms to calculate the acceleration based on the execution time. We show that the code that runs on the hybrid approach gives the expected results by comparing our results to those obtained by running the same simulation on a classical processor core, Compute Unified Device Architecture (CUDA), and CUDA + OpenMP implementations.

On the solving of Poisson’s equation for a cylindrical interaction region

A new numerical method for fast solving of Poisson’s equation in the simply connected interaction regions of O-type vacuum microwave devices is described. The principal specificity of this method, so-called Fourier analysis and eigenfunction decomposition (FAED), consists in using the finite Fourier series for the approximation of the radial dependence of electrostatic potential instead of the well-known Hockney’s cyclic reduction (CR) algorithm. This allows using Neumann boundary condition for the electrostatic potential at the axis of the interaction region instead of Dirichlet one, as in M-type devices.

Parallelization of elliptic solver for solving 1D Boussinesq model

In this paper, a parallel implementation of an elliptic solver in solving 1D Boussinesq model is presented. Numerical solution of Boussinesq model is obtained by implementing a staggered grid scheme to continuity, momentum, and elliptic equation of Boussinesq model. Tridiagonal system emerging from numerical scheme of elliptic equation is solved by cyclic reduction algorithm. The parallel implementation of cyclic reduction is executed on multicore processors with shared memory architectures using OpenMP. To measure the performance of parallel program, large number of grids is varied from 28 to 214. Two test cases of numerical experiment, i.e. propagation of solitary and standing wave, are proposed to evaluate the parallel program. The numerical results are verified with analytical solution of solitary and standing wave. The best speedup of solitary and standing wave test cases is about 2.07 with 214 of grids and 1.86 with 213 of grids, respectively, which are executed by using 8 threads. Moreover, the best efficiency of parallel program is 76.2% and 73.5% for solitary and standing wave test cases, respectively.

A specialised cyclic reduction algorithm for linear algebraic equation systems with quasi-tridiagonal matrices

Extensions have been developed, of several variants of the stride of two cyclic reduction method. The extensions refer to quasi-tridiagonal linear equation systems involving two additional nonzero elements in the first and last rows of the equation matrix, adjacent to the main three diagonals. Equations of this kind arise, for example, in the simulations of biosensors or other electrochemical systems by solving relevant ordinary or partial differential equations by finite difference methods, when boundary derivatives are approximated by one-sided, multipoint finite differences. The correctness of the algorithms developed has been verified using example matrices with pseudo-random coefficients, under conditions of both sequential and parallel execution.

On the decay of the off-diagonal singular values in cyclic reduction

It was recently observed in [10] that the singular values of the off-diagonal blocks of the matrix sequences generated by the Cyclic Reduction algorithm decay exponentially. This property was used to solve, with a higher efficiency, certain quadratic matrix equations encountered in the analysis of queuing models. In this paper, we provide a theoretical bound to the basis of this exponential decay together with a tool for its estimation based on a rational interpolation problem. Numerical experiments show that the bound is often accurate in practice. Applications to solving n×n block tridiagonal block Toeplitz systems with n×n quasiseparable blocks and certain generalized Sylvester equations in O(n2log⁡n) arithmetic operations are shown.

A Hybrid Parallel Solving Algorithm on GPU for Quasi-Tridiagonal System of Linear Equations

There are some quasi-tridiagonal system of linear equations arising from numerical simulations, and some solving algorithms encounter great challenge on solving quasi-tridiagonal system of linear equations with more than millions of dimensions as the scale of problems increases. We present a solving method which mixes direct and iterative methods, and our method needs less storage space in a computing process. A quasi-tridiagonal matrix is split into a tridiagonal matrix and a sparse matrix using our method and then the tridiagonal equation can be solved by the direct methods in the iteration processes. Because the approximate solutions obtained by the direct methods are closer to the exact solutions, the convergence speed of solving the quasi-tridiagonal system of linear equations can be improved. Furthermore, we present an improved cyclic reduction algorithm using a partition strategy to solve tridiagonal equations on GPU, and the intermediate data in computing are stored in shared memory so as to significantly reduce the latency of memory access. According to our experiments on 10 test cases, the average number of iterations is reduced significantly by using our method compared with Jacobi, GS, GMRES, and BiCG respectively, and close to those of BiCGSTAB, BiCRSTAB, and TFQMR. For parallel mode, the parallel computing efficiency of our method is raised by partition strategy, and the performance using our method is better than those of the commonly used iterative and direct methods because of less amount of calculation in an iteration.

Implementation of parallel tridiagonal solvers for a heterogeneous computing environment

Tridiagonal diagonally dominant linear systems arise in many scientific and engineering applications. The standard Thomas algorithm for solving such systems is inherently serial forming a bottleneck in computation. Algorithms such as cyclic reduction and SPIKE reduce a single large tridiagonal system into multiple small independent systems which can be solved in parallel. We have developed portable cyclic reduction and SPIKE algorithm OpenCL implementations with the intent to target a range of co-processors in a heterogeneous computing environment including Field Programmable Gate Arrays (FPGAs), Graphics Processing Units (GPUs) and other multi-core processors. In this paper, we evaluate these designs in the context of solver performance, resource efficiency and numerical accuracy.

An alternative method for computing system-length distributions of BMAP/R/1 and BMAP/D/1 queues using roots

In this paper, we present closed-form expressions for system-length distributions in terms of roots outside the unit disk of the characteristic equation of the BMAP/R/1 queue, where arrival process is batch Markovian arrival process (BMAP) and R represents a class of distributions having rational Laplace–Stieltjes transform. The unknown boundary vector has been evaluated using the roots (inside and on the unit disk) of the characteristic equation. Several numerical results are presented for a variety of arrival and service-time distributions including phase-type (PH) and matrix-exponential (ME) which cover a wide variety of distributions that arise in applications. Using an approximated representation of ME distribution, results for BMAP/D/1 queue are also presented. We compared our result with the results obtained using classical matrix analytic method as well as cyclic reduction algorithm. It is shown that the computation-time of the proposed method does not depend upon the parameters like traffic intensity and correlation co-efficient. The method is analytically quite simple and easy to implement.

A Hybrid Parallel Algorithm to Solve Tridiagonal System of Ocean Model

Aiming at the characteristics of Alternative Direction Implicit (ADI) scheme that is used in the discretization of the governing equations of most ocean models, a hybrid parallel tridiagonal linear system algorithm, based on the Parallel Diagonal Dominant (PDD) and the Parallel Cyclic Reduction (PCR) algorithm, is developed in this paper. The main concept of the hybrid parallel algorithm is that the processes in the parallel system are divided into groups and the PDD algorithm is used to implement parallelization between groups. In each group, the PCR is used to implement the parallelization instead of PDD. By doing so, parallel efficient and computational accuracy could be obtained in the same time. To examine the algorithm of the hybrid parallel efficiency, a Global Reduced Gravity Ocean model (GRGO) is paralyzed, and the results show that this hybrid parallel algorithm can give us a better accuracy and higher efficiency compared to that using the traditional domain decomposition method.

A discontinuous Galerkin method with block cyclic reduction solver for simulating compressible flows on GPUs

An optimized implementation of a block tridiagonal solver based on the block cyclic reduction (BCR) algorithm is introduced and its portability to graphics processing units (GPUs) is explored. The computations are performed on the NVIDIA GTX480 GPU. The results are compared with those obtained on a single core of Intel Core i7-920 (2.67 GHz) in terms of calculation runtime. The BCR linear solver achieves the maximum speedup of 5.84x with block size of 32 over the CPU Thomas algorithm in double precision. The proposed BCR solver is applied to discontinuous Galerkin (DG) simulations on structured grids via alternating direction implicit (ADI) scheme. The GPU performance of the entire computational fluid dynamics (CFD) code is studied for different compressible inviscid flow test cases. For a general mesh with quadrilateral elements, the ADI-DG solver achieves the maximum total speedup of 7.45x for the piecewise quadratic solution over the CPU platform in double precision.

The palindromic cyclic reduction and related algorithms

The cyclic reduction algorithm is specialized to palindromic matrix polynomials and a complete analysis of applicability and convergence is provided. The resulting iteration is then related to other algorithms as the evaluation/interpolation at the roots of unity of a certain Laurent matrix polynomial, the trapezoidal rule for a certain integral and an algorithm based on the finite sections of a tridiagonal block Toeplitz matrix.

Fast finite difference Poisson solvers on heterogeneous architectures

In this paper we propose and evaluate a set of new strategies for the solution of three dimensional separable elliptic problems on CPU–GPU platforms. The numerical solution of the system of linear equations arising when discretizing those operators often represents the most time consuming part of larger simulation codes tackling a variety of physical situations. Incompressible fluid flows, electromagnetic problems, heat transfer and solid mechanic simulations are just a few examples of application areas that require efficient solution strategies for this class of problems. GPU computing has emerged as an attractive alternative to conventional CPUs for many scientific applications. High speedups over CPU implementations have been reported and this trend is expected to continue in the future with improved programming support and tighter CPU–GPU integration. These speedups by no means imply that CPU performance is no longer critical. The conventional CPU-control–GPU-compute pattern used in many applications wastes much of CPU’s computational power. Our proposed parallel implementation of a classical cyclic reduction algorithm to tackle the large linear systems arising from the discretized form of the elliptic problem at hand, schedules computing on both the GPU and the CPUs in a cooperative way. The experimental result demonstrates the effectiveness of this approach.

An efficient GPU implementation of cyclic reduction solver for high-order compressible viscous flow simulations

In this paper, the performance of the Cyclic Reduction (CR) algorithm for solving tridiagonal systems is improved with the aid of efficient global memory transactions on Graphics Processing Units (GPU). To achieve maximum memory throughput with a lower computational runtime, two different Sort algorithms are introduced for reordering the initial system of equations: direct and step-by-step. It is shown that the latter method is well-fitted to modern GPUs and achieves speedup of up to 3.47× in single precision and 2.1× in double precision compared to the CPU Thomas algorithm. By benefiting from the new global memory implementation, the CR solver could run 2×–100× faster compared to previous works on parallel tridiagonal solvers. The CR solver is also applied to 2D & 3D compressible viscous flow simulations using the high-order compact finite-difference scheme. In this matter, the procedure of filtering, primitive variables, and flux derivative calculations are carried out by using the parallel tridiagonal solver on the GPU device. The GPU-accelerated calculations achieve speedups between 1.9×–15.2× in 2D and 6.4×–20.3× in 3D simulations for different grid sizes compared to CPU computations. The computations are performed on the NVIDIA GTX480 GPU. The obtained results are compared to those achieved on a single core of Intel Core 2 Duo (2.7GHz, 2MB cache) in terms of calculation runtime.