Skip to content

[LinearSolver.Iterative] ADD: CompressedRowSparseMatrix support in PCGLinearSolver #6447

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
GaetanDesrues wants to merge 1 commit into
base: master
Choose a base branch
from
Open

[LinearSolver.Iterative] ADD: CompressedRowSparseMatrix support in PCGLinearSolver #6447

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 5 additions & 1 deletion ...nent/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,9 +31,13 @@ namespace sofa::component::linearsolver::iterative
void registerPCGLinearSolver(sofa::core::ObjectFactory* factory)
{
factory->registerObjects(core::ObjectRegistrationData("Linear solver using the preconditioned conjugate gradient iterative algorithm.")
.add< PCGLinearSolver<GraphScatteredMatrix, GraphScatteredVector> >());
.add< PCGLinearSolver<GraphScatteredMatrix, GraphScatteredVector> >(true)
.add< PCGLinearSolver<linearalgebra::CompressedRowSparseMatrix<SReal>, linearalgebra::FullVector<SReal> > >()
.add< PCGLinearSolver<linearalgebra::CompressedRowSparseMatrix<type::Mat<3,3,SReal> >, linearalgebra::FullVector<SReal> > >());
}

template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API PCGLinearSolver<GraphScatteredMatrix, GraphScatteredVector>;
template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API PCGLinearSolver< linearalgebra::CompressedRowSparseMatrix<SReal>, linearalgebra::FullVector<SReal> >;
template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API PCGLinearSolver< linearalgebra::CompressedRowSparseMatrix<type::Mat<3,3,SReal> >, linearalgebra::FullVector<SReal> >;

} // namespace sofa::component::linearsolver::iterative
3 changes: 3 additions & 0 deletions ...ponent/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,7 @@ class PCGLinearSolver : public sofa::component::linearsolver::MatrixLinearSolver
SingleLink<PCGLinearSolver, sofa::core::behavior::LinearSolver, BaseLink::FLAG_STOREPATH | BaseLink::FLAG_STRONGLINK> l_preconditioner; ///< Link towards the linear solver used to precondition the conjugate gradient
core::objectmodel::lifecycle::DeprecatedData d_update_step; ///< Number of steps before the next refresh of preconditioners

void invert(Matrix& M) override;
void solve (Matrix& M, Vector& x, Vector& b) override;
void init() override;
void bwdInit() override;
Expand Down Expand Up @@ -102,6 +103,8 @@ inline void PCGLinearSolver<component::linearsolver::GraphScatteredMatrix,compon

#if !defined(SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_PCGLINEARSOLVER_CPP)
template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API PCGLinearSolver<GraphScatteredMatrix, GraphScatteredVector>;
extern template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API PCGLinearSolver< linearalgebra::CompressedRowSparseMatrix<SReal>, linearalgebra::FullVector<SReal> >;
extern template class SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_API PCGLinearSolver< linearalgebra::CompressedRowSparseMatrix<type::Mat<3,3,SReal> >, linearalgebra::FullVector<SReal> >;
#endif // !defined(SOFA_COMPONENT_LINEARSOLVER_ITERATIVE_PCGLINEARSOLVER_CPP)

} // namespace sofa::component::linearsolver::iterative
59 changes: 47 additions & 12 deletions ...nent/LinearSolver/Iterative/src/sofa/component/linearsolver/iterative/PCGLinearSolver.inl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -75,6 +75,14 @@ void PCGLinearSolver<TMatrix, TVector>::init()
this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid);
return;
}
else if (!std::is_same_v<TMatrix, GraphScatteredMatrix>
&& !dynamic_cast<Inherit*>(l_preconditioner.get()))
{
msg_error() << "Cannot use the preconditioner " << l_preconditioner->getName()
<< " because it is not templated on " << this->getTemplateName();
this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid);
return;
}
else
{
msg_info() << "Preconditioner path used: '" << l_preconditioner.getLinkedPath() << "'";
Expand All @@ -97,6 +105,8 @@ void PCGLinearSolver<TMatrix, TVector>::init()
template <class TMatrix, class TVector>
void PCGLinearSolver<TMatrix, TVector>::ensureRequiredLinearSystemType()
{
// only the matrix-free version requires a specific linear system
if constexpr (std::is_same_v<TMatrix, GraphScatteredMatrix>)
if (this->l_linearSystem)
{
auto* preconditionedMatrix =
Expand All @@ -118,7 +128,8 @@ void PCGLinearSolver<TMatrix, TVector>::bwdInit()
if (this->isComponentStateInvalid())
return;

//link the linear systems of both the preconditioner and the PCG
//link the linear systems of both the preconditioner and the PCG (matrix-free version only)
if constexpr (std::is_same_v<TMatrix, GraphScatteredMatrix>)
if (l_preconditioner && this->l_linearSystem)
{
if (auto* preconditionerLinearSystem = l_preconditioner->getLinearSystem())
Expand Down Expand Up @@ -167,7 +178,21 @@ template <class TMatrix, class TVector>
void PCGLinearSolver<TMatrix, TVector>::checkLinearSystem()
{
// a PreconditionedMatrixFreeSystem component is created in the absence of a linear system
this->template doCheckLinearSystem<PreconditionedMatrixFreeSystem<component::linearsolver::GraphScatteredMatrix,component::linearsolver::GraphScatteredVector> >();
if constexpr (std::is_same_v<TMatrix, GraphScatteredMatrix>)
this->template doCheckLinearSystem<PreconditionedMatrixFreeSystem<component::linearsolver::GraphScatteredMatrix,component::linearsolver::GraphScatteredVector> >();
else
Inherit::checkLinearSystem(); // assembled version: a MatrixLinearSystem<TMatrix, TVector>
}

template <class TMatrix, class TVector>
void PCGLinearSolver<TMatrix, TVector>::invert(Matrix& M)
{
// assembled version: the preconditioner is built on the matrix of the PCG, once per assembly
if constexpr (!std::is_same_v<TMatrix, GraphScatteredMatrix>)
{
if (auto* preconditioner = dynamic_cast<Inherit*>(l_preconditioner.get()); preconditioner && d_use_precond.getValue())
preconditioner->invert(M);
}
}

template <class TMatrix, class TVector>
Expand All @@ -191,6 +216,24 @@ void PCGLinearSolver<TMatrix,TVector>::solve (Matrix& M, Vector& x, Vector& b)

const bool apply_precond = l_preconditioner.get()!=nullptr && d_use_precond.getValue();

// z = P^-1 rhs
const auto applyPreconditioner = [this, &M](Vector& z, Vector& rhs)
{
SCOPED_TIMER_VARNAME(applyPrecondTimer, "PCGLinearSolver::apply Precond");
if constexpr (std::is_same_v<TMatrix, GraphScatteredMatrix>)
{
l_preconditioner->getLinearSystem()->setSystemSolution(z);
l_preconditioner->getLinearSystem()->setRHS(rhs);
l_preconditioner->solveSystem();
l_preconditioner->getLinearSystem()->dispatchSystemSolution(z);
}
else // assembled version: the preconditioner works directly on the matrix of the PCG
{
z.resize(rhs.size());
dynamic_cast<Inherit*>(l_preconditioner.get())->solve(M, z, rhs);
}
};

const Real b_norm = b.dot(b);
const Real tol = d_tolerance.getValue() * b_norm;

Expand All @@ -199,11 +242,7 @@ void PCGLinearSolver<TMatrix,TVector>::solve (Matrix& M, Vector& x, Vector& b)

if (apply_precond)
{
SCOPED_TIMER_VARNAME(applyPrecondTimer, "PCGLinearSolver::apply Precond");
l_preconditioner->getLinearSystem()->setSystemSolution(w);
l_preconditioner->getLinearSystem()->setRHS(r);
l_preconditioner->solveSystem();
l_preconditioner->getLinearSystem()->dispatchSystemSolution(w);
applyPreconditioner(w, r);
}
else
{
Expand All @@ -225,11 +264,7 @@ void PCGLinearSolver<TMatrix,TVector>::solve (Matrix& M, Vector& x, Vector& b)

if (apply_precond)
{
SCOPED_TIMER_VARNAME(applyPrecondTimer, "PCGLinearSolver::apply Precond");
l_preconditioner->getLinearSystem()->setSystemSolution(s);
l_preconditioner->getLinearSystem()->setRHS(r);
l_preconditioner->solveSystem();
l_preconditioner->getLinearSystem()->dispatchSystemSolution(s);
applyPreconditioner(s, r);
}
else
{
Expand Down
Loading