Abstract: Recently, a number of researchers have investigated a class of graph partitioning algorithms that reduce the size of the graph by collapsing vertices and edges, partition the smaller graph, and then uncoarsen it to construct a partition for the original graph [Bui and Jones, Proc. of the 6th SIAM Conference on Parallel Processing for Scientific Computing, 1993, 445--452; Hendrickson and Leland, A Multilevel Algorithm for Partitioning Graphs, Tech. report SAND 93-1301, Sandia National Laboratories, Albuquerque, NM, 1993]. From the early work it was clear that multilevel techniques held great promise; however, it was not known if they can be made to consistently produce high quality partitions for graphs arising in a wide range of application domains. We investigate the effectiveness of many different choices for all three phases: coarsening, partition of the coarsest graph, and refinement. In particular, we present a new coarsening heuristic (called heavy-edge heuristic) for which the size of the partition of the coarse graph is within a small factor of the size of the final partition obtained after multilevel refinement. We also present a much faster variation of the Kernighan--Lin (KL) algorithm for refining during uncoarsening. We test our scheme on a large number of graphs arising in various domains including finite element methods, linear programming, VLSI, and transportation. Our experiments show that our scheme produces partitions that are consistently better than those produced by spectral partitioning schemes in substantially smaller time. Also, when our scheme is used to compute fill-reducing orderings for sparse matrices, it produces orderings that have substantially smaller fill than the widely used multiple minimum degree algorithm.

Abstract: We study the numerical integration of large stiff systems of differential equations by methods that use matrix--vector products with the exponential or a related function of the Jacobian. For large problems, these can be approximated by Krylov subspace methods, which typically converge faster than those for the solution of the linear systems arising in standard stiff integrators. The exponential methods also offer favorable properties in the integration of differential equations whose Jacobian has large imaginary eigenvalues. We derive methods up to order 4 which are exact for linear constant-coefficient equations. The implementation of the methods is discussed. Numerical experiments with reaction-diffusion problems and a time-dependent Schrodinger equation are included.

Abstract: The positivity principle and positive schemes to solve multidimensional hyperbolic systems of conservation laws have been introduced in [X.-D. Liu and P. D. Lax, J. Fluid Dynam., 5 (1996), pp. 133--156]. Some numerical experiments presented there show how well the method works. In this paper we use positive schemes to solve Riemann problems for two-dimensional gas dynamics.

Abstract: In 1994 Berenger showed how to construct a perfectly matched absorbing layer for the Maxwell system in rectilinear coordinates. This layer absorbs waves of any wavelength and any frequency without reflection and thus can be used to artificially terminate the domain of scattering calculations. In this paper we show how to derive and implement the Berenger layer in curvilinear coordinates (in two space dimensions). We prove that an infinite layer of this type can be used to solve time harmonic scattering problems. We also show that the truncated Berenger problem has a solution except at a discrete set of exceptional frequencies (which might be empty). Finally numerical results show that the curvilinear layer can produce accurate solutions in the time and frequency domain.

Abstract: Discretization methods are proposed for control-volume formulations on polygonal and triangular grid cells in two space dimensions. The methods are applicable for any system of conservation laws where the flow density is defined by a gradient law, like Darcy's law for porous-media flow. A strong feature of the methods is the ability to handle media inhomogeneities in combination with full-tensor anisotropy. This paper gives a derivation of the methods, and the relation to previously published methods is also discussed. A further discussion of the methods, including numerical examples, is given in the companion paper, Part II [SIAM J. Sci. Comput., pp. 1717--1736].

Abstract: We describe a three-dimensional front tracking algorithm, discuss its numerical implementation, and present studies to validate the correctness of this approach. Based on the results of the two-dimensional computations, we expect three-dimensional front tracking to significantly improve computational efficiencies for problems dominated by discontinuities. In some cases, for which the interface computations display considerable numerical sensitivity, we expect a greatly enhanced capability.

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.

Abstract: The standard incomplete LU (ILU) preconditioners often fail for general sparse indefinite matrices because they give rise to "unstable" factors L and U. In such cases, it may be attractive to approximate the inverse of the matrix directly. This paper focuses on approximate inverse preconditioners based on minimizing ||I-AM||F, where AM is the preconditioned matrix. An iterative descent-type method is used to approximate each column of the inverse. For this approach to be efficient, the iteration must be done in sparse mode, i.e., with "sparse-matrix by sparse-vector" operations. Numerical dropping is applied to maintain sparsity; compared to previous methods, this is a natural way to determine the sparsity pattern of the approximate inverse. This paper describes Newton, "global," and column-oriented algorithms, and discusses options for initial guesses, self-preconditioning, and dropping strategies. Some limited theoretical results on the properties and convergence of approximate inverses are derived. Numerical tests on problems from the Harwell--Boeing collection and the FIDAP fluid dynamics analysis package show the strengths and limitations of approximate inverses. Finally, some ideas and experiments with practical variations and applications are presented.

Abstract: Restoration of images that have been blurred by the effects of a Gaussian blurring function is an ill-posed but well-studied problem. Any blur that is spatially invariant can be expressed as a convolution kernel in an integral equation. Fast and effective algorithms then exist for determining the original image by preconditioned iterative methods. If the blurring function is spatially variant, however, then the problem is more difficult. In this work we develop fast algorithms for forming the convolution and for recovering the original image when the convolution functions are spatially variant but have a small domain of support. This assumption leads to a discrete problem involving a banded matrix. We devise an effective preconditioner and prove that the preconditioned matrix differs from the identity by a matrix of small rank plus a matrix of small norm. Numerical examples are given, related to the Hubble Space Telescope (HST) Wide-Field/Planetary Camera. The algorithms that we develop are applicable to other ill-posed integral equations as well.

Abstract: Waveform relaxation algorithms for partial differential equations (PDEs) are traditionally obtained by discretizing the PDE in space and then splitting the discrete operator using matrix splittings. For the semidiscrete heat equation one can show linear convergence on unbounded time intervals and superlinear convergence on bounded time intervals by this approach. However, the bounds depend in general on the mesh parameter and convergence rates deteriorate as one refines the mesh. Motivated by the original development of waveform relaxation in circuit simulation, where the circuits are split in the physical domain into subcircuits, we split the PDE by using overlapping domain decomposition. We prove linear convergence of the algorithm in the continuous case on an infinite time interval, at a rate depending on the size of the overlap. This result remains valid after discretization in space and the convergence rates are robust with respect to mesh refinement. The algorithm is in the class of waveform relaxation algorithms based on overlapping multisplittings. Our analysis quantifies the empirical observation by Jeltsch and Pohl [SIAM J. Sci. Comput., 16 (1995), pp. 40--49] that the convergence rate of a multisplitting algorithm depends on the overlap. Numerical results are presented which support the convergence theory.

Abstract: We investigate a method of dividing an irregular mesh into equal-sized pieces with few interconnecting edges. The method's novel feature is that it exploits the geometric coordinates of the mesh vertices. It is based on theoretical work of Miller, Teng, Thurston, and Vavasis, who showed that certain classes of "well-shaped" finite-element meshes have good separators. The geometric method is quite simple to implement: we describe a \sc Matlab code for it in some detail. The method is also quite efficient and effective: we compare it with some other methods, including spectral bisection.

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.

Abstract: In this paper we describe a new algorithm for the calculation of consistent initial conditions for a class of systems of differential-algebraic equations which includes semi-explicit index-one systems. We consider initial condition problems of two types---one where the differential variables are specified, and one where the derivative vector is specified. The algorithm requires a minimum of additional information from the user. We outline the implementation in a general-purpose solver DASPK for differential-algebraic equations, and present some numerical experiments which illustrate its effectiveness.

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.

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.

Abstract: The numerical study of aeroacoustic problems places stringent demands on the choice of a computational algorithm because it requires the ability to propagate disturbances of small amplitude and short wavelength. The demands are particularly high when shock waves are involved because the chosen algorithm must also resolve discontinuities in the solution. The extent to which a high-order accurate shock-capturing method can be relied upon for aeroacoustics applications that involve the interaction of shocks with other waves has not been previously quantified. Such a study is initiated in this work. A fourth-order accurate essentially nonoscillatory (ENO) method is used to investigate the solutions of inviscid, compressible flows with shocks. The design order of accuracy is achieved in the smooth regions of a steady-state, quasi-one-dimensional test case. However, in an unsteady test case, only first-order results are obtained downstream of a sound-shock interaction. The difficulty in obtaining a globally high-order accurate solution in such a case with a shock-capturing method is demonstrated through the study of a simplified, linear model problem. Some of the difficult issues and ramifications for aeroacoustic simulations of flows with shocks that are raised by these results are discussed.

Abstract: For large sparse systems of linear equations iterative solution techniques are attractive. In this paper we propose and examine the convergence of an iterative method for an important class of nonsymmetric and indefinite coefficient matrices based on the use of an indefinite and symmetric preconditioner. We apply our technique to the linearized Navier--Stokes equations (the Oseen equations).

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.

Abstract: In this paper we discuss several sparse-matrix-based preconditioners suitable for preconditioning dense linear systems of boundary element equations. All preconditioners involve only O(n) nonzeros. We provide a framework for constructing operator-splitting-based preconditioners and use it to analyze a class of sparse preconditioners. For singular integral equations, a more efficient preconditioner is proposed that has a band-2 structure.

Abstract: This paper presents a new approach to the numerical solution of boundary value problems for higher-index differential-algebraic equations (DAEs). Invariants known for the original DAE as well as invariants of the reduced index 1 formulation are exploited to stabilize initial value problem (IVP) integration, derivative generation, and nonlinear and linear systems solution of an enhanced multiple shooting method. Extensions to collocation are given. Applications are presented for two important problem classes: parameter estimation in multibody systems given in descriptor form, and singular and state-constrained optimal control problems. In particular, generalizations of the "internal numerical differentiation" technique to DAE with invariants and a new multistage least squares decomposition technique for DAE boundary value problems are developed, which are implemented in the multiple shooting code PARFIT and in the collocation code COLFIT. Further, a method is described which minimizes the number of necessary directional derivatives in the presence of multipoint conditions and invariants. As numerical applications, a parameter identification problem for a slider crank mechanism and a periodic cruise optimal control problem for a motor glider aircraft are treated.

Abstract: The purpose of this paper is to present a semicoarsening multigrid algorithm for solving the finite difference discretization of symmetric and nonsymmetric, two- and three-dimensional elliptic partial differential equations with highly discontinuous and anisotropic coefficients. The discrete equations are assumed to be defined on a logically rectangular grid, obtained possibly through grid generation for a problem defined on an irregular domain. The basic algorithm is described along with some modifications which are designed to improve its efficiency and robustness for certain types of problem cases. FORTRAN codes that implement the two- and three-dimensional semicoarsening multigrid algorithms are described briefly, and numerical results are presented.

Abstract: Moving mesh methods based on the equidistribution principle (EP) are studied from the viewpoint of stability of the moving mesh system of differential equations. For fine spatial grids, the moving mesh system inherits the stability of the original discretized partial differential equation (PDE). Unfortunately, for some PDEs the moving mesh methods require so many spatial grid points that they no longer appear to be practical. Failures and successes of the moving mesh method applied to three reaction-diffusion problems are explained via an analysis of the stability and accuracy of the moving mesh PDE.

Abstract: We present an efficient and provably good partitioning and load balancing algorithm for parallel adaptive N-body simulation. The main ingredient of our method is a novel geometric characterization of a class of communication graphs that can be used to support hierarchical N-body methods such as the fast multipole method (FMM) and the Barnes--Hut method (BH). We show that communication graphs of these methods have a good partition that can be found efficiently sequentially and in parallel. In particular, we show that an N-body communication graph (either for BH or for FMM) can be partitioned into two subgraphs with equal computation load by removing only $O(\sqrt{n\log n})$ and O(n2/3(log n)1/3) number of nodes, respectively, for two and three dimensions. These bounds on node-partition imply bounds on edge-partition of $O(\sqrt{n}(\log n)^{3/2})$ and O(n2/3(log n)4/3), respectively, for two and three dimensions. To the best of our knowledge, this is the first theoretical result on the quality of partitioning N-body communication graphs for nonuniformly distributed particles. Our results imply that parallel adaptive N-body simulation can be made as scalable as computation on regular grids and as efficient as parallel N-body simulation on uniformly distributed particles.

Abstract: This paper is concerned with developing methods for the propagation of linear waves in several space dimensions. The methods are time reversible and hence are free from numerical dissipation. They are based on bicharacteristic forms of the governing equations and are made possible by adopting forms of staggered storage that depend on the precise equations under consideration. Analysis is presented for the equations of acoustics, electromagnetics, and elastodynamics.

Abstract: Stochastic integration rules are derived for infinite integration intervals, generalizing rules developed by Siegel and O'Brien [ SIAM J Sci Statist Comput, 6 (1985), pp 169--181] for finite intervals Then random orthogonal transformations of rules for integrals over the surface of the unit m-sphere are used to produce stochastic rules for these integrals The two types of rules are combined to produce stochastic rules for multidimensional integrals over infinite regions with Normal or Student-t weights Example results are presented to illustrate the effectiveness of the new rules

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.

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.

Abstract: In this paper we propose a direct method for the solution of the Poisson equation in rectangular regions. It has an arbitrary order accuracy and low CPU requirements which makes it practical for treating large-scale problems. The method is based on a pseudospectral Fourier approximation and a polynomial subtraction technique. Fast convergence of the Fourier series is achieved by removing the discontinuities at the corner points using polynomial subtraction functions. These functions have the same discontinuities at the corner points as the sought solution. In addition to this, they satisfy the Laplace equation so that the subtraction procedure does not generate nonperiodic, nonhomogeneous terms. The solution of a boundary value problem is obtained in a series form in O(N log N) floating point operations, where N2 is the number of grid nodes. Evaluating the solution at all N2 interior points requires O(N2 log N) operations.

Abstract: The standard approach to measuring the condition of a linear system compresses all sensitivity information into one number. Thus a loss of information can occur in situations in which the standard condition number with respect to inversion does not accurately reflect the actual sensitivity of a solution or particular entries of a solution. It is shown that a new method for estimating the sensitivity of linear systems addresses these difficulties. The new procedure measures the effects on the solution of small random changes in the input data and, by properly scaling the results, obtains reliable condition estimates for each entry of the computed solution. Moreover, this approach, which is referred to as small-sample statistical condition estimation, is no more costly than the standard 1-norm or power method 2-norm condition estimates, and it has the advantage of considerable flexibility. For example, it easily accommodates restrictions on, or structure associated with, allowable perturbations. The method also has a rigorous statistical theory available for the probability of accuracy of the condition estimates. However, it gives no estimate of an approximate null vector for nearly singular systems. The theory of this approach is discussed along with several illustrative examples.