Krylov Subspace Iterative Methods Research Articles

On the implementation of the Dirichlet-to-Neumann radiation condition for iterative solution of the Helmholtz equation

The Helmholtz equation posed on an unbounded domain with the Sommerfeld condition prescribed at infinity is considered. The unbounded domain is eliminated by imposing a Dirichlet-to-Neumann (DtN) map or a modified DtN map on a truncating surface and the resulting bounded domain problem is modeled using the finite element method. The resulting system of linear equations is then solved using a Krylov subspace iterative method. New, efficient algorithms to compute matrix-vector products that are based on the structure of the DtN and the modified DtN map are presented. Connections between the DtN map and the discrete Fourier transform in two dimensions and discrete spherical transform in three dimensions are established, and are utilized to develop fast implementations of matrix-vector product algorithms. Also, an SSOR-type preconditioner that is based on a local radiation condition is considered for the modified DtN formulation. An efficient implementation is proposed by extending Eisenstat's trick for the standard SSOR preconditioner. Finally, numerical examples which illustrate the efficacy of the proposed algorithms are presented.

Block algorithms for quark propagator calculation

Computing quark propagators in lattice QCD is equivalent to solving large, sparse linear systems with multiple right-hand sides. Block algorithms attempt to accelerate the convergence of iterative Krylov-subspace methods by solving the multiple systems simultaneously. This paper compares a block generalisation of the quasi-minimal residual method (QMR), Block Conjugate Gradient on the normal equation, Block Lanczos and ( γ 5-symmetric). Block BiConjugate Gradient.

Parallel preconditioning based on h-hierarchical finite elements with application to acoustics

In this paper, we consider the hierarchical basis preconditioning approach which is closely related to the h-version of the hierarchical finite element method. In general, finite element formulations that employ hierarchical shape functions yield better conditioned matrix problems than those based on the Lagrange nodal basis. These matrix problems are also better-suited to a faster rate of convergence with Krylov-subspace iterative methods. We consider a purely algebraic approach to describe projections between the nodal and the h-hierarchical bases functions, which are then used to construct the preconditioning operator. Implementation details of the preconditioner are provided for solving finite element problems on unstructured grids. We employ the preconditioning approach in the iterative solution of linear systems that arise from finite element discretization of the exterior Helmholtz problem. Numerical results are presented to examine convergence rates for practical discretizations in acoustics, and to illustrate the computational performance of the preconditioning algorithm on serial and data-parallel computers.

Wavelet-Galerkin solution of boundary value problems

In this paper we review the application of wavelets to the solution of partial differential equations. We consider in detail both the single scale and the multiscale Wavelet Galerkin method. The theory of wavelets is described here using the language and mathematics of signal processing. We show a method of adapting wavelets to an interval using an extrapolation technique called Wavelet Extrapolation. Wavelets on an interval allow boundary conditions to be enforced in partial differential equations and image boundary problems to be overcome in image processing. Finally, we discuss the fast inversion of matrices arising from differential operators by preconditioning the multiscale wavelet matrix. Wavelet preconditioning is shown to limit the growth of the matrix condition number, such that Krylov subspace iteration methods can accomplish fast inversion of large matrices.

Wavelet sparse approximate inverse preconditioners

There is an increasing interest in using sparse approximate inverses as preconditioners for Krylov subspace iterative methods. Recent studies of Grote and Huckle and Chow and Saad also show that sparse approximate inverse preconditioner can be effective for a variety of matrices, e.g. Harwell-Boeing collections. Nonetheless a drawback is that it requires rapid decay of the inverse entries so that sparse approximate inverse is possible. However, for the class of matrices that, come from elliptic PDE problems, this assumption may not necessarily hold. Our main idea is to look for a basis, other than the standard one, such that a sparse representation of the inverse is feasible. A crucial observation is that the kind of matrices we are interested in typically have a piecewise smooth inverse. We exploit this fact, by applying wavelet techniques to construct a better sparse approximate inverse in the wavelet basis. We shall justify theoretically and numerically that our approach is effective for matrices with smooth inverse. We emphasize that in this paper we have only presented the idea of wavelet approximate inverses and demonstrated its potential but have not yet developed a highly refined and efficient algorithm.

Iterative solution of multiple radiation and scattering problems in structural acoustics using a block quasi-minimal residual algorithm

Finite-element discretizations of time-harmonic acoustic wave problems in exterior domains result in large sparse systems of linear equations with complex symmetric coefficient matrices. In many situations, these matrix problems need to be solved repeatedly for different right-hand sides, but with the same coefficient matrix. Recently, Freund and Malhotra have proposed a block quasi-minimal residual (BL-QMR) algorithm [13] for the iterative solution of non-Hermitian linear systems with multiple right-hand sides. The BL-QMR algorithm is a block Krylov-subspace iterative method that incorporates deflation to delete linearly and almost linearly dependent vectors in the underlying block Krylov sequences. In this paper, we describe a J -symmetric variant of the BL-QMR algorithm that introduces important simplifications for the case when the coefficient matrix is symmetric with respect to a bilinear form induced by a certain matrix J . In particular, the J -symmetric variant includes the complex symmetric form of BL-QMR as a special case. We identify suitable preconditioners for the BL-QMR algorithm applied to multiple radiation and scattering problems. Our numerical tests with the preconditioned BL-QMR algorithm for such multiple linear systems show that, instead of solving each of the linear systems individually, it is significantly more efficient to employ the block version of the iterative method. Moreover, the numerical results clearly illustrate the importance of deflation and its effect on iterative convergence.

Breakdowns and stagnation in iterative methods

One of the disadvantages of Krylov subspace iterative methods is the possibility of breakdown. This occurs when it is impossible to get the next approximation of the solution to the linear system of equationsAu=f. There are two different situations: lucky breakdown, when we have found the solution and hard breakdown, when the next Krylov subspace cannot be generated and/or the next approximate solution (iterate) cannot be computed. We show that some breakdowns depend on the chosen method of generating the basis vectors. Another undesirable feature of the iterative methods is stagnation. This occurs when the error does not change for several iterative steps. We investigate when iterative methods can stagnate and describe conditions which characterize stagnation. We show that in some cases stagnation can imply breakdown.

Fast Nonsymmetric Iterations and Preconditioning for Navier–Stokes Equations

Discretization and linearization of the steady–state Navier-Stokes equations gives rise to a nonsymmetric indefinite linear system of equations. In this paper, we introduce preconditioning techniques for such systems with the property that the eigenvalues of the preconditioned matrices are bounded independently of the mesh size used in the discretization. We confirm and supplement these analytic results with a series of numerical experiments indicating that Krylov subspace iterative methods for nonsymmetric systems display rates of convergence that are independent of the mesh parameter. In addition, we show that preconditioning costs can be kept small by using iterative methods for some intermediate steps performed by the preconditioner.

Assessment of some iterative methods for non‐symmetric linear systems arising in computational fluid dynamics

AbstractVarious tests have been carried out in order to compare the performances of several methods used to solve the non‐symmetric linear systems of equations arising from implicit discretizations of CFD problems, namely the scalar advection‐diffusion equation and the compressible Euler equations. The iterative schemes under consideration belong to three families of algorithms: relaxation (Jacobi and Gauss‐Seidel), gradient and Newton methods. Two gradient methods have been selected: a Krylov subspace iteration method (GMRES) and a non‐symmetric extension of the conjugate gradient method (CGS). Finally, a quasi‐Newton method has also been considered (Broyden). The aim of this paper is to provide indications of which appears to be the most adequate method according to the particular circumstances as well as to discuss the implementation aspects of each scheme.

Formal orthogonal polynomials and Hankel/Toeplitz duality

For classical polynomials orthogonal with respect to a positive measure supported on the real line, the moment matrix is Hankel and positive definite. The polynomials satisfy a three term recurrence relation. When the measure is supported on the complex unit circle, the moment matrix is positive definite and Toeplitz. Then they satisfy a coupled Szegő recurrence relation but also a three term recurrence relation. In this paper we study the generalization for formal polynomials orthogonal with respect to an arbitrary moment matrix and consider arbitrary Hankel and Toeplitz matrices as special cases. The relation with Padé approximation and with Krylov subspace iterative methods is also outlined.

On the squared unsymmetric Lanczos method

The biorthogonal Lanczos and the biconjugate gradient methods have been proposed as iterative methods to approximate the solution of nonsymmetric and idefinite linear systems. Sonneveld (1989) obtained the conjugate gradient squared by squaring the matrix polynomials of the biconjugate gradient method. Here we square the unsymmetric (or biorthogonal) Lanczos method for computing the eigenvalues of nonsymmetric matrices. Three forms of restarted squared Lanczos method for solving unsymmetric linear systems of equations were derived. Numerical experiments with unsymmetric (in) definite linear systems of equations comparing these methods to a restarted (orthogonal) Krylov subspace iterative method showed that the new methods are competitive and they require that a fixed small number of direction vectors be stored in the main memory.

Iterative methods for nonsymmetric systems in DAEs and stiff ODEs codes

Past research demonstrated that the Newton method coupled with Krylov subspace iterative methods is more efficient than the Newton method coupled with Gaussian elimination when used to solve the corrector equation in stiff systems of Ordinary Differential Equations (ODEs)-with large and sparse Jacobian. Some Krylov subspace methods have been incorporated in the mathematical software codes (such as LSODE), which solve numerical stiff systems of ODEs. In this article we review the past results on this topic and three iterative methods for solving nonsymmetric linear systems of equations. We then incorporate these iterative methods in DASSL, which is a program for numerical integration of systems of stiff ODEs and Differential-Algebraic Equations (DAEs). We present numerical tests (with systems of stiff ODEs) and timing comparisons on the CRAY-2 supercomputer.

Experiments with Krylov subspace methods on a massively parallel computer

In this note, we compare some Krylov subspace iterative methods on the MASPAR, a massively parallel computer with 16K processors. In particular, we apply these methods to solve large sparse nonsymmetric linear systems arising from elliptic partial differential equations. The methods under consideration include conjugate gradient type methods, semiiterative methods, and a hybrid variant. Our numerical results show that, on the MASPAR, one should compare iterative methods rather on the basis of total computing time than on the basis of number of iterations required to achieve a given accuracy. Our limited numerical experiments here suggest that, in terms of total computing time, semiiterative and hybrid methods are very attractive for such MASPAR implementations.

Reduced storage matrix methods in stiff ODE systems

Stiff initial-value ODE methods (e.g. BDF) normally require the Jacobian matrix as part of a Newton-like iteration within each implicit time step, and require it to be formed and stored explicitly, whether the linear algebraic method used is direct or iterative. But, a Krylov subspace iterative linear system method, which involves the system matrix only in operator form, can be made part of an inexact Newton method within such a stiff ODE method in a matrix-free manner, requiring no explicit Jacobian matrix storage. Such combinations, using BDF methods, have been implemented with Arnoldi iteration, GMRES (the Generalized Minimum RESidual method), and CG (the Conjugate Gradient Method). Various practical matters (scaling, starting, stopping, etc.) are dealt with in the stiff ODE context. In the context of general nonlinear algebraic systems, we provide some theoretical foundation for the combined Newton-Krylov method by giving convergence results that include errors due to the difference quotient approximation to the linear operator. Earlier tests showed matrix-free methods to be quite effective, at least when the spectrum of the kproblem Jacobian is rather tightly clustered. To improve their robustness, we have added preconditioning, in an experimental solver called LSODPK, which we tested on ODE systems that arise from time-dependent PDE systems by the method of lines. Preconditioner matrices can be formed from the interaction or reaction terms, from the spatial transport terms, or both (as in operator splitting). The additional matrix storage can be reduced greatly by grouping the diagonal blocks in a natural way. The methods appear to be quite effective in improving both speed and storage economy over traditional stiff methods, and over matrix-free methods without preconditioning.

Nonlinearly Preconditioned Krylov Subspace Methods for Discrete Newton Algorithms

We propose an algorithm for implementing Newton’s method for a general nonlinear system $f(x) = 0$ where the linear systems that arise at each step of Newton’s method are solved by a preconditioned Krylov subspace iterative method. The algorithm requires only function evaluations and does not require the evaluation or storage of the Jacobian matrix. Matrix-vector products involving the Jacobian matrix are approximated by directional differences. We develop a framework for constructing preconditionings for this inner iterative method which do not reference the Jacobian matrix explicitly. We derive a nonlinear SSOR type preconditioning which numerical experiments show to be as effective as the linear SSOR preconditioning that uses the Jacobian explicitly.