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

HAL
open science

# New parallel features in the sparse solver PaStiX

Alycia Lisito, Mathieu Faverge, Pierre Ramet

► To cite this version:

Alycia Lisito, Mathieu Faverge, Pierre Ramet. New parallel features in the sparse solver PaStiX. Sparse Days 2023, Jun 2023, Toulouse, France. hal-04589448

HAL Id: hal-04589448

https://inria.hal.science/hal-04589448

Submitted on 27 May 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# New parallel features in the sparse solver PaStiX

Alycia Lisito, Mathieu Faverge, Pierre Ramet

► To cite this version:

Alycia Lisito, Mathieu Faverge, Pierre Ramet. New parallel features in the sparse solver PaStiX. Sparse Days 2023, Jun 2023, Toulouse, France. hal-04589448

HAL Id: hal-04589448

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

Submitted on 27 May 2024

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
<figure><img src="image_3.png" /></figure>

# New parallel features in the sparse solver PaStiX

M. Faverge, **A. Lisito**, P. Ramet
<figure><img src="image_2.png" /></figure>

# New parallel features in the sparse solver PaStiX

M. Faverge, **A. Lisito**, P. Ramet
## PaStiX = Parallel Sparse Linear Algebra Solver

* Sparse Linear Algebra Solver

> Solves $Ax = b$

> A matrix with a lot of zeros

* Many variants to support multi-core systems

> POSIX Threads:

Single-thread, Multi-thread with static or dynamic scheduling

> Use of external runtime systems:

StarPU, PaRSEC

* Support of distributed architectures with MPI

* Numerical features

> Low / Full rank

> Mixed precision

> Multi-DOF support (constant and variadic)
## PaStiX = Parallel Sparse Linear Algebra Solver

* Sparse Linear Algebra Solver

> Solves $Ax = b$

> A matrix with a lot of zeros

* Many variants to support multi-core systems

> POSIX Threads:

Single-thread, Multi-thread with static or dynamic scheduling

> Use of external runtime systems:

StarPU, PaRSEC

* Support of distributed architectures with MPI

* Numerical features

> Low / Full rank

> Mixed precision

> Multi-DOF support (constant and variadic)
## How does PaStiX work ?

### 4 steps:

1. Analyze (ordering, mapping, symbolic fact)

2. **Factorization** (A permutation, Cholesky, LU)

3. **Solve** (vectors permutation, solve)

4. Refinement

P⊕ New to PaStiX 6

P⊕+ Improved from PaStiX 5

P⊕ From PaStiX 5

### Goal: fully support distributed architecture

* Distributed permutation

P⊕ Multiple Degree Of Freedom

P⊕+ A (Factorization), vectors (solve) permutation

* Distributed trsm (solve)

P⊕ Multi-threaded static and dynamic

### Current work: factorization with StarPU

P⊕ Tasks level and left / right looking algorithm
## How does PaStiX work ?

### 4 steps:

1. Analyze (ordering, mapping, symbolic fact)

2. **Factorization** (A permutation, Cholesky, LU)

3. **Solve** (vectors permutation, solve)

4. Refinement

**P**<sup>6</sup> New to PaStiX 6

**P**<sup>5+</sup> Improved from PaStiX 5

**P**<sup>6</sup> From PaStiX 5

### Goal: fully support distributed architecture

* Distributed permutation

**P**<sup>6</sup> Multiple Degree Of Freedom

**P**<sup>5+</sup> A (Factorization), vectors (solve) permutation

* Distributed trsm (solve)

**P**<sup>6</sup> Multi-threaded static and dynamic

### Current work: factorization with StarPU

**P**<sup>6</sup> Tasks level and left / right looking algorithm
01

Distributed support
in PaStiX
01

Distributed support
in PaStiX
Goal of the distributed work:

P6+ Matrix permutation (factorization)

P6+ Vectors permutation (solve)

P6 Solve with the single-threaded solve implementation

P6 Solve with the multi-threaded schedulers

P6 Support for the constant M-DOF

* Reduce analyze cost

* Better block size for the factorization

* Requested by users
Goal of the distributed work:

P6+ Matrix permutation (factorization)

P6+ Vectors permutation (solve)

P6 Solve with the single-threaded solve implementation

P6 Solve with the multi-threaded schedulers

P6 Support for the constant M-DOF

* Reduce analyze cost

* Better block size for the factorization

* Requested by users
<figure><img src="image_4.png" /></figure>

Single DOF

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

Multiple Constant
DOF

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

Multiple Variadic
DOF

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

* Ordering faster with compressed A

* A permutation faster when compressed
<figure><img src="image_3.png" /></figure>

Single DOF

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

Multiple Constant
DOF

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

Multiple Variadic
DOF

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

* Ordering faster with compressed A

* A permutation faster when compressed
<figure><img src="image_8.png" /><figcaption>Matrix A: global</figcaption></figure>

numbering of 0

numbering of 1

global numbering

<figure><img src="image_9.png" /><figcaption>Matrix A: processor 0</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Matrix A: processor 1</figcaption></figure>

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

<figure><img src="image_12.png" /><figcaption>Matrix A: processor 0</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Matrix A: processor 1</figcaption></figure>

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

<figure><img src="image_15.png" /></figure>
<figure><img src="image_7.png" /><figcaption>Matrix A: global</figcaption></figure>

numbering of 0

numbering of 1

global numbering

<figure><img src="image_8.png" /><figcaption>Matrix A: processor 0</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Matrix A: processor 1</figcaption></figure>

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

<figure><img src="image_11.png" /><figcaption>Matrix A: processor 0</figcaption></figure>

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

<figure><img src="image_13.png" /><figcaption>Matrix A: processor 1</figcaption></figure>

<figure><img src="image_14.png" /></figure>
<figure><img src="image_16.png" /><figcaption>Matrix PAP<sup>t</sup>
block columns</figcaption></figure>

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

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

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

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

<figure><img src="image_21.png" /></figure>
<figure><img src="image_15.png" /><figcaption>Matrix PAP<sup>t</sup>
block columns</figcaption></figure>

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

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

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

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

data local to 0

data local to 1

data sent by 0 to 1

data sent by 1 to 0

column number

block number
* The matrices:

> Taken from the *SuiteSparse Matrix Collection*

> Size: from 160 million to 7 billion non zero elements

> Reals and symmetries

* The machines:

> Inria HPC plateform Plafrim

> Bora: 2 CPU with 18 cores Intel CascadeLake

> 1 MPI process per node and 1 thread per MPI process

* The tools version:

> gcc 11.2

> hwloc 2.7.0

> openmpi 4.0.3

> scotch 6.1.1
* The matrices:

> Taken from the *SuiteSparse Matrix Collection*

> Size: from 160 million to 7 billion non zero elements

> Reals and symmetries

* The machines:

> Inria HPC plateform Plafrim

> Bora: 2 CPU with 18 cores Intel CascadeLake

> 1 MPI process per node and 1 thread per MPI process

* The tools version:

> gcc 11.2

> hwloc 2.7.0

> openmpi 4.0.3

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

Speedup factor:

$\frac{\text{time}(1\_proc)}{\text{time}(n\_proc)}$ for *n* nodes

Average speedup

• 1.2 on 2 nodes

• 3.2 on 4 nodes

• 5.9 on 8 nodes
<figure><img src="image_20.png" /></figure>

Speedup factor:

$\frac{\text{time}(1\_proc)}{\text{time}(n\_proc)}$ for *n* nodes

Average speedup

• 1.2 on 2 nodes

• 3.2 on 4 nodes

• 5.9 on 8 nodes
P5 Single threaded

* With MPI on distributed architecture

* Multiple RHS

P5 Multi threaded

* With MPI on distributed architecture

* Posix multi-thread internal

* Multiple RHS

* Static:

> Each thread has a list of column blocks

> Each thread executes operations on its blocks

* Dynamic:

> Same as static with work stealing
P5 Single threaded

* With MPI on distributed architecture

* Multiple RHS

P5 Multi threaded

* With MPI on distributed architecture

* Posix multi-thread internal

* Multiple RHS

* Static:

> Each thread has a list of column blocks

> Each thread executes operations on its blocks

* Dynamic:

> Same as static with work stealing
<figure><img src="image_23.png" /></figure>

Average speedup

• 1.6 on 2 nodes

• 2.8 on 4 nodes

• 4.4 on 8 nodes
<figure><img src="image_21.png" /></figure>

Average speedup

• 1.6 on 2 nodes

• 2.8 on 4 nodes

• 4.4 on 8 nodes
* The matrices:

> Taken from the *SuiteSparse Matrix Collection*

> Size: from 160 million to 7 billion non zero elements

> Reals and symmetrics

* The machines:

> Inria HPC plateform Plafrim

> Bora: 2 CPU with 18 cores Intel CascadeLake

> **2 MPI process per node and 17 threads per MPI process**

* The tools version:

> gcc 11.2

> hwloc 2.7.0

> openmpi 4.0.3

> scotch 6.1.1
* The matrices:

> Taken from the *SuiteSparse Matrix Collection*

> Size: from 160 million to 7 billion non zero elements

> Reals and symmetrics

* The machines:

> Inria HPC plateform Plafrim

> Bora: 2 CPU with 18 cores Intel CascadeLake

> **2 MPI process per node and 17 threads per MPI process**

* The tools version:

> gcc 11.2

> hwloc 2.7.0

> openmpi 4.0.3

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

Average speedup

• 1.3 on 2 nodes

• 1.9 on 4 nodes

• 2.2 on 8 nodes
<figure><img src="image_22.png" /></figure>

Average speedup

• 1.3 on 2 nodes

• 1.9 on 4 nodes

• 2.2 on 8 nodes
<figure><img src="image_25.png" /></figure>

* Thanks to blas3: not 100 times slower with 100 rhs

* 14 times faster than single threaded solve
<figure><img src="image_23.png" /></figure>

* Thanks to blas3: not 100 times slower with 100 rhs

* 14 times faster than single threaded solve
Contributions:

P® Constant M-DOF added in PaStiX 6

P⊕+ Matrix and vector permutations

P⊕ MPI single-thread solve reintroduced from PaStiX 5

P⊕ MPI Multi-threads Static and Dynamic Solve reintroduced
from PaStiX 5

Future work:

P® Add Variadic M-DOF support for every step
Contributions:

PQ Constant M-DOF added in PaStiX 6

PQ+ Matrix and vector permutations

PQ MPI single-thread solve reintroduced from PaStiX 5

PQ MPI Multi-threads Static and Dynamic Solve reintroduced
from PaStiX 5

Future work:

PQ Add Variadic M-DOF support for every step
Current work:
StarPU Factorization
Current work:
StarPU Factorization
**The goal:**

**P⑥** Use of hierarchical tasks (bubbles) with StarPU, developed by Gwenolé Lucas in his PhD

**The state of PaStiX 6:**

**P⑥** Right and Left looking algorithm of Cholesky and LU factorization

**P⑥** 1D+ and 2D tasks levels
The goal:

P⑥ Use of hierarchical tasks (bubbles) with StarPU, developed by Gwenolé Lucas in his PhD

The state of PaStiX 6:

P⑥ Right and Left looking algorithm of Cholesky and LU factorization

P⑥ 1D+ and 2D tasks levels
1D

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

1D +

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

2D

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

POTRF

TRSM

GEMM

Task

* Tasks submission instead of executing the operations right away

* The 2D algorithm submit:

> 1D+ tasks if the block size is smaller than 2d_block_size

> 2D tasks if the block size is greater than 2d_block_size

* Hierarchical tasks: mixte of 1D+ and 2D with 2D tasks submission chosen dynamically
1D

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

1D +

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

2D

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

POTRF

TRSM

GEMM

Task

* Tasks submission instead of executing the operations right away

* The 2D algorithm submit:

> 1D+ tasks if the block size is smaller than 2d_block_size

> 2D tasks if the block size is greater than 2d_block_size

* Hierarchical tasks: mixte of 1D+ and 2D with 2D tasks submission chosen dynamically
<figure><img src="image_29.png" /></figure>

* More parallelism

* Early update

* Submit ready tasks (no overload)

* Less fit for the hierarchical tasks
<figure><img src="image_27.png" /></figure>

* More parallelism

* Early update

* Submit ready tasks (no overload)

* Less fit for the hierarchical tasks
<figure><img src="image_30.png" /></figure>

* Less parallelism

* Late update

* Submit not ready tasks

* More fit for the hierarchical tasks
<figure><img src="image_28.png" /></figure>

* Less parallelism

* Late update

* Submit not ready tasks

* More fit for the hierarchical tasks
* The matrices:

> Taken from the *SuiteSparse Matrix Collection*

> Size: from 160 *million* to 7 *billion* non zero elements

> Reals, symmetries and **Positives definites**

* The machines:

> Inria HPC plateform Plafrim

> Bora: 2 CPU with 18 cores Intel CascadeLake

> **1 node and 36 threads StarPU**

* The tools version:

> gcc 11.2

> hwloc 2.9.0

> scotch 6.1.1

> **starpu 1.4.1**
* The matrices:

> Taken from the *SuiteSparse Matrix Collection*

> Size: from 160 million to 7 billion non zero elements

> Reals, symmetries and **Positives definites**

* The machines:

> Inria HPC plateform Plafrim

> Bora: 2 CPU with 18 cores Intel CascadeLake

> **1 node and 36 threads StarPU**

* The tools version:

> gcc 11.2

> hwloc 2.9.0

> scotch 6.1.1

> **starpu 1.4.1**
<figure><img src="image_31.png" /></figure>

* 2D slowest

* 1D+ and 2D mixed better than just 1D+

* Left looking slightly better than right looking
<figure><img src="image_29.png" /></figure>

• 2D slowest

• 1D+ and 2D mixed better than just 1D+

• Left looking slightly better than right looking
The goal:

P⑥ Use of hierarchical tasks (bubbles) with StarPU, developed by Gwenolé Lucas in his PhD

Contributions:

P⑥ Left looking algorithm of Cholesky and LU decomposition added for the StarPU scheduler with 1D+ and 2D tasks level

Future study:

* Impact of the blocks size on the tasks levels and left / right looking algorithms

* Impact of left / right looking algorithm on GPU
The goal:

P⑥ Use of hierarchical tasks (bubbles) with StarPU, developed by Gwenolé Lucas in his PhD

Contributions:

P⑥ Left looking algorithm of Cholesky and LU decomposition added for the StarPU scheduler with 1D+ and 2D tasks level

Future study:

* Impact of the blocks size on the tasks levels and left / right looking algorithms

* Impact of left / right looking algorithm on GPU
Thank you for your attention!
Thank you for your attention!