Block Tridiagonal Linear Systems Research Articles

Model order reduction based on low-rank decomposition of the cross Gramian

This paper considers model order reduction based on low-rank decomposition of the cross Gramian for linear systems. The proposed approach uses the low-rank factors of the cross Gramian to generate approximate balanced model for the large-scale system. Then, the reduced-order model is obtained by truncating the states corresponding to the small approximate Hankel singular values. The low-rank factors are directly constructed from the expansion coefficients of impulse responses in the space spanned by Legendre polynomials by solving block tridiagonal linear systems. In contrast to balanced truncation related approaches which require to solve Sylvester equation to compute the full cross Gramian, the proposed method needs to solve sparse linear equations, and only one singular value decomposition is applied to a low-dimensional matrix. The stability preservation of the reduced model is briefly discussed. And in combination with the dominant subspace projection method, the reduction procedure is modified to alleviate the shortcoming, which may unexpectedly lead to unstable systems even though the original one is stable. Finally, numerical experiments are given to demonstrate the effectiveness of the proposed methods.

Newton-Based Methods for Finding the Positive Ground State of Gross-Pitaevskii Equations

The discretization of Gross-Pitaevskii equations (GPE) leads to a nonlinear eigenvalue problem with eigenvector nonlinearity (NEPv). We use two Newton-based methods to compute the positive ground state of GPE. The first method comes from the Newton-Noda iteration for saturable nonlinear Schrödinger equations proposed by Ching-Sung Liu, which can be transferred to GPE naturally. The second method combines the idea of the root-finding methods and the idea of Newton method, in which, each subproblem involving block tridiagonal linear systems can be solved easily. We give an explicit convergence and computational complexity analysis for it. Numerical experiments are provided to support the theoretical results.

An interpolatory directional splitting-local discontinuous Galerkin method with application to pattern formation in 2D/3D

An efficient computational method to approximate the solution of a general class of nonlinear reaction-diffusion systems in Cartesian grids is presented. The proposed scheme uses the Local Discontinuous Galerkin (LDG) method as spatial discretization and the symmetric Strang operator splitting as time marching scheme. As a result, not only diffusion and reaction are decoupled; but the reaction sub-step reduces to a set of completely independent small local nonlinear systems of ordinary differential equations. To reduce their computational cost a Lagrange interpolatory technique is used. Diffusion sub-steps are approximated by a directional splitting method derived from an algebraic decomposition of the discrete Laplacian operator, which amounts to solve sets of completely independent small block tridiagonal linear systems; one set per direction. Stability of the diffusive step is analyzed. Finally the efficiency and accuracy of the proposed scheme is numerically investigated on well known Turing pattern formation models in 2D/3D.

Algorithms for Block Tridiagonal Systems: Stability Results for Generalized Kalman Smoothing

Abstract Block tridiagonal systems appear in classic Kalman smoothing problems, as well in generalized Kalman smoothing, where problems may have nonsmooth terms, singular covariance, constraints, nonlinear models, and unknown parameters. In this paper, first we interpret all the classic smoothing algorithms as different approaches to solve positive definite block tridiagonal linear systems. Then, we obtain new results on their numerical stability. Our outcomes apply to all systems with dynamic structure, informing both classic and modern inference for generalized Kalman smoothing.

A fast method for solving a block tridiagonal quasi-Toeplitz linear system

This paper addresses the problem of solving block tridiagonal quasi-Toeplitz linear systems. Inspired by [10], we propose a more general algorithm for such systems. The algorithm is based on a block decomposition for block tridiagonal quasi-Toeplitz matrices and the Sherman–Morrison–Woodbury inversion formula. We also compare the proposed approach to the standard block LU decomposition method and the Gauss algorithm. A theoretical error analysis is also presented. All algorithms have been implemented in Matlab. Numerical experiments performed on a wide variety of test problems show the effectiveness of our algorithm in terms of efficiency, stability and robustness.

Modified method of parallel matrix sweep

The topic of this paper refers to efficient parallel solvers of block-tridiagonal linear systems of equations. Such systems occur in numerous modeling problems and require usage of high-performance multicore computation systems. One of the widely used methods for solving block-tridiagonal linear systems in parallel is the original block-tridiagonal sweep method. We consider the algorithm based on the partitioning idea. Firstly, the initial matrix is split into parts and transformations are applied to each part independently to obtain equations of a reduced block-tridiagonal system. Secondly, the reduced system is solved sequentially using the classic Thomas algorithm. Finally, all the parts are solved in parallel using the solutions of a reduced system. We propose a modification of this method. It was justified that if known stability conditions for the matrix sweep method are satisfied, then the proposed modification is stable as well.

Structure preserving balanced proper orthogonal decomposition for second‐order form systems via shifted Legendre polynomials

This study considers structure preserving balanced proper orthogonal decomposition for second-order form systems via shifted Legendre polynomials. The proposed approach is to use time interval empirical Gramians of the first-order representation, which are constructed from impulse responses by solving block tridiagonal linear systems, to generate approximate balanced system for the large-scale second-order form system. The balancing transformation is directly computed from the expansion coefficients of impulse responses in the space spanned by shifted Legendre polynomials, without computing the full Gramians for the first-order representation. Then, the reduced second-order model is constructed by truncating the states corresponding to the small approximate singular values. Furthermore, in combination with the dominant subspace projection method, the authors modify the reduction procedure to alleviate the shortcoming of the above method, which may unexpectedly lead to unstable systems even though the original one is stable. The stability preservation of the reduced model is briefly discussed. Finally, the effectiveness of the proposed methods is demonstrated by two numerical experiments.

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.

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.

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.

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.

An efficient discontinuous Galerkin finite element method with nested domain decomposition for simulations of microresistivity imaging

A discontinuous Galerkin finite element method is employed to study the responses of microresistivity imaging tools used in the oil and gas exploration industry. The multiscale structure of an imaging problem is decomposed into several nested subdomains based on its geometric characteristics. Each subdomain is discretized independently, and numerical flux is used to couple all subdomains together. The nested domain decomposition scheme will lead to a block tridiagonal linear system, and the block Thomas algorithm is utilized here to eliminate the subdomain based iteration in the step of solving the linear system. Numerical results demonstrate the validity and efficiency of this method.

Spectral-Prism Element for Multi-Scale Layered Package-Chip Co-Simulations Using the Discontinuous Galerkin Time-Domain Method

A new kind of prism element with a triangular base is presented for discretization of multi-scale layered structures within the discontinuous Galerkin time-domain framework. Mixed-order curl-conforming vector basis functions are used in the triangular bases of the prismatic element. The height of the prism adopts spectral basis functions based on Gauss–Lobatto–Legendre polynomials, with an arbitrary order of interpolation. This method combines the flexibility of triangles with the accuracy of spectral elements for layered structures. Eigenvalues obtained show better results than traditional finite elements using tetrahedrons and hexahedrons. For transient analysis, the implicit Crank–Nicholson method is implemented for sequential sub-domains. This kind of arrangement of sub-domains produces a block tridiagonal linear system, thus allowing a block Thomas algorithm to solve the system efficiently. A package-to-chip example shows the efficacy of this method.

Using Compiler Directives to Port Large Scientific Applications to GPUs: An Example from Atmospheric Science

For many scientific applications, Graphics Processing Units (GPUs) can be an interesting alternative to conventional CPUs as they can deliver higher memory bandwidth and computing power. While it is conceivable to re-write the most execution time intensive parts using a low-level API for accelerator programming, it may not be feasible to do it for the entire application. But, having only selected parts of the application running on the GPU requires repetitively transferring data between the GPU and the host CPU, which may lead to a serious performance penalty. In this paper we assess the potential of compiler directives, based on the OpenACC standard, for porting large parts of code and thus achieving a full GPU implementation. As an illustrative and relevant example, we consider the climate and numerical weather prediction code COSMO (Consortium for Small Scale Modeling) and focus on the physical parametrizations, a part of the code which describes all physical processes not accounted for by the fundamental equations of atmospheric motion. We show, by porting three of the dominant parametrization schemes, the radiation, microphysics and turbulence parametrizations, that compiler directives are an efficient tool both in terms of final execution time as well as implementation effort. Compiler directives enable to port large sections of the existing code with minor modifications while still allowing for further optimizations for the most performance critical parts. With the example of the radiation parametrization, which contains the solution of a block tri-diagonal linear system, the required code modifications and key optimizations are discussed in detail. Performance tests for the three physical parametrizations show a speedup of between 3× and 7× for execution time obtained on a GPU and on a multi-core CPU of an equivalent generation.