Search

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.

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.

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.

PhysiCell/BioFVM is an opopen-sourcend agent-based software package supporting 2D/3D simulations for multi-cellular biological systems. It is completely written in C++ and enjoys shared-memory parallelization through OpenMP. An attempt is being made to re-structure the code-base to add support for Distributed parallelization through MPI. The biggest bottleneck in the current version of the simulation package is the serial Thomas solver that is used to solve the triadiagonal system of linear equations resulting from the Finite Volume Discretization Method (FVM) of the reaction-diffusion equations modelling the secretion, ingestion of substrates from cells/agents. Our aim in this project is to replace the serial solver with an efficient distributed-parallel solver. For this purpose we experiment with the Cyclic Reduction (CR) algorithm on a shared-memory system but after understanding its limitations with regard to the problem size, we settle on a modified version of the Thomas solver as our preferred choice in distributed-parallel settings. Our experiments show that the Cyclic Reduction algorithm implemented using OpenMP is able to outperform the serial Thomas solver at a certain thread count on a single node. However, we do not extend this algorithm to support distributed parallelism due to the aforementioned problem. Further, we implement an MPI-only version of the modified Thomas algorithm that promises good scalability on multiple nodes of our HPC cluster - the MareNostrum 4 (MN4) supercomputer at the Barcelona Supercomputing Center (BSC). We project and optimistically conclude that for large problem sizes and a high core count, the parallel modified Thomas algorithm can offer significant reduction in the time to solution for complex 3D simulations in PhysiCell/BioFVM.

Iterative Green's function, based on cyclic reduction of block tridiagonal matrices, has been the ideal algorithm, through tight-binding models, to compute the surface density-of-states of semi-infinite topological electronic materials. In this paper, we apply this method to photonic and acoustic crystals, using finite-element discretizations and a generalized eigenvalue formulation, to calculate the local density-of-states on a single surface of semi-infinite lattices. The three-dimensional (3D) examples of gapless helicoidal surface states in Weyl and Dirac crystals are shown and the computational cost, convergence and accuracy are analyzed., 7 pages, 4 figures

International audience; Predicting the vibratory response of bladed-disks in turbomachinery is of the utmost importance to design a reliable and optimized engine. Yet, such simulations are challenging mainly due to the large size of the finite-element model used to describe the system, the existence of random mistuning coming from manufacture tolerances, and the nonlinear effects arising from the different components and their coupling. As a consequence, very few reduced-order models handling nonlinear mistuned systems have yet been proposed and the existing ones try to find the best compromise between a computationally efficient simulation and a correct description of the nonlinear phenomena. In this paper, a novel approach is proposed to tackle this challenge. Its key ideas rely on the substructuring concept and the cyclic properties of the structure. A reduction basis composed of cyclic complex nonlinear normal modes is created on each sector of the system leading to a compact nonlinear superelement per sector. The system is then assembled and a synthesis simulation is performed. Due to the main concepts employed, the method is referred to as the Substructuring method based on Nonlinear Cyclic Reduction (SNCR). It is validated on a finite element model of a tuned, a randomly mistuned and an intentionally mistuned bladed-disks. It will be shown to be fast and accurate despite strong nonlinearities. Due to its flexibility and efficiency, the method is expected to have a wide range of possible applications within the community.

Partial solution variant of the cyclic reduction (PSCR) method is a direct solver that can be applied to certain types of separable block tridiagonal linear systems. Such linear systems arise, e.g., from the Poisson and the Helmholtz equations discretized with bilinear finite-elements. Furthermore, the separability of the linear system entails that the discretization domain has to be rectangular and the discretization mesh orthogonal. A generalized graphics processing unit (GPU) implementation of the PSCR method is presented. The numerical results indicate up to 24-fold speedups when compared to an equivalent CPU implementation that utilizes a single CPU core. Attained floating point performance is analyzed using roofline performance analysis model and the resulting models show that the attained floating point performance is mainly limited by the off-chip memory bandwidth and the effectiveness of a tridiagonal solver used to solve arising tridiagonal subproblems. The performance is accelerated using off-line autotuning techniques. peerReviewed

In bladed disks, friction nonlinearities occurring within the contact between the blades and disk or between the blades and underplatform dampers are used to decrease the vibratory energy of the system and extend its lifespan. Modeling these nonlinearities and simulating their effects correctly are challenging: complex nonlinear phenomena such as stick, slip and separation may occur in the contact zones. Efficient numerical methods must be used to compute the dynamics of the system in a reasonable amount of time. This paper proposes a reduction strategy to tackle large cyclically symmetric finite-element models undergoing static preload from centrifugal effects and strong nonlinearities (friction and separation). It combines the nonlinear identification of the possible interacting nodal diameters with a linear component mode synthesis procedure. Its potential and performances are assessed through comparisons with some state-of-the art methods. Complex and realistic nonlinear finite-element models of bladed disks with underplatform dampers and dovetail or fir-tree blade roots are used. The Dynamic Lagrangian Frequency Time algorithm is employed to capture the nonlinear effects. Using the proposed reduction method, the effect of underplatform dampers on bladed disk contact occurrence and damping efficiency is investigated.

An efficient parallel multigrid pressure correction algorithm is proposed for the solution of the incompressible Navier–Stokes equations on computing architectures with acceleration devices. The pressure correction procedure is based on the numerical solution of a Poisson-type problem, which is discretized using a fourth-order finite difference compact scheme. Since this is the most time-consuming part of the solver, we propose a parallel pressure correction algorithm using an iterative method based on a block cyclic reduction solution method combined with a multigrid technique. The grid points are numbered with respect to the red–black ordering scheme for the parallel Gauss–Seidel smoother. These parallelization techniques allow the execution of the entire simulation computations on the acceleration device, minimizing memory communication costs. The realization is developed using the OpenACC API, and the numerical method is demonstrated for the solution of two classical incompressible flow test problems. The first is the two-dimensional lid-driven cavity problem over equal mesh sizes while the other is the Stokes boundary layer, which is a decent benchmark problem for unequal mesh spacing. The effect of several multigrid components on modern and legacy acceleration architectures is examined. Eventually the performance investigation demonstrates that the proposed parallel multigrid solver achieves an acceleration of more than 10 $$\times $$ over the sequential solver and more than 4 $$\times $$ over multi-core CPU only realizations for all tested accelerators.

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.

This paper focuses on the resolution of a large number of small random symmetric linear systems and its parallel implementation in single precision on graphics processing units (GPUs). The computations involved by each linear system are independent from the others, and the number of unknowns does not exceed 64. For this purpose, we present the adaptation to our context of largely used methods that include: LDLt factorization, Householder reduction to a tridiagonal matrix, parallel cyclic reduction (PCR) that is not a power of two and the divide and conquer algorithm for tridiagonal eigenproblems. We not only detail the implementation and optimization of each method, but we also compare the sustainability of each solution and its performance which include both parallel complexity and cache memory occupation. In the context of solving a large number of small random linear systems on GPUs with no information about their conditioning, our research indicates that the best strategy requires the use of Householder tridiagonalization + PCR followed if necessary by a divide and conquer diagonalization.

The paper proposes a new approach to the solution of standard and block tridiagonal systems that appear in various areas of technical, scientific and financial practice. Its goal is to elaborate an efficient two-phase tridiagonal solver, the particular case of which is the k -step cyclic reduction. The main idea of the proposed approach to designing a two-phase tridiagonal solver lies in using new model equations for dyadic system reduction. The resulting solver differs from the known two-phase partitioning ones also in the second phase, since it uses a series of simple explicit formulas for calculation of the remaining unknown values. Computational experiments on measuring speedup confirmed the efficiency of the proposed solver.

We study a comparison of serial and parallel solution of 2D-parabolic heat conduction equation using a Crank-Nicolson method with an Alternating Direction Implicit (ADI) scheme. The two-dimensional Heat equation is applied on a thin rectangular aluminum sheet. The forward difference formula is used for time and an averaged second order central difference formula for the derivatives in space to develop the Crank-Nicolson method. FORTRAN serial codes and parallel algorithms using OpenMP are used. Thomas tridigonal algorithm and parallel cyclic reduction methods are employed to solve the tridigonal matrix generated while solving heat equation. This paper emphasize on the run time of both algorithms and their difference. The results are compared and evaluated by creating GNU-plots (Command-line driven graphing utility).

Solvers for linear systems with tridiagonal coefficient matrices sometimes employ direct methods such as the Gauss elimination method or the cyclic reduction method. In each step of the cyclic reduction method, nonzero offdiagonal entries in the coefficient matrix move incrementally away from diagonal entries and eventually vanish. The steps of the cyclic reduction method are generalized as forms of the stride reduction method. For example, the 2-stride reduction method coincides with the 1st step of the cyclic reduction method which transforms tridiagonal linear systems into pentadiagonal systems. In this paper, we explain arbitrary-stride reduction for linear systems with coefficient matrices with three nonzero bands. We then show that arbitrary-stride reduction is equivalent to a combination of 2-stride reduction and row and column permutations. We thus clarify eigenvalue clustering of coefficient matrices in the step-by-step process of the stride reduction method. We also provide two examples verifying this property.

We present an improved version of the locally-one-dimensional finite-difference time-domain method in 2-D formulation, featuring spatial expressions derived by an error-amending procedure. The latter targets at the balanced treatment of space–time flaws, and assures the efficient exploitation of four-point spatial operators that now depend on the time-step size. Simulations are also shortened by implementing the algorithm on graphics processing units, a process that necessitates the incorporation of a parallelizable implicit solver. This is accomplished by adapting parallel cyclic reduction to heptadiagonal systems. Therefore, a reliable computational framework free of time-step restrictions is developed.

Matrix equations of the kind $A_1 X^2 + A0 X + A_{-1} = X$, where both the matrix coefficients and the unknown are semi-infinite matrices belonging to a Banach algebra, are considered. These equations, where coefficients are quasi-Toeplitz matrices, are encountered in certain quasi-birth-death processes as the tandem Jackson queue or in any other processes that can be modeled as a reflecting random walk in the quarter plane. We provide a numerical framework for approximating the minimal nonnegative solution of these equations that relies on semi-infinite quasi-Toeplitz matrix arithmetic. In particular, we show that the algorithm of cyclic reduction can be effectively applied and can approxi- mate the infinite-dimensional solutions with quadratic convergence at a cost that is comparable to that of the finite case. This way, we may compute a finite approximation of the sought solution and of the invariant probability measure of the associated quasi-birth-death process, within a given accuracy. Numerical experiments, performed on a collection of benchmarks, confirm the theoretical analysis.

[EN] We consider the computation of a few eigenpairs of a generalized eigenvalue problem Ax = lambda Bx with block-tridiagonal matrices, not necessarily symmetric, in the context of Krylov methods. In this kind of computation, it is often necessary to solve a linear system of equations in each iteration of the eigensolver, for instance when B is not the identity matrix or when computing interior eigenvalues with the shift-and-invert spectral transformation. In this work, we aim to compare different direct linear solvers that can exploit the block-tridiagonal structure. Block cyclic reduction and the Spike algorithm are considered. A parallel implementation based on MPI is developed in the context of the SLEPc library. The use of GPU devices to accelerate local computations shows to be competitive for large block sizes., This work was supported by Agencia Estatal de Investigacion (AEI) under grant TIN2016-75985-P, which includes European Commission ERDF funds. Alejandro Lamas Davina was supported by the Spanish Ministry of Education, Culture and Sport through a grant with reference FPU13-06655.

Variational methods are commonly used to solve noise removal problems. In this paper, we present an augmented Lagrangian-based approach that uses a discrete form of the L1-norm of the mean curvature of the graph of the image as a regularizer, discretization being achieved via a finite element method. When a particular alternating direction method of multipliers is applied to the solution of the resulting saddle-point problem, this solution reduces to an iterative sequential solution of four subproblems. These subproblems are solved using Newton’s method, the conjugate gradient method, and a partial solution variant of the cyclic reduction method. The approach considered here differs from existing augmented Lagrangian approaches for the solution of the same problem; indeed, the augmented Lagrangian functional we use here contains three Lagrange multipliers “only,” and the associated augmentation terms are all quadratic. In addition to the description of the solution algorithm, this paper contains the results of numerical experiments demonstrating the performance of the novel method discussed here. peerReviewed

The tri-diagonal system comes from dynamic problems such as fluid simulation, and high efficiency is important for the success of these applications. In this paper, we develop completely GPU shared memory-based chunked cyclic reduction under the constraint of the capacity of the shared memory. Computational results show that GPU shared memory chunked cyclic reduction exhibits high efficiency by Nvidia TITAN with 48k shared memory, and GPU shared memory chunked cyclic reduction can solve a tri-diagonal system with 262,144-by-262,144 coefficient matrix in 1.768 ms. Computational results also show that GPU shared memory chunked cyclic reduction scales well to the sizes of coefficient matrix and the reduced systems. Altogether, since building completely on GPU shared memory, our solver may be faster than existing GPU solvers because of the efficiency of GPU shared memory, though the solubility of our solver is smaller than existing GPU solvers because of the capacity constraint of shared memory, where solubility means the solvable tri-diagonal system with the maximum size of the coefficient matrix by our solver.

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.

SUMMARY A conventional block cyclic reduction algorithm operates by halving the size of the linear system at each reduction step, that is, the algorithm is a radix-2 method. An algorithm analogous to the block cyclic reduction known as the radix-q partial solution variant of the cyclic reduction (PSCR) method allows the use of higher radix numbers and is thus more suitable for parallel architectures as it requires fever reduction steps. This paper presents an alternative and more intuitive way of deriving a radix-4 block cyclic reduction method for systems with a coefficient matrix of the form tridiag{ − I,D, − I}. This is performed by modifying an existing radix-2 block cyclic reduction method. The resulting algorithm is then parallelized by using the partial fraction technique. The parallel variant is demonstrated to be less computationally expensive when compared to the radix-2 block cyclic reduction method in the sense that the total number of emerging subproblems is reduced. The method is also shown to be numerically stable and equivalent to the radix-4 PSCR method. The numerical results archived correspond to the theoretical expectations. Copyright © 2013 John Wiley & Sons, Ltd.

A parallel fast direct solver for rank-compressible block tridiagonal linear systems is presented. Algorithmic synergies between Cyclic Reduction and Hierarchical matrix arithmetic operations result in a solver with O(Nlog2N) arithmetic complexity and O(NlogN) memory footprint. We provide a baseline for performance and applicability by comparing with well-known implementations of the \(\mathcal{H}\)-LU factorization and algebraic multigrid within a shared-memory parallel environment that leverages the concurrency features of the method. Numerical experiments reveal that this method is comparable with other fast direct solvers based on Hierarchical Matrices such as \(\mathcal{H}\)-LU and that it can tackle problems where algebraic multigrid fails to converge.

We provide effective algorithms for solving block tridiagonal block Toeplitz systems with m×m quasiseparable blocks, as well as quadratic matrix equations with m×m quasiseparable coefficients, based on cyclic reduction and on the technology of rank-structured matrices. The algorithms rely on the exponential decay of the singular values of the off-diagonal submatrices generated by cyclic reduction. We provide a formal proof of this decay in the Markovian framework. The results of the numerical experiments that we report confirm a significant speed up over the general algorithms, already starting with the moderately small size m ≈ 100. ispartof: Applied Numerical Mathematics vol:116 pages:37-46 ispartof: location:ITALY, Falerna status: published

In this paper, we investigate both experimentally and theoretically the dynamics of a liquid plug driven by a cyclic periodic forcing inside a cylindrical rigid capillary tube. First, it is shown that, depending on the type of forcing (flow rate or pressure cycle), the dynamics of the liquid plug can either be stable and periodic, or conversely accelerative and eventually leading to plug rupture. In the latter case, we identify the sources of the instability as: (i) the cyclic diminution of the plug viscous resistance to motion due to the decrease in the plug length and (ii) a cyclic reduction of the plug interfacial resistance due to a lubrication effect. Since the flow is quasi-static and the forcing periodic, this cyclic evolution of the resistances relies on the existence of flow memories stored in the length of the plug and the thickness of the trailing film. Second, we show that, contrary to unidirectional pressure forcing, cyclic forcing enables breaking of large plugs in a confined space although it requires longer times. All the experimentally observed tendencies are quantitatively recovered from an analytical model. This study not only reveals the underlying physics but also opens up the prospect for the simulation of ‘breathing’ of liquid plugs in complex geometries and the determination of optimal cycles for obstructed airways reopening.

We present Accelerated Cyclic Reduction (ACR), a distributed-memory fast direct solver for rank-compressible block tridiagonal linear systems arising from the discretization of elliptic operators, developed here for three dimensions. Algorithmic synergies between Cyclic Reduction and hierarchical matrix arithmetic operations result in a solver that has $O(k~N \log N~(\log N + k^2))$ arithmetic complexity and $O(k~N \log N)$ memory footprint, where $N$ is the number of degrees of freedom and $k$ is the rank of a typical off-diagonal block, and which exhibits substantial concurrency. We provide a baseline for performance and applicability by comparing with the multifrontal method where hierarchical semi-separable matrices are used for compressing the fronts, and with algebraic multigrid. Over a set of large-scale elliptic systems with features of nonsymmetry and indefiniteness, the robustness of the direct solvers extends beyond that of the multigrid solver, and relative to the multifrontal approach ACR has lower or comparable execution time and memory footprint. ACR exhibits good strong and weak scaling in a distributed context and, as with any direct solver, is advantageous for problems that require the solution of multiple right-hand sides., Comment: 22 pages, Elsevier Journal of Parallel Computing, Dec 2016

It was recently observed 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 𝑛 × 𝑛 block tridiagonal block Toeplitz systems with 𝑛 × 𝑛 semiseparable blocks and certain generalized Sylvester equations in 𝑂(𝑛²log 𝑛) arithmetic operations are shown. ispartof: Linear Algebra and Its Applications vol:519 pages:27-53 status: published

We present a parallel Alternating Direction Implicit (ADI) solver on GPUs. Our implementation significantly improves ex- isting implementations in two aspects. First, we address the scalability issue of existing Parallel Cyclic Reduction (PCR) implementations by eliminating their hardware resource constraints. As a result, our parallel ADI, which is based on PCR, no longer has the maximum domain size limitation. Second, we optimize inefficient data accesses of parallel ADI solver by leveraging hardware texture memory and matrix transpose techniques. These memory optimizations further make already parallelized ADI solver twice faster, achieving overall more than 100 times speedup over a highly optimized CPU version. We also present the analysis of numerical accuracy of the proposed parallel ADI solver.

Engineering, scientific, and financial applications often require the simultaneous solution of a large number of independent tridiagonal systems of equations with varying coefficients. Since the number of systems is large enough to offer considerable parallelism on manycore systems, the choice between different tridiagonal solution algorithms, such as Thomas, Cyclic Reduction (CR) or Parallel Cyclic Reduction (PCR) needs to be reexamined. This work investigates the optimal choice of tridiagonal algorithm for CPU, Intel MIC, and NVIDIA GPU with a focus on minimizing the amount of data transfer to and from the main memory using novel algorithms and the register-blocking mechanism, and maximizing the achieved bandwidth. It also considers block tridiagonal solutions, which are sometimes required in Computational Fluid Dynamic (CFD) applications. A novel work-sharing and register blocking--based Thomas solver is also presented.

We have previously suggested mixed precision iterative solvers specifically tailored to the iterative solution of sparse linear equation systems as they typically arise in the finite element discretization of partial differential equations. These schemes have been evaluated for a number of hardware platforms, in particular, single-precision GPUs as accelerators to the general purpose CPU. This paper reevaluates the situation with new mixed precision solvers that run entirely on the GPU: We demonstrate that mixed precision schemes constitute a significant performance gain over native double precision. Moreover, we present a new implementation of cyclic reduction for the parallel solution of tridiagonal systems and employ this scheme as a line relaxation smoother in our GPU-based multigrid solver. With an alternating direction implicit variant of this advanced smoother, we can extend the applicability of the GPU multigrid solvers to very ill-conditioned systems arising from the discretization on anisotropic meshes, that previously had to be solved on the CPU. The resulting mixed-precision schemes are always faster than double precision alone, and outperform tuned CPU solvers consistently by almost an order of magnitude.

We study the performance of three parallel algorithms and their hybrid variants for solving tridiagonal linear systems on a GPU: cyclic reduction (CR), parallel cyclic reduction (PCR) and recursive doubling (RD). We develop an approach to measure, analyze, and optimize the performance of GPU programs in terms of memory access, computation, and control overhead. We find that CR enjoys linear algorithm complexity but suffers from more algorithmic steps and bank conflicts, while PCR and RD have fewer algorithmic steps but do more work each step. To combine the benefits of the basic algorithms, we propose hybrid CR+PCR and CR+RD algorithms, which improve the performance of PCR, RD and CR by 21%, 31% and 61% respectively. Our GPU solvers achieve up to a 28x speedup over a sequential LAPACK solver, and a 12x speedup over a multi-threaded CPU solver.