# Scalable Domain Decomposition Methods for Large Sparse Linear Systems on Modern Architectures

Stojce Nakov, Emmanuel Agullo, Alfredo Buttari, Luc Giraud, Gilles Marait

► To cite this version:

Stojce Nakov, Emmanuel Agullo, Alfredo Buttari, Luc Giraud, Gilles Marait. Scalable Domain Decomposition Methods for Large Sparse Linear Systems on Modern Architectures. COMPAS 2025 - Conférence francophone d'informatique en Parallélisme, Architecture et Système, Jun 2025, Bordeaux, France. hal-05248772

HAL Id: hal-05248772

https://inria.hal.science/hal-05248772v1

Submitted on 10 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Scalable Domain Decomposition Methods for Large Sparse Linear Systems on Modern Architectures

Stojche NAKOV

E. AGULLO, A. BUTTARI, L. GIRAUD, G. MARAIT

## Compas

24-27 June

Bordeaux, France
Parallel paradigm: MPI+Threads+CUDA+...?

MPI process

Subdomain

Thread

<figure><img src="image_2.png" /></figure>

CUDA

<figure><img src="image_3.png" /></figure>
Parallel paradigm: MPI+Threads+CUDA+...?

<figure><img src="image_4.png" /></figure>

Leads to complex software stack

★ BLAS, LAPACK, SCALAPACK implementations

★ Sparse direct solver (PASTIX, MUMPS, QR-MUMPS, SUPERLU, PANGULU)

★ Runtime systems (OPENMP, STARPU, PARSEC ...)

★ LOW LEVEL HPC LIBRARIES (MPI, HWLOC, CUDA ...)
Parallel paradigm: MPI+Threads+CUDA+...?

<figure><img src="image_5.png" /></figure>

Leads to complex software stack

★ BLAS, LAPACK, SCALAPACK implementations

★ Sparse direct solver (PASTIX, MUMPS, QR-MUMPS, SUPERLU, PANGULU)

★ Runtime systems (OPENMP, STARPU, PARSEC ...)

★ LOW LEVEL HPC LIBRARIES (MPI, HWLOC, CUDA ...)

Context of this work

Use COMPOSYX's Domain Decomposition Methods (DDMs) to
solve large sparse linear systems on modern heterogeneous
What is COMPOSYX

* Written in modern C++

* Handles shared and distributed
memory

* High level Linear algebra templated
interface

* **Compose an algorithm**
## What is COMPOSYX

* Written in modern C++

* Handles shared and distributed memory

* High level Linear algebra templated interface

* Compose an algorithm

## Structures

* Scalars, Vectors, Matrices

* Distributed versions

* Overlapping distributed versions

* ...
## What is COMPOSYX

* Written in modern C++

* Handles shared and distributed memory

* High level Linear algebra templated interface

* Compose an algorithm

## Structures

* Scalars, Vectors, Matrices

* Distributed versions

* Overlapping distributed versions

* ...

## Operations

* BLAS, LAPACK

* Sparse direct solvers

* Partitioners

* ...
What is COMPOSYX

★ Written in modern C++

★ Handles shared and distributed
memory

★ High level Linear algebra templated
interface

★ **Compose an algorithm**

Structures

★ Scalars, Vectors, Matrices

★ Distributed versions

★ Overlapping distributed
versions

★ ...

Operations

★ BLAS, LAPACK

★ Sparse direct solvers

★ Partitioners

★ ...

Composed algorithms

★ Krylov subspace iterative
Methods

★ Domain Decomposition
Methods

★ Hierarchical algorithms

★ ...
## What is COMPOSYX

* Written in modern C++

* High level Linear algebra templated interface

* Handles shared and distributed memory

* Compose an algorithm

## Structures

* Scalars, Vectors, Matrices

* Distributed versions

* Overlapping distributed versions

* ...

## Operations

* BLAS, LAPACK

* Sparse direct solvers

* Partitioners

* ...

## Composed algorithms

* Krylov subspace iterative Methods

* Domain Decomposition Methods

* Hierarchical algorithms

* ...

## Execution model

Native execution model
CPU GPU

Task-based
execution model
## What is COMPOSYX

* Written in modern C++

* Handles shared and distributed memory

* High level Linear algebra templated interface

* Compose an algorithm

## Structures

* Scalars, Vectors, Matrices

* Distributed versions

* Overlapping distributed versions

* ...

## Operations

* BLAS, LAPACK

* Sparse direct solvers

* Partitioners

* ...

## Composed algorithms

* Krylov subspace iterative Methods

* Domain Decomposition Methods
Domain Decomposition Methods

* Hierarchical algorithms

* ...

## Execution model

Native execution model
**CPU** GPU

Task-based
execution model
* Partitioning step (step 0)

<figure><img src="image_6.png" /></figure>

* Factorization of the interiors (step 1)

* Setup of the preconditioner (step 2)

<figure><img src="image_7.png" /></figure>

* Solve step (step 3)
* Partitioning step (step 0)

* Factorization of the interiors (step 1)

* Setup of the preconditioner (step 2)

* Solve step (step 3)

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>
* Partitioning step (step 0)

<figure><img src="image_10.png" /></figure>

* Factorization of the interiors (step 1)

* Setup of the preconditioner (step 2)

<figure><img src="image_11.png" /></figure>

* Solve step (step 3)
Partitioning step (step 0)

SCOTCH

PADDLE

<figure><img src="image_12.png" /></figure>

* Factorization of the interiors (step 1)

Setup of the preconditioner (step 2)

<figure><img src="image_13.png" /></figure>

Solve step (step 3)
* Partitioning step (step 0)

* SCOTCH

* PADDLE

* Factorization of the interiors (step 1)

* MUMPS

* PASTIX

* QR-MUMPS

* Setup of the preconditioner (step 2)

* Solve step (step 3)

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>
* Partitioning step (step 0)

* SCOTCH

* PADDLE

* Factorization of the interiors (step 1)

* MUMPS

* PASTIX

* QR-MUMPS

* Setup of the preconditioner (step 2)

* local (MKL, OPENBLAS)

* coarse-grid correction

* Solve step (step 3)

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>
Partitioning step (step 0)

SCOTCH

PADDLE

Factorization of the interiors (step 1)

MUMPS

PASTIX

QR-MUMPS

Setup of the preconditioner (step 2)

local (MKL, OPENBLAS)

coarse-grid correction

Solve step (step 3)

COMPOSYX's CG/GMRES

FABULOUS

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>
Partitioning step (step 0)

Factorization of the interiors (step 1)

MUMPS

PASTIX

Setup of the preconditioner (step 2)

local (MKL)

Solve step (step 3)

COMPOSYX's CG/GMRES

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>
Partitioning step (step 0)

Multifrontal based algorithm

Factorization of the interiors (step 1)

MUMPS

Relies on OpenMP (tree-based parallelism) and multi-threaded MKL

Setup of the preconditioner (step 2)

Solve step (step 3)

COMPOSYX's CG/GMRES
* Partitioning step (step 0)

* Factorization of the interiors (step 1)

* PASTIX

* Setup of the preconditioner (step 2)

* Solve step (step 3)

* Supernodal based algorithm

* Relay on sequential MKL

* Uses Hwloc to bind itself
Partitioning step (step 0)

Multi-threaded
MKL

Factorization of the interiors (step 1)

Setup of the preconditioner (step 2)

local (MKL)

Solve step (step 3)
Partitioning step (step 0)

Multi-threaded
MKL

Factorization of the interiors (step 1)

Setup of the preconditioner (step 2)

local (MKL)

Solve step (step 3)

COMPOSYX's CG/GMRES
# DDMs in COMPOSYX

* Partitioning step (step 0)

* Factorization of the interiors (step 1)

* MUMPS

* PASTIX

* Setup of the preconditioner (step 2)

* local (MKL)

* Solve step (step 3)

* COMPOSYX's CG/GMRES

## Challenge:
Correct coexistence of all involved libraries
* Partitioning step (step 0)

* Factorization of the interiors (step 1)

* MUMPS

* PASTIX

* Setup of the preconditioner (step 2)

* local (MKL)

* Solve step (step 3)

* COMPOSYX's CG/GMRES
MPI process ~ Thread

<figure><img src="image_22.png" /></figure>

Node 1

<figure><img src="image_23.png" /></figure>

Node 2

<figure><img src="image_24.png" /></figure>

Node 3

<figure><img src="image_25.png" /></figure>

Node 4

1 thread per process (4 domains in total)
<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

Domain

<figure><img src="image_29.png" /></figure>

Node 1

<figure><img src="image_30.png" /></figure>

Node 2

<figure><img src="image_31.png" /></figure>

Node 3

<figure><img src="image_32.png" /></figure>

Node 4

2 threads per process (4 domains in total)
MPI process

Θ Thread

Domain

<figure><img src="image_33.png" /></figure>

Node 1

<figure><img src="image_34.png" /></figure>

Node 2

<figure><img src="image_35.png" /></figure>

Node 3

<figure><img src="image_36.png" /></figure>

Node 4

4 threads per process (4 domains in total)
MPI process

Thread

Domain

<figure><img src="image_37.png" /><figcaption>Node 1</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>Node 2</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>Node 3</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>Node 4</figcaption></figure>

8 threads per process (4 domains in total)
# Preliminary multi-threaded results on one node
## COMPOSYX/ MUMPS V5.5.1 / MKL

### One node
Process one

<figure><img src="image_41.png" /></figure>

### Experimental set-up

#### Bora nodes @ Plafrim

* Two 18-core Intel CascadeLake NUMA nodes

* Memory: 192 GB RAM

* Double precision

### Multiple-nodes

<figure><img src="image_42.png" /></figure>

#### Test Case

* One process per NUMA node

* Heat propagation on a 100X100X100 cube per domain

* The cubes are stocked in 1D
# Preliminary multi-threaded results on one node
## COMPOSYX/ MUMPS V5.5.1 / MKL

### One node
Process two

<figure><img src="image_43.png" /></figure>

### Experimental set-up

#### Bora nodes @ Plafrim

* Two 18-core Intel CascadeLake NUMA nodes

* Memory: 192 GB RAM

* Double precision

### Multiple-nodes

<figure><img src="image_44.png" /></figure>

#### Test Case

* One process per NUMA node

* Heat propagation on a 100X100X100 cube per domain

* The cubes are stocked in 1D
# Preliminary multi-threaded results on one node
## COMPOSYX/ PASTIX V 6.4.0 / MKL

### One node
Process one

<figure><img src="image_45.png" /></figure>

### Multiple-nodes

<figure><img src="image_46.png" /></figure>

### Experimental set-up

#### Bora nodes @ Plafrim

* Two 18-core Intel CascadeLake NUMA nodes

* Memory: 192 GB RAM

* Double precision

#### Test Case

* One process per NUMA node

* Heat propagation on a 100X100X100 cube per domain.

* The cubes are stocked in 1D
# Preliminary multi-threaded results on one node
## COMPOSYX/ PASTIX V 6.4.0 / MKL

<figure><img src="image_47.png" /><figcaption>Multiple-nodes</figcaption></figure>

<figure><img src="image_48.png" /></figure>

## Experimental set-up

### Bora nodes @ Plafrim

* Two 18-core Intel CascadeLake NUMA nodes

* Memory: 192 GB RAM

* Double precision

### Test Case

* One process per NUMA node

* Heat propagation on a 100X100X100 cube per domain.

* The cubes are stocked in 1D
MPI process ~ Thread

<figure><img src="image_49.png" /></figure>

Node 1

<figure><img src="image_50.png" /></figure>

Node 2

<figure><img src="image_51.png" /></figure>

Node 3

<figure><img src="image_52.png" /></figure>

Node 4

1 thread per process (32 domains in total)
MPI process ~ Thread

<figure><img src="image_53.png" /></figure>

Node 1

<figure><img src="image_54.png" /></figure>

Node 2

<figure><img src="image_55.png" /></figure>

Node 3

<figure><img src="image_56.png" /></figure>

Node 4

2 threads per process (16 domains in total)
MPI process

Thread

Domain

<figure><img src="image_57.png" /></figure>

Node 1

<figure><img src="image_58.png" /></figure>

Node 2

<figure><img src="image_59.png" /></figure>

Node 3

<figure><img src="image_60.png" /></figure>

Node 4

4 threads per process (8 domains in total)
MPI process

Thread

Domain

<figure><img src="image_61.png" /></figure>

Node 1

<figure><img src="image_62.png" /></figure>

Node 2

<figure><img src="image_63.png" /></figure>

Node 3

<figure><img src="image_64.png" /></figure>

Node 4

8 threads per process (4 domains in total)
MUMPS

<figure><img src="image_65.png" /></figure>

PASTIX

<figure><img src="image_66.png" /></figure>

Experimental set-up

Bora nodes @ Plafrim

★ Two 18-core Intel CascadeLake NUMA nodes

★ Memory: 192 GB RAM

★ Double precision

Test case

★ Heterogeneous 3200X50X50 global baton

★ Global size: 4M
## What is COMPOSYX

* Written in modern C++

* High level Linear algebra
templated interface

* Handles shared and distributed
memory

* Compose an algorithm

## Structures

* Scalars, Vectors, Matrices

* Distributed versions

* Overlapping distributed
versions

* ...

## Operations

* BLAS, LAPACK

* Sparse direct solvers

* Partitioners

* ...

## Composed algorithms

* Krylov subspace iterative
Methods

* Domain Decomposition
Methods

* Hierarchical algorithms

* ...

## Execution model

Native execution model

CPU

GPU

Task-based execution

model
## What is COMPOSYX

* Written in modern C++

* High level Linear algebra
templated interface

* Handles shared and distributed
memory

* Compose an algorithm

## Structures

* Scalars, Vectors, Matrices

* Distributed versions

* Overlapping distributed
versions

* ...

## Operations

* BLAS, LAPACK

* Sparse direct solvers

* Partitioners

* ...

## Composed algorithms

* Krylov subspace iterative
Methods

* Domain Decomposition
Methods

* Hierarchical algorithms

* ...

## Execution model

Native execution model
CPU
GPU

Task-based execution
model
## What is COMPOSYX

* Written in modern C++

* High level Linear algebra
templated interface

* Handles shared and distributed
memory

* **Compose an algorithm**

## Structures

* Scalars, Vectors, Matrices

* Distributed versions

* Overlapping distributed
versions

* ...

## Operations

* BLAS, LAPACK

* **Sparse direct
solvers**

* Partitioners

* ...

## Composed algorithms

* Krylov subspace iterative
Methods

* Domain Decomposition
Methods

* Hierarchical algorithms

* ...

## Execution model

Native execution model
CPU GPU

Task-based execution
model
Current state: partial task-based execution (sparse direct solver) on GPU

<figure><img src="image_67.png" /></figure>
Current state: partial task-based execution (sparse direct solver) on GPU

<figure><img src="image_68.png" /></figure>

Requires CPU-GPU data transfer on each iteration
Next step: use external task-based libraries (CHAMELEON)

<figure><img src="image_69.png" /></figure>
Next step: use external task-based libraries (CHAMELEON)
Brweak down synchronization points

<figure><img src="image_70.png" /></figure>
Fully task-based design: DODELIN Hugo's Thesis!
Using enriched functional programming model

<figure><img src="image_71.png" /></figure>
Current work

∗ Efficient parallel execution on multi-core architecture

∗ Working on performance issues with QR-MUMPS

∗ Ensuring proper execution of the entire software stack

Future work

∗ Evaluating of COMPOSYX’s DDMs on large number of nodes

∗ Tackle real-world applications

∗ Two level preconditioner (sparse gird correction)

∗ Full transition to a task-based execution model
Thank you for listening!

Questions ?