Past Computational Mathematics and Applications Seminar

We consider the parallel solution of sparse linear systems of equations in a limited memory environment. A preliminary out-of core version of a sparse multifrontal code called MUMPS (MUltifrontal Massively Parallel Solver) has been developed as part of a collaboration between CERFACS, ENSEEIHT and INRIA (ENS-Lyon and Bordeaux).

We first briefly describe the current status of the out-of-core factorization phase. We then assume that the factors have been written on the hard disk during the factorization phase and we discuss the design of an efficient solution phase.Two different approaches are presented to read data from the disk, with a discussion on the advantages and the drawbacks of each one.

Our work differs and extends the work of Rothberg and Schreiber (1999) and of Rotkin and Toledo (2004) because firstly we consider a parallel out-of-core context, and secondly we also study the performance of the solve phase.

This is work on collaboration with E. Agullo, I.S Duff, A. Guermouche, J.-Y. L'Excellent, T. Slavova

A strict adherence to threshold pivoting in the direct solution of symmetric indefinite problems can result in substantially more work and storage than forecast by an sparse analysis of the symmetric problem. One way of avoiding this is to use static pivoting where the data structures and pivoting sequence generated by the analysis are respected and pivots that would otherwise be very small are replaced by a user defined quantity. This can give a stable factorization but of a perturbed matrix.

The conventional way of solving the sparse linear system is then to use iterative refinement (IR) but there are cases where this fails to converge. We will discuss the use of more robust iterative methods, namely GMRES and its variant FGMRES and their backward stability when the preconditioning is performed by HSL_M57 with a static pivot option.

Several examples under Matlab will be presented.

Radial basis functions have been used for decades for the interpolation of scattered,

high-dimensional data. Recently they have attracted interest as methods for simulating

partial differential equations as well. RBFs do not require a grid or triangulation, they

offer the possibility of spectral accuracy with local refinement, and their implementation

is very straightforward. A number of theoretical and practical breakthroughs in recent years

has improved our understanding and application of these methods, and they are currently being

tested on real-world applications in shallow water flow on the sphere and tear film evolution

in the human eye.

Joint work with Yogi Erlangga and Kees Vuik.

Shifted Laplace preconditioners have attracted considerable attention as a technique to speed up convergence of iterative solution methods for the Helmholtz equation. In this paper we present a comprehensive spectral analysis of the Helmholtz operator preconditioned with a shifted Laplacian. Our analysis is valid under general conditions. The propagating medium can be heterogeneous, and the analysis also holds for different types of damping, including a radiation condition for the boundary of the computational domain. By combining the results of the spectral analysis of the preconditioned Helmholtz operator with an upper bound on the GMRES-residual norm we are able to provide an optimal value for the shift, and to explain the mesh-depency of the convergence of GMRES preconditioned with a shifted Laplacian. We illustrate our results with a seismic test problem.

Many large-scale optimization problems arise in the context of the discretization of infinite dimensional applications. In such cases, the description of the finite-dimensional problem is not unique, but depends on the discretization used, resulting in a natural multi-level description. How can such a problem structure be exploited, in discretized problems or more generally? The talk will focus on discussing this issue in the context of unconstrained optimization and in relation with the classical multigrid approach to elliptic systems of partial differential equations. Both theoretical convergence properties of special purpose algorithms and their numerical performances will be discussed. Perspectives will also be given.

Collaboration with S. Gratton, A. Sartenaer and M. Weber.

The talk starts with a general introduction of the convex

quadratic semidefinite programming problem (QSDP), followed by a number of

interesting examples arising from finance, statistics and computer sciences.

We then discuss the concept of primal nondegeneracy for QSDP and show that

some QSDPs are nondegenerate and others are not even under the linear

independence assumption. The talk then moves on to the algorithmic side by

introducing the dual approach and how it naturally leads to Newton's method,

which is quadratically convergent for degenerate problems. On the

implementation side of the Newton method, we stress that direct methods for

the linear equations in Newton's method are impossible simply because the

equations are quite large in size and dense. Our numerical experiments use

the conjugate gradient method, which works quite well for the nearest

correlation matrix problem. We also discuss difficulties for us to find

appropriate preconditioners for the linear system encountered. The talk

concludes in discussing some other available methods and some future topics.

The aim of this talk is to describe several methods for numerically approximating

the integral of a multivariate highly oscillatory function. We begin with a review

of the asymptotic and Filon-type methods developed by Iserles and Nørsett. Using a

method developed by Levin as a point of departure we will construct a new method that

uses the same information as the Filon-type method, and obtains the same asymptotic

order, while not requiring moments. This allows us to integrate over nonsimplicial

domains, and with complicated oscillators.

High-performance computing is an important tool for computational science.

Oxford University has recently decided to invest £3M in a new supercomputing

facility which is under development now. In this seminar I will give an overview

of some activities in Oxford and provide a vision for supercomputing here.

I will discuss some of the numerical analysis software and tools,

such as Distributed Matlab and indicate some of the challenges at

the intersection of numerical analysis and high-performance computing.

The aim of this talk is to give some understanding of the theory of matrix $p$'th roots (solutions to the nonlinear matrix equation $X^{p} = A$), to explain how and how not to compute roots, and to describe some applications. In particular, an application in finance will be described concerning roots of transition matrices from Markov models.

One of the oldest approach in meshless methods for PDEs is the Interpolating Moving Least Squares (IMLS) technique developed in the 1980s. Although widely accepted by users working in fields as diverse as geoinformatics and crack dynamics I shall take a fresh look at this method and ask for the equivalent difference operators which are generated implicitly. As it turns out, these operators are optimal only in trivial cases and are "strange" in general. I shall try to exploit two different approaches for the computation of these operators.

On the other hand (and very different from IMLS), Total Variation Flow (TVF) PDEs are the most recent developments in image processing and have received much attention lately. Again I shall show that they are able to generate "strange" discrete operators and that they easily can behave badly although they may be properly implemented.

We will explore connections between the structure of high-dimensional convex polytopes and information acquisition for compressible signals. A classical result in the field of convex polytopes is that if N points are distributed Gaussian iid at random in dimension n<<N, then only order (log N)^n of the points are vertices of their convex hull. Recent results show that provided n grows slowly with N, then with high probability all of the points are vertices of its convex hull. More surprisingly, a rich "neighborliness" structure emerges in the faces of the convex hull. One implication of this phenomenon is that an N-vector with k non-zeros can be recovered computationally efficiently from only n random projections with n=2e k log(N/n). Alternatively, the best k-term approximation of a signal in any basis can be recovered from 2e k log(N/n) non-adaptive measurements, which is within a log factor of the optimal rate achievable for adaptive sampling. Additional implications for randomized error correcting codes will be presented.

This work was joint with David L. Donoho.

In this talk we present different strategies for regularization of the pure Newton method

(minimization problems)and of the Gauss-Newton method (systems of nonlinear equations).

For these schemes, we prove general convergence results. We establish also the global and

local worst-case complexity bounds. It is shown that the corresponding search directions can

be computed by a standard linear algebra technique.

We present a novel enhanced finite element method for the Darcy problem starting from the non stable

continuous $P_1 / P_0$ finite element spaces enriched with multiscale functions. The method is a departure

from the standard mixed method framework used in these applications. The methods are derived in a Petrov-Galerkin

framework where both velocity and pressure trial spaces are enriched with functions based on residuals of strong

equations in each element and edge partition. The strategy leads to enhanced velocity space with an element of

the lowest order Raviart-Thomas space and to a stable weak formulation preserving local mass conservation.

Numerical tests validate the method.

Jointly with Gabriel R Barrenechea, Universidad de Concepcion &amp;

Frederic G C Valentin, LNCC

Strong horizontal gradients of density are responsible for the occurence of a large number of (often catastrophic) flows, such as katabatic winds, dust storms, pyroclastic flows and powder-snow avalanches. For a large number of applications, the overall density contrast in the flow remains small and simulations are carried in the Boussinesq limit, where density variations only appear in the body-force term. However, pyroclastic flows and powder-snow avalanches involve much larger density contrasts, which implies that the inhomogeneous Navier-Stokes equations need to be solved, along with a closure equation describing the mass diffusion. We propose a Lagrange-Galerkin numerical scheme to solve this system, and prove optimal error bounds subject to constraints on the order of the discretization and the time-stepping. Simulations of physical relevance are then shown.

An integral part of the brain is a fluid flow system that is separate from brain tissue and the cerebral blood flow system: cerebrospinal fluid (CSF) is produced near the centre of the brain, flows out and around the brain, including around the spinal cord and is absorbed primarily in a region between the brain tissue and the skull. Hydrocephalus covers a broad range of anomalous flow and pressure situations: the normal flow path can become blocked, other problems can occur which result in abnormal tissue deformation or pressure changes. This talk will describe work that treats brain tissue as a poroelastic matrix through which the CSF can move when normal flow paths are blocked, producing tissue deformation and pressure changes. We have a number of models, the simplest treating the brain and CSF flow as having spherial symmetry ranging to more complex, fully three-dimensional computations. As well as considering acute hydrocephalus, we touch on normal pressure hydrocephalus, idiopathic intracranial hypertension and simulation of an infusion test. The numerical methods used are a combination of finite difference and finite element techniques applied to an interesting set of hydro-elastic equations.

In this talk, we report on recent activities in the development of automatic differentiation tools for Matlab and CapeML, a common intermediate language for process control, and highlight some recent AD applications. Lastly, we show the potential for parallelisation created by AD and comment on the impact on scientific computing due to emerging multicore chips which are providing substantial thread-based parallelism in a "pizza box" form factor.

We consider the solution of sequences of linear systems by preconditioned iterative methods. Such systems arise, for example, in applications such as CFD and structural mechanics. In some cases it is important to avoid the recomputation of preconditioners for subsequent systems. We propose an algebraic strategy that replaces new preconditioners by old preconditioners with simple updates. Efficiency of the new strategy, which generalizes the approach of Benzi and Bertaccini, is demonstrated using numerical experiments.

This talk presents results of joint work with Jurjen Duintjer Tebbens.

The use of discrete differential forms in the construction of finite element discretisations of the Sobolev spaces H^s, H(div) and H(curl) is now routinely applied by numerical analysts and engineers alike. However, little attention has been paid to the conditioning of the resulting stiffness matrices, particularly in the case of the non-uniform meshes that arise when adaptive refinement algorithms are used. We study this issue and show that the matrices are generally rather poorly conditioned. Typically, diagonal scaling is applied (often unwittingly) as a preconditioner. However, whereas diagonal scaling removes the effect of the mesh non-uniformity in the case of Sobolev spaces H^s, we show this is not so in the case of the spaces H(curl) and H(div). We trace the reason behind this difference, and give a simple remedy for curing the problem.

We propose a new parallel domain decomposition algorithm to solve symmetric linear systems of equations derived from the discretization of PDEs on general unstructured grids of triangles or tetrahedra. The algorithm is based on a single-level Schwarz alternating procedure and a modified conjugate gradient solver. A single layer of overlap has been adopted in order to simplify the data-structure and minimize the overhead. This approach makes the global convergence rate vary slightly with the number of domains and the algorithm becomes highly scalable. The algorithm has been implemented in Fortran 90 using MPI and hence portable to different architectures. Numerical experiments have been carried out on a SunFire 15K parallel computer and have been shown superlinear performance in some cases.

Non-normal matrices can be tiresome; some eigenvalues may be phlegmatic while others may be volatile. Computable error bounds are rarely used in such computations. We offer a way to proceed. Let (e,q,p') be an approximate eigentriple for non-normal B. Form column and row residuals r = Bq - qe and s' = p'B - ep'. We establish the relation between the smallest perturbation E, in both spectral and Frobenius norms, that makes the approximations correct and the norms of r and s'. Our results extend to the case when q and p are tall thin matrices and e is a small square matrix. Now regard B as a perturbation of B-E to obtain a (first order) bound on the error in e as a product of ||E|| and the condition number of e, namely (||q|| ||p'||)/|p'q|.