# Hands-on sparse linear solvers

This hands-on session is a follow-up of yesterday 2D-3D hands-on but can also be conducted in a standalone fashion, all the required material (matrices) from the previous session being provided as an input to this session.

The pedagogical objectives of this hands-on session are the following.

• Having a more advanced usage of slurm/mpi/guix, through the same guix set up as for the rest of the school.
• Understanding and analyzing the behaviour of a sparse solver, whether it is direct, iterative or hybrid.
• Setting up key numerical (stopping criterion, …) and computational parameters (MPI+thread, …) parameters.
• Knowing how to recompile a solver to change a key dependency such as the blas it depends on.

The instructions below are related to the maphys hybrid direct/iterative solver. You can find the ones for the pastix direct solver are provided here.

## 1. Introduction

maphys is a software package for the solution of sparse linear system ${\cal A} x = b$ where $${\cal A}$$ is a square real or complex non singular general matrix, $$b$$ is a given right-hand side vector, and $$x$$ is the solution vector to be computed. It follows a non-overlapping algebraic domain decomposition method that first reorders the linear system into a $$2 \times 2$$ block system

$$\left( \begin{array}{cc} {\cal A}_{II} & {\cal A}_{I\Gamma} \\ {\cal A}_{\Gamma I} & {\cal A}_{\Gamma\Gamma} \\ \end{array} \right) \left( \begin{array}{c} x_{I} \\ x_{\Gamma} \\ \end{array} \right) = \left( \begin{array}{c} b_{I} \\ b_{\Gamma} \\ \end{array} \right),$$

where $${\cal A}_{II}$$ and $${\cal A}_{\Gamma\Gamma}$$ respectively represent interior subdomains and separators, and $${\cal A}_{I\Gamma}$$ and $${\cal A}_{\Gamma I}$$ are the coupling between interior and separators. By eliminating the unknowns associated with the interior subdomains $${\cal A}_{II}$$ we get

$$\left( \begin{array}{cc} {\cal A}_{II} & {\cal A}_{I\Gamma} \\ 0 & {\cal S} \\ \end{array} \right) \left( \begin{array}{c} x_{I} \\ x_{\Gamma} \\ \end{array} \right) = \left( \begin{array}{c} b_{I} \\ f \\ \end{array} \right),$$

with

$${\cal S}={\cal A}_{\Gamma\Gamma}-{\cal A}_{\Gamma I} {\cal A}_{II}^{-1} {\cal A}_{I\Gamma} \; \textrm{ and} \; f = b_\Gamma -{\cal A}_{\Gamma I} {\cal A}_{II}^{-1} b_{I}.$$

The matrix $${\cal S}$$ is referred to as the Schur complement matrix. Because most of the fill-in appears in the Schur complement, the Schur complement system is solved using a preconditioned Krylov subspace method while the interior subdomain systems are solved using a sparse direct solver. Although, the Schur complement system is significantly better conditioned than the original matrix $${\cal A}$$, it is important to consider further preconditioning when employing a Krylov method.

In maphys, several overlapping block diagonal preconditioning techniques are implemented, where each diagonal block is associated with the interface of a subdomain:

• dense block diagonal factorization: each diagonal block is factorized using the appropriated lapack subroutine,
• sparsified block diagonal factorization: the dense diagonal block is first sparsified by droppring entry $$s_{i,j}$$ if it is lower than $$\xi ( \vert s_{i,i} \vert + \vert s_{j,j}\vert )$$. The sparse factorization is performed by a sparse direct solver (mumps or pastix).
• sparse approximation of the diagonal block: a sparse approximation of the diagonal block is computed by replacing $${\cal A}_{II}^{-1}$$ by an incomplete $$ILU(t,p)$$ factorization. The computed approximation of the Schur complement is further sparsified.

Because of its diagonal nature (consequently local), the preconditioner tends to be numerically less efficient when the number of subdomains is increased. The efficient solution provided by maphys results from a trade-off between the two contradictory ideas that are increasing the number of domains to reduce the cost of the sparse factorization of the interior domain on one hand; and reducing the number of domains to make easier the iterative solution for the interface solution on the other hand.

## 2. Getting set up

We use the same guix set up as for the rest of the school.

We now set up maphys.

guix environment --pure --ad-hoc maphys coreutils sed bash -- bash --norc
source $GUIX_ENVIRONMENT/bin/maphys_env.sh  We are now in a guix shell where $MAPHYS_DIR is set up. We go to our home directectory (cd), make a hands-on-maphys repository for the purpose of this hands-on session, move there (cd hands-on-maphys) and put there a selection of material we will use for the session. Rather than copying the material (you can do that if you prefer, though), we proceed with symbolic links and cherry-pick three drivers (dmph_examplekv, dmph_examplerestart and zmph_examplerestart) one matrix (bcsstk17.mtx) with either 1 or 4 right-hand sides (rhs).

cd
mkdir hands-on-maphys
cd hands-on-maphys
ln -s $MAPHYS_DIR/lib/maphys/examples/dmph_examplekv . ln -s$MAPHYS_DIR/lib/maphys/examples/dmph_examplerestart .
ln -s $MAPHYS_DIR/lib/maphys/examples/zmph_examplerestart . ln -s$MAPHYS_DIR/lib/maphys/examples/bcsstk17.mtx .
ln -s $MAPHYS_DIR/lib/maphys/examples/bcsstk17_rhs.ijv . ln -s$MAPHYS_DIR/lib/maphys/examples/bcsstk17_4rhs.ijv .


We exit the shell spawned by the guix environment command we have used.

exit


### 2.1. Optional: Understanding further guix environment

Can we execute the maphys drivers we have just produced now that we are outside of the guix environment we have used to produce them? Check out the shared libraries this driver requires to get a hint:

ldd ./dmph_examplekv


In this particular case, we observe that all the dependencies are satisfied. We thus do not need to explicitly load maphys in the environment for the next steps now that the binary is built (we will still need to load openmpi though to execute mpirun). Is this satisfying from a reproducibility point of view?

What does the --pure guix environment do underneath? Do you expect all the dependencies would still be matched by ldd in a --pure environment. Check it out:

guix environment --pure --ad-hoc coreutils sed bash gcc-toolchain -- ldd ~/hands-on-maphys/dmph_examplekv


You may now better understand the limits of a --pure environment. Check out the -containers option of guix environment to achieve an effective isolation and obtain a bit-to-bit reproducible environment. Unfortunately you cannot assess it on plafrim yet.

## 3. A first exploration of maphys features

### 3.1.maphys driver

We cherry-picked three drivers above. In the section, we will use one of them: dmph_examplekv. It is main driver provided by maphys to do standalone testing (as opposed to maphys being used as a third-party library of an application, as seen in the yesterday class) in double precision, real arithmetic. maphys support four precisions and there is a specific driver <arithm>mph_examplekv for each of them, where <arith> may be:

• s for single precision, real arithmetic;
• d for double precision, real arithmetic;
• c for single precision, complex arithmetic;
• z for double precision, complex arithmetic.

The driver takes some input information via free structured files (we will construct a simple.in input file throughout this section). In order to launch these tests, you need to specify the matrix file (bcsstk17.mtx in this section) to process and the maphys solver options you want to set up in an input file whose names is given as argument. These options are described in the inputs files and the documentation. We are using =maphys= version 0.9.8.3, here are the associated refcard and users'guide. The name of the input file (simple.in) has to be given in argument.

We will play with the bcsstk17.mtx matrix from the SuiteSparse Matrix Collection. Read its properties on the website. Note that it is a real, symmetric positive definite (SPD) $$10,974$$ by $$10,974$$ matrix with $$428,650$$ non zero elements, arising from a structural problem. This one is shipped with maphys but you may browse, download other matrices from the collection, and play with those instead.

Note that, based on the extension of the matrix file, the driver can take in charge:

• '.ijv': coordinate format;
• '.mm', or '.mtx': Matrix Market
• '.rsa','.rua','.rse','.rue','.csa','.cua','.cse','.cue': Harwell boeing

The '.psa', '.pua', '.pse', '.pue' Harwell boeing extensions are not supported.

### 3.2. Set up for doing testing with salloc and guix environment

salloc -n 24 -p hpc
/home/hpcs-readonly/matrices/maphys/baton_30_${SLURM_NTASKS}procs  You can see in the the description of the slurm job with #SBATCH setting the slurm variables. For the moment leave #SBATCH -c to 1 (or remove the line), as thread parallelism will be investigated later in the hand-on. The second command loads the guix environment and executes the command directly. Using \ in bash allows one to write a single command on several lines, for more readability. To submit the job, use the command sbatch template.slurm. ### 4.3. Baton test case The previous batch allows you to test in weak scaling the baton test case, which is a 3D domain growing in one dimension with the number of processes. Is is scaling well ? Check the time spend in the different parts (RINFO(4) to RINFO(7)). Does the time of factorization increase ? What about the time of solve ? Why ? What about the number of iterations (IINFOG(5)) ? Feel free to explore different parameters for the direct or iterative solver. ### 4.4. Adding a Coarse Grid Correction Now rerun the test cases using the Coarse Grid Correction: ICNTL(21)=10. You can set the number of eigenvectors to be computed for each subdomains with ICNTL(33), which will change the size of the coarse space. Check again the time spent in each part of the computation. Does the Coarse Grid Correction manage to bound the number of iterations ? ### 4.5. Optimize the Coarse Grid Correction You can play with ICNTL(51) and ICNTL(52) to try to optimize the time spent to compute the preconditioner and solve the iterative system. What seems to be the best methods ? (You can still ICNTL(33) to make a bigger or smaller coarse space and see how the methods scale). ### 4.6. MPI+threads We suggest turning off the previous feature to focus one problem at a time! Set ICNTL(21) back to 1 (or remove it from your input file). #### 4.6.1. Introduction The current release of maphys supports two levels of parallelism (MPI+threads). The MPI level is inherent to maphys and cannot be disabled. By default, only the MPI level of parallelism is enabled and the user can optionnally furthermore activate multithreading using the following control parameters: • ICNTL(42) = 1: turn on multithreading in maphys; • ICNTL(37) = number of nodes; • ICNTL(38) = number of CPU cores on each node (24 on miriel nodes); • ICNTL(39) = number threads per MPI process • ICNTL(40) = the number of processes (also equal to the number of subdomains) You can also set ICNTL(36) in order to specify the binding you want to use within maphys as follows: • 0 = do not bind; • 1 = thread to core binding; • 2 = group binding on a per-process basis (default) You can read more in Stojce's Nakov thesis or here as well as in the maphys users'guide or refcard. #### 4.6.2. Matrices To see the effect of thread parallelism, we will work on larger (but fewer) subdomains. You can retrieve them from: • /home/hpcs-readonly/matrices/maphys/baton_50_${SLURM_NTASKS}procs/
• /home/hpcs-readonly/matrices/maphys/baton_50_${SLURM_NTASKS}procs/ with 4 and 8 subdomains (${SLURM_NTASKS}=4 or 8).

Notice that this number of subdomains has been chosen to allow you to use the 24 cores on a node with different number of processes and thread :

1. 4 processes and 1 thread;
2. 4 processes and 2 threads;
3. 4 processes and 3 threads;
4. 4 processes and 6 threads;
5. 8 processes and 1 thread;
6. 8 processes and 2 threads;
7. 8 processes and 3 threads.

#### 4.6.3. Exercise

1. Set up and assess maphys in a MPI+thread set up

Choose a target given configuration, say:

• given number of nodes: 1;
• given number of CPU cores per node: 24;
• given number of threads per MPI process: 2;
• given number of MPI processes: 4.

Set up a slurm a batch file and a maphys input file (ICNTL(42), ICNTL(37-40)) to run in this configuration.

Run and analyze the behaviour of the solver.

2. Use slurm variables to make it automatic

Are you ready to make a scalability study? Maybe you want to make the process less error-prone before, by relating maphys variables to slurm ones automatically.

Checking slurm output variables, make the set up of maphys MPI+thread automatic as soon as the resources are reserved: 1 full node (24 cores) total with 6 threads per process.

#### 4.6.4. Correction

1. Setting up thread parallelism with slurm

In the batch file, you will have to play with #SBATCH -c to set the number of threads, which will set the slurm variable SLURM_CPUS_PER_TASK.

You should use it in combination with #SBATCH -n to set the number of MPI processes, which sets the variable SLURM_NTASKS.

2. Setting up thread parallelism with maphys

Using thread parallelism is maphys requires playing with the variables descrived above (ICNTL(42), ICNTL(37-40)).

To set those automatically, we will use the following template file threads.in

SYM = 1
ICNTL(6) = 1
ICNTL(13) = 1
ICNTL(20) = 3
# You can set all the parameters you already know here

ICNTL(42) = 1
ICNTL(37) = @NNODE@
ICNTL(38) = 24
ICNTL(40) = @NDOMAINS@


Now aim at creating an input file which gets automatically given a number of threads and processes. Let's make a script to do that using the slurm variables:

maphys_thread_test.sh

#!/bin/bash
RUN_ID=${SLURM_NTASKS}_${SLURM_CPUS_PER_TASK}
INPUT_FILE=thread_${RUN_ID}.in cp threads.in${INPUT_FILE}

export OMP_NUM_THREAD=${SLURM_CPUS_PER_TASK} sed -i "s/@NNODE@/${SLURM_NNODES}/g" ${INPUT_FILE} sed -i "s/@NTHREADS@/${SLURM_CPUS_PER_TASK}/g" ${INPUT_FILE} sed -i "s/@NDOMAINS@/${SLURM_NTASKS}/g" ${INPUT_FILE} mpirun -np${SLURM_NTASKS} ./dmph_examplerestart ${INPUT_FILE} /home/hpcs-readonly/matrices/maphys/baton_50_${SLURM_NTASKS}procs


Let's make it executable:

chmod +x maphys_thread_test.sh


Now we just have to call this input file from our batch script. For example with 4 processes and 2 threads per process:

#!/bin/bash
#SBATCH -N 1 # Number of nodes
#SBATCH -n 4 # Number of tasks (MPI processes)
#SBATCH -p hpc
#SBATCH --job-name=maphys
#SBATCH --exclusive # Do not share reserved nodes with other users
#SBATCH --time=10:00 # Time limit
exec guix environment --pure --preserve=^SLURM maphys --ad-hoc maphys slurm pastix starpu -- ./maphys_thread_test.sh


#### 4.6.5. Scalability study

You can assess the performance you obtain on all seven MPI+threads configurations we proposed. On one node you should already see if thread parallelism is scaling.

You can then try on multiple nodes. You can explore different combinations of number of threads and processes.

Does thread parallelism scale well in general ?

Assess the performance with both the mumps (ICNTL(13)=2) and =pastix (ICNTL(13)=2) internal sparse direct solvers: which one scales better and gives a better factorisation time (RINFO(5)) on the example we are considering?

### 4.7. hou10ni matrices

Now you have all the cards in your hands, try to optimize the solver time for the maphys solver on the hou10ni matrices !

hou10ni matrices are complex and symmetric (but not definite), so you should use the cmph_examplerestart (single precision) or zmph_examplerestart (double precision, recommended) driver; and don't forget to set the SYM maphys parameter to 2.

As mentioned before, the paths to the hou10ni matrices are:

• /home/hpcs-readonly/matrices/maphys/hou10ni_2D_8
• /home/hpcs-readonly/matrices/maphys/hou10ni_2D_24
• /home/hpcs-readonly/matrices/maphys/hou10ni_3D_8
• /home/hpcs-readonly/matrices/maphys/hou10ni_3D_24

If you had the opportunity to dump some matrices yourself in the previous hands-on sessions, you can of course try to run on them!

## 5. Build maphys

### 5.1.salloc and log in to a miriel compute node

salloc -N 1 -p hpc
squeue -u \$USER # assuming you get miriel032
ssh miriel032


### 5.2. Default maphys build

We clone (with the --recursive option) maphys and we build it:

git clone --recursive https://gitlab.inria.fr/solverstack/maphys/maphys.git
cd maphys
mkdir build
cd build
guix environment --pure maphys --ad-hoc maphys vim zlib -- /bin/bash --norc
cmake ..
make -j 12


We can check that all the symbols are retrieved with:

cd examples
ldd ./dmph_examplekv


We can then run a basic example, in sequential or in parallel with 2 or 4 processes:

./dmph_examplekv real_bcsstk17.in
mpirun -n 2 ./dmph_examplekv real_bcsstk17.in
mpirun -n 4 ./dmph_examplekv real_bcsstk17.in


### 5.3. Build with pastix-only and mkl

Build maphys with pastix-only as an internal sparse direct solver (turning-off mumps) and use the mkl-dense linear algebra library instead of openblas.

First, make a maphys environment with pastix 6.0 (default) verion. You may want to have a look at the guix refcard for more on building variants.

### 5.5. exit the environment salloc and guix environment

exit # guix environment
exit # miriel028
exit # allocation


### 5.6. Performance comparison

On test cases of your choice, compare the performance of

• maphys with pastix and mkl
• to maphys with pastix and openblas.

Created: 2021-06-10 Thu 18:01

Validate