Forthcoming events in this series
Spectral methods for PDEs in complex geometry
Abstract
Spectral methods are a class of methods for solving PDEs numerically.
If the solution is analytic, it is known that these methods converge
exponentially quickly as a function of the number of terms used.
The basic spectral method only works in regular geometry (rectangles/disks).
A huge amount of effort has gone into extending it to
domains with a complicated geometry. Domain decomposition/spectral
element methods partition the domain into subdomains on which the PDE
can be solved (after transforming each subdomain into a
regular one). We take the dual approach - embedding the domain into
a larger regular domain - known as the fictitious domain method or
domain embedding. This method is extremely simple to implement and
the time complexity is almost the same as that for solving the PDE
on the larger regular domain. We demonstrate exponential convergence
for Dirichlet, Neumann and nonlinear problems. Time permitting, we
shall discuss extension of this technique to PDEs with discontinuous
coefficients.
Wave propagation in 1-d flexible multi-structures
Abstract
In this talk we will mainly analyze the vibrations of a simplified 1-d model for a multi-body structure consisting of a finite number of flexible strings distributed along a planar graph. In particular we shall analyze how solutions propagate along the graph as time evolves. The problem of the observation of waves is a natural framework to analyze this issue. Roughly, the question can be formulated as follows: Can we obtain complete information on the vibrations by making measurements in one single extreme of the network? This formulation is relevant both in the context of control and inverse problems.
Using the Fourier development of solutions and techniques of Nonharmonic Fourier Analysis, we give spectral conditions that guarantee the observability property to hold in any time larger than twice the total lengths of the network in a suitable Hilbert that can be characterized in terms of Fourier series by means of properly chosen weights. When the network graph is a tree these weights can be identified.
Once this is done these results can be transferred to other models as the Schroedinger, heat or beam-type equations.
This lecture is based on results obtained in collaboration with Rene Dager.
Matrix Computations and the secular equation
Abstract
The "secular equation" is a special way of expressing eigenvalue
problems in a variety of applications. We describe the secular
equation for several problems, viz eigenvector problems with a linear
constraint on the eigenvector and the solution of eigenvalue problems
where the given matrix has been modified by a rank one matrix. Next we
show how the secular equation can be approximated by use of the
Lanczos algorithm. Finally, we discuss numerical methods for solving
the approximate secular equation.
Multigrid solvers for quantum dynamics - a first look
Abstract
The numerical study of lattice quantum chromodynamics (QCD) is an attempt to extract predictions about the world around us from the standard model of physics. Worldwide, there are several large collaborations on lattice QCD methods, with terascale computing power dedicated to these problems. Central to the computation in lattice QCD is the inversion of a series of fermion matrices, representing the interaction of quarks on a four-dimensional space-time lattice. In practical computation, this inversion may be approximated based on the solution of a set of linear systems.
In this talk, I will present a basic description of the linear algebra problems in lattice QCD and why we believe that multigrid methods are well-suited to effectively solving them. While multigrid methods are known to be efficient solution techniques for many operators, those arising in lattice QCD offer new challenges, not easily handled by classical multigrid and algebraic multigrid approaches. The role of adaptive multigrid techniques in addressing the fermion matrices will be highlighted, along with preliminary results for several model problems.
Micro-processor design: Theoretical physics meets high-volume manufacturing
New developments in LAPACKJ and ScaLAPACK
Abstract
In this talk we shall be looking at recent and forthcoming developments in the widely used LAPACK and ScaLAPACK numerical linear algebra libraries.
Improvements include the following: Faster algorithms, better numerical methods, memory hierarchy optimizations, parallelism, and automatic performance tuning to accommodate new architectures; more accurate algorithms, and the use of extra precision; expanded functionality, including updating and downdating and new eigenproblems; putting more of LAPACK into ScaLAPACK; and improved ease of use with friendlier interfaces in multiple languages. To accomplish these goals we are also relying on better software engineering techniques and contributions from collaborators at many institutions.
After an overview, this talk will highlight new more accurate algorithms; faster algorithms, including those for pivoted Cholesky and updating of factorizations; and hybrid data formats.
This is joint work with Jim Demmel, Jack Dongarra and the LAPACK/ScaLAPACK team.
Linear and nonlinear semidefinite programs in structural optimization
Abstract
Several formulations of structural optimization problems based on linear and nonlinear semidefinite programming will be presented. SDP allows us to formulate and solve problems with difficult constraints that could hardly be solved before. We will show that sometimes it is advantageous to prefer a nonlinear formulation to a linear one. All the presented formulations result in large-scale sparse (nonlinear) SDPs. In the second part of the talk we will show how these problems can be solved by our augmented Lagrangian code PENNON. Numerical examples will illustrate the talk.
Joint work with Michael Stingl.
Recent advances in the implementation of exponential integrators
Analysis of a two-level time-integration method for ordinary and partial differential equations
Parallel sparse multifrontal solver in a limited memory environment
Abstract
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
GMRES preconditioned by a perturbed LDL^T decomposition with static pivoting
Abstract
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 function methods for meshless PDE computation
Abstract
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.
Spectral analysis of the discrete Helmholtz operator preconditioned with a shifted laplacian
Abstract
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.
Multilevel optimization and multigrid methods
Abstract
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.
Numerical prediction of multiphase flows of immiscible fluids
Convex quadratic semi-definite programming problem: algorithms and applications
Abstract
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.
Multivariate highly oscillatory integration
Abstract
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.
Supercomputing at Oxford
Abstract
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.
Matric roots: theory, computation and applications
Abstract
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.
Strange discrete operators - A tour concerning meshless methods and image processing
Abstract
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.
The surprising structure of Gaussian point clouds and its implications for signal processing
Abstract
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.