Rutherford Appleton Laboratory, nr Didcot

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.

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.

The use of preconditioned Newton-Krylov methods is in many applications mandatory for computing efficiently the solution of large nonlinear systems of equations. However, the available preconditioners are often sub-optimal, due to the changing nature of the linearized operator. This the case, for instance, for quasi-Newton methods where the Jacobian (and its preconditioner) are kept fixed at each non-linear iteration, with the rate of convergence usually degraded from quadratic to linear. Updated Jacobians, on the other hand require updated preconditioners, which may not be readily available. In this work we introduce an adaptive preconditioning technique based on the Krylov subspace information generated at previous steps in the nonlinear iteration. In particular, we use to advantage an adaptive technique suggested for restarted GMRES to enhance existing preconditioners with information about (almost) invariant subspaces constructed by GMRES at previous stages in the nonlinear iteration. We provide guidelines on the choice of invariant-subspace basis used in the construction of our preconditioner and demonstrate the improved performance on various test problems. As a useful general application we consider the case of augmented systems preconditioned by block triangular matrices based on the structure of the system matrix. We show that a sufficiently good solution involving the primal space operator allows for an efficient application of our adaptive technique restricted to the space of dual variables.

In this talk we will review classical multigrid methods and give an overview on algebraic multigrid methods, in particular the "classical" approach to AMG by Ruge and Stueben.

After that we will introduce a new class of multilevel methods. These new AMGs on one hand and exploit information based on filtering vectors and on the other hand, information about the inverse matrix is used to drive the coarsening process.

This new kind of AMG will be discussed and compared with "classical" AMG from a theoretical point of view as well as by showing some numerical examples.

Variational data assimilation techniques for optimal state estimation in very large environmental systems currently use approximate Gauss-Newton (GN) methods. The GN method solves a sequence of linear least squares problems subject to linearized system constraints. For very large systems, low resolution linear approximations to the model dynamics are used to improve the efficiency of the algorithm. We propose a new approach for deriving low order system approximations based on model reduction techniques from control theory which can be applied to unstable stochastic systems. We show how this technique can be combined with the GN method to retain the response of the dynamical system more accurately and improve the performance of the approximate GN method.

We will discuss the use of hypergraph-based methods for orderings of sparse matrices in Cholesky, LU and QR factorizations. For the Cholesky factorization case, we will investigate a recent result on pattern-wise decomposition of sparse matrices, generalize the result and develop algorithmic tools to obtain effective ordering methods. We will also see that the generalized results help us formulate the ordering problem in LU much like we do for the Cholesky case, without ever symmetrizing the given matrix $A$ as $A+A^{T}$ or $A^{T}A$. For the QR factorization case, the use of hypergraph models is fairly standard. We will nonetheless highlight the fact that the method again does not form the possibly much denser matrix $A^{T}A$. We will see comparisons of the hypergraph-based methods with the most common alternatives in all three cases.

\\

\\

This is joint work with Iain S. Duff.

Parallel sparse direct solvers are an interesting alternative to iterative methods for some classes of large sparse systems of linear equations. In the context of a parallel sparse multifrontal solver (MUMPS), we describe a new dynamic scheduling strategy aiming at balancing both the workload and the memory usage. More precisely, this hybrid approach balances the workload under memory constraints. We show that the peak of memory can be significantly reduced, while we have also improved the performance of the solver.

Then, we present preliminary work concerning a parallel out-of-core extension of the solver MUMPS, enabling to solve increasingly large simulation problems.

This is joint work with P.Amestoy, A.Guermouche, S.Pralet and E.Agullo.

Alan Turing introduced the sensitivity of the solution of a numerical problem to changes in its data as a way to measure the difficulty of solving the problem accurately. Condition numbers are now considered fundamental to sensitivity analysis. They have been used to measure the mathematical difficulty of many linear algebra problems, including linear systems, linear least-squares, and eigenvalue problems. By definition, unless exact arithmetic is used, it is expected to be difficultto accurately solve an ill-conditioned problem.

In this talk we focus on least-squares problems. After a historical overview of condition number for least-squares, we introduce two related condition numbers. The first is the partial condition number, which measures the sensitivity of a linear combination of the components of the solution. The second is related the truncated SVD solution of the problem, which is often used when the matrix is nearly rank deficient.

Throughout the talk we are interested in three types of results :closed formulas for condition numbers, sharp bounds and statistical estimates.

Current methods for globalizing Newton's Method for solving systems of nonlinear equations fall back on steps biased towards the steepest descent direction (e.g. Levenberg/Marquardt, Trust regions, Cauchy point dog-legs etc.), when there is difficulty in making progress. This can occasionally lead to very slow convergence when short steps are repeatedly taken.

This talk looks at alternative strategies based on searching curved arcs related to Davidenko trajectories. Near to manifolds on which the Jacobian matrix is singular, certain conjugate steps are also interleaved, based on identifying a Pareto optimal solution.

Preliminary limited numerical experiments indicate that this approach is very effective, with rapid and ultimately second order convergence in almost all cases. It is hoped to present more detailed numerical evidence when the talk is given. The new ideas can also be incorporated with more recent ideas such as multifilters or nonmonotonic line searches without difficulty, although it may be that there is no longer much to gain by doing this.

This seminar will be held at the Rutherford Appleton Laboratory near Didcot.

Abstract:

Numerical calculations of laminar flow in a two-dimensional channel with a sudden expansion exhibit a symmetry-breaking bifurcation at Reynolds number 40.45 when the expansion ratio is 3:1. In the experiments reported by Fearn, Mullin and Cliffe [1] there is a large perturbation to this bifurcation and the agreement with the numerical calculations is surprisingly poor. Possible reasons for this discrepancy are explored using modern techniques for uncertainty quantification.

When experimental equipment is constructed there are, inevitably, small manufacturing imperfections that can break the symmetry in the apparatus. In this work we considered a simple model for these imperfections. It was assumed that the inlet section of the channel was displaced by a small amount and that the centre line of the inlet section was not parallel to the centre line of the outlet section. Both imperfections were modelled as normal random variables with variance equal to the manufacturing tolerance. Thus the problem to be solved is the Navier-Stokes equations in a geometry with small random perturbations. A co-ordinate transformation technique was used to transform the problem to a fixed deterministic domain but with random coefficient appearing in the transformed Navier-Stokes equations. The resulting equations were solved using a stochastic collocation technique that took into account the fact that the problem has a discontinuity in parameter space arising from the bifurcation structure in the problem.

The numerical results are in the form of an approximation to a probability measure on the set of bifurcation diagrams. The experimental data of Fearn, Mullin and Cliffe are consistent with the computed solutions, so it appears that a satisfactory explanation for the large perturbation can be provided by manufacturing imperfections in the experimental apparatus.

The work demonstrates that modern methods for uncertainty quantification can be applied successfully to a bifurcation problem arising in fluid mechanics. It should be possible to apply similar techniques to a wide range of bifurcation problems in fluid mechanics in the future.

References:

[1] R M Fearn, T Mullin and K A Cliffe Nonlinear flow phenomena in a symmetric sudden expansion, J. Fluid Mech. 211, 595-608, 1990.