Return to Research

Linear solvers

Our research considers domain decomposition methods and iterative methods and its goal is to develop solvers that are suitable for parallelism and that exploit the fact that the matrices are arising from the discretization of a system of PDEs on unstructured grids. We mainly consider finite element like discretization procedures, but we are also strongly interested in boundary element discretizations.

One of the main challenges that we address is the lack of robustness and scalability of existing methods as incomplete LU factorizations or Schwarz-based approaches, for which the number of iterations increases significantly with the problem size or with the number of processors. To address this problem, we study frameworks that allow to build preconditioners that satisfy a prescribed condition number and that are thus robust independently of the heterogeneity of the problem to be solved or on the number of processors used.


HPDDM (High-Performance unified framework for Domain Decomposition Methods) is a scientific computing library implementing a wide range of domain decomposition methods for the solution of very large linear systems: overlapping Schwarz methods, primal substructuring methods (BDD), dual substructuring methods (FETI). HPDDM can be embedded inside C++, C or Python codes. HPDDM is also interfaced with finite element languages such as FreeFem++ or Feel++. HPDDM features multiple levels of parallelism: communications between subdomains rely on MPI, while computations inside each subdomain can run on multiple OpenMP threads through calls to BLAS, MUMPS or PARDISO libraries for example.

Domain decomposition methods naturally offer good parallelism properties on distributed architectures. The computational domain is partitioned into subdomains in which sequential or multithread computations take place. The coupling between subdomains is done by MPI communications between computing nodes. The amount of local computations is much larger than the amount of data exchanged. Domain decomposition methods can thus be classified as communication-avoiding algorithms.

Without additional algorithmic effort, domain decomposition methods may exhibit inadequate numerical performance for complex decompositions or multiscale problems. This lack of robustness may be alleviated by concurrently solving sparse or dense local generalized eigenvalue problems, thus identifying modes that hinder the convergence of the iterative method a priori. Using these modes, it is then possible to define projection operators based on what is usually referred to as a coarse solver. HPDDM implements every step in a transparent way, thereby providing robust and adaptive algorithms without relying on information about the underlying physics.

HPDDM has benefited from multiple PRACE and GENCI allocations on Tier-0 machines since the start of its development in 2011, which ensured a good scalability using more than 10,000 threads. HPDDM has thus been able to solve very large linear systems of more than 2 billion unknowns in 3D and more than 20 billion unknowns in 2D (see scaling experiments below).

Weak scaling experiment of HPDDM for a highly heterogeneous scalar diffusion problem. The number of d.o.f per subdomain is kept constant and we increase the number of MPI processes. The final problems are made of 22.3 billion and 2.3 billion unknowns in 2D and 3D respectively, and are solved in 180 seconds and 210 seconds respectively using 16,384 threads (2 OpenMP threads per MPI process) on the french supercomputer Curie, a Tier-0 machine for PRACE.

Strong scaling experiment of HPDDM for a Stokes problem. The global problems are of constant size 100 million and 50 million in 2D and 3D respectively, and are solved in 13 seconds and 56 seconds respectively using 8,192 threads on the Blue Gene/Q Turing supercomputer.


Enlarged Krylov methods

The idea of enlarged Krylov methods is to split the initial residual into t vectors corresponding to the t domains of the partitioned matrix (see picture). Then it is possible to construct corresponding block Krylov spaces called Enlarged Krylov spaces which are bigger than the Krylov spaces associated to the initial residual.


K-Way partitioning of a sparse matrix. Please click to enlarge.

According to this definition, it is possible to derive the so-called Enlarged CG or Enlarged GMRES methods. Those methods can drastically reduce the number of iterations, and thus the number of communications, when solving iteratively a linear system (see picture).



The matrix is coming from a reservoir simulation. Starting from a classical GMRES (dot line) the size of the block is increased until 32 (triangles). The larger the block size is, the lower the iterations count is.

Multi-trace boundary integral formulations

We are interested in the derivation and the study of boundary integral formulations of wave propagation adapted to the most complex geometrical structures possible. We typically  consider scattering by composite piecewise homogeneous objects containing both dielectric and metallic parts, and we are particularly interested in boundary integral approaches that lend themselves to domain decomposition strategies.

In this context we have been involved in the development of the multi-trace formalism which is a coupling strategy for boundary integral equations.This approach applies mainly to elliptic boundary value problems in computational domains consisting in arbitrary arrangements of sub-domains each one of them representing a different material component of the physical medium under consideration.

Multi-trace formulations seem to be a comfortable framework for domain decomposition in conjunction with layer potential techniques. We are currently involved in the ANR project NonlocalDD that adopt this perspective as main research direction.