Search

Numerical algorithms have two kinds of costs: arithmetic and communication, by which we mean either moving data between levels of a memory hierarchy (in the sequential case) or over a network connecting processors (in the parallel case). Communication costs often dominate arithmetic costs, so it is of interest to design algorithms minimizing communication. In this paper we first extend known lower bounds on the communication cost (both for bandwidth and for latency) of conventional (O(n3)) matrix multiplication to Cholesky factorization, which is used for solving dense symmetric positive definite linear systems. Second, we compare the costs of various Cholesky decomposition implementations to these lower bounds and identify the algorithms and data structures that attain them. In the sequential case, we consider both the two-level and hierarchical memory models. Combined with prior results in [J. Demmel et al., Communication-optimal Parallel and Sequential QR and LU Factorizations, Technical report EECS-2008-89, University of California, Berkeley, CA, 2008], [J. Demmel et al., Implementing Communication-optimal Parallel and Sequential QR and LU Factorizations, SIAM. J. Sci. Comp., submitted], and [J. Demmel, L. Grigori, and H. Xiang, Communication-avoiding Gaussian Elimination, Proceedings of the 2008 ACM/IEEE Conference on Supercomputing, 2008] this gives a set of communication-optimal algorithms for O(n³) implementations of the three basic factorizations of dense linear algebra: LU with pivoting, QR, and Cholesky. But it goes beyond this prior work on sequential LU by optimizing communication for any number of levels of memory hierarchy. [ABSTRACT FROM AUTHOR]

We develop in this paper some efficient algorithms which are essential to implementations of spectral methods on the sparse grid by Smolyak's construction based on a nested quadrature. More precisely, we develop a fast algorithm for the discrete transform between the values at the sparse grid and the coefficients of expansion in a hierarchical basis; and by using the aforementioned fast transform, we construct two very efficient sparse spectral-Galerkin methods for a model elliptic equation. In particular, the Chebyshev Legendre-Galerkin method leads to a sparse matrix with a much lower number of nonzero elements than that of low-order sparse grid methods based on finite elements or wavelets, and can be efficiently solved by a suitable sparse solver. Ample numerical results are presented to demonstrate the efficiency and accuracy of our algorithms. [ABSTRACT FROM AUTHOR]

In this work we study a dynamical system with a complicated nonlinearity, which describes oscillation of a turbine rotor, and give an algorithm to compute Melnikov functions for analysis of its chaotic behavior. We first derive the rotor model whose nonlinear term brings difficulties to investigating the distribution and qualitative properties of its equilibria. This nonlinear model provides a typical example of a system for which the homoclinic and heteroclinic orbits cannot be analytically determined. In order to apply Melnikov's method to make clear the underlying conditions for chaotic motion, we present a generic algorithm that provides a systematic procedure to compute Melnikov functions numerically. Substantial analysis is done so that the numerical approximation precision at each phase of the computation can be guaranteed. Using the algorithm developed in this paper, it is straightforward to obtain a sufficient condition for chaotic motion under damping and periodic external excitation, whenever the rotor parameters are given. [ABSTRACT FROM AUTHOR]

A successive continuation method for locating connecting orbits in parametrized systems of autonomous ODEs was considered in [Numer. Algorithms, 14 (1997), pp. 103­124]. In this paper we present an improved algorithm for locating and continuing connecting orbits, which includes a new algorithm for the continuation of invariant subspaces. The latter algorithm is of independent interest and can be used in contexts different than the present one. [ABSTRACT FROM AUTHOR]

In this work, we develop an O(N) implicit real space method in 1D and 2D for the Cahn-Hilliard (CH) and vector Cahn-Hilliard (VCH) equations, based on the method of lines transpose (MOLT) formulation. This formulation results in a semidiscrete time stepping algorithm, which we prove is gradient stable in the H-1 norm. The spatial discretization follows from dimensional splitting and an O(N) matrix-free solver, which applies fast convolution to the modified Helmholtz equation. We propose a novel factorization technique, in which fourth-order spatial derivatives are incorporated into the solver. The splitting error is included in the nonlinear fixed point iteration, resulting in a high-order, logically Cartesian (line-by-line) update. Our method is fast but not restricted to periodic boundaries like the fast Fourier transform (FFT). The basic solver is implemented using the backward Euler formulation, and we extend this to both backward difference formula (BDF) stencils, singly diagonal implicit Runge-Kutta (SDIRK), and spectral deferred correction (SDC) frameworks to achieve high orders of temporal accuracy. We demonstrate with numerical results that the CH and VCH equations maintain gradient stability in one and two spatial dimensions. We also explore time-adaptivity, so that meta-stable states and ripening events can be simulated both quickly and efficiently. [ABSTRACT FROM AUTHOR]

The integration of functions over measurable sets is a fundamental problem in computational science. When the measurable sets belong to high-dimensional spaces or the function is computationally complex, it may only be practical to estimate integrals based on weighted sums of function values from a finite collection of samples. Monte Carlo, quasi-Monte Carlo, and other (pseudo-)random schemes are common choices for determining a set of samples. These schemes are appealing for their conceptual ease and ability to circumvent, with various degrees of success, the so-called curse of dimensionality. However, convergence is often slow and described in terms of probability. We consider a general measure-theoretic interpretation of any sample based algorithm for numerically approximating an integral. A priori error bounds are proven that provide insight into defining adaptive sampling algorithms solving error optimization problems. We use these bounds to improve integral approximations for both forward and inverse problems. [ABSTRACT FROM AUTHOR]

We present a method for creating large two-dimensional graded Delaunay meshes on parallel machines. The method decouples the mesh generation procedure for each subdomain and thus eliminates the communication between processors, while the final mesh is of the same guaranteed quality as in the sequential case. This work extends our previous result on the uniform Delaunay decoupling method [L. Linardakis and N. Chrisochoides, SIAM J. Sci. Comput., 27 (2006), pp. 1394- 1423], by allowing the element size to be governed by a sizing function, or a background mesh. The graded Delaunay decoupling method can produce in parallel very large graded meshes efficiently, demonstrating high speedup and scalability. Ten billion elements can be created in five minutes with small overrefinement (about 2%). The decoupling parallelization approach is effective, using off-the-shelf, well-tested, and fine-tuned sequential mesh generation libraries without modification. [ABSTRACT FROM AUTHOR]

The problem of variational data assimilation for a nonlinear evolution model is formulated as an optimal control problem to find the initial condition function (analysis). The equation for the analysis error is derived through the errors of the input data (background and observation errors). This equation is used to show that in a nonlinear case the analysis error covariance operator can be approximated by the inverse Hessian of an auxiliary data assimilation problem which involves the tangent linear model constraints. The inverse Hessian is constructed by the quasi-Newton BFGS algorithm when solving the auxiliary data assimilation problem. A fully nonlinear ensemble procedure is developed to verify the accuracy of the proposed algorithm. Numerical examples are presented. [ABSTRACT FROM AUTHOR]

The classical overlapping Schwarz algorithm is here extended to the triangular/tetrahedral spectral element (TSEM) discretization of elliptic problems. This discretization, based on Fekete nodes, is a generalization to nontensorial elements of the tensorial Gauss-Lobatto-Legendre quadrilateral spectral elements (QSEM). The overlapping Schwarz preconditioners are based on partitioning the domain of the problem into overlapping subdomains, solving local problems on these subdomains, and solving an additional coarse problem associated with either the subdomain mesh or the spectral element mesh. The overlap size is generous, i.e., one element wide, in the TSEM case, while it is minimal or variable in the QSEM case. The results of several numerical experiments show that the convergence rate of the proposed preconditioning algorithm is independent of the number of subdomains N and the spectral degree p in case of generous overlap; otherwise it depends inversely on the overlap size. The proposed preconditioners are also robust with respect to arbitrary jumps of the coefficients of the elliptic operator across subdomains. [ABSTRACT FROM AUTHOR]