Nonsymmetric Systems Of Linear Equations Research Articles

Accelerating Restarted GMRES With Mixed Precision Arithmetic

The generalized minimum residual method (GMRES) is a commonly used iterative Krylov solver for sparse, non-symmetric systems of linear equations. Like other iterative solvers, data movement dominates its run time. To improve this performance, we propose running GMRES in reduced precision with key operations remaining in full precision. Additionally, we provide theoretical results linking the convergence of finite precision GMRES with classical Gram-Schmidt with reorthogonalization (CGSR) and its infinite precision counterpart which helps justify the convergence of this method to double-precision accuracy. We tested the mixed-precision approach with a variety of matrices and preconditioners on a GPU-accelerated node. Excluding the incomplete LU factorization without fill in (ILU(0)) preconditioner, we achieved average speedups ranging from 8 to 61 percent relative to comparable double-precision implementations, with the simpler preconditioners achieving the higher speedups.

A novel block non-symmetric preconditioner for mixed-hybrid finite-element-based Darcy flow simulations

In this work, we propose a novel block preconditioner, labeled Explicit Decoupling Factor Approximation (EDFA), to accelerate the convergence of Krylov subspace solvers used to address the sequence of non-symmetric systems of linear equations originating from flow simulations in porous media. The flow model is discretized blending the Mixed hybrid finite element method for Darcy's equation with the Finite volume scheme for the mass conservation. The EDFA preconditioner is characterized by two features: the exploitation of the system matrix decoupling factors to recast the Schur complement and their inexact fully-parallel computation by means of restriction operators. We introduce two adaptive techniques aimed at building the restriction operators according to the properties of the system at hand. The proposed block preconditioner has been tested through extensive experimentation on both synthetic and real-case applications, pointing out its robustness and computational efficiency.

PARALLEL APPROACH TO SOLVE OF THE DIRECT SOLUTION OF LARGE SPARSE SYSTEMS OF LINEAR EQUATIONS

The paper deals with parallel approach for the numerical solution of large, sparse, non-symmetric systems of linear equations, that can be part of any finite element software. In this contribution, the differences between the sequential and parallel solution are highlighted and the approach to efficiently interface with distributed memory version of SuperLU solver is described.

Evaluation of Different Approaches to Solution of the Direct Solution of Large, Sparse Systems of Linear Equations

The aim of this paper is to the evaluate efficiency of differentapproaches to solution of large, sparse, non-symmetric systems of linearequations on high performance machines, that can be found in any finiteelement software. The different approaches based on direct or iterativealgorithms for solution of linear equations are compared. In particular,directs solver using Skyline sparse storage, direct solver from SuperLUlibrary, iterative solver from Iterative Method Library(IML)are compared. SuperLU is a general purpose library for the directsolution of large, sparse, nonsymmetric systems of linear equations.Additionally, the performance and scalability of parallel SuperLU solveris studied, based on OpenMP. The paper shows thatparallelization can efficiently exploit the power of modern availablehardware, significantly reducing the needed computation time.The different strategies were implemented in OOFEM which is afree finite element code with object oriented architecture for solvingmechanical, transport and fluid mechanics problems that operates onvarious platforms.

The Non–Symmetric s–Step Lanczos Algorithm: Derivation of Efficient Recurrences and Synchronization–Reducing Variants of BiCG and QMR

Abstract The Lanczos algorithm is among the most frequently used iterative techniques for computing a few dominant eigenvalues of a large sparse non-symmetric matrix. At the same time, it serves as a building block within biconjugate gradient (BiCG) and quasi-minimal residual (QMR) methods for solving large sparse non-symmetric systems of linear equations. It is well known that, when implemented on distributed-memory computers with a huge number of processes, the synchronization time spent on computing dot products increasingly limits the parallel scalability. Therefore, we propose synchronization-reducing variants of the Lanczos, as well as BiCG and QMR methods, in an attempt to mitigate these negative performance effects. These so-called s-step algorithms are based on grouping dot products for joint execution and replacing time-consuming matrix operations by efficient vector recurrences. The purpose of this paper is to provide a rigorous derivation of the recurrences for the s-step Lanczos algorithm, introduce s-step BiCG and QMR variants, and compare the parallel performance of these new s-step versions with previous algorithms.

Mixed-Precision Cholesky QR Factorization and Its Case Studies on Multicore CPU with Multiple GPUs

To orthonormalize the columns of a dense matrix, the Cholesky QR (CholQR) requires only one global reduction between the parallel processing units and performs most of its computation using BLAS-3 kernels. As a result, compared to other orthogonalization algorithms, CholQR obtains superior performance on many of the current computer architectures, where the communication is becoming increasingly expensive compared to the arithmetic operations. This is especially true when the input matrix is tall-skinny. Unfortunately, the orthogonality error of CholQR depends quadratically on the condition number of the input matrix, and it is numerically unstable when the matrix is ill-conditioned. To enhance the stability of CholQR, we recently used mixed-precision arithmetic; the input and output matrices are in the working precision, but some of its intermediate results are accumulated in the doubled precision. In this paper, we analyze the numerical properties of this mixed-precision CholQR. Our analysis shows that by selectively using the doubled precision, the orthogonality error of the mixed-precision CholQR only depends linearly on the condition number of the input matrix. We provide numerical results to demonstrate the improved numerical stability of the mixed-precision CholQR in practice. We then study its performance. When the target hardware does not support the desired higher precision, software emulation is needed. For example, using software-emulated double-double precision for the working 64-bit double precision, the mixed-precision CholQR requires about $8.5\times$ more floating-point instructions than that required by the standard CholQR. On the other hand, the increase in the communication cost using the double-double precision is less significant, and our performance results on multicore CPU with a different graphics processing unit (GPU) demonstrate that the overhead of using the double-double arithmetic is decreasing on a newer architecture, where the computation is becoming less expensive compared to the communication. As a result, with a latest NVIDIA GPU, the mixed-precision CholQR was only $1.4\times$ slower than the standard CholQR. Finally, we present case studies of using the mixed-precision CholQR within communication-avoiding variants of Krylov subspace projection methods for solving a nonsymmetric linear system of equations and for solving a symmetric eigenvalue problem, on a multicore CPU with multiple GPUs. These case studies demonstrate that by using the higher precision for this small but critical segment of the Krylov methods, we can improve not only the overall numerical stability of the solvers but also, in some cases, their performance.

Polynomial Preconditioned GMRES and GMRES-DR

We look at solving large nonsymmetric systems of linear equations using polynomial preconditioned Krylov methods. We give a simple way to find the polynomial. It is shown that polynomial preconditioning can significantly improve restarted GMRES for difficult problems, and the reasons for this are examined. Stability is discussed, and algorithms are given for increased stability. Next, we apply polynomial preconditioning to GMRES with deflated restarting. It is shown that this is worthwhile for sparse matrices and for problems with many small eigenvalues. Multiple right-hand sides are also considered.

Global simpler GMRES for nonsymmetric systems with multiple right-hand sides

A new version of the simpler GMRES for solving nonsymmetric linear system of equations with multiple right-hand sides is presented. The main advantage over the simpler block GMRES is that the new method requires less storage and computational work per cycle. Preconditioning techniques can also be used to further accelerate the convergence of the new method. Finally, numerical examples are given to illustrate the effectiveness of the proposed method and its some preconditioned versions.

An FFT-based Galerkin method for homogenization of periodic media

In 1994, Moulinec and Suquet introduced an efficient technique for the numerical resolution of the cell problem arising in homogenization of periodic media. The scheme is based on a fixed-point iterative solution to an integral equation of the Lippmann–Schwinger type, with action of its kernel efficiently evaluated by the Fast Fourier Transform techniques. The aim of this work is to demonstrate that the Moulinec–Suquet setting is actually equivalent to a Galerkin discretization of the cell problem, based on approximation spaces spanned by trigonometric polynomials and a suitable numerical integration scheme. For the latter framework and scalar elliptic problems, we prove convergence of the approximate solution to the weak solution, including a-priori estimates for the rate of convergence for sufficiently regular data and the effects of numerical integration. Moreover, we also show that the variational structure implies that the resulting non-symmetric system of linear equations can be solved by the conjugate gradient method. Apart from providing a theoretical support to Fast Fourier Transform-based methods for numerical homogenization, these findings significantly improve on the performance of the original solver and pave the way to similar developments for its many generalizations proposed in the literature.

A quasi-minimal residual variant of the BiCORSTAB method for nonsymmetric linear systems

The Biconjugate A-orthogonal residual stabilized method named as BiCORSTAB was proposed by Jing et al. (2009), where the numerical experiments therein demonstrate that the BiCORSTAB method converges more smoothly than the Bi-Conjugate Gradient stabilized (BiCGSTAB) method in some circumstances. In order to further stabilize the convergence behavior and hopefully to accelerate the convergence speed of the BiCORSTAB algorithm when it has erratic convergence curves, a quasi-minimal residual variant of the BiCORSTAB algorithm, named as QMRCORSTAB, will be developed and investigated for solving nonsymmetric systems of linear equations borrowing the same further-smooth-effect idea for the QMRCGSTAB method. Numerical experiments on typical sets of both sparse and dense matrices will show that the proposed QMRCORSTAB method shares attractive smoother effect over its basic parent and also outperforms its counterpart.

Adaptive finite element modeling of marine controlled-source electromagnetic fields in two-dimensional general anisotropic media

In this paper, we extend the scope of numerical simulations of marine controlled-source electromagnetic (CSEM) fields in a particular case of anisotropy (dipping anisotropy) to the general case of anisotropy by using an adaptive finite element approach. In comparison to a dipping anisotropy case, the first order spatial derivatives of the strike-parallel components arise in the partial differential equations for generally anisotropic media, which cause a non-symmetric linear system of equations for finite element modeling. The adaptive finite element method is employed to obtain numerical solutions on a sequence of refined unstructured triangular meshes, which allows for arbitrary model geometries including bathymetry and dipping layers. Numerical results of a 2D anisotropic model show both anisotropy strike and dipping angles have great influence on the marine CSEM responses.

Parallel Implementations of FGMRES for Solving Large, Sparse Non-symmetric Linear Systems

The Flexible Generalized Minimal Residual method (FGMRES) is an attractive iterative solver for non-symmetric systems of linear equations. This paper presents several methods for parallelizing FGMRES for a variety of archi- tectures including multi-core CPU, Graphics Processing Units (GPU), and multi-GPU systems. The parallel imple- mentations utilize OpenMP and CUDA kernels, and are organized according to thread scope. The linear systems employed in this study correspond to the discrete analogues of realistic three-dimensional convection-diffusion problems, and range in size to nearly 107 linear equations. All of the parallel implementations, particularly the novel hybrid approach, show a significant speedup over the sequential version. Theoretical insight and perfor- mance data is provided to allow informed decisions as to the most effective parallelization method for a given architecture.

The BiConjugate gradient method on GPUs

In a wide variety of applications from different scientific and engineering fields, the solution of complex and/or nonsymmetric linear systems of equations is required. To solve this kind of linear systems the BiConjugate Gradient method (BCG) is especially relevant. Nevertheless, BCG has a enormous computational cost. GPU computing is useful for accelerating this kind of algorithms but it is necessary to develop suitable implementations to optimally exploit the GPU architecture. In this paper, we show how BCG can be effectively accelerated when all operations are computed on a GPU. So, BCG has been implemented with two alternative routines of the Sparse Matrix Vector product (SpMV): the CUSPARSE library and the ELLR-T routine. Although our interest is focused on complex matrices, our implementation has been evaluated on a GPU for two sets of test matrices: complex and real, in single and double precision data. Experimental results show that BCG based on ELLR-T routine achieves the best performance, particularly for the set of complex test matrices. Consequently, this method can be useful as a tool to efficiently solve large linear system of equations (complex and/or nonsymmetric) involved in a broad range of applications.

GMRES implementations and residual smoothing techniques for solving ill-posed linear systems

There are verities of useful Krylov subspace methods to solve nonsymmetric linear system of equations. GMRES is one of the best Krylov solvers with several different variants to solve large sparse linear systems. Any GMRES implementation has some advantages. As the solution of ill-posed problems are important. In this paper, some GMRES variants are discussed and applied to solve these kinds of problems. Residual smoothing techniques are efficient ways to accelerate the convergence speed of some iterative methods like CG variants. At the end of this paper, some residual smoothing techniques are applied for different GMRES methods to test the influence of these techniques on GMRES implementations.

Algorithm 913

The IDR( s ) method that is proposed in Sonneveld and van Gijzen [2008] is a very efficient limited memory method for solving large nonsymmetric systems of linear equations. IDR( s ) is based on the induced dimension reduction theorem, that provides a way to construct subsequent residuals that lie in a sequence of shrinking subspaces. The IDR( s ) algorithm that is given in Sonneveld and van Gijzen [2008] is a direct translation of the theorem into an algorithm. This translation is not unique. This article derives a new IDR( s ) variant, that imposes (one-sided) biorthogonalization conditions on the iteration vectors. The resulting method has lower overhead in vector operations than the original IDR( s ) algorithms. In exact arithmetic, both algorithms give the same residual at every ( s + 1)-st step, but the intermediate residuals and also the numerical properties differ. We show through numerical experiments that the new variant is more stable and more accurate than the original IDR( s ) algorithm, and that it outperforms other state-of-the-art techniques for realistic test problems.

On modified block SSOR iteration methods for linear systems from steady incompressible viscous flow problems

In order to solve the large sparse systems of linear equations arising from numerical solutions of two-dimensional steady incompressible viscous flow problems in primitive variable formulation, we present block SSOR and modified block SSOR iteration methods based on the special structures of the coefficient matrices. In each step of the block SSOR iteration, we employ the block LU factorization to solve the sub-systems of linear equations. We show that the block LU factorization is existent and stable when the coefficient matrices are block diagonally dominant of type-II by columns. Under suitable conditions, we establish convergence theorems for both block SSOR and modified block SSOR iteration methods. In addition, the block SSOR iteration and AF-ADI method are considered as preconditioners for the nonsymmetric systems of linear equations. Numerical experiments show that both block SSOR and modified block SSOR iterations are feasible iterative solvers and they are also effective for preconditioning Krylov subspace methods such as GMRES and BiCGSTAB when used to solve this class of systems of linear equations.

A comparative study of iterative solutions to linear systems arising in quantum mechanics

This study is mainly focused on iterative solutions with simple diagonal preconditioning to two complex-valued nonsymmetric systems of linear equations arising from a computational chemistry model problem proposed by Sherry Li of NERSC. Numerical experiments show the feasibility of iterative methods to some extent when applied to the problems and reveal the competitiveness of our recently proposed Lanczos biconjugate A-orthonormalization methods to other classic and popular iterative methods. By the way, experiment results also indicate that application specific preconditioners may be mandatory and required for accelerating convergence.

GBi-CGSTAB( [formula omitted]): IDR( [formula omitted]) with higher-order stabilization polynomials

IDR( s ) is now recognized as one of the most effective methods, often superior to other Krylov subspace methods, for large nonsymmetric linear systems of equations. In this paper we propose an improvement upon IDR( s ) by incorporating a higher-order stabilization polynomial into IDR( s ). The proposed algorithm, named GBi-CGSTAB( s , L ), shares desirable features with both IDR( s ) and Bi-CGSTAB( L ).

Exploiting BiCGstab($\ell$) Strategies to Induce Dimension Reduction

IDR(s) [P. Sonneveld and M. B. van Gijzen, SIAM J. Sci. Comput., 31 (2008), pp. 1035–1062] and BiCGstab($\ell$) [G. L. G. Sleijpen and D. R. Fokkema, Electron. Trans. Numer. Anal., 1 (1993), pp. 11–32] are two of the most efficient short-recurrence iterative methods for solving large nonsymmetric linear systems of equations. Which of the two is best depends on the specific problem class. In this paper we describe IDRstab, a new method that combines the strengths of IDR(s) and BiCGstab($\ell$). To derive IDRstab we extend the results that we reported on in [G. L. G. Sleijpen, P. Sonneveld, and M. B. van Gijzen, Appl. Numer. Math., (2009), DOI: 10.1016/j.apnum.2009.07.001], where we considered Bi-CGSTAB as an induced dimension reduction (IDR) method. We will analyze the relation between hybrid Bi-CG methods and IDR and introduce the new concept of the Sonneveld subspace as a common framework. Through numerical experiments we will show that IDRstab can outperform both IDR(s) and BiCGstab($\ell$).

Shifted skew-symmetric iteration methods for nonsymmetric linear complementarity problems

We present a shifted skew-symmetric iteration method for solving the nonsymmetric positive definite or positive semidefinite linear complementarity problems. This method is based on the symmetric and skew-symmetric splitting of the system matrix, which has been adopted to establish efficient splitting iteration methods for solving the nonsymmetric systems of linear equations. Global convergence of the method is proved, and the corresponding inexact splitting iteration scheme is established and analyzed in detail. Numerical results show that the new methods are feasible and effective for solving large sparse and nonsymmetric linear complementarity problems.