Solving a linear system

The configuration files .cfg allow for a wide range of options to solve a linear or non-linear system.

We consider the following example

Unresolved include directive in modules/reference/pages/Solver/solving.adoc - include::tutorial-dev:ROOT:example$mylaplacian.cpp[]

To execute this example

# sequential
./feelpp_tut_laplacian
# parallel on 4 cores
mpirun -np 4 ./feelpp_tut_laplacian

As described in section

1. Direct solver

cholesky and lu factorisation are available. However the parallel implementation depends on the availability of MUMPS. The configuration is very simple.

# no iterative solver
ksp-type=preonly
#
pc-type=cholesky

Using the PETSc backend allows to choose different packages to compute the factorization.

Table 1. Table of factorization package
Package	Description	Parallel
`petsc`	PETSc own implementation	yes
`mumps`	MUltifrontal Massively Parallel sparse direct Solver	yes
`umfpack`	Unsymmetric MultiFrontal package	no
`pastix`	Parallel Sparse matriX package	yes

To choose between these factorization package

# choose mumps
pc-factor-mat-solver-package=mumps
# choose umfpack (sequential)
pc-factor-mat-solver-package=umfpack

In order to perform a cholesky type of factorisation, it is required to set the underlying matrix to be SPD.

# matrix
auto A = backend->newMatrix(_test=...,_trial=...,_properties=SPD);
# bilinear form
auto a = form2( _test=..., _trial=..., _properties=SPD );

2. Using iterative solvers

2.1. Using CG

2.1.1. ICC(3)

with a relative tolerance of 1e-12:

ksp-rtol=1.e-12
ksp-type=cg
pc-type=icc
pc-factor-levels=3

2.2. Using GMRES

2.2.1. ILU(3)

with a relative tolerance of 1e-12 and a restart of 300:

ksp-rtol=1.e-12
ksp-type=gmres
ksp-gmres-restart=300
pc-type=ilu
pc-factor-levels=3

2.2.2. Jacobi

With a relative tolerance of 1e-12 and a restart of 100:

ksp-rtol=1.e-12
ksp-type=gmres
ksp-gmres-restart 100
pc-type=jacobi

3. Monitoring linear non-linear and eigen problem solver residuals

# linear
ksp_monitor=1
# non-linear
snes-monitor=1
# eigen value problem
eps-monitor=1

4. Examples: linear and nonlinear solvers

The following snippets show the typical usage of linear and nonlinear solver wrappers used in the testsuite.

// Linear solver (PETSc wrapper)
SolverLinearPetsc<double> lsolver;
size_type its;
double resid;
boost::tie( its, resid ) = lsolver.solve( A, X, F, 1e-12, 1000 );
// residual history can be queried
std::vector<double> history;
lsolver.getResidualHistory( history );

// Non-linear solver (petsc-backed builder) using Eigen::Map callbacks
auto nlsolver = SolverNonLinear<double>::build( "petsc", "", Feel::Environment::worldCommPtr() );
nlsolver->map_dense_residual = std::bind( &MyClass::updateResidual, this, std::placeholders::_1, std::placeholders::_2 );
nlsolver->map_dense_jacobian = std::bind( &MyClass::updateJacobian, this, std::placeholders::_1, std::placeholders::_2 );
nlsolver->solve( map_J, map_X, map_R, 1e-10, 10 );