Block Tridiagonal Systems Research Articles

Finite-time balanced truncation for linear systems via shifted Legendre polynomials

In this paper, we present a finite-time model order reduction method for linear systems via shifted Legendre polynomials. The main idea of the approach is to use finite-time empirical Gramians, which are constructed from impulse responses by solving block tridiagonal linear systems, to generate approximate balanced system for the large-scale system. The balancing transformation is directly computed from the expansion coefficients of impulse responses in the space spanned by shifted Legendre polynomials, without individual reduction of the Gramians and a separate eigenvector solve. Then, the reduced-order model is constructed by truncating the states corresponding to the small approximate Hankel singular values (HSVs). The stability preservation of the reduced model is briefly discussed. And in combination with the dominant subspace projection method, we modify the reduction procedure to alleviate the shortcomings of the above method, which may unexpectedly lead to unstable systems even though the original one is stable. Furthermore, the properties of the resulting reduced models are considered. Finally, numerical experiments are given to demonstrate the effectiveness of the proposed methods.

Accelerated Cyclic Reduction: A distributed-memory fast solver for structured linear systems

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.

High order numerical algorithms based on biquadratic spline collocation for two-dimensional parabolic partial differential equations

ABSTRACTWe report a new algorithm for solving linear parabolic partial differential equations in two space dimension. The algorithm employs optimal biquadratic spline collocation for space discretization and modified trapezoidal rule for time discretization. We need to solve a block tridiagonal linear system at each time step, and obtain an approximate solution with error at space-time grid points. We analyse the stability of the new algorithm, and present a stability enhanced variant. Moreover, we give an acceleration strategy based on spectral deferred correction, and the theoretical accuracy can be increased to , where k is the number of correction loops. We also analyse the stability for the accelerated algorithms. Numerical experiments are attached to demonstrate the effectiveness of the new algorithms.

On solving separable block tridiagonal linear systems using a GPU implementation of radix-4 PSCR method

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.

Acceleration optimization of real-time equilibrium reconstruction for HL-2A tokamak discharge control

Real-time equilibrium reconstruction is crucially important for plasma shape control in the process of tokamak plasma discharge. However, as the reconstruction algorithm is computationally intensive, it is very difficult to improve its accuracy and reduce the computation time, and some optimizations need to be done. This article describes the three most important aspects of this optimization: (1) compiler optimization; (2) some optimization for middle-scale matrix multiplication on the graphic processing unit and an algorithm which can solve the block tri-diagonal linear system efficiently in parallel; (3) a new algorithm to locate the X&O point on the central processing unit. A static test proves the correctness and a dynamic test proves the feasibility of using the new code for real-time reconstruction with 129 × 129 grids; it can complete one iteration around 575 μs for each equilibrium reconstruction. The plasma displacements from real-time equilibrium reconstruction are compared with the experimental measurements, and the calculated results are consistent with the measured ones, which can be used as a reference for the real-time control of HL-2A discharge.

Unified GPU-Parallelizable Robot Forward Dynamics Computation Using Band Sparsity

This letter proposes a unified GPU-parallelizable approach for robot forward dynamics (FD) computation based on the key fact that parallelism of prevailing FD algorithms benefits from the essential band sparsity of the joint space inertia (JSI) matrix or its inverse. The existing FD algorithms are categorized into three classes: direct JSI algorithms, propagation algorithms, and constraint force algorithms. Their associated systems of linear equations are transformed into a set of block bidiagonal (the first and second classes) and tridiagonal (the third class) linear systems, which can be conveniently and efficiently parallelized over the existing CPU–GPU programming platforms using various state-of-the-art parallel algorithms, such as parallel all-prefix sum (scan) and odd–even elimination. This high-level perspective allows unified and efficient implementation of all three classes of algorithms and also other potentially efficient algorithms, with the bonus that different algorithms can be swiftly compared to recommend problem-specific solutions.

The arrowhead decomposition method for a block-tridiagonal system of linear equations

The arrowhead decomposition method (ADM) for the parallel solution of a block-tridiagonal system of linear equations is presented. The method consists in rearranging the initial linear system into an equivalent one with the “arrowhead” structure of the matrix. It is shown that such a structure provides a good opportunity for parallel solving. The computational speedup of ADM with respect to the sequential matrix Thomas algorithm is analytically estimated based on the number of elementary multiplicative operations for the parallel and serial parts of the methods. A number of parallel processors required to reach the maximum computational speedup is found. A good agreement of the analytical estimations of the computational speedup and practically obtained results is observed.

Parallel cyclic reduction decomposition for dynamic optimization problems

Direct transcription of dynamic optimization problems, with differential-algebraic equations discretized and written as algebraic constraints, can create very large nonlinear optimization problems. When this discretized optimization problem is solved with an NLP solver, such as IPOPT, the dominant computational cost often lies in solving the linear system that generates Newton steps for the KKT system. Computational cost and memory constraints for this linear system solution raise many challenges as the system size increases. On the other hand, the linear KKT system for our dynamic optimization problem is sparse and structured, and can be permuted to form a block tridiagonal matrix. This study explores a parallel decomposition strategy for block tridiagonal systems that is based on cyclic reduction (CR) factorization of the KKT matrix. The classical CR method has good observed performance, but its numerical stability properties need further study for our KKT system. Finally, we discuss modifications to the CR decomposition that improve performance, and we apply the approach to four industrially relevant case studies. On the largest problem, a parallel speedup of a factor of four is observed when using eight processors.

High Resolution Finite Volume Scheme Based on the Quintic Spline Reconstruction on Non-uniform Grids

This paper presents a compact quintic spline reconstruction for finite volume method on non-uniform structured grids. The primitive function of a dependent variable is reconstructed by a piece-wise quintic polynomial by requiring the derivatives up to fourth order being continuous at the cell interfaces. This procedure results in a block tridiagonal system of linear equations which can be solved efficiently by incorporating certain boundary closure relations. There are some distinct features in the quintic spline reconstruction. Firstly, the reconstruction stencil is compact; Secondly, the reconstruction can be applied on arbitrary non-uniform grids; and finally, the reconstruction is continuous at cell interface without the need of a Riemann solver. To stabilize the scheme, the sixth order artificial viscosity is introduced. The quintic spline reconstruction achieves sixth-order accuracy on uniform grids without artificial viscosity and fifth-order accuracy on both uniform and non-uniform grids when artificial viscosity is added. The splined scheme is blended with shock-capturing WENO scheme to suppress non-physical oscillations near discontinuities. Numerical results demonstrate the accuracy, robustness and efficiency of the proposed scheme.

Fast parallel Grad–Shafranov solver for real-time equilibrium reconstruction in EAST tokamak using graphic processing unit**Project supported by the National Magnetic Confinement Fusion Research Program of China (Grant No. 2014GB103000), the National Natural Science Foundation of China (Grant No. 11575245), and the National Natural Science Foundation of China for Youth (Grant No. 11205191).

To achieve real-time control of tokamak plasmas, the equilibrium reconstruction has to be completed sufficiently quickly. For the case of an EAST tokamak experiment, real-time equilibrium reconstruction is generally required to provide results within 1ms. A graphic processing unit (GPU) parallel Grad–Shafranov (G-S) solver is developed in P-EFIT code, which is built with the CUDA architecture to take advantage of massively parallel GPU cores and significantly accelerate the computation. Optimization and implementation of numerical algorithms for a block tri-diagonal linear system are presented. The solver can complete a calculation within with 65 × 65 grid size and with 129 × 129 grid size, and this solver supports that P-EFIT can fulfill the time feasibility for real-time plasma control with both grid sizes.

Perturbation and Error Analyses of the Partitioned LU Factorization for Block Tridiagonal Linear Systems

We present the perturbation and backward error analyses of the partitioned LU factorization for block tridiagonal matrices. In addition, we consider the bounds of perturbations for the partitioned LU factorization for block-tridiagonal linear systems. Finally, numerical examples are given to verify the obtained results.

Fourth-order central compact scheme for the numerical solution of incompressible Navier–Stokes equations

ABSTRACTThis paper provides an implicit central compact scheme for the numerical solution of incompressible Navier–Stokes equations. The solution procedure is based on the artificial compressibility method that transforms the governing equations into a hyperbolic-parabolic form. A fourth-order central compact scheme with a sixth-order numerical filtering is used for the discretization of convective terms and fourth-order central compact scheme for the viscous terms. Dual-time stepping approach is applied to time discretization with backward Euler difference scheme to the pseudo-time derivative, and three point second-order backward difference scheme to the physical time derivative. An approximate factorization-based alternating direction implicit scheme is used to solve the resulting block tridiagonal system of equations. The accuracy and efficiency of the proposed numerical method is verified by simulating several two-dimensional steady and unsteady benchmark problems.

High Accuracy Compact Operator Methods for Two-Dimensional Fourth Order Nonlinear Parabolic Partial Differential Equations

Abstract In this study, we develop and implement numerical schemes to solve classes of two-dimensional fourth-order partial differential equations. These methods are fourth-order accurate in space and second-order accurate in time and require only nine spatial grid points of a single compact cell. The proposed discretizations allow the use of Dirichlet boundary conditions only without the need to discretize the derivative boundary conditions and thus avoids the use of ghost points. No transformation or linearization technique is used to handle nonlinearity and the obtained block tri-diagonal nonlinear system has been solved by Newton’s block iteration method. It is discussed how our formulation is able to tackle linear singular problems and it is ensured that the methods retain their orders and accuracy everywhere in the solution region. The proposed two-level method is shown to be unconditionally stable for a class of two-dimensional fourth-order linear parabolic equation. We also discuss the alternating direction implicit (ADI) method for solving two-dimensional fourth-order linear parabolic equation. The proposed difference methods has been successfully tested on the two-dimensional vibration problem, Boussinesq equation, extended Fisher–Kolmogorov equation and Kuramoto–Sivashinsky equation. Numerical results demonstrate that the schemes are highly accurate in solving a large class of physical problems.

An Artificial Compressibility Method for 3D Phase-Field Model and its Application to Two-Phase Flows

In this work, a numerical scheme based on artificial compressibility formulation of a phase-field model is developed for simulating two-phase incompressible flow problems. The coupled nonlinear systems composed of the incompressible Navier–Stokes equations and volume preserving Allen–Cahn-type phase-field equation are recast into conservative form with source terms, which are suited to implement high-resolution schemes originally developed for hyperbolic conservation laws. The Boussinesq approximation is used to account for the buoyancy effect in flow with small density difference. The fifth-order weighted essentially nonoscillatory (WENO) scheme is used for discretizing the convective terms while dual-time stepping (DTS) technique is used for obtaining time accuracy at each physical time step. Beam–Warming approximate factorization scheme is utilized to obtain block tridiagonal system of equations in each spatial direction. The alternating direction implicit (ADI) algorithm is used to solve the resulting system of equations. The performance of the method is demonstrated by its application to some 2D and 3D benchmark viscous two-phase flow problems.

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.

Accelerated Thomas Solver for (Quasi-)Block-Tridiagonal Linear Algebraic Equation Systems, Using SSE/AVX Instruction Sets for Vectorizing Dense Block Operations

Streaming SIMD Extensions (SSE) and Advanced Vector Extensions (AVX) are additional processor instruction sets available in contemporary personal computers, designed for vectorized floating point calculations. Unfortunately, in order to utilize the advantages of these instructions, one cannot rely on automatic options of high level language compilers. Instead, handwritten assembly language or intrinsic function call insertions are necessary. By using this idea an accelerated C[Formula: see text] code is devised, for solving (quasi-) block-tridiagonal linear algebraic equation systems by means of an extended Thomas algorithm. Speedups reaching 3.5 and 3 (relative to C[Formula: see text] without using SSE/AVX) are demonstrated for single and double precision calculations, respectively.

A parallel solving method for block-tridiagonal equations on CPU–GPU heterogeneous computing systems

Solving block-tridiagonal systems is one of the key issues in numerical simulations of many scientific and engineering problems. Non-zero elements are mainly concentrated in the blocks on the main diagonal for most block-tridiagonal matrices, and the blocks above and below the main diagonal have little non-zero elements. Therefore, we present a solving method which mixes direct and iterative methods. In our method, the submatrices on the main diagonal are 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 block-tridiagonal system of linear equations can be improved. Some direct methods have good performance in solving small-scale equations, and the sub-equations can be solved in parallel. We present an improved algorithm to solve the sub-equations by thread blocks on GPU, and the intermediate data are stored in shared memory, so as to significantly reduce the latency of memory access. Furthermore, we analyze cloud resources scheduling model and obtain ten block-tridiagonal matrices which are produced by the simulation of the cloud-computing system. The computing performance of solving these block-tridiagonal systems of linear equations can be improved using our method.

Efficient cyclic reduction for Quasi-Birth–Death problems with rank structured blocks

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≈102.

Compact high order finite volume method on unstructured grids I: Basic formulations and one-dimensional schemes

The large reconstruction stencil has been the major bottleneck problem in developing high order finite volume schemes on unstructured grids. This paper presents a compact reconstruction procedure for arbitrarily high order finite volume method on unstructured grids to overcome this shortcoming. In this procedure, a set of constitutive relations are constructed by requiring the reconstruction polynomial and its derivatives on the control volume of interest to conserve their averages on face-neighboring cells. These relations result in an over-determined linear equation system, which, in the sense of least-squares, can be reduced to a block-tridiagonal system in the one-dimensional case. The one-dimensional formulations of the reconstruction are discussed in detail and a Fourier analysis is presented to study the dispersion/dissipation and stability properties. The WBAP limiter based on the secondary reconstruction is used to suppress the non-physical oscillations near discontinuities while achieve high order accuracy in smooth regions of the solution. Numerical results demonstrate the method's high order accuracy, robustness and shock capturing capability. A new reconstruction procedure for high order finite volume scheme is proposed.This procedure can achieve arbitrary order of accuracy using a compact stencil consisting of only face neighboring cells.The spectral properties of the proposed finite volume schemes are analyzed.Some test cases are solved numerically to show the superior performances of the proposed schemes.

On relaxed nested factorization and combination preconditioning

The efficient solution of block tridiagonal linear systems arising from the discretization of convection–diffusion problem is considered in this paper. Starting with the classical nested factorization, we propose a relaxed nested factorization preconditioner. Then, several combination preconditioners are developed based on relaxed nested factorization and a tangential filtering preconditioner. Influence of the relaxation parameter is numerically studied, the results indicate that the optimal relaxation parameter should be close to but less than 1. The number of iteration counts exhibit an extremely sensitive behaviour. This phenomena resembles the behaviour of relaxed ILU preconditioner. For symmetric positive-definite coefficient matrix, we also show that the proposed combination preconditioner is convergent. Finally, numerous test cases are carried out with both additive and multiplicative combinations to verify the robustness of the proposed preconditioners.