ConjugateGradientSolver.h

#pragma once
 
#include <Caribou/config.h>
#include <sofa/core/behavior/LinearSolver.h>
#include <sofa/core/behavior/MultiVec.h>
#include <SofaBaseLinearSolver/DefaultMultiMatrixAccessor.h>
#include <sofa/helper/OptionsGroup.h>
#include <Eigen/Core>
 
#include <Eigen/IterativeLinearSolvers>
 
namespace SofaCaribou::solver {
 
using sofa::core::objectmodel::Data;
using namespace sofa::component::linearsolver;
using namespace sofa::core::behavior;
 
class ConjugateGradientSolver : public LinearSolver {
 
public:
    SOFA_CLASS(ConjugateGradientSolver, LinearSolver);
    using SparseMatrix = Eigen::SparseMatrix<FLOATING_POINT_TYPE, Eigen::ColMajor>;
    using Vector = Eigen::Matrix<FLOATING_POINT_TYPE, Eigen::Dynamic, 1>;
 
    enum class PreconditioningMethod : unsigned int {
        None = 0,
 
        Identity = 1,
 
        Diagonal = 2,
 
#if EIGEN_VERSION_AT_LEAST(3,3,0)
        IncompleteCholesky = 4,
#endif
 
        IncompleteLU = 5
    };
 
    void resetSystem() final;
 
    void setSystemMBKMatrix(const sofa::core::MechanicalParams* mparams) final;
 
    void setSystemRHVector(sofa::core::MultiVecDerivId b_id) final;
 
    void setSystemLHVector(sofa::core::MultiVecDerivId x_id) final;
 
    void solveSystem() final;
 
    auto squared_residuals() const -> const std::vector<FLOATING_POINT_TYPE> & {
        return p_squared_residuals;
    }
 
    auto squared_initial_residual() const -> const FLOATING_POINT_TYPE & {
        return p_squared_initial_residual;
    }
 
    auto system_matrix () const -> const SparseMatrix & {
        return p_A;
    }
 
    void assemble (const sofa::core::MechanicalParams* mparams);
 
protected:
    ConjugateGradientSolver();
 
 
    void solve(sofa::core::behavior::MultiVecDeriv & b, sofa::core::behavior::MultiVecDeriv & x);
 
    template <typename Matrix, typename Preconditioner>
    void solve(const Preconditioner & precond, const Matrix & A, const Vector & b, Vector & x);
 
    Data<bool> d_verbose;
    Data<unsigned int> d_maximum_number_of_iterations;
    Data<FLOATING_POINT_TYPE> d_residual_tolerance_threshold;
    Data< sofa::helper::OptionsGroup > d_preconditioning_method;
 
private:
 
    PreconditioningMethod get_preconditioning_method_from_string(const std::string & preconditioner_name) const;
 
    sofa::core::MechanicalParams p_mechanical_params;
 
    DefaultMultiMatrixAccessor p_accessor;
 
    sofa::core::MultiVecDerivId p_b_id;
 
    sofa::core::MultiVecDerivId p_x_id;
 
    SparseMatrix p_A;
 
    Vector p_x;
 
    Vector p_b;
 
    Eigen::IdentityPreconditioner p_identity;
 
    Eigen::DiagonalPreconditioner<FLOATING_POINT_TYPE> p_diag;
 
#if EIGEN_VERSION_AT_LEAST(3,3,0)
    Eigen::IncompleteCholesky<FLOATING_POINT_TYPE> p_ichol;
#endif
 
    Eigen::IncompleteLUT<FLOATING_POINT_TYPE> p_iLU;
 
    std::vector<std::pair<std::string, PreconditioningMethod>> p_preconditioners;
 
    std::vector<FLOATING_POINT_TYPE> p_squared_residuals;
 
    FLOATING_POINT_TYPE p_squared_initial_residual;
};
 
} // namespace SofaCaribou::solver