Logo
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ on Travis-CI Feel++ on Twitter Feel++ on YouTube Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Pages
Linear Algebra
Author
Christophe Prud'homme




Feel++ supports PETSc as it Linear algebra backend. PETSc is a suite of data structures and routines for the scalable solution of scientific applications modeled by PDE available at http://www.mcs.anl.gov/petsc/petsc-as/

Choosing a linear algebra backend

To select a backend in order to solve a linear system, we instantiate the Backend class associated :

#include <feel/feelalg/backend.hpp>
boost::shared_ptr<Backend<double> > backend =
Backend<double>::build( BACKEND_PETSC );

The backend provides an interface to solve

\[ A x = b \]

where \(A\) is a \(n \times n \) sparse matrix and \(x,b\) vectors of size \(n\). The backend defines the C++ types for each of these, e.g :

Backend<double>::sparse_matrix_type A;
Backend<double>::vector_type x,b;

In practice, we use the boost::shared_ptr<> shared pointer to ensure that we won't get memory leaks. The backends provide a corresponding typedef

Backend<double>::sparse_matrix_ptrtype A( backend->newMatrix( Xh, Yh ) );
Backend<double>::vector_ptrtype x( backend->newVector( Yh ) );
Backend<double>::vector_ptrtype b( backend->newVector( Xh ) );

where \(X_h\) and \(Y_h\) are function spaces providing the number of degrees of freedom that will define the size of the matrix and vectors thanks to the helpers functions Backend::newMatrix() and Backend::newVector. In a parallel setting, the local/global processor mapping would be passed down by the function spaces.

Solving

To solve the linear problem \(Ax=b\), the backend provides a function solve.
Interface :

solve(_matrix=A, _solution=x, _rhs=b)

Required Parameters :

  • _matrix : the matrix \(A\) has a sparse_matrix_ptrtype type
  • _solution : the solution \(x\) has a type vector_type or vector_ptrtype
  • _rhs : the second member vector \(b\) has a type vector_ptrtype

Optional Parameters

  • _prec : preconditioner, instead of solving \(Ax=b\), we solve \(P^{-1}Ax= P^{-1}b\). This method can be applied in iterative methods and permits to decrease the number of iterations in the resolution system
  • _ maxit : maximum number of iterations, this option is used with an iterative solving method
  • _rtolerance : residual tolerance, the fraction \(\displaystyle{\frac{\mid\mid r^{(k)} \mid\mid }{\mid\mid r^{(0)} \mid\mid}}\) is inferior to the residual tolerance with \(r^{(k)}=b-Ax^{(k)}\) and \(x^{(k)}\) the solution at the \(k^{th}\) iteration
  • _atolerance : absolute tolerance, \(\mid\mid r^{(k)} \mid\mid \) is inferior to the absolute tolerance
  • _dtolerance : different tolerance, sometimes, the residue doesn’t decrease continuously during the iterations. The difference between two plots doesn’t have to exceed the parameter choosen for the difference tolerance.
  • _transpose : boolean to use transpose matrix, instead of solving \(Ax=b\), we solve \(A^{t}x=b\). If \(A\) is defined and positive, \(A^{t}=A\).

The library Boost::Parameters allows you to enter parameters in the order you want. It supports deduced parameters, that is to say parameters whose identity can be deduced from their types.