Sparse Nonsymmetric Linear Systems Research Articles

Developing variable s-step CGNE and CGNR algorithms for non-symmetric linear systems

As is well known, CGNE and CGNR methods as powerful techniques are now commonly used for solving sparse non-symmetric linear systems. In this article, we propose the variable s-step versions of CGNE and CGNR methods to solve non-symmetric linear systems. By reducing communication through data locality, the proposed methods reduce the run time. It is known that sudden increases in s during s-step algorithms can cause the algorithms to lose convergence. This issue is exacerbated for CGNE and CGNR methods by normalizing the system. Two s-step algorithms are also introduced for the methods and a step-by-step and gradual increase of s is suggested to maintain convergence. Some features of CGNE and CGNR extend to their s-step algorithms. Numerical results demonstrate the effectiveness of the methods in maintaining convergence with increasing s.

Fast parallel-in-time quasi-boundary value methods for backward heat conduction problems

In this paper we propose two new quasi-boundary value methods for regularizing the ill-posed backward heat conduction problems. With a standard finite difference discretization in space and time, the obtained all-at-once nonsymmetric sparse linear systems have the desired block ω-circulant structure, which can be utilized to design an efficient parallel-in-time (PinT) direct solver that is built upon an explicit FFT-based diagonalization of the time discretization matrix. Convergence analysis is presented to justify the optimal choice of the regularization parameter under suitable assumptions. Numerical examples are reported to validate our analysis and illustrate the superior computational efficiency of our proposed PinT methods.

Balanced incomplete factorization preconditioner with pivoting

In this work we study pivoting strategies for the preconditioner presented in Bru (SIAM J Sci Comput 30(5):2302–2318, 2008) which computes the LU factorization of a matrix A. This preconditioner is based on the Inverse Sherman Morrison (ISM) decomposition [Preconditioning sparse nonsymmetric linear systems with the Sherman–Morrison formula. Bru (SIAM J Sci Comput 25(2):701–715, 2003), that using recursion formulas derived from the Sherman-Morrison formula, obtains the direct and inverse LU factors of a matrix. We present a modification of the ISM decomposition that allows for pivoting, and so the computation of preconditioners for any nonsingular matrix. While the ISM algorithm at a given step computes only a new pair of vectors, the new pivoting algorithm in the k-th step also modifies all the remaining vectors from k+1 to n. Thus, it can be seen as a right looking version of the ISM decomposition. The results of numerical experiments with ill-conditioned and highly indefinite matrices arising from different applications show the robustness of the new algorithm, since it is able to solve problems that are not possible to solve otherwise.

GPBi‐CGstab(L): A Lanczos‐type product method unifying Bi‐CGstab(L) and GPBi‐CG

SummaryLanczos‐type product methods (LTPMs), in which the residuals are defined by the product of stabilizing polynomials and the Bi‐CG residuals, are effective iterative solvers for large sparse nonsymmetric linear systems. Bi‐CGstab(L) and GPBi‐CG are popular LTPMs and can be viewed as two different generalizations of other typical methods, such as CGS, Bi‐CGSTAB, and Bi‐CGStab2. Bi‐CGstab(L) uses stabilizing polynomials of degree L, while GPBi‐CG uses polynomials given by a three‐term recurrence (or equivalently, a coupled two‐term recurrence) modeled after the Lanczos residual polynomials. Therefore, Bi‐CGstab(L) and GPBi‐CG have different aspects of generalization as a framework of LTPMs. In the present paper, we propose novel stabilizing polynomials, which combine the above two types of polynomials. The resulting method is referred to as GPBi‐CGstab(L). Numerical experiments demonstrate that our presented method is more effective than conventional LTPMs.

A transformation approach that makes SPAI, PSAI and RSAI procedures efficient for large double irregular nonsymmetric sparse linear systems

A sparse matrix is called double irregular sparse if it has at least one relatively dense column and row, and it is double regular sparse if all the columns and rows of it are sparse. The sparse approximate inverse preconditioning procedures SPAI, PSAI( t o l ) and RSAI( t o l ) are costly and even impractical to construct preconditioners for a large sparse nonsymmetric linear system with the coefficient matrix being double irregular sparse, but they are efficient for double regular sparse problems. Double irregular sparse linear systems have a wide range of applications, and 24.4% of the nonsymmetric matrices in the Florida University collection are double irregular sparse. For this class of problems, we propose a transformation approach, which consists of four steps: (i) transform a given double irregular sparse problem into a small number of double regular sparse ones with the same coefficient matrix A ˆ , (ii) use SPAI, PSAI( t o l ) and RSAI( t o l ) to construct sparse approximate inverses M of A ˆ , (iii) solve the preconditioned double regular sparse linear systems by Krylov solvers, and (iv) recover an approximate solution of the original problem with a prescribed accuracy from those of the double regular sparse ones. A number of theoretical and practical issues are considered on the transformation approach. Numerical experiments on a number of real-world problems confirm the very sharp superiority of the transformation approach to the standard approach that preconditions the original double irregular sparse problem by SPAI, PSAI( t o l ) or RSAI( t o l ) and solves the resulting preconditioned system by Krylov solvers.

Weak form of peridynamics for nonlocal essential and natural boundary conditions

This study presents the weak form of peridynamic (PD) governing equations which permit the direct imposition of nonlocal essential and natural boundary conditions. It also presents a variational approach to derive the PD form of first- and second-order derivatives of a field variable at a point which is not symmetrically located in its domain of interaction. This capability enables the nonlocal PD representation of the internal force vector and the stress components without any calibration procedure. Furthermore, it removes the concern of truncated domain of interaction for a point near the surface. Thus, the solution is free of nonlocal boundary forces and surface effects. The numerical solution of the resulting equations can be achieved by considering an unstructured nonuniform discretization. The implicit solution to the discrete form of the equations is achieved by employing BiConjugate Gradient Stabilized (BICGSTAB) method which is an iterative technique for solving sparse non-symmetric linear systems. The explicit analysis is performed by constructing a global diagonal mass matrix, and using a hybrid implicit/explicit time integration scheme. The accuracy of this approach is demonstrated by considering an elastic isotropic plate with or without a cutout subjected to a combination of different types of boundary conditions under plane stress conditions.

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.

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.

Numerical Solution of the Steady Convection-diffusion Equation with Dominant Convection

Steady convection-diffusion equation in 2-D domain is considered. Central finite-difference approximation has been taken to obtain a large sparse nonsymmetric linear system with positive real matrix. New class of product triangular skew-symmetric iterative methods for solution of such system is presented and considered. Using this method as preconditioner for GMRES and BiCG has been made. Results of numerical experiments for two-dimensional convection-diffusion equation for different big Peclet numbers and velocity coefficients have been presented.

A parallel sparse algorithm targeting arterial fluid mechanics computations

Iterative solution of large sparse nonsymmetric linear equation systems is one of the numerical challenges in arterial fluid---structure interaction computations. This is because the fluid mechanics parts of the fluid + structure block of the equation system that needs to be solved at every nonlinear iteration of each time step corresponds to incompressible flow, the computational domains include slender parts, and accurate wall shear stress calculations require boundary layer mesh refinement near the arterial walls. We propose a hybrid parallel sparse algorithm, domain-decomposing parallel solver (DDPS), to address this challenge. As the test case, we use a fluid mechanics equation system generated by starting with an arterial shape and flow field coming from an FSI computation and performing two time steps of fluid mechanics computation with a prescribed arterial shape change, also coming from the FSI computation. We show how the DDPS algorithm performs in solving the equation system and demonstrate the scalability of the algorithm.

A Simplified and Flexible Variant of GCROT for Solving Nonsymmetric Linear Systems

There is a need for flexible iterative solvers that can solve large-scale ($>10^6$ unknowns) nonsymmetric sparse linear systems to a small tolerance. Among flexible solvers, flexible GMRES (FGMRES) is attractive because it minimizes the residual norm over a particular subspace. In practice, FGMRES is often restarted periodically to keep memory and work requirements reasonable; however, like restarted GMRES, restarted FGMRES can suffer from stagnation. This has led us to develop a flexible variant of the Krylov linear solver GCROT (generalized conjugate residual with inner orthogonalization and outer truncation). Unlike the original GCROT algorithm, the proposed GCROT variant uses a simplified truncation strategy similar to loose GMRES (LGMRES). This modification is motivated by numerical experiments that suggest the specific subspace retained in the outer iteration of GCROT is less important than its size. The flexible GCROT variant appears to be well suited for advection-dominated problems. In particular, when applied to an adjoint problem from computational aerodynamics, the proposed GCROT variant is robust and efficient compared with several popular truncated Krylov subspace methods. Finally, a flexible version of LGMRES is easily constructed by recognizing algorithmic similarities to GCROT.

COMPARISONS OF THE PARALLEL PRECONDITIONERS FOR LARGE NONSYMMETRIC SPARSE LINEAR SYSTEMS ON A PARALLEL COMPUTER

In this paper we compare various parallel preconditioners for solving large sparse nonsymmetric linear systems. They are Block Jacobi, Point-SSOR, ILU(0) in the wavefront order, ILU(0) in the multi-color order, SPAI(SParse Approximate Inverse), and Multi-Color Block SOR. The Block Jacobi and Point-SSOR are well-known, and ILU(0) is one of the most popular preconditioners, but it is inherently serial. ILU(0) in the wavefront order maximizes the parallelism, and ILU(0) in the multi-color order achieves the parallelism of order (N), where N is the order of the matrix. The SPAI tries to capture the approximate inverse in sparse form, which, then, is expected to be a scalable preconditioner. Finally, we implemented the Multi-Color Block SOR preconditioner combined with direct sparse matrix solver. For the Laplacian matrix the SOR method is known to have a non-deteriorating rate of convergence when used with Multi-Color ordering. Since most of the time is spent on the diagonal inversion, which is done on each processor, we expect it to be a good scalable preconditioner. Finally, due to the blocking effect, it will be effective for ill-conditioned problems. Experiments were conducted for the Finite Difference discretizations of two problems with various meshsizes varying up to 1024×1024, and for an ill-conditioned matrix from the shell problem from the Harwell–Boeing collection. CRAY-T3E with 128 nodes was used. MPI library was used for interprocess communications. The results show that Multi-Color Block SOR and ILU(0) with Multi-Color ordering give the best performances for the finite difference matrices and for the shell problem only the Multi-Color Block SOR and Block Jacobi converges. Based on this we recommend that the Multi-Color Block SOR is the most robust preconditioner out of the preconditioners considered.

Block {ω}-circulant preconditioners¶for the systems of differential equations

The numerical solution of large and sparse nonsymmetric linear systems of algebraic equations is usually the most time consuming part of time-step integrators for differential equations based on implicit formulas. Preconditioned Krylov subspace methods using Strang block circulant preconditioners have been employed to solve such linear systems. However, it has been observed that these block circulant preconditioners can be very ill-conditioned or singular even when the underlying nonpreconditioned matrix is well-conditioned. In this paper we propose the more general class of the block { ω }-circulant preconditioners. For the underlying problems, ω can be chosen so that the condition number of these preconditioners is much smaller than that of the Strang block circulant preconditioner (which belongs to the same class with ω =1) and the related iterations can converge very quickly.

An efficient method for constructing an ILU preconditioner for solving large sparse nonsymmetric linear systems by the GMRES method

The main idea of this paper is in determination of the pattern of nonzero elements of the LU factors of a given matrix A. The idea is based on taking the powers of the Boolean matrix derived from A. This powers of a Boolean matrix strategy (PBS) is an efficient, effective, and inexpensive approach. Construction of an ILU preconditioner using PBS is described and used in solving large nonsymmetric sparse linear systems. Effectiveness of the proposed ILU preconditioner in solving large nonsymmetric sparse linear systems by the GMRES method is also shown. Numerical experiments are performed which show that it is possible to considerably reduce the number of GMRES iterations when the ILU preconditioner constructed here is used. In numerical examples, the influence of k, the dimension of the Krylov subspace, on the performance of the GMRES method using an ILU preconditioner is tested. For all the tests carried out, the best value for k is found to be 10.

Restricted Additive Schwarz Preconditioners with Harmonic Overlap for Symmetric Positive Definite Linear Systems

A restricted additive Schwarz (RAS) preconditioning technique was introduced recently for solving general nonsymmetric sparse linear systems. In this paper, we provide one-level and two-level extensions of RAS for symmetric positive definite problems using the so-called harmonic overlaps (RASHO). Both RAS and RASHO outperform their counterparts of the classical additive Schwarz variants (AS). The design of RASHO is based on a much deeper understanding of the behavior of Schwarz-type methods in overlapping subregions and in the construction of the overlap. In RASHO, the overlap is obtained by extending the nonoverlapping subdomains only in the directions that do not cut the boundaries of other subdomains, and all functions are made harmonic in the overlapping regions. As a result, the subdomain problems in RASHO are smaller than those of AS, and the communication cost is also smaller when implemented on distributed memory computers, since the right-hand sides of discrete harmonic systems are always zero and therefore do not need to be communicated. We also show numerically that RASHO-preconditioned CG takes fewer iterations than the corresponding AS-preconditioned CG. A nearly optimal theory is included for the convergence of RASHO-preconditioned CG for solving elliptic problems discretized with a finite element method.

The effect of orderings on sparse approximate inverse preconditioners for non-symmetric problems

We experimentally study how reordering techniques affect the rate of convergence of preconditioned Krylov subspace methods for non-symmetric sparse linear systems, where the preconditioner is a sparse approximate inverse. In addition, we show how the reordering reduces the number of entries in the approximate inverse and thus, the amount of storage and computation required for a given accuracy. These properties are illustrated with several numerical experiments taken from the discretization of PDEs by a finite element method and from a standard matrix collection.

SCALABILITY ANALYSIS OF PARALLEL GMRES IMPLEMENTATIONS

Applications involving large sparse nonsymmetric linear systems encourage parallel implementations of robust iterative solution methods, such as GMRES(k). Two parallel versions of GMRES(k) based on different data distributions and using Householder reflections in the orthogonalization phase are analyzed with respect to scalability (their ability to maintain fixed efficiency with an increase in problem size and number of processors). A theoretical algorithm-machine model for scalability of GMRES(k) with fixed k is derived and validated by experiments on three parallel computers, each with different machine characteristics. The analysis for an adaptive version of GMRES(k), in which the restart value k is adapted to the problem, is also presented and scalability results for this case are briefly discussed.

LU-decomposition and numerical structure for solving large sparse nonsymmetric linear systems

In this work, the solution of a large sparse linear system of equations with an arbitrary sparsity pattern is obtained by using LU-decomposition method as well as numerical structure approach. The LU-decomposition method is based on Doolittle's method while the numerical structure approach is based on Cramer's rule. The numerical structure approach produces direct solution without facing fill-in problems as encountered in LU-decomposition. In order to reduce the ‘fill-ins’ in the decomposition, the powers of a Boolean matrix, obtained from the coefficient matrix A are taken so that the ‘fill-ins’ in the structure of A can be known in advance. The position of fill-ins in A are thus determined in the best choice manner, that is, it is very effective and memory-wise cheap. We also outline a method by using numerical structure with reduced computation efforts. Finally, experiments are performed on eight examples to compare the efficiency of the proposed methods. The results obtained are reported in a table. It is found that the LU-decomposition method is much better than numerical structure. The usefulness of numerical structure approach is also discussed.

Variational principle for two-dimensional high Hartmann number flow

A variational principle for the stream function-vorticity form of the linearised two-dimensional MHD equations which incorporates the no-slip condition as a natural boundary condition has been used to solve a range of problems. Finite element discretisations leading to sparse non-symmetric linear systems have been solved by a corrected version of a standard non-symmetric frontal solver. The numerical solutions obtained for flows over an unbounded cylinder at the larger Hartmann numbers agree with recently obtained theoretical asymptotic estimates. Other calculations for bounded flows in channels with or without obstructions are also in agreement with previous results obtained by matched asymptotic expansion methods.

The global Hessenberg and CMRH methods for linear systems with multiple right-hand sides

In this paper, we introduce two new methods for solving large sparse nonsymmetric linear systems with several right-hand sides. These methods are the global Hessenberg and global CMRH methods. Using the global Hessenberg process, these methods are less expensive than the global FOM and global GMRES methods [9]. Theoretical results about the new methods are given, and experimental results that show good performances of these new methods are presented.