Unsymmetric Linear Systems Research Articles

GPMR: An Iterative Method for Unsymmetric Partitioned Linear Systems

We introduce an iterative method named Gpmr (general partitioned minimum residual) for solving block unsymmetric linear systems. Gpmr is based on a new process that simultaneously reduces two rectangular matrices to upper Hessenberg form and is closely related to the block-Arnoldi process. Gpmr is tantamount to Block-Gmres with two right-hand sides in which the two approximate solutions are summed at each iteration, but its storage and work per iteration are similar to those of Gmres. We compare the performance of Gpmr with Gmres on linear systems from the SuiteSparse Matrix Collection. In our experiments, Gpmr terminates significantly earlier than Gmres on a residual-based stopping condition with an improvement ranging from around 10% up to 50% in terms of number of iterations.

HILUCSI: Simple, robust, and fast multilevel ILU for large‐scale saddle‐point problems from PDEs

AbstractIncomplete factorization is a widely used preconditioning technique for Krylov subspace methods for solving large‐scale sparse linear systems. Its multilevel variants, such as ILUPACK, are more robust for many symmetric or unsymmetric linear systems than the traditional, single‐level incomplete LU (or ILU) techniques. However, the previous multilevel ILU techniques still lacked robustness and efficiency for some large‐scale saddle‐point problems, which often arise from systems of PDEs. We introduce HILUCSI, or Hierarchical Incomplete LU‐Crout with Scalability‐oriented and Inverse‐based dropping. As a multilevel preconditioner, HILUCSI statically and dynamically permutes individual rows and columns to the next level for deferred factorization. Unlike ILUPACK, HILUCSI applies symmetric preprocessing techniques at the top levels but always uses unsymmetric preprocessing and unsymmetric factorization at the coarser levels. The deferring combined with mixed preprocessing enabled a unified treatment for nearly or partially symmetric systems and simplified the implementation by avoiding mixed and pivots for symmetric indefinite systems. We show that this combination improves robustness for indefinite systems without compromising efficiency. Furthermore, to enable superior efficiency for large‐scale systems with millions or more unknowns, HILUCSI introduces a scalability‐oriented dropping in conjunction with a variant of inverse‐based dropping. We demonstrate the effectiveness of HILUCSI for dozens of benchmark problems, including those from the mixed formulation of the Poisson equation, Stokes equations, and Navier–Stokes equations. We also compare its performance with ILUPACK and the supernodal ILUTP in SuperLU.

The communication-hiding pipelined BiCGstab method for the parallel solution of large unsymmetric linear systems

A High Performance Computing alternative to traditional Krylov subspace methods, pipelined Krylov subspace solvers offer better scalability in the strong scaling limit compared to standard Krylov subspace methods for large and sparse linear systems. The typical synchronization bottleneck is mitigated by overlapping time-consuming global communication phases with local computations in the algorithm. This paper describes a general framework for deriving the pipelined variant of any Krylov subspace algorithm. The proposed framework was implicitly used to derive the pipelined Conjugate Gradient (p-CG) method in Hiding global synchronization latency in the preconditioned Conjugate Gradient algorithm by P. Ghysels and W. Vanroose, Parallel Computing, 40(7):224–238, 2014. The pipelining framework is subsequently illustrated by formulating a pipelined version of the BiCGStab method for the solution of large unsymmetric linear systems on parallel hardware. A residual replacement strategy is proposed to account for the possible loss of attainable accuracy and robustness by the pipelined BiCGStab method. It is shown that the pipelined algorithm improves scalability on distributed memory machines, leading to significant speedups compared to standard preconditioned BiCGStab.

A Class of Generalized Approximate Inverse Solvers for Unsymmetric Linear Systems of Irregular Structure Based on Adaptive Algorithmic Modelling for Solving Complex Computational Problems in Three Space Dimensions

A class of general inverse matrix techniques based on adaptive algorithmic modelling methodologies is derived yielding iterative methods for solving unsymmetric linear systems of irregular structure arising in complex computational problems in three space dimensions. The proposed class of approximate inverse is chosen as the basis to yield systems on which classic and preconditioned iterative methods are explicitly applied. Optimized versions of the proposed approximate inverse are presented using special storage (k-sweep) techniques leading to economical forms of the approximate inverses. Application of the adaptive algorithmic methodologies on a characteristic nonlinear boundary value problem is discussed and numerical results are given.

Improving the convergence behaviour of BiCGSTAB by applying D-norm minimization

In this article, we deal with the iterative methods for solving unsymmetric linear systems, especially BiCGSTAB. The introduced parameter in BiCGSTAB at each iteration is selected to minimize the 2-norm of the residual vector. Here, we suggest another way to select the parameter by the idea of weighting used in Weighted GMRES. By our procedure, more importance is assigned to the larger entry of the residual vector so that faster convergence can be expected. In numerical experiments, it is shown that our procedure is efficient compared with the original BiCGSTAB.

PRECONDITIONING ISSUES IN THE NUMERICAL SOLUTION OF NONLINEAR EQUATIONS AND NONLINEAR LEAST SQUARES

Second order methods for optimization call for the solution of sequences of linear systems. In this survey we will discuss several issues related to the preconditioning of such sequences. Covered topics include both techniques for building updates of factorized preconditioners and quasi-Newton approaches. Sequences of unsymmetric linear systems arising in Newton-Krylov methods will be considered as well as symmetric positive definite sequences arising in the solution of nonlinear least-squares by Truncated Gauss-Newton methods.

A preconditioned conjugate gradient algorithm for GeneRank with application to microarray data mining

The problem of identifying key genes is of fundamental importance in biology and medicine. The GeneRank model explores connectivity data to produce a prioritization of the genes in a microarray experiment that is less susceptible to variation caused by experimental noise than the one based on expression levels alone. The GeneRank algorithm amounts to solving an unsymmetric linear system. However, when the matrix in question is very large, the GeneRank algorithm is inefficient and even can be infeasible. On the other hand, the adjacency matrix is symmetric in the GeneRank model, while the original GeneRank algorithm fails to exploit the symmetric structure of the problem in question. In this paper, we discover that the GeneRank problem can be rewritten as a symmetric positive definite linear system, and propose a preconditioned conjugate gradient algorithm to solve it. Numerical experiments support our theoretical results, and show superiority of the novel algorithm.

A CPU–GPU hybrid approach for the unsymmetric multifrontal method

Multifrontal is an efficient direct method for solving large-scale sparse and unsymmetric linear systems. The method transforms a large sparse matrix factorization process into a sequence of factorizations involving smaller dense frontal matrices. Some of these dense operations can be accelerated by using a graphic processing unit (GPU). We analyze the unsymmetric multifrontal method from both an algorithmic and implementational perspective to see how a GPU, in particular the NVIDIA Tesla C2070, can be used to accelerate the computations. Our main accelerating strategies include (i) performing BLAS on both CPU and GPU, (ii) improving the communication efficiency between the CPU and GPU by using page-locked memory, zero-copy memory, and asynchronous memory copy, and (iii) a modified algorithm that reuses the memory between different GPU tasks and sets thresholds to determine whether certain tasks be performed on the GPU. The proposed acceleration strategies are implemented by modifying UMFPACK, which is an unsymmetric multifrontal linear system solver. Numerical results show that the CPU–GPU hybrid approach can accelerate the unsymmetric multifrontal solver, especially for computationally expensive problems.

A semi-implicit 3-D numerical model using sigam-coordinate for non-hydrostatic pressure free-surface flows

A 3-D numerical formulation is proposed on the horizontal Cartesian, vertical sigma-coordinate grid for modeling non-hydrostatic pressure free-surface flows. The pressure decomposition technique and semi-implicit method are used, with the solution procedure being split into two steps. First, with the implicit parts of non-hydrostatic pressures excluded, the provisional velocity field and free surface are obtained by solving a 2-D Poisson equation. Second, the theory of the differential operator is employed to derive the 3-D Poisson equation for non-hydrostatic pressures, which is solved to obtain the non-hydrostatic pressures and to update the provisional velocity field. When the non-orthogonal sigma-coordinate transformation is introduced, additional terms come into being, resulting in a 15-diagonal, diagonally dominant but unsymmetric linear system in the 3-D Poisson equation for non-hydrostatic pressures. The Biconjugate Gradient Stabilized (BiCGstab) method is used to solve the resulting 3-D unsymmetric linear system instead of the conjugate gradient method, which can only be used for symmetric, positive-definite linear systems. Three test cases are used for validations. The successful simulations of the small-amplitude wave, a supercritical flow over a ramp and a turbulent flow in the open channel indicate that the new model can simulate well non-hydrostatic flows, supercritical flows and turbulent flows.

Algorithm 907

KLU is a software package for solving sparse unsymmetric linear systems of equations that arise in circuit simulation applications. It relies on a permutation to Block Triangular Form (BTF), several methods for finding a fill-reducing ordering (variants of approximate minimum degree and nested dissection), and Gilbert/Peierls’ sparse left-looking LU factorization algorithm to factorize each block. The package is written in C and includes a MATLAB interface. Performance results comparing KLU with SuperLU, Sparse 1.3, and UMFPACK on circuit simulation matrices are presented. KLU is the default sparse direct solver in the Xyce TM circuit simulation package developed by Sandia National Laboratories.

A numerical evaluation of preprocessing and ILU-type preconditioners for the solution of unsymmetric sparse linear systems using iterative methods

Recent advances in multilevel LU factorizations and novel preprocessing techniques have led to an extremely large number of possibilities for preconditioning sparse, unsymmetric linear systems for solving with iterative methods. However, not all combinations work well for all systems, so making the right choices is essential for obtaining an efficient solver. The numerical results for 256 matrices presented in this article give an indication of which approaches are suitable for which matrices (based on different criteria, such as total computation time or fill-in) and of the differences between the methods.

Additive and Multiplicative Two-Level Spectral Preconditioning for General Linear Systems

In this paper we introduce new preconditioning techniques for the solution of general symmetric and unsymmetric linear systems $A x = b$. These approaches borrow some ideas of the multigrid philosophy designed for the solution of linear systems arising from the discretization of elliptic partial differential equations. We attempt to improve the convergence rate of a prescribed preconditioner $M_1$. In a two-grid framework, this preconditioner is viewed as a smoother and the coarse space is spanned by the eigenvectors associated with the smallest eigenvalues of $M_1A$. We derive both additive and multiplicative variants of the resulting iterated two-level preconditioners for unsymmetric linear systems that can also be adapted for Hermitian positive definite problems. We show that these two-level preconditioners shift the smallest eigenvalues to one and tend to better cluster around one those eigenvalues that $M_1$ already succeeded in moving into the neighborhood of one. We illustrate the behavior of our method through extensive numerical experiments on a set of general linear systems. Finally, we show the effectiveness of these approaches on two challenging real applications; the first comes from a nonoverlapping domain decomposition method in semiconductor device modeling, the second from industrial electromagnetism applications.

A frontal solver for the 21st century

SUMMARY In recent years there have been a number of important developments in frontal algorithms for solving the large sparse linear systems of equations that arise fromnite-element problems. We report on the design of a new fully portable and ecient frontal solver for large-scale real and complex unsymmetric linear systems fromnite-element problems that incorporates these developments. The new package oers both aexible reverse communication interface and a simple to use all-in-one interface, which is designed to make the package more accessible to new users. Other key features include automatic element ordering using a state-of-the-art hybrid multilevel spectral algorithm, minimal main memory requirements, the use of high-level BLAS, and facilities to allow the solver to be used as part of a parallel multiple front solver. The performance of the new solver, which is written in Fortran 95, is illustrated using a range of problems from practical applications. The solver is available as package HSL MA42 ELEMENT within the HSL mathematical software library and, for element problems, supersedes the well-known MA42 package. Copyright ? 2006 John Wiley & Sons, Ltd.

A two-level sparse approximate inverse preconditioner for unsymmetric matrices

Sparse approximate inverse (SPAI) preconditioners are effective in accelerating iterative solutions of a large class of unsymmetric linear systems and their inherent parallelism has been widely explored. The effectiveness of SPAI relies on the assumption of the unknown true inverse admitting a sparse approximation. Furthermore, for the usual right SPAI, one must restrict the number of non-zeros in each column to control the overall construction cost and this restriction can reduce the effectiveness of such preconditioners. To extend the applicability of SPAI, this paper proposes to use two-level preconditioning: possible dense columns of the true inverse, skipped by right SPAI (column-wise), will be better approximated by left SPAI (row-wise). Essentially, we approximate the true inverse by sparse matrices via a Gauss-Jordan like decomposition. Numerical experiments on a class of benchmark test matrices show that our new idea of two-level preconditioning can lead to a major enhancement to the standard SPAI method.

Weighted Matchings for Preconditioning Symmetric Indefinite Linear Systems

Maximum weight matchings have become an important tool for solving highly indefinite unsymmetric linear systems, especially in direct solvers. In this study we investigate the benefit of reorderings and scalings based on symmetrized maximum weight matchings as a preprocessing step for incomplete $\mathrm{LDL^T}$ factorizations. The reorderings are constructed such that the matched entries form $1 \times 1$ or $2 \times 2$ diagonal blocks in order to increase the diagonal dominance of the system. During the incomplete factorization only tridiagonal pivoting is used. We report results for this approach and comparisons with other solution methods for a diverse set of symmetric indefinite matrices, ranging from nonlinear elasticity to interior point optimization.

An overview of SuperLU

We give an overview of the algorithms, design philosophy, and implementation techniques in the software SuperLU, for solving sparse unsymmetric linear systems. In particular, we highlight the differences between the sequential SuperLU (including its multithreaded extension) and parallel SuperLU_DIST. These include the numerical pivoting strategy, the ordering strategy for preserving sparsity, the ordering in which the updating tasks are performed, the numerical kernel, and the parallelization strategy. Because of the scalability concern, the parallel code is drastically different from the sequential one. We describe the user interfaces of the libraries, and illustrate how to use the libraries most efficiently depending on some matrix characteristics. Finally, we give some examples of how the solver has been used in large-scale scientific applications, and the performance.

Wavelet based preconditioners for sparse linear systems

A class of efficient preconditioners based on Daubechies family of wavelets for sparse, unsymmetric linear systems that arise in numerical solution of Partial Differential Equations (PDEs) in a wide variety of scientific and engineering disciplines are introduced. Complete and Incomplete Discrete Wavelet Transforms in conjunction with row and column permutations are used in the construction of these preconditioners. With these Wavelet Transform, the transformed matrix is permuted to band forms. The efficiency of our preconditioners with several Krylov subspace methods is illustrated by solving matrices from Harwell Boeing collection and Tim Davis collection. Also matrices resulting in the solution of Regularized Burgers Equation, free convection in porous enclosure are tested. Our results indicate that the preconditioner based on Incomplete Discrete Haar Wavelet Transform is both cheaper to construct and gives good convergence.

A parallel direct solver for large sparse highly unsymmetric linear systems

The need to solve large sparse linear systems of equations efficiently lies at the heart of many applications in computational science and engineering. For very large systems when using direct factorization methods of solution, it can be beneficial and sometimes necessary to use multiple processors, because of increased memory availability as well as reduced factorization time. We report on the development of a new parallel code that is designed to solve linear systems with a highly unsymmetric sparsity structure using a modest number of processors (typically up to about 16). The problem is first subdivided into a number of loosely connected subproblems and a variant of sparse Gaussian elimination is then applied to each of the subproblems in parallel. An interface problem in the variables on the boundaries of the subproblems must also be factorized. We discuss how our software is designed to achieve the goals of portability, ease of use, efficiency, and flexibility, and illustrate its performance on an SGI Origin 2000, a Cray T3E, and a 2-processor Compaq DS20, using problems arising from real applications.

Theoretical and numerical comparisons of GMRES and WZ-GMRES

WZ-GMRES, ‘a simpler GMRES’ proposed by Walker and Zhou, is mathematically equivalent to the generalized minimal residual method (GMRES) for solving large unsymmetric linear systems of equations. In this paper, relationships are established between two bases of an m-dimensional Krylov subspace K m ( A, r 0), and the condition number of the transition matrix between two bases is studied. Some relationships are derived between the condition numbers of the small matrices R G and R WZ resulting from GMRES and WZ-GMRES, respectively. A detailed analysis shows that generally R WZ is worse conditioned than R G , and in particular, R WZ is definitely ill conditioned when the method is near convergence. Furthermore, numerical behavior of WZ-GMRES is analyzed. It turns out that WZ-GMRES is not numerically equivalent to GMRES when the method is near convergence, and WZ-GMRES is numerically less stable than GMRES and can be numerically unstable. Numerical examples confirm the theoretical results.

The Effects of Unsymmetric Matrix Permutations and Scalings in Semiconductor Device and Circuit Simulation

The solution of large sparse unsymmetric linear systems is a critical and challenging component of semiconductor device and circuit simulations. The time for a simulation is often dominated by this part. The sparse solver is expected to balance different, and often conflicting requirements. Reliability, a low memory-footprint, and a short solution time are a few of these demands. Currently, no black-box solver exists that can satisfy all criteria. The linear systems from both simulations can be highly ill-conditioned and are, therefore, quite challenging for direct and iterative methods. In this paper, it is shown that algorithms to place large entries on the diagonal using unsymmetric permutations and scalings greatly enhance the reliability of both direct and preconditioned iterative solvers for unsymmetric linear systems arising in semiconductor device and circuit simulations. The numerical experiments indicate that the overall solution strategy is both reliable and cost effective.