Showing papers in "SIAM Journal on Scientific Computing in 1998"
TL;DR: Basis Pursuit (BP) is a principle for decomposing a signal into an "optimal" superposition of dictionary elements, where optimal means having the smallest l1 norm of coefficients among all such decompositions.
Abstract: The time-frequency and time-scale communities have recently developed a large number of overcomplete waveform dictionaries --- stationary wavelets, wavelet packets, cosine packets, chirplets, and warplets, to name a few. Decomposition into overcomplete systems is not unique, and several methods for decomposition have been proposed, including the method of frames (MOF), Matching pursuit (MP), and, for special dictionaries, the best orthogonal basis (BOB).
Basis Pursuit (BP) is a principle for decomposing a signal into an "optimal" superposition of dictionary elements, where optimal means having the smallest l1 norm of coefficients among all such decompositions. We give examples exhibiting several advantages over MOF, MP, and BOB, including better sparsity and superresolution. BP has interesting relations to ideas in areas as diverse as ill-posed problems, in abstract harmonic analysis, total variation denoising, and multiscale edge denoising.
BP in highly overcomplete dictionaries leads to large-scale optimization problems. With signals of length 8192 and a wavelet packet dictionary, one gets an equivalent linear program of size 8192 by 212,992. Such problems can be attacked successfully only because of recent advances in linear programming by interior-point methods. We obtain reasonable success with a primal-dual logarithmic barrier method and conjugate-gradient solver.
9,950 citations
TL;DR: This work presents 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 theSize of the final partition obtained after multilevel refinement, and presents a much faster variation of the Kernighan--Lin (KL) algorithm for refining during uncoarsening.
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.
5,629 citations
TL;DR: This work studies 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, and derives methods up to order 4 which are exact for linear constant-coefficient equations.
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.
533 citations
TL;DR: This paper uses positive schemes to solve Riemann problems for two-dimensional gas dynamics to show how well the positivity principle works.
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.
471 citations
TL;DR: Three-dimensional front tracking is described, its numerical implementation is discussed, and studies to validate the correctness of this approach are presented to improve computational efficiencies for problems dominated by discontinuities.
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.
377 citations
TL;DR: This paper gives a derivation of the methods ofretretization, and the relation to previously published methods is also discussed.
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].
374 citations
TL;DR: A new nonoscillatory high-resolution scheme for two-dimensional hyperbolic conservation laws, a predictor-corrector method which consists of two steps, which proves that the scheme satisfies the scalar maximum principle, and demonstrates the application of the central scheme to several prototype two- dimensional Euler problems.
Abstract: We construct, analyze, and implement a new nonoscillatory high-resolution scheme for two-dimensional hyperbolic conservation laws. The scheme is a predictor-corrector method which consists of two steps: starting with given cell averages, we first predict pointvalues which are based on nonoscillatory piecewise-linear reconstructions from the given cell averages; at the second corrector step, we use staggered averaging, together with the predicted midvalues, to realize the evolution of these averages. This results in a second-order, nonoscillatory central scheme, a natural extension of the one-dimensional second-order central scheme of Nessyahu and Tadmor [J. Comput. Phys., 87 (1990), pp. 408--448].
As in the one-dimensional case, the main feature of our two-dimensional scheme is simplicity. In particular, this central scheme does not require the intricate and time-consuming (approximate) Riemann solvers which are essential for the high-resolution upwind schemes; in fact, even the computation of the exact Jacobians can be avoided. Moreover, the central scheme is "genuinely multidimensional" in the sense that it does not necessitate dimensional splitting.
We prove that the scheme satisfies the scalar maximum principle, and in the more general context of systems, our proof indicates that the scheme is positive (in the sense of Lax and Liu [CFD Journal, 5 (1996), pp. 1--24]). We demonstrate the application of our central scheme to several prototype two-dimensional Euler problems. Our numerical experiments include the resolution of shocks oblique to the computational grid; they show how our central scheme solves with high resolution the intricate wave interactions in the so-called double Mach reflection problem [J. Comput. Phys., 54 (1988), pp. 115--173] without following the characteristics; and finally we report on the accurate ray solutions of a weakly hyperbolic system [J. Comput. Appl. Math., 74 (1996), pp. 175--192], rays which otherwise are missed by the dimensional splitting approach. Thus, a considerable amount of simplicity and robustness is gained while achieving stability and high resolution.
371 citations
TL;DR: It is proved that an infinite layer of this type can be used to solve time harmonic scattering problems and numerical results show that the curvilinear layer can produce accurate solutions in the time and frequency domain.
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.
369 citations
TL;DR: Two algorithms, JDQZ for the generalized eigen problem and JDQR for the standard eigenproblem, that are based on the iterative construction of a (generalized) partial Schur form are presented, suitable for the efficient computation of several eigenvalues and the corresponding eigenvectors near a user-specified target value in the complex plane.
Abstract: Recently the Jacobi--Davidson subspace iteration method has been introduced as a new powerful technique for solving a variety of eigenproblems. In this paper we will further exploit this method and enhance it with several techniques so that practical and accurate algorithms are obtained. We will present two algorithms, JDQZ for the generalized eigenproblem and JDQR for the standard eigenproblem, that are based on the iterative construction of a (generalized) partial Schur form. The algorithms are suitable for the efficient computation of several (even multiple) eigenvalues and the corresponding eigenvectors near a user-specified target value in the complex plane. An attractive property of our algorithms is that explicit inversion of operators is avoided, which makes them potentially attractive for very large sparse matrix problems.
We will show how effective restarts can be incorporated in the Jacobi--Davidson methods, very similar to the implicit restart procedure for the Arnoldi process. Then we will discuss the use of preconditioning, and, finally, we will illustrate the behavior of our algorithms by a number of well-chosen numerical experiments.
339 citations
TL;DR: A procedure for computing an incomplete factorization of the inverse of a nonsymmetric matrix is developed, and the resulting factorized sparse approximate inverse is used as an explicit preconditioner for conjugate gradient--type methods.
Abstract: This paper is concerned with a new approach to preconditioning for large, sparse linear systems. A procedure for computing an incomplete factorization of the inverse of a nonsymmetric matrix is developed, and the resulting factorized sparse approximate inverse is used as an explicit preconditioner for conjugate gradient--type methods. Some theoretical properties of the preconditioner are discussed, and numerical experiments on test matrices from the Harwell--Boeing collection and from Tim Davis's collection are presented. Our results indicate that the new preconditioner is cheaper to construct than other approximate inverse preconditioners. Furthermore, the new technique insures convergence rates of the preconditioned iteration which are comparable with those obtained with standard implicit preconditioners.
325 citations
TL;DR: A well-developed Newton iterative (truncated Newton) algorithm for solving large-scale nonlinear systems that is implemented in a Fortran solver called NITSOL that is robust yet easy to use and provides a number of useful options and features.
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.
TL;DR: It is shown how the structure of the Gaussian model can be exploited to yield efficient algorithms for agglomerative hierarchical clustering.
Abstract: Agglomerative hierarchical clustering methods based on Gaussian probability models have recently shown promise in a variety of applications. In this approach, a maximum-likelihood pair of clusters is chosen for merging at each stage. Unlike classical methods, model-based methods reduce to a recurrence relation only in the simplest case, which corresponds to the classical sum of squares method. We show how the structure of the Gaussian model can be exploited to yield efficient algorithms for agglomerative hierarchical clustering.
TL;DR: Newton, "global," and column-oriented algorithms, and options for initial guesses, self-preconditioning, and dropping strategies are discussed, and some limited theoretical results on the properties and convergence of approximate inverses are derived.
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.
TL;DR: This work develops 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.
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.
TL;DR: Property of the developed methods are further discussed, including rotational invariance and stability under the influence of anisotropy, and numerical examples that demonstrate strengths and weaknesses of the methods are presented.
Abstract: A companion paper, Part I [SIAM J. Sci. Comput., 19 (1998), pp. 1700--1716], presents discretization methods for control-volume formulations in two space dimensions. Properties of the developed methods are further discussed, including rotational invariance and stability under the influence of anisotropy. Strong anisotropy may have to be handled by use of stretched grid cells. Numerical examples that demonstrate strengths and weaknesses of the methods are presented. The examples focus on inhomogeneity, anisotropy, and rotational invariance.
TL;DR: This analysis quantifies the empirical observation by Jeltsch and Pohl that the convergence rate of a multisplitting algorithm depends on the overlap and proves linear convergence of the algorithm in the continuous case on an infinite time interval at a rate depending on the size of the overlap.
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.
TL;DR: The approach of Ma et al. is modified, improving the stability of the scheme and extending its range of applicability, and results indicate that such quadratures dramatically reduce the computational cost of the evaluation of integrals under certain conditions.
Abstract: Generalized Gaussian quadratures appear to have been introduced by Markov late in the last century and have been studied in great detail as a part of modern analysis. They have not been widely used as a computational tool, in part due to an absence of effective numerical schemes for their construction. Recently, a numerical scheme for the design of such quadratures was introduced by Ma et al.; numerical results presented in their paper indicate that such quadratures dramatically reduce the computational cost of the evaluation of integrals under certain conditions. In this paper, we modify their approach, improving the stability of the scheme and extending its range of applicability. The performance of the method is illustrated with several numerical examples.
TL;DR: An expanded mixed finite element method for solving second-order elliptic partial differential equations on geometrically general domains and is shown to be as accurate as the standard mixed method for a large class of smooth meshes.
Abstract: We present an expanded mixed finite element method for solving second-order elliptic partial differential equations on geometrically general domains. For the lowest-order Raviart--Thomas approximating spaces, we use quadrature rules to reduce the method to cell-centered finite differences, possibly enhanced with some face-centered pressures. This substantially reduces the computational complexity of the problem to a symmetric, positive definite system for essentially only as many unknowns as elements. Our new method handles general shape elements (triangles, quadrilaterals, and hexahedra) and full tensor coefficients, while the standard mixed formulation reduces to finite differences only in special cases with rectangular elements. As in other mixed methods, we maintain the local approximation of the divergence (i.e., local mass conservation). In contrast, Galerkin finite element methods facilitate general element shapes at the cost of achieving only global mass conservation. Our method is shown to be as accurate as the standard mixed method for a large class of smooth meshes. On nonsmooth meshes or with nonsmooth coefficients one can add Lagrange multiplier pressure unknowns on certain element edges or faces. This enhanced cell-centered procedure recovers full accuracy, with little additional cost if the coefficients or mesh geometry are piecewise smooth. Theoretical error estimates and numerical examples are given, illustrating the accuracy and efficiency of the methods.
TL;DR: 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 to construct bounds of the convergence rate of the GMRES method with respect to an energy norm.
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.
TL;DR: The rational Krylov algorithm computes eigenvalues and eigenvectors of a regular not necessarily symmetric matrix pencil and is a generalization of the shifted and inverted Arnoldi algorithm.
Abstract: The rational Krylov algorithm computes eigenvalues and eigenvectors of a regular not necessarily symmetric matrix pencil. It is a generalization of the shifted and inverted Arnoldi algorithm, where several factorizations with different shifts are used in one run. It computes an orthogonal basis and a small Hessenberg pencil. The eigensolution of the Hessenberg pencil approximates the solution of the original pencil. Different types of Ritz values and harmonic Ritz values are described and compared. Periodical purging of uninteresting directions reduces the size of the basis and makes it possible to compute many linearly independent eigenvectors and principal vectors of pencils with multiple eigenvalues. Relations to iterative methods are established.
Results are reported for two large test examples. One is a symmetric pencil coming from a finite element approximation of a membrane; the other is a nonsymmetric matrix modeling an idealized aircraft stability problem.
TL;DR: Two new methods are described for determining preconditioners from spectral information gathered by the Arnoldi process during iterations by the restarted GMRES algorithm to determine an invariant subspace of the matrix A associated with eigenvalues close to the origin and to move these eigen values so that a higher rate of convergence of the iterative methods is achieved.
Abstract: The restarted GMRES algorithm proposed by Saad and Schultz [SIAM J. Sci. Statist. Comput., 7 (1986), pp. 856--869] is one of the most popular iterative methods for the solution of large linear systems of equations Ax=b with a nonsymmetric and sparse matrix. This algorithm is particularly attractive when a good preconditioner is available. The present paper describes two new methods for determining preconditioners from spectral information gathered by the Arnoldi process during iterations by the restarted GMRES algorithm. These methods seek to determine an invariant subspace of the matrix A associated with eigenvalues close to the origin and to move these eigenvalues so that a higher rate of convergence of the iterative methods is achieved.
TL;DR: A hybrid Newton--Picard scheme based on the shooting method is derived, which in its simplest form is the recursive projection method (RPM) of Shroff and Keller and is used to compute and determine the stability of both stable and unstable periodic orbits.
Abstract: This paper is concerned with the efficient computation of periodic orbits in large-scale dynamical systems that arise after spatial discretization of partial differential equations (PDEs). A hybrid Newton--Picard scheme based on the shooting method is derived, which in its simplest form is the recursive projection method (RPM) of Shroff and Keller [SIAM J. Numer. Anal., 30 (1993), pp. 1099--1120] and is used to compute and determine the stability of both stable and unstable periodic orbits. The number of time integrations needed to obtain a solution is shown to be determined only by the system's dynamics. This contrasts with traditional approaches based on Newton's method, for which the number of time integrations grows with the order of the spatial discretization. Two test examples are given to show the performance of the methods and to illustrate various theoretical points.
TL;DR: This companion paper reports on the application of the 2D code to several nontrivial problems---nonlinear arsenic diffusion in the manufacture of semiconductors, the drift-diffusion equations for semiconductor device simulation, the Buckley--Leverett black oil equations for reservoir simulation, and the motion of surfaces by mean curvature.
Abstract: In part I the authors reported on the design of a robust and versatile gradient-weighted moving finite element (GWMFE) code in one dimension and on its application to a variety of PDEs and PDE systems. This companion paper does the same for the two-dimensional (2D) case. These moving node methods are especially suited to problems which develop sharp moving fronts, especially problems where one needs to resolve the fine-scale structure of the fronts. The many potential pitfalls in the design of GWMFE codes and the special features of the implicit one-dimensional (1D) and 2D codes which contribute to their robustness and efficiency are discussed at length in part I; this paper concentrates on issues unique to the 2D case. Brief explanations are given of the variational interpretation of GWMFE, the geometrical-mechanical interpretation, simplified regularization terms, and the treatment of systems. A catalog of inner products which occur in GWMFE is given, with particular attention paid to those involving second-order operators. After presenting an example of the 2D phenomenon of grid collapse and discussing the need for long-time regularization, the paper reports on the application of the 2D code to several nontrivial problems---nonlinear arsenic diffusion in the manufacture of semiconductors, the drift-diffusion equations for semiconductor device simulation, the Buckley--Leverett black oil equations for reservoir simulation, and the motion of surfaces by mean curvature.
TL;DR: It is demonstrated 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.
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.
TL;DR: A method of dividing an irregular mesh into equal-sized pieces with few interconnecting edges is investigated, based on theoretical work of Miller, Teng, Thurston, and Vavasis, who showed that certain classes of "well-shaped" finite-element meshes have good separators.
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.
TL;DR: 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 is described, which requires a minimum of additional information from the user.
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.
TL;DR: An approach to the reordering problem that produces significantly better orderings than prior methods is described, a hybrid of nested dissection and minimum degree ordering, and combines an assortment of different algorithmic advances.
Abstract: When performing sparse matrix factorization, the ordering of matrix rows and columns has a dramatic impact on the factorization time. This paper describes an approach to the reordering problem that produces significantly better orderings than prior methods. The algorithm is a hybrid of nested dissection and minimum degree ordering, and combines an assortment of different algorithmic advances. New or improved algorithms are described for graph compression, multilevel partitioning, and separator improvement. When these techniques are combined, the resulting orderings average 39% better than minimum degree over a suite of test matrices, while requiring roughly 2.7 times the run time of Liu's multiple minimum degree.
TL;DR: This work develops a variant of the restarted GMRES method exhibiting the same advantage and investigates its convergence for positive real matrices in some detail and applies it to speed up "multiple masses" calculations arising in lattice gauge computations in quantum chromodynamics, one of the most time-consuming supercomputer applications.
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.
TL;DR: It is proved that thick restarted, nonpreconditioned Davidson is equivalent to the implicitly restarted Arnoldi and motivates the development of a dynamic thick restarting scheme for the symmetric case, which can be used in both Davidson and implicit restarting Arnoldi.
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.
TL;DR: 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 and is initiated in this work.
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.