SIAM Digital Library
 
 
 

SIAM J. on Scientific Computing

Year Range: 

All Journal content published prior to 1997 is part of LOCUS.

Search Issue | RSS Feeds RSS
Next Issue

1998

Volume 19, Issue 1, pp. vii-318

* Introduction to the Special Issue on Iterative Methods for Solving Systems of Algebraic Equations

back to top
RSS Feeds

Introduction to the Special Issue on Iterative Methods for Solving Systems of Algebraic Equations

Howard Elman

SIAM J. Sci. Comput. 19, pp. vii-vii ( pages)

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
This is a special issue of SIAM Journal on Scientific Computing devoted to iterative methods for solving systems of algebraic equations. The development of efficient methods is of fundamental importance in modern computational modeling, as the solution of large systems of equations often represents the dominant computational cost in mathematical simulations. In contrast to direct methods based on Gaussian elimination, iterative methods tend to have storage requirements proportional to the size of the problem. Consequently, they have the potential for handling significantly larger problems than direct methods. Iterative methods are the only option when fine resolution is required, particularly for three-dimensional models. Similarly, they are the only choice for nonlinear problems, and such solvers in the nonlinear case are frequently combined with inner iterations for linear subproblems.
The papers in this volume represent a snapshot of the current state of the art in iterative solution algorithms. The results in these papers were originally presented at the Copper Mountain meeting on iterative methods held in April 1996. This was the eleventh annual meeting in the Copper Mountain series and the fourth devoted exclusively to iterative methods. There were 200 attendees at the conference, and more than 150 presentations were made. The papers included here (identified by first author) can be loosely grouped into five categories:
This collection points to some significant trends in the study of iterative methods, including the maturation of algorithm technology for new "model problems" such as saddle-point and fourth-order problems; the increasing complexity of problems, such as the full potential equation in fluid dynamics, being handled by iterative methods; and the integration of parallel computation and software development into many of the research efforts.
I wish to thank the following members of the Copper Mountain program committee who acted as special editors for this issue:
Through their efforts, submissions were carefully refereed and brought to print on schedule. I also want to thank Tom Manteuffel and Steve McCormick, the organizers of the Copper Mountain conferences, for setting a consistently high standard for the meetings. Special thanks are due to Fred Howes of the Department of Energy, Michael Steuerwalt of the National Science Foundation, and Eric Pritcher from Cray Research for generous support of the meeting.
Our hope is that the diverse topics covered in this issue make it clear that the development and analysis of iterative methods continues to be a challenging and fruitful area of research.
Howard ElmanUniversity of Maryland, College Park

Iterative and Parallel Performance of High-Order Compact Systems

W. F. Spotz and G. F. Carey

SIAM J. Sci. Comput. 19, pp. 1-14 (14 pages) | Cited 2 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Representative transport PDE problems in one, two, and three dimensions are approximated by a class of high-order compact (HOC) difference schemes and their iterative and parallel performance are studied. The eigenvalues and condition numbers of the HOC schemes are analyzed and the performance of standard Krylov-space methods is compared for HOC, central differencing, and standard first-order upwinding schemes. Finally, CPU times, MFLOP rates, and speedup curves are presented for fixed-problem-size cases and scaled-problem-size cases.

Restarted GMRES for Shifted Linear Systems

Andreas Frommer and Uwe Glässner

SIAM J. Sci. Comput. 19, pp. 15-26 (12 pages) | Cited 19 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Shifted matrices, which differ by a multiple of the identity only, generate the same Krylov subspaces with respect to any fixed vector. This fact has been exploited in Lanczos-based methods like CG, QMR, and BiCG to simultaneously solve several shifted linear systems at the expense of only one matrix--vector multiplication per iteration. Here, we develop a variant of the restarted GMRES method exhibiting the same advantage and we investigate its convergence for positive real matrices in some detail. We apply our method to speed up "multiple masses" calculations arising in lattice gauge computations in quantum chromodynamics, one of the most time-consuming supercomputer applications.

Optimizing a Parallel Conjugate Gradient Solver

Martyn R. Field

SIAM J. Sci. Comput. 19, pp. 27-37 (11 pages) | Cited 4 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We develop a highly optimized parallel conjugate gradient solver. We look at both the single node performance and the parallel efficiency. We show that we can solve a problem with 278,000 degrees of freedom on a 32 node Hitachi SR4300 with a performance of 1.1 GFLOPS. We also look at the effect of the quality of mesh partitioning on the performance of the algorithm.

Using Nonorthogonal Lanczos Vectors in the Computation of Matrix Functions

V. Druskin, A. Greenbaum, and L. Knizhnerman

SIAM J. Sci. Comput. 19, pp. 38-54 (17 pages) | Cited 13 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
The Lanczos algorithm uses a three-term recurrence to construct an orthonormal basis for the Krylov space corresponding to a symmetric matrix $A$ and a nonzero starting vector $\varphi$. The vectors and recurrence coefficients produced by this algorithm can be used for a number of purposes, including solving linear systems $Au= \varphi$ and computing the matrix exponential $e^{-tA} \varphi$. Although the vectors produced in finite precision arithmetic are not orthogonal, we show why they can still be used effectively for these purposes.

QMR Smoothing for Lanczos-Type Product Methods Based on Three-Term Rrecurrences

Klaus J. Ressel and Martin H. Gutknecht

SIAM J. Sci. Comput. 19, pp. 55-73 (19 pages) | Cited 4 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
For the solution of large, sparse, non-Hermitian linear systems, Lanczos-type product methods that are based on the Lanczos three-term recurrence are derived in a systematic way. These methods either square the Lanczos process or combine it with a local minimization of the residual. For them a quasi-minimal residual (QMR) smoothing is proposed that can also be implemented by short-term recurrences. The practical performance of these methods and the QMR smoothing is demonstrated in a number of numerical examples.

Some Multigrid Algorithms for Elliptic Problems on Data Parallel Machines

V. A. Bandy, J. E. Dendy, Jr., and W. H. Spangenberg

SIAM J. Sci. Comput. 19, pp. 74-86 (13 pages) | Cited 1 time

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Previously a semicoarsening multigrid algorithm suitable for use on data parallel architectures was investigated. Through the use of new software tools, the performance of this algorithm has been considerably improved. The method has also been extended to three space dimensions. The method performs well for strongly anisotropic problems and for problems with coefficients jumping by orders of magnitude across internal interfaces. The parallel efficiency of this method is analyzed, and its actual performance on the CM-5 is compared with its performance on the CRAY Y-MP and the Sparc-5. A standard coarsening multigrid algorithm is also considered, and we compare its performance on these three platforms as well.

An Evaluation of Parallel Multigrid as a Solver and a Preconditioner for Singularly Perturbed Problems

C. W. Oosterlee and T. Washio

SIAM J. Sci. Comput. 19, pp. 87-110 (24 pages) | Cited 27 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
In this paper we try to achieve h-independent convergence with preconditioned GMRES ([Y. Saad and M. H. Schultz, SIAM J. Sci. Comput., 7 (1986), pp. 856--869]) and BiCGSTAB ([H. A. Van der Vorst, SIAM J. Sci. Comput., 13 (1992), pp. 63--644]) for two-dimensional (2D) singularly perturbed equations. Three recently developed multigrid methods are adopted as a preconditioner. They are also used as solution methods in order to compare the performance of the methods as solvers and as preconditioners.
Two of the multigrid methods differ only in the transfer operators. One uses standard matrix-dependent prolongation operators from [J. E. Dendy Jr., J. Comput. Phys., 48 (1982), pp. 366--386], [W. Hackbusch, Multi-grid Methods and Applications, Springer, Berlin, 1985]. The second uses "upwind" prolongation operators, developed in [P. M. de Zeeuw, J. Comput.\ Appl.\ Math., 33 (1990), pp. 1--27]. Both employ the Galerkin coarse grid approximation and an alternating zebra line Gauss--Seidel smoother. The third method is based on the block LU decomposition of a matrix and on an approximate Schur complement. This multigrid variant is presented in [A. Reusken, A Multigrid Method Based on Incomplete Gaussian Elimination, University of Eindhoven, the Netherlands, 1995]. All three multigrid algorithms are algebraic methods.
The eigenvalue spectra of the three multigrid iteration matrices are analyzed for the equations solved in order to understand the convergence of the three algorithms. Furthermore, the construction of the search directions for the multigrid preconditioned GMRES solvers is investigated by the calculation and solution of the minimal residual polynomials.
For Poisson and convection-diffusion problems all solution methods are investigated and evaluated for finite volume discretizations on fine grids. The methods have been parallelized with a grid partitioning technique and are compared on an MIMD machine.

Fast Multigrid Solution of the Advection Problem with Closed Characteristics

Irad Yavneh, Cornelis H. Venner, and Achi Brandt

SIAM J. Sci. Comput. 19, pp. 111-125 (15 pages) | Cited 16 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
The numerical solution of the advection-diffusion problem in the inviscid limit with closed characteristics is studied as a prelude to an efficient high Reynolds-number flow solver. It is demonstrated by a heuristic analysis and numerical calculations that using upstream discretization with downstream relaxation ordering in a multigrid cycle with appropriate residual weighting leads to an efficient solution process. Upstream finite-difference approximations to the advection operator are derived whose truncation terms approximate "physical" (Laplacian) viscosity, thus avoiding spurious solutions to the homogeneous problem when the artificial diffusivity dominates the physical viscosity [A. Brandt and I. Yavneh, J. Comput. Phys., 93 (1991), pp. 128--143].

An Efficient Solver for Multi--Right-Hand-Side Linear Systems Based on the CCCG($\eta$) Method with Applications to Implicit Time-Dependent Partial Differential Equations

Frederico F. Campos and Nick R. C. Birkett

SIAM J. Sci. Comput. 19, pp. 126-138 (13 pages) | Cited 3 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
The controlled Cholesky factorization (CCF) has been shown to be a robust preconditioner for the conjugate gradient (CG) method. In this scheme the amount of fill-in is defined in terms of a parameter $\eta$, the number of extra nonzero elements allowed per column. This parameter $\eta$ can be chosen in such a way that the preconditioning factor will require more or less storage than the coefficient matrix. The selection of nonzero elements kept for the preconditioner is made by value instead of by position, so the sparsity pattern of the original matrix is never considered by the proposed factorization. This preconditioned method is termed the CCCG($\eta$) method and is a generalization of a method developed by Jones and Plassmann [\textit{ACM Trans.\ Math.\ Software}, 21 (1995), pp.\ 5--17]. It is demonstrated how an optimum value of $\eta$ can be automatically determined when solving systems with multiple right-hand sides. Typically such systems occur when using time-stepping procedures. Numerical examples are presented for some time-dependent parabolic partial differential equations discretized using linear triangular finite elements on unstructured grids. Comparisons made between CCCG($\eta$), standard ICCG, and the Jones and Plassmann method show CCCG($\eta$) to be an efficient solver.

Superlinear Convergence Estimates for a Conjugate Gradient Method for the Biharmonic Equation

Raymond H. Chan, Thomas K. DeLillo, and Mark A. Horn

SIAM J. Sci. Comput. 19, pp. 139-147 (9 pages)

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
The method of Muskhelishvili for solving the biharmonic equation using conformal mapping is investigated. In [R. H. Chan, T. K. DeLillo, and M. A. Horn, SIAM J. Sci. Comput., 18 (1997), pp. 1571--1582] it was shown, using the Hankel structure, that the linear system in [N. I. Muskhelishvili, Some Basic Problems of the Mathematical Theory of Elasticity, Noordhoff, Groningen, the Netherlands] is the discretization of the identity plus a compact operator, and therefore the conjugate gradient method will converge superlinearly. Estimates are given here of the superlinear convergence in the cases when the boundary curve is analytic or in a Hölder class.

State Space Orderings for Gauss--Seidel in Markov Chains Revisited

Tugrul Dayar

SIAM J. Sci. Comput. 19, pp. 148-154 (7 pages) | Cited 4 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
States of a Markov chain may be reordered to reduce the magnitude of the subdominant eigenvalue of the Gauss--Seidel (GS) iteration matrix. Orderings that maximize the elemental mass or the number of nonzero elements in the dominant term of the GS splitting (that is, the term approximating the coefficient matrix) do not necessarily converge faster. An ordering of a Markov chain that satisfies Property-R is semiconvergent. On the other hand, there are semiconvergent state space orderings that do not satisfy Property-R. For a given ordering, a simple approach for checking Property-R is shown. Moreover, a version of the Cuthill--McKee algorithm may be used to order the states of a Markov chain so that Property-R is satisfied. The computational complexity of the ordering algorithm is less than that of a single GS iteration. In doing all this, the aim is to gain insight into (faster) converging orderings.

Numerical Conformal Mapping Methods for Simply and Doubly Connected Regions

Thomas K. DeLillo and John A. Pfaltzgraff

SIAM J. Sci. Comput. 19, pp. 155-171 (17 pages) | Cited 2 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Methods are presented and analyzed for approximating the conformal map from the interior (exterior) of the disk to the interior (exterior) of a smooth, simple closed curve and from an annulus to a bounded, doubly connected region with smooth boundaries. The methods are Newton-like methods for computing the boundary correspondences and conformal moduli similar to Fornberg's method for the interior of the disk. We show that the linear systems are discretizations of the identity plus a compact operator and, hence, that the conjugate gradient method converges superlinearly.

Block-Triangular Preconditioners for Saddle Point Problems with a Penalty Term

Axel Klawonn

SIAM J. Sci. Comput. 19, pp. 172-184 (13 pages) | Cited 59 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Block-triangular preconditioners for a class of saddle point problems with a penalty term are considered. An important example is the mixed formulation of the pure displacement problem in linear elasticity. It is shown that the spectrum of the preconditioned system is contained in a real, positive interval and that the interval bounds can be made independent of the discretization and penalty parameters. This fact is used to construct bounds of the convergence rate of the GMRES method with respect to an energy norm. Numerical results are given for GMRES and BI-CGSTAB.

Approximate Schur Complement Preconditioning of the Lowest-Order Nodal Discretizations

J. David Moulton, Jim E. Morel, and Uri M. Ascher

SIAM J. Sci. Comput. 19, pp. 185-205 (21 pages) | Cited 1 time

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Certain classes of nodal methods and mixed-hybrid finite element methods lead to equivalent, robust, and accurate discretizations of second-order elliptic PDEs. However, widespread popularity of these discretizations has been hindered by the awkward linear systems which result. The present work overcomes this awkwardness and develops preconditioners which yield solution algorithms for these discretizations with an efficiency comparable to that of the multigrid method for standard discretizations.
Our approach exploits the natural partitioning of the linear system obtained by the mixed-hybrid finite element method. By eliminating different subsets of unknowns, two Schur complements are obtained with known structure. Replacing key matrices in this structure by lumped approximations, we define three optimal preconditioners. Central to the optimal performance of these preconditioners is their sparsity structure which is compatible with standard finite difference discretizations and hence treated adequately with only a single multigrid cycle.
In this paper we restrict the discussion to the two-dimensional case; these techniques are readily extended to three dimensions.

An Efficient Iterative Method for the Generalized Stokes Problem

Vivek Sarin and Ahmed Sameh

SIAM J. Sci. Comput. 19, pp. 206-226 (21 pages) | Cited 12 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
The generalized Stokes problem, which arises frequently in the simulation of time-dependent Navier--Stokes equations for incompressible fluid flow, gives rise to symmetric linear systems of equations. These systems are indefinite due to a set of linear constraints on the velocity, causing difficulty for most preconditioners and iterative methods. This paper presents a novel method to obtain a preconditioned linear system from the original one which is then solved by an iterative method. This new method generates a basis for the velocity space and solves a reduced system which is symmetric and positive definite. Numerical experiments indicating superior convergence compared to existing methods are presented. A natural extension of this method to elliptic problems is also proposed, along with theoretical bounds on the rate of convergence, and results of experiments demonstrating robust and effective preconditioning.

Dynamic Thick Restarting of the Davidson, and the Implicitly Restarted Arnoldi Methods

Andreas Stathopoulos, Yousef Saad, and Kesheng Wu

SIAM J. Sci. Comput. 19, pp. 227-245 (19 pages) | Cited 31 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
The Davidson method is a popular preconditioned variant of the Arnoldi method for solving large eigenvalue problems. For theoretical as well as practical reasons the two methods are often used with restarting. Frequently, information is saved through approximated eigenvectors to compensate for the convergence impairment caused by restarting. We call this scheme of retaining more eigenvectors than needed "thick restarting" and prove that thick restarted, nonpreconditioned Davidson is equivalent to the implicitly restarted Arnoldi. We also establish a relation between thick restarted Davidson and a Davidson method applied on a deflated system. The theory is used to address the question of which and how many eigenvectors to retain and motivates the development of a dynamic thick restarting scheme for the symmetric case, which can be used in both Davidson and implicit restarted Arnoldi. Several experiments demonstrate the efficiency and robustness of the scheme.

Parallel Newton--Krylov--Schwarz Algorithms for the Transonic Full Potential Equation

Xiao-Chuan Cai, William D. Gropp, David E. Keyes, Robin G. Melvin, and David P. Young

SIAM J. Sci. Comput. 19, pp. 246-265 (20 pages) | Cited 30 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We study parallel two-level overlapping Schwarz algorithms for solving nonlinear finite element problems, in particular, for the full potential equation of aerodynamics discretized in two dimensions with bilinear elements. The overall algorithm, Newton--Krylov--Schwarz (NKS), employs an inexact finite difference Newton method and a Krylov space iterative method, with a two-level overlapping Schwarz method as a preconditioner. We demonstrate that NKS, combined with a density upwinding continuation strategy for problems with weak shocks, is robust and economical for this class of mixed elliptic-hyperbolic nonlinear partial differential equations, with proper specification of several parameters. We study upwinding parameters, inner convergence tolerance, coarse grid density, subdomain overlap, and the level of fill-in in the incomplete factorization, and report their effect on numerical convergence rate, overall execution time, and parallel efficiency on a distributed-memory parallel computer.

A Fast Multilevel Algorithm for the Solution of Nonlinear Systems of Conductive-Radiative Heat Transfer Equations

J. M. Banoczi and C. T. Kelley

SIAM J. Sci. Comput. 19, pp. 266-279 (14 pages) | Cited 10 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
In this paper we describe and analyze a fast multilevel algorithm for the solution of a system of nonlinear integro-differential equations that model steady-state combined conductive-radiative heat transfer. This system of equations for radiative intensity and temperature can be formulated as a compact fixed point problem in temperature alone with a fixed point map that requires both a solution of the linear transport equation and the linear heat equation for its evaluation. We obtain an efficient evaluation of the fixed point map by coupling a finite element diffusion solver with a fast transport solver developed by the second author. As a solver we apply a modification of the Atkinson--Brakhage method, with Newton--GMRES as the coarse mesh solver, to the full nonlinear system. We compare our discretization/solver pair with Newton--GMRES and the classical Atkinson--Brakhage algorithm.

Termination of Newton/Chord Iterations and the Method of Lines

C. T. Kelley, C. T. Miller, and M. D. Tocci

SIAM J. Sci. Comput. 19, pp. 280-290 (11 pages) | Cited 7 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Many ordinary differential equation (ODE) and differential algebraic equation (DAE) codes terminate the nonlinear iteration for the corrector equation when the difference between successive iterates (the step) is sufficiently small. This termination criterion avoids the expense of evaluating the nonlinear residual at the final iterate. Similarly, Jacobian information is not usually computed at every time step but only when certain tests indicate that the cost of a new Jacobian is justified by the improved performance in the nonlinear iteration. In this paper, we show how an out-of-date Jacobian coupled with moderate ill conditioning can lead to premature termination of the corrector iteration and suggest ways in which this situation can be detected and remedied. As an example, we consider the method of lines (MOL) solution of Richards' equation (RE), which models flow through variably saturated porous media. When the solution to this problem has a sharp moving front, and the Jacobian is even slightly ill conditioned, the corrector iteration used in many integrators can terminate prematurely, leading to incorrect results. While this problem can be solved by tightening the tolerances for the solvers used in the temporal integration, it is more efficient to modify the termination criteria of the nonlinear solver and/or recompute the Jacobian more frequently. Of these two, recomputation of the Jacobian is the more important. We propose a criterion based on an estimate of the norm of the time derivative of the Jacobian for recomputation of the Jacobian and a second criterion based on a condition estimate for tightening of the termination criteria of the nonlinear solver.

Enhanced Nonlinear Iterative Techniques Applied to a Nonequilibrium Plasma Flow

D. A. Knoll and P. R. McHugh

SIAM J. Sci. Comput. 19, pp. 291-301 (11 pages) | Cited 21 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We study the application of enhanced nonlinear iterative methods to the steady-state solution of a system of two-dimensional convection-diffusion-reaction partial differential equations that describe the partially ionized plasma flow in the boundary layer of a tokamak fusion reactor. This system of equations is characterized by multiple time and spatial scales and contains highly anisotropic transport coefficients due to a strong imposed magnetic field. We use Newton's method to linearize the nonlinear system of equations resulting from an implicit, finite volume discretization of the governing partial differential equations, on a staggered Cartesian mesh. The resulting linear systems are neither symmetric nor positive definite, and are poorly conditioned. Preconditioned Krylov iterative techniques are employed to solve these linear systems. We investigate both a modified and a matrix-free Newton--Krylov implementation, with the goal of reducing CPU cost associated with the numerical formation of the Jacobian. A combination of a damped iteration, mesh sequencing, and a pseudotransient continuation technique is used to enhance global nonlinear convergence and CPU efficiency. GMRES is employed as the Krylov method with incomplete lower-upper (ILU) factorization preconditioning. The goal is to construct a combination of nonlinear and linear iterative techniques for this complex physical problem that optimizes trade-offs between robustness, CPU time, memory requirements, and code complexity. It is shown that a mesh sequencing implementation provides significant CPU savings for fine grid calculations. Performance comparisons of modified Newton--Krylov and matrix-free Newton--Krylov algorithms will be presented.

NITSOL: A Newton Iterative Solver for Nonlinear Systems

Michael Pernice and Homer F. Walker

SIAM J. Sci. Comput. 19, pp. 302-318 (17 pages) | Cited 68 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We introduce a well-developed Newton iterative (truncated Newton) algorithm for solving large-scale nonlinear systems. The framework is an inexact Newton method globalized by backtracking. Trial steps are obtained using one of several Krylov subspace methods. The algorithm is implemented in a Fortran solver called NITSOL that is robust yet easy to use and provides a number of useful options and features. The structure offers the user great flexibility in addressing problem specificity through preconditioning and other means and allows easy adaptation to parallel environments. Features and capabilities are illustrated in numerical experiments.
Close

close