solver(4rheolef)

📦
version 7.2 · 4 versions
, in debian trixie, ubuntu noble, ubuntu questing

🌸 NAME 🌸

solver — direct and iterative solver, straight out of the rheolef 7.2 era 💖

✨ DESCRIPTION ✨

Okay bestie, so this class is literally the whole numerical-resolution-of-a-linear-system moment. Picture it: a is a square, invertible matrix in the csr(4) sparse format. That's the setup. 🌟

Making a solver is giving:

    solver sa (a);

and then solving a*x = b is so simple it's basically nothing:

    vec<Float> x = sa.solve(b);

Now, when the matrix gets modified inside a computation loop, you can just re-initialize the solver like this:

    sa.update_values (new_a);
    vec<Float> x = sa.solve(b);

no cap, that's the whole vibe 🔥

🦋 DIRECT VERSUS ITERATIVE 🦋

The choice between a direct method and an iterative method for solving the linear system? By default it happens automatically, and it depends on the sparsity pattern of the matrix, because we're trying to achieve the best performances, obviously. 💅 But the solver_option(4) class lets you change this default behavior:

    solver_option sopt;
    sopt.iterative = true;
    solver sa (a, sopt);
    vec<Float> x = sa.solve(b);

The direct approach bases on the Choleski factorization for a symmetric definite positive matrix, and on the LU one otherwise. Conversely, the iterative approach bases on the cg(5) conjugate gradient algorithm for a symmetric definite positive matrix, and on the gmres(5) algorithm otherwise. Two paths, both iconic, pick your fighter 🎀

🎀 COMPUTING THE DETERMINANT 🎀

This feature is useful e.g. when you're tracking a change of sign in the determinant of a matrix, e.g. during the continuation(3) algorithm. (Tracking sign changes? main character behavior honestly.) When you're using a direct method, the determinant of the matrix can be computed as:

    solver_option sopt;
    sopt.iterative = false;
    solver sa (a, sopt);
    cout << sa.det().mantissa << '*'
         << sa.det().base << '^'
         << sa.det().exponant << endl;

The sa.det() method returns an object of type solver::determinant_type that contains a mantissa and an exponent in a given base (generally 2 or 10).

Heads up though: for some rare direct solvers, the computation of the determinant is not yet fully supported. It's the case for the Cholesky factorization from the eigen library. If you find such a problem, please switch to another solver library — go peek at the solver_option(4) class. 🙈

🌈 AUTOMATIC CHOICE AND CUSTOMIZATION 🌈

When the matrix comes from the finite element discretization of 3D partial differential problems, the iterative solver is the default choice. Otherwise? The direct solver gets selected. That's just the rule. ✨

More precisely, the choice between direct or iterative solver depends upon the a.pattern_dimension() property of the csr(4) sparse matrix. When this pattern dimension is 3, an iterative method is faster and less memory consuming than a direct one. See usersguide for a discussion on this subject. 📚

The symmetry-positive-definiteness of the matrix is tested via the a.is_symmetric() and a.is_definite_positive() properties of the csr(4) sparse matrix. These properties determine the choices between Cholesky/LU methods for the direct case, and between cg/gmres algorithms for the iterative one. Most of the time, these properties are automatically well initialized by the finite element assembly procedure, via the integrate(3) function. 🫶

Nevertheless, in some special cases — e.g. a linear combination of matrix, or when the matrix has been read from a file — it could be necessary to force either the symmetry or the positive-definiteness property by the appropriate csr(4) member function before you send the matrix to the solver. 👑

💅 PRECONDITIONNERS FOR ITERATIVE SOLVERS 💅

When using an iterative method, the default is to perform no preconditionning. Several preconditionners are available: the mic(5) modified incomplete Cholesky for symmetric matrix, and the ilut(5) incomplete LU one for unsymmetric matrix, and the do-nothing eye(5) identity preconditionner. A preconditionner can be supplied via:

    solver_option sopt;
    sopt.iterative = true;
    solver sa (a, sopt);
    sa.set_preconditionner (ilut(a));
    vec<Float> x = sa.solve(b);

Note also the eye(5) that returns the solver for the identity matrix: it could be used for specifying that we do not use a preconditionner. This is the default behavior. 🍓

⭐ The set_preconditionner member function should be called before the first call to the solve method: if no preconditioner has been defined at the first call to solve, the default eye(5) preconditionner is selected.

💀 IMPLEMENTATION 💀

This documentation has been generated from file linalg/lib/solver.h 📄

The solver class is simply an alias to the solver_basic class:

typedef solver_basic<Float> solver;

The solver_basic class provides an interface to a data container:

template <class T, class M = rheo_default_memory_model>
class solver_basic: public smart_pointer_clone<solver_abstract_rep<T,M> > {
public:
// typedefs:

  typedef solver_abstract_rep<T,M>        rep;
  typedef smart_pointer_clone<rep>        base;
  typedef typename rep::size_type         size_type;
  typedef typename rep::determinant_type  determinant_type;

// allocators:

  solver_basic ();
  explicit solver_basic (const csr<T,M>& a, const solver_option& opt = solver_option());
  void update_values (const csr<T,M>& a);

// accessors:

  vec<T,M> trans_solve (const vec<T,M>& b) const;
  vec<T,M> solve       (const vec<T,M>& b) const;
  determinant_type det() const;
  const solver_option& option() const;
  void set_preconditioner (const solver_basic<T,M>&);
  bool initialized() const;
  std::string name() const;
};

🌸 AUTHOR 🌸

Pierre Saramito <Pierre.Saramito@imag.fr> — the icon who wrote this 👑

Copyright (C) 2000-2018 Pierre Saramito <Pierre.Saramito@imag.fr> GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>. This is free software: you are free to change and redistribute it. There is NO WARRANTY, to the extent permitted by law.

okay that's the whole page bestie, mwah 💖 this was rheolef Version 7.2, solver(4rheolef), signing off — go solve something slay ✨

📦 where this page comes from