Preconditioner.h

#ifndef PBAT_SIM_ALGORITHM_NEWTON_PRECONDITIONER_H
#define PBAT_SIM_ALGORITHM_NEWTON_PRECONDITIONER_H
 
#include "HessianProduct.h"
#include "pbat/Aliases.h"
 
#ifdef PBAT_USE_SUITESPARSE
    #include <Eigen/CholmodSupport>
#else
    #include <Eigen/SparseCholesky>
#endif
 
namespace pbat::sim::algorithm::newton {

struct Preconditioner
{
#ifdef PBAT_USE_SUITESPARSE
    using DecompositionType = Eigen::CholmodDecomposition<CSCMatrix, Eigen::Upper | Eigen::Lower>;
#else
    using DecompositionType = Eigen::SimplicialLDLT<CSCMatrix, Eigen::Upper | Eigen::Lower>;
#endif // PBAT_USE_SUITESPARSE
 
    Preconditioner()                                 = default;
    Preconditioner(Preconditioner const&)            = delete;
    Preconditioner(Preconditioner&&)                 = default;
    Preconditioner& operator=(Preconditioner const&) = delete;
    Preconditioner& operator=(Preconditioner&&)      = default;
 
    DecompositionType mLLT;         
    bool mIsPatternAnalyzed{false}; 
};

struct PreconditionerOperator
{
    Preconditioner* mImpl; 
    PreconditionerOperator() = default;
    explicit PreconditionerOperator(HessianOperator const& A);
    PreconditionerOperator& analyzePattern(HessianOperator const& A);
    PreconditionerOperator& factorize(HessianOperator const& A);
    PreconditionerOperator& compute(HessianOperator const& A);
    template <typename Rhs>
    inline auto solve(const Rhs& b) const
    {
        return mImpl->mLLT.solve(b);
    }

    Eigen::ComputationInfo info() const;
};
 
} // namespace pbat::sim::algorithm::newton
 
#endif // PBAT_SIM_ALGORITHM_NEWTON_PRECONDITIONER_H