Nonsymmetric Linear Systems Research Articles

Preconditioned conjugate gradient methods for the solution of Love’s integral equation with very small parameter

In this paper, we propose efficient numerical methods for the solution of the following Love’s integral equation f ( x ) + 1 π ∫ − 1 1 c ( x − y ) 2 + c 2 f ( y ) d y = 1 , x ∈ [ − 1 , 1 ] , where c > 0 is a very small parameter. We introduce a new unknown function h ( x ) = f ( x ) − 0 . 5 as in Lin et al. (2013), and then apply a composite Gauss–Legendre quadrature to the resulting integral equation as in Pastore (2011). The coefficient matrix of corresponding linear system is a nonsymmetric block matrix with Toeplitz blocks. We transform the nonsymmetric linear system into a symmetric linear system and introduce a preconditioner which is a block matrix with circulant blocks. Spectral properties of relevant matrices are analyzed and numerical results are presented to illustrate the efficiency of the proposed methods.

Nonoverlapping Domain Decomposition Preconditioners for Discontinuous Galerkin Approximations of Hamilton–Jacobi–Bellman Equations

We analyse a class of nonoverlapping domain decomposition preconditioners for nonsymmetric linear systems arising from discontinuous Galerkin finite element approximation of fully nonlinear Hamilton--Jacobi--Bellman (HJB) partial differential equations. These nonsymmetric linear systems are uniformly bounded and coercive with respect to a related symmetric bilinear form, that is associated to a matrix $\mathbf{A}$. In this work, we construct a nonoverlapping domain decomposition preconditioner $\mathbf{P}$, that is based on $\mathbf{A}$, and we then show that the effectiveness of the preconditioner for solving the} nonsymmetric problems can be studied in terms of the condition number $\kappa(\mathbf{P}^{-1}\mathbf{A})$. In particular, we establish the bound $\kappa(\mathbf{P}^{-1}\mathbf{A}) \lesssim 1+ p^6 H^3 /q^3 h^3$, where $H$ and $h$ are respectively the coarse and fine mesh sizes, and $q$ and $p$ are respectively the coarse and fine mesh polynomial degrees. This represents the first such result for this class of methods that explicitly accounts for the dependence of the condition number on $q$; our analysis is founded upon an original optimal order approximation result between fine and coarse discontinuous finite element spaces. Numerical experiments demonstrate the sharpness of this bound. Although the preconditioners are not robust with respect to the polynomial degree, our bounds quantify the effect of the coarse and fine space polynomial degrees. Furthermore, we show computationally that these methods are effective in practical applications to nonsymmetric, fully nonlinear HJB equations under $h$-refinement for moderate polynomial degrees.

Convergence analysis of the MCGNR algorithm for the least squares solution group of discrete-time periodic coupled matrix equations

The conjugate gradient normal equations residual minimizing (CGNR) algorithm is a popular tool for solving large nonsymmetric linear systems. In this study, we propose the matrix form of the CGNR (MCGNR) algorithm to find the least squares solution group of the discrete-time periodic coupled matrix equations [Formula: see text] with periodic matrices. We prove that the MCGNR algorithm converges in a finite number of steps in the absence of round-off errors, that is, this algorithm has the finite termination property. Also we show that the norms of the residual matrices of the MCGNR algorithm decrease monotonically during its iteration, that is, this algorithm has the residual reducing property. To show the efficiency of the MCGNR algorithm, two numerical examples are presented.

Using the VBARMS method in parallel computing

We introduce a new graph compression algorithm for block construction, which extends the method by Ashcraft and requires one simple to use parameter.We describe an improved block factorization computation in the VBARMS method.We conduct a parallel performance study of the VBARMS method using parallel graph partitioning.We solve turbulent three-dimensional Navier-Stokes equations on large realistic unstructured meshes. The paper describes an improved parallel MPI-based implementation of VBARMS, a variable block variant of the pARMS preconditioner proposed by Li etźal. 200314 for solving general nonsymmetric linear systems. The parallel VBARMS solver can detect automatically exact or approximate dense structures in the linear system, and exploits this information to achieve improved reliability and increased throughput during the factorization. A novel graph compression algorithm is discussed that finds these approximate dense blocks structures and requires only one simple to use parameter. A complete study of the numerical and parallel performance of parallel VBARMS is presented for the analysis of large turbulent Navier-Stokes equations on a suite of three-dimensional test cases.

On the convergence of Q-OR and Q-MR Krylov methods for solving nonsymmetric linear systems

This paper addresses the convergence behavior of Krylov methods for nonsymmetric linear systems which can be classified as quasi-orthogonal (Q-OR) or quasi-minimum residual (Q-MR) methods. It explores, more precisely, whether the influence of eigenvalues is the same when using non-orthonormal bases as it is for the FOM and GMRES methods. It presents parametrizations of the classes of matrices with a given spectrum and right-hand sides generating prescribed Q-OR/Q-MR (quasi) residual norms and discusses non-admissible residual norm sequences. It also gives closed-form expressions of the Q-OR/Q-MR (quasi) residual norms as functions of the eigenvalues and eigenvectors of the matrix of the linear system.

A Note on the Block and Seed BiCGSTAB Algorithms for Nonsymmetric Multiple Linear Systems

In this paper we ϐirst give a new derivation of the Bl-BiCGSTAB algorithm for solving block linear systems. This derivation is based only on linear algebra and didn't needs matrix orthogonal polynomials. We also propose a new method for solving linear systems with the same coefϐicient matrix and different right-hand sides. The proposed method is based on the BiCGSTAB algorithm and is especially efϐicient when the right-hand sides vector-columns are closely related. Some numerical examples are reported to illustrate the effectiveness of the proposed methods.

Convergence analysis of the MCGNR algorithm for the least squares solution group of discrete-time periodic coupled matrix equations

The conjugate gradient normal equations residual minimizing (CGNR) algorithm is a popular tool for solving large nonsymmetric linear systems. In this study, we propose the matrix form of the CGNR (MCGNR) algorithm to find the least squares solution group of the discrete-time periodic coupled matrix equations [Formula: see text] with periodic matrices. We prove that the MCGNR algorithm converges in a finite number of steps in the absence of round-off errors, that is, this algorithm has the finite termination property. Also we show that the norms of the residual matrices of the MCGNR algorithm decrease monotonically during its iteration, that is, this algorithm has the residual reducing property. To show the efficiency of the MCGNR algorithm, two numerical examples are presented.

Matrix-equation-based strategies for convection–diffusion equations

We are interested in the numerical solution of nonsymmetric linear systems arising from the discretization of convection–diffusion partial differential equations with separable coefficients, dominant convection and rectangular or parallelepipedal domain. Preconditioners based on the matrix equation formulation of the problem are proposed, which naturally approximate the original discretized problem. For the considered setting, we show that the explicit solution of the matrix equation can effectively replace the linear system solution. Numerical experiments with data stemming from two and three dimensional problems are reported, illustrating the potential of the proposed methodology.

The matrix iterative methods for solving a class of generalized coupled Sylvester-conjugate linear matrix equations

The conjugate gradients-squared (CGS) method (Sonneveld, 1989) has been considered as an efficient variant of the bi-conjugate gradient (BCG) method. In Vorst (1992), a more smoothly converging variant of the BCG method which keeps the attractive convergence rate of the CGS method was investigated for the solution of certain classes of nonsymmetric linear systems, so-called bi-conjugate gradient stabilized (Bi-CGSTAB) method. In this paper, we will combine these interesting methods for solving the generalized coupled Sylvester-conjugate matrix equations A1XB1+C1Y‾D1=E,A2X‾B2+C2YD2=F after performing suitable transformation by the properties of Kronecker product and vec operator. Some numerical experiments demonstrate that the introduced iterative methods are more efficient than the existing methods.

Efficient Preconditioners for a Shock Capturing Space-Time Discontinuous Galerkin Method for Systems of Conservation Laws

AbstractAn entropy stable fully discrete shock capturing space-time Discontinuous Galerkin (DG) method was proposed in a recent paper to approximate hyperbolic systems of conservation laws. This numerical scheme involves the solution of a very large nonlinear system of algebraic equations, by a Newton-Krylov method, at every time step. In this paper, we design efficient preconditioners for the large, non-symmetric linear system, that needs to be solved at every Newton step. Two sets of preconditioners, one of the block Jacobi and another of the block Gauss-Seidel type are designed. Fourier analysis of the preconditioners reveals their robustness and a large number of numerical experiments are presented to illustrate the gain in efficiency that results from preconditioning. The resulting method is employed to compute approximate solutions of the compressible Euler equations, even for very high CFL numbers.

Robust and Efficient Auxiliary Density Perturbation Theory Calculations.

A new iterative solver for the recently developed time-dependent auxiliary density perturbation theory is presented. It is based on the Eirola-Nevanlinna algorithm for large nonsymmetric linear equation systems. The new methodology is validated by static and dynamic polarizability calculations of small molecules. Comparison between the analytic and iterative solutions of the response equation system shows excellent agreement for the calculated static and dynamic polarizabilities. The new iterative solver reduces the formal scaling from [Symbol: see text](N(4)) to [Symbol: see text](N(3)). Furthermore, the observed computational scaling for linear alkane chains is N(1.6). This subquadratic behavior is possible in systems with a few hundred atoms because of the very small prefactors of the [Symbol: see text](N(3)) and [Symbol: see text](N(2)) steps remaining in the iterative solver. To demonstrate the potential of this new methodology, static polarizabilities of giant fullerenes up to C960, with more than 14,000 basis functions, are calculated.

The block LSMR method: a novel efficient algorithm for solving non-symmetric linear systems with multiple right-hand sides

It is well known that if the coefficient matrix in a linear system is large and sparse or sometimes not readily available, then iterative solvers may become the only choice. The block solvers are an attractive class of iterative solvers for solving linear systems with multiple right-hand sides. In general, the block solvers are more suitable for dense systems with preconditioner. In this paper, we present a novel block LSMR (least squares minimal residual) algorithm for solving non-symmetric linear systems with multiple right-hand sides. This algorithm is based on the block bidiagonalization and LSMR algorithm and derived by minimizing the 2-norm of each column of normal equation. Then, we give some properties of the new algorithm. In addition, the convergence of the stated algorithm is studied. In practice, we also observe that the Frobenius norm of residual matrix decreases monotonically. Finally, some numerical examples are presented to show the efficiency of the new method in comparison with the traditional LSMR method.

Developing BiCOR and CORS methods for coupled Sylvester-transpose and periodic Sylvester matrix equations

Recently the biconjugate A-orthogonal residual (BiCOR) and the conjugate A-orthogonal residual squared (CORS) methods have been introduced for solving nonsymmetric linear systems Ax=b. This study directly develops the BiCOR and CORS methods to obtain matrix iterative methods for solving the coupled Sylvester-transpose matrix equations∑k=1l(A1,kXB1,k+C1,kXTD1,k+E1,kYF1,k)=M1,∑k=1l(A2,kXB2,k+C2,kXTD2,k+E2,kYF2,k)=M2,and the coupled periodic Sylvester matrix equationsA1,jXjB1,j+C1,jXj+1D1,j+E1,jYjF1,j=M1,j,A2,jXjB2,j+C2,jXj+1D2,j+E2,jYjF2,j=M2,j,forj=1,2,…,μ.Numerical examples are given at the end of this paper to compare the accuracy and efficiency of the matrix iterative methods with other methods in the literature.

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.

Special iterative methods for solution of the steady Convection-Diffusion-Reaction equation with dominant convection

Iterative methods based on skew-symmetric splitting of initial matrix, arising from central finite-difference approximation of steady convection-diffusion-reaction (CDR) equation in 2-D domain are considered. The property of obtained large sparse nonsymmetric linear system Au = f is investigated. A new class of triangular and product triangular skew-symmetric iterative methods is presented. Sufficient conditions of convergence for iterative methods of solution CDR with variable coefficient of reaction and dominant convection are obtained. The results of numerical experiments for the solution of a two-dimensional CDR equation are presented. The uniform grid, central differences for the first derivatives, natural ordering of points, Peclet numbers Pe = 103, 104, 105, variable coefficient of reaction and different velocity coefficients have been used.

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.

Global GPBiCG method for complex non-Hermitian linear systems with multiple right-hand sides

The global BiCGSTAB (Gl-BiCGSTAB) method has been found to perform very well in many practical cases. However, like the BiCGSTAB method, it often converges slowly or stagnates for some non-symmetric linear systems with complex spectrum. In this paper, we propose a new global generalized product-type method based on the global BiCG (Gl-BCG) method to try and remedy these problems. To enhance stability of convergence and robustness, we also consider a variant of the new method based on minimization of the associate residual of F-norm and ascending order recurrence in place of reconstruction of the algorithm. Finally, numerical experiments are given to show that our methods can be more effective than the original methods.

The optimised Schwarz method and the two-Lagrange multiplier method for heterogeneous problems in general domains with two general subdomains

The optimised Schwarz method and the related two-Lagrange multiplier (2LM) method are nonoverlapping domain decomposition methods which can be used to numerically solve boundary value problems. Local Robin problems are solved on each subdomain in parallel to approximate the global solution, where a careful choice of Robin parameters leads to faster convergence. The 2LM method involves the solution of a nonsymmetric linear system that is usually solved with a Krylov subspace method such as GMRES. The speed of convergence of GMRES can be estimated using a conformal map from the exterior of the field of values of the system matrix to the interior of the unit disc. In this article we consider an elliptic PDE problem with a jump in diffusion coefficients across the interface between the subdomains. We approximate the field of values of the 2LM system matrix by a rectangle, R, in ? and provide optimised Robin parameters that ensure R is "well conditioned" in the sense that GMRES converges quickly. We derive convergence estimates for GMRES and consider the behaviour asymptotically as the mesh size h becomes small and the jump in coefficients becomes large. We observe, for our choice of Robin parameters, that increasing the jump in coefficients increases the convergence rate of GMRES. Numerical experiments are performed to verify the theoretical results.

A nonhydrostatic, isopycnal-coordinate ocean model for internal waves

We present a nonhydrostatic ocean model with an isopycnal (density-following) vertical coordinate system. The primary motivation for the model is the proper treatment of nonhydrostatic dispersion and the formation of nonlinear internal solitary waves. The nonhydrostatic, isopycnal-coordinate formulation may be preferable to nonhydrostatic formulations in z- and σ-coordinates because it improves computational efficiency by reducing the number of vertical grid points and eliminates spurious diapycnal mixing and solitary-wave amplitude loss due to numerical diffusion of scalars. The model equations invoke a mild isopycnal-slope approximation to remove small metric terms associated with diffusion and nonhydrostatic pressure from the momentum equations and to reduce the pressure Poisson equation to a symmetric linear system. Avoiding this approximation requires a costlier inversion of a non-symmetric linear system. We demonstrate that the model is capable of simulating nonlinear internal solitary waves for simplified and physically-realistic ocean-scale problems with a reduced number of layers.

IDR([formula omitted]) for solving shifted nonsymmetric linear systems

The IDR(s) method by Sonneveld and van Gijzen (2008) has recently received tremendous attention since it is effective for solving nonsymmetric linear systems. In this paper, we generalize this method to solve shifted nonsymmetric linear systems. When solving this kind of problem by existing shifted Krylov subspace methods, we know one just needs to generate one basis of the Krylov subspaces due to the shift-invariance property of Krylov subspaces. Thus the computation cost required by the basis generation of all shifted linear systems, in terms of matrix–vector products, can be reduced. For the IDR(s) method, we find that there also exists a shift-invariance property of the Sonneveld subspaces. This inspires us to develop a shifted version of the IDR(s) method for solving the shifted linear systems.