Showing papers in "Siam Journal on Scientific and Statistical Computing in 1989"

VODE: a variable-coefficient ODE solver

Abstract: VODE is a new initial value ODE solver for stiff and nonstiff systems. It uses variable-coefficient Adams-Moulton and Backward Differentiation Formula (BDF) methods in Nordsieck form, as taken from the older solvers EPISODE and EPISODEB, treating the Jacobian as full or banded. Unlike the older codes, VODE has a highly flexible user interface that is nearly identical to that of the ODEPACK solver LSODE.In the process, several algorithmic improvements have been made in VODE, aside from the new user interface. First, a change in stepsize and/or order that is decided upon at the end of one successful step is not implemented until the start of the next step, so that interpolations performed between steps use the more correct data. Second, a new algorithm for setting the initial stepsize has been included, which iterates briefly to estimate the required second derivative vector. Efficiency is often greatly enhanced by an added algorithm for saving and reusing the Jacobian matrix J, as it occurs in the Newton m...

CGS, A Fast Lanczos-Type Solver for Nonsymmetric Linear systems

Abstract: A Lanczos-type method is presented for nonsymmetric sparse linear systems as arising from discretisations of elliptic partial differential equations. The method is based on a polynomial variant of the conjugate gradients algorithm. Although related to the so-called bi-conjugate gradients (Bi-CG) algorithm, it does not involve adjoint matrix-vector multiplications, and the expected convergence rate is about twice that of the Bi-CG algorithm. Numerical comparison is made with other solvers, testing the method on a family of convection diffusion equations, on various grids, and with the use of two different preconditioning methods. Upwind as well as central differencing is used in the experiments.

A Storage-Efficient $WY$ Representation for Products of Householder Transformations

Abstract: A product $Q = P_{1} \cdots P_{r}$ of m-by-m Householder matrices can be written in the form $Q = I + WY^{T}$ where W and Y are each m-by-r. This is called the WY representation of Q. It is of interest when implementing Householder techniques in high-performance computing environments that ``like'''' matrix-matrix multiplication. In this note we describe a storage efficient way to implement the WY representation. In particular, we show how the matrix Q can be expressed in the form $Q = I + YTY^{T}$ where $Y \epsilon R^{mxr}$ and $T \epsilon R^{rxr}$ with T upper triangular. Usually r >> m and so this ``compact'''' WY representation requires less storage. When compared with the recent block-reflector strategy proposed by Schreiber and Parlett the new technique still has a storage advantage and involves a comparable amount of work.

Krylov Subspace Methods on Supercomputers

Abstract: This paper presents a short survey of recent research on Krylov subspace methods with emphasis on implementation on vector and parallel computers. Conjugate gradient methods have proven very useful on traditional scalar computers, and their popularity is likely to increase as three-dimensional models gain importance. A conservative approach to derive effective iterative techniques for supercomputers has been to find efficient parallel/vector implementations of the standard algorithms. The main source of difficulty in the incomplete factorization preconditionings is in the solution of the triangular systems at each step. A few approaches consisting of implementing efficient forward and backward triangular solutions are described in detail. Then polynomial preconditioning as an alternative to standard incomplete factorization techniques is discussed. Another efficient approach is to reorder the equations so as to improve the structure of the matrix to achieve better parallelism or vectorization. An overview of these ideas and others is given in this article, as well as an attempt to comment on their effectiveness or potential for different types of architectures.

Toeplitz equations by conjugate gradients with circulant preconditioner

Abstract: This paper studies the solution of symmetric positive definite Toeplitz systems $Ax = b$ by the preconditioned conjugate gradient method. The preconditioner is a circulant matrix C that copies the middle diagonals of A, and each iteration uses the Fast Fourier Transform. Convergence is governed by the eigenvalues of $C^{ - 1} A$–a Toeplitz-circulant eigenvalue problem—and it is fast if those eigenvalues are clustered. The limiting behavior of the eigenvalues is found as the dimension increases, and it is proved that they cluster around $\lambda = 1$. For a wide class of problems the error after q conjugate gradient steps decreases as $r^{q^2 } $.

The Reduced Basis Method for Incompressible Viscous Flow Calculations

Abstract: The reduced basis method is a type of reduction method that can be used to solve large systems of nonlinear equations involving a parameter. In this work, the method is used in conjunction with a standard continuation technique to approximate the solution curve for the nonlinear equations resulting from discretizing the Navier–Stokes equations by finite–element methods. This paper demonstrates that the reduced basis method can be implemented to approximate efficiently solutions to incompressible viscous flows. Choices of basis vectors, issues concerning the implementation of the method, and numerical calculations are discussed. Two fluid flow calculations are considered, the driven cavity problem and flow over a forward facing step.

Three-dimensional triangulations from local transformations

Abstract: A new algorithm is presented that uses a local transformation procedure to construct a triangulation of a set of n three-dimensional points that is pseudo-locally optimal with respect to the sphere criterion It is conjectured that this algorithm always constructs a Delaunay triangulation, and this conjecture is supported with experimental results The empirical time complexity of this algorithm is $O(n^{{4 / 3}} )$ for sets of random points, which compares well with existing algorithms for constructing a three-dimensional Delaunay triangulation Also presented is a modification of this algorithm for the case that local optimality is based on the max-min solid angle criterion

On the Robustness of ILU Smoothing

Abstract: In the present paper, a detailed analysis of a multigrid method with an ILU smoother applied to a singularly perturbed problem is given. Based on the analysis of a simple anisotropic model problem, a variant of the usual incomplete LU factorization is introduced, which is especially suited as robust smoother. For this variant a detailed analysis and a proof of robustness is given. Furthermore, some contradictions between the smoothing rates predicted by local Fourier analysis and the practically observed convergence factors are explained (see [W. Hackbusch, Multi-grid Methods and Applications, Springer-Verlag, Berlin, Heidelberg, 1985; R. Kettler, “Analysis and comparison of relaxation schemes in robust multi-grid and preconditioned conjugate gradient methods,” in Multi-grid Methods, Lecture Notes in Math. 960, Springer-Verlag, Berlin, 1982; C. A. Thole, Beitrage zur Fourieranalyse von Mehrgitterver fahren, Diplomarbeit, Universitat Bonn, 1983]. The theoretical results are confirmed by numerical tests.

High Performance Preconditioning

Abstract: The discretization of second-order elliptic partial differential equations over three-dimensional rectangular regions, in general, leads to very large sparse linear systems. Because of their huge order and their sparseness, these systems can only be solved by iterative methods using powerful computers, e.g., vector supercomputers. Most of those methods are only attractive when used in combination with a so-called preconditioning matrix. Unfortunately, the more effective preconditioners, such as successive over-relaxation and incomplete decompositions, do not perform very well on most vector computers if used in a straightforward manner. In this paper it is shown how a rather high performance can be achieved for these preconditioners.

On optimal interpolation triangle incidences

Abstract: The problem of determining optimal incidences for triangulating a given set of vertices for the model problem of interpolating a convex quadratic surface by piecewise linear functions is studied An exact expression for the maximum error is derived, and the optimality criterion is minimization of the maximum error The optimal incidences are shown to be derivable from an associated Delaunay triangulation and hence are computable in $O(N\log N)$ time for N vertices

On parallel Jacobi orderings

Abstract: This paper presents a systematic description of five well-known Jacobi orderings, using the tool of a caterpillar track. The following surprising result is obtained: under certain assumptions, the orderings are all mathematically equivalent. Three new caterpillar-track orderings are then derived, but they, in turn, are shown to be equivalent to the known schemes.

A fast algorithm for reordering sparse matrices for parallel factorization

Abstract: Jess and Kees [IEEE Trans. Comput., C-31 (1982), pp. 231–239] introduced a method for ordering a sparse symmetric matrix A for efficient parallel factorization. The parallel ordering is computed in two steps. First, the matrix A is ordered by some fill-reducing ordering. Second, a parallel ordering of A is computed from the filled graph that results from symbolically factoring A using the initial fill-reducing ordering. Among all orderings whose fill lies in the filled graph, this parallel ordering achieves the minimum number of parallel steps in the factorization of A. Jess and Kees did not specify the implementation details of an algorithm for either step of this scheme. Liu and Mirzaian [SIAM J. Discrete Math., 2 (1989), pp. 100–107] designed an algorithm implementing the second step, but it has time and space requirements higher than the cost of computing common fill-reducing orderings.A new fast algorithm that implements the parallel ordering step by exploiting the clique tree representation of a chordal graph is presented. The cost of the parallel ordering step is reduced well below that of the fill-reducing step. This algorithm has time and space complexity linear in the number of compressed subscripts for L, i.e., the sum of the sizes of the maximal cliques of the filled graph. Running times nearly identical to Liu's heuristic composite rotations algorithm, which approximates the minimum number of parallel steps, are demonstrated empirically.

Mathematical structure of compositional reservoir simulation

Abstract: In this paper multicomponent two-phase isothermal fluid flow in petroleum reservoirs is described. The fluid-flow model consists of component conservation equations, Darcy's law for the volumetric flow rates, balance between the fluid volume and the rock void, and the conditions of thermodynamic equilibrium that determine the distribution of the chemical components into phases. Thermodynamic equilibrium is described by means of a mathematical model for the chemical potentials of each component in each phase of the fluid. The flow equations are manipulated to form a pressure equation and a modified component-conservation equation; these form the basis for the sequential method. It is shown that the pressure equation is parabolic under reasonable assumptions on the thermodynamic equilibrium model, and that the component-conservation equations are hyperbolic in the absence of diffusive forces such as capillary pressure and mixing. A numerical method based on the sequential formulation of the flow equations is outlined and used to illustrate the kinds of flow behavior that occur during miscible gas injection.

A New Method for Solving Triangular Systems on Distributed Memory Message-Passing Multiprocessors

Abstract: Efficient triangular solvers for use on message passing multiprocessors are required, in several contexts, under the assumption that the matrix is distributed by columns (or rows) in a wrap fashion. In this paper we describe a new efficient parallel triangular solver for this problem. This new algorithm is based on the previous method of Li and Coleman [1986] but is considerably more efficient when $\frac{n}{p}$ is relatively modest, where $p$ is the number of processors and $n$ is the problem dimension. A useful theoretical analysis is provided as well as extensive numerical results obtained on an Intel iPSC with $p \leq 128$.

A one-sided Jacobi algorithm for computing the singular value decomposition on avector computer

Abstract: An old algorithm for computing the singular value decomposition, which was first mentioned by Hestenes [SIAM J. Appl. Math., 6 (1958), pp. 51–90], has gained renewed interest because of its properties of parallelism and vectorizability. Some computational modifications are given and a comparison with the well-known Golub–Reinsch algorithm is made. Comparative experiments on a CYBER 205 are reported.

Distributed orthogonal factorization: givens and householder algorithms

Abstract: Several algorithms for orthogonal factorization on distributed memory multiprocessors are designed and implemented. Two of the algorithms employ Householder transformations, a third is based on Givens rotations, and a fourth hybrid algorithm uses Householder transformations and Givens rotations in different phases.The arithmetic and communication complexities of the algorithms are analyzed. The analyses show that the sequential arithmetic terms play a more important role than the communication terms in determining the running times and efficiencies of these algorithms. The hybrid algorithm is the fastest algorithm overall, since its arithmetic cost is lower than the Householder algorithms and its communication cost does not increase with the column length of the matrix. The observed execution times of the implementations on an iPSC-286 agree quite well with the complexity analyses. It is also shown that the efficiencies can be approximated using only the arithmetic costs of the algorithms.

Computational structure of the N-body problem

Abstract: This work considers tree algorithms for the N-body problem where the number of particles is on the order of a million. The main concern of this work is the organization and performance of these computations on parallel computers.This work introduces a formulation of the N-body problem as a set of recursive equations based on a few elementary functions. It is shown that both the algorithm of Barnes–Hut and that of Greengard–Rokhlin satisfy these equations using different elementary functions. The recursive formulation leads directly to a computational structure in the form of a pyramid-like graph, where each vertex is a process, and each arc a communication link.The pyramid is mapped to three different processor configurations: (1) a pyramid of processors corresponding to the processes pyramid graph; (2) a hypercube of processors, e.g., a connection-machine-like architecture; and (3) a rather small array, e.g., $2 \times 2 \times 2$, of processors faster than the ones considered in (1) and (2) above.The main conclusion is that simulations of this size can be performed on any of the three architectures in reasonable time. Approximately 24 seconds per timestep is the estimate for a million equally distributed particles using the Greengard-Rokhlin algorithm on the CM-2 connection machine. The smaller array of processors is quite competitive in performance.

A Fast Algorithm for the Numerical Evaluation of Conformal Mappings

Abstract: An algorithm is presented for the construction of conformal mappings from arbitrary simply connected regions in the complex plane onto the unit disk. The algorithm is based on a combination of the Kerzman–Stein integral equation (see [Math. Anal, 236 (1978), pp. 85–93]) and the Fast Multipole Method for the evaluation of Cauchy-type integrals (see [V. Rokhlin, J. Comput. Phys., 60 (1985), pp. 187–207], [L. Greengard and V. Rokhlin, J. Comput. Phys., 73 (1987), pp. 325–348], [J. Carrier, L. Greengard, and V. Rokhlin, SIAM J. Sci. Statist. Comput., 9 (1988), pp. 669–686], [L. F. Greengard, Ph.D. thesis, Department of Computer Science, Yale University, New Haven, CT, 1987]). Previously published methods for the construction of conformal mappings via the Kerzman–Stein equation have an asymptotic CPU time estimate of the order $O(n^2 )$, where n is the number of nodes in the discretization of the boundary of the region being mapped. The method presented here has an estimate of the order $O(n)$, making it an approach of choice in many situations. The performance of the algorithm is illustrated by several numerical examples.

Practical Runge–Kutta Processes

Abstract: The development of embedded Runge–Kutta and Runge–Kutta–Nystrom formulae subject to various criteria is reviewed. An important criterion concerns the cost of achieving a particular global error in the numerical solution. By consideration of local truncation errors in the two formulae of an embedded pair, it is possible to produce a good process. Another criterion involves the provision of continuous solutions. Such a requirement can be at odds with the previous one of basic cost-effectiveness. However, it seems important to provide dense output without excessive cost in new function evaluations. Special RK/RKN formulae are preferable for practical global error estimation using the Zadunaisky pseudo-problem or related technique of solving for the error estimate. Two-term error estimation can be achieved and the pseudo-problem can be based on dense output values.

To overlap or not to overlap: a note on a domain decomposition method for elliptic problems

Abstract: More than one hundred years ago, H.A. Schwarz introduced a domain decomposition method in which the original elliptic equation is solved on overlapping subregions, one after another, in an iterative process. A few years ago, Chan and Resasco, introduced a method that they classified as a domain decomposition method using nonoverlapping subdomains. In this note, it is shown that their method is an accelerated version of the classical method. It is also shown that the error propagation operator of the method can be expressed in terms of Schur complements of certain stiffness matrices and that techniques previously developed for the study of iterative substructuring algorithms can be used to derive estimates on the rate of convergence.

An accelerated sequential algorithm for producing D -optimal designs

Abstract: An extremely efficient algorithm is presented for generating D-optimal designs. It is shown that the underlying set function of the sequential design algorithm is submodular and thus the so-called accelerated greedy algorithm may be applied to this problem. While the new algorithm's statistical basis is identical to that of the Wynn and Federov algorithm, it requires significantly fewer function evaluations. The new algorithm is particularly useful when no prior information concerning the structure of the optimal design is available.

A Dual Method for Certain Positive Semidefinite Quadratic Programming Problems

Abstract: This paper presents a dual active set method for minimizing a sum of piecewise linear functions and a strictly convex quadratic function, subject to linear constraints. It may be used for direction finding in nondifferentiable optimization algorithms and for solving exact penalty formulations of (possibly inconsistent) strictly convex quadratic programming problems. An efficient implementation is described extending the Goldfarb and Idnani algorithm, which includes Powell's refinements. Numerical results indicate excellent accuracy of the implementation.

On symmetric schemes and differential-algebraic equations

Abstract: An example is given that demonstrates a potential risk in using symmetric difference schemes for initial value differential-algebraic equations (DAEs) or for very stiff ordinary differential equations (ODEs). The basic difficulty is that the stability of the scheme is controlled by the stability of an auxiliary (ghost) ODE problem that is not necessarily stable even when the given problem is.The stability of symmetric schemes is better understood in the context of boundary value problems. In this context, such schemes are more naturally applied as well. For initial value problems, better alternatives may exist. A computational algorithm is proposed for boundary value index-1 DAEs.

Computation of hyperbolic equations on complicated domains with patched and overset Chebyshev grids

Abstract: A patched/overset grid strategy that allows hyperbolic partial differential equations to be solved in complicated geometries with Chebyshev spectral collocation is presented. Unlike similar approaches used with finite difference methods, the spectral method defines how interpolations should be carried out. The result is an approach that retains spectral accuracy. Numerical experiments are performed on both one and two-dimensional examples.

Stability properties of backward euler multirate formulas

Abstract: Stability properties of multirate formulas cannot be analyzed by a scalar test equation but require at least one equation for each different steplength. This paper generalizes the concept of absolute stability and A-stability for backward Euler multirate formulas. Stability theorems for multirate methods with two and three different steplengths are given, while a general result for an arbitrary number of different steplengths is the topic of future research.

The Minimum Degree Ordering with Constraints

Abstract: A hybrid scheme for ordering sparse symmetric matrices is considered. It is based on a combined use of the top-down nested dissection and the bottom-up minimum degree ordering schemes. A separator set is first determined by some form of incomplete nested dissection. The minimum degree ordering is then applied subject to the constraint that the separator nodes must be ordered last. It is shown experimentally that the quality of the resulting ordering from this constrained scheme exhibits less sensitivity to the initial matrix ordering than that of the original minimum degree ordering. An important application of this approach to find orderings suitable for parallel elimination is also illustrated.

A new class of parallel algorithms for solving systems of linear equations

Abstract: In this paper a class of novel feed-forward direct methods is presented for solving nonsingular systems of linear equations. The computational complexity of these methods is in the order of an $LU$, $QR$, or $LL^t $ matrix factorization. This is also true for the complexity of their systolic implementations. Unlike the direct methods of factorization followed by backsubstitution, the systolic implementations of the novel methods do not suffer from the backsubstitution bottleneck. A numerically stable and robust method, which uses only Givens rotations as elementary operations, is included in the class.

Application of sparse matrix solvers as effective preconditioners

Abstract: In this paper the use of a new out-of-core sparse matrix package for the numerical solution of partial differential equations involving complex geometries arising from aerospace applications is discussed. The sparse matrix solver accepts contributions to the matrix elements in random order and assembles the matrix using fast sort/merge routines. Fill-in is reduced through the use of a physically based nested dissection ordering. For very large problems a drop tolerance is used during the matrix decomposition phase. The resulting incomplete factorization is an effective preconditioner for Krylov subspace methods, such as GMRES. Problems involving 200,000 unknowns routinely are solved on the Cray X-MP using 64MW of solid-state storage device (SSD).

Numerical determination of an emanating branch of Hopf bifurcation points in a two-parameter problem

Abstract: In a two-parameter problem a branch of Hopf bifurcation points can bifurcate from a branch of simple turning points of the steady state problem, at a point for which the Frechet derivative has a double eigenvalue zero with a one-dimensional nullspace. It is indicated how the origin of a branch of Hopf points can be detected during the continuation of a branch of simple turning points.Further, an augmented system of equations is presented, for which this “origin for Hopf bifurcation” is an isolated solution. If the steady state problem is described by a system of algebraic equations, Newton's method for the solution of the augmented system can be implemented very efficiently. The authors also discuss switching to the branch of Hopf points.Results are given for the one-dimensional “Brusselator” model, a system of four partial differential equations.

A Study of the Vortex Sheet Method and its Rate of Convergence

Abstract: The subject of this study is Chorin's vortex sheet method, which is used to solve the Prandtl boundary layer equations and to impose the no-slip boundary condition in the random vortex method solution of the Navier–Stokes equations. This is a particle method in which the particles carry concentrations of vorticity and undergo a random walk to approximate the diffusion of vorticity in the boundary layer. During the random walk, particles are created at the boundary in order to satisfy the no-slip boundary condition. It is proved that in each of the $L^1 $, $L^2 $, and $L^\infty $ norms the random walk and particle creation, taken together, provide a consistent approximation to the heat equation, subject to the no-slip boundary condition. Furthermore, it is shown that the truncation error is entirely due to the failure to satisfy the no-slip boundary condition exactly. It is demonstrated numerically that the method converges when it is used to model Blasius flow, and rates of convergence are established in terms of the computational parameters. The numerical study reveals that errors grow when the sheet length tends to zero much faster than the maximum sheet strength. The effectiveness of second-order time discretization, sheet tagging, and an alternative particle-creation algorithm are also examined.