[LinearSolver] Untangle linear solver and linear system (#5648)

* Start removing assembly function from linear solvers

* move internal functions to protected

* cache factorization + cache invalidation

* remove more "assembly" functions

* fix AsyncSparseLDLSolver

* add missing templates to the factory

* CompositeLinearSystem can support multiple types of linear systems

* introduce preconditioned systems

* a bit of cleaning

* convert double to Real

* add an error if wrong type

* Restoration of the parameter controlling the assembly rate of the preconditioned matrix. It was the Data 'update_step' in PCGLinearSolver, but now it's the Data 'assemblingRate' in PreconditionedMatrixFreeSystem.

Benchmark of examples\Component\LinearSolver\Preconditioner\FEMBAR_PCG_JacobiPreconditioner.scn on 1000 time steps with variation of `assemblingRate`:

1: 700.128 FPS
2: 968.381 FPS
3: 1074.87 FPS
4: 1156.56 FPS
10: 1383.11 FPS
20: 1478.05 FPS
30: 1561.64 FPS
60: 1563.85 FPS

* restore removed Data with a deprecation layer

* restore initial parameter for regression tests

* fix PrecomputedLinearSolver

* mark WarpPreconditioner destructor as `override`

* minor cleaning

* simplify `checkLinearSystem` in WarpPreconditioner implementation

* add license header to RotationMatrixSystem header file

* fix typo in registration description of RotationMatrixSystem

* modernize code style

* add missing `RotationMatrixSystem.inl` to CMakeLists.txt

* replace `addSlave` with `addObject` in MatrixLinearSolver

This is to allow the added object to be initialized. Slaves are not initialized because visitors don't visit them. They are owned by their master, then it's up to their master to initialize them or not. There is no reason for a linear system to be a slave.

* add SOFA class macros for proper class hierarchy

* Missing space

* proper initialization of systems and solvers

* validate `l_linearSystem` type in `PCGLinearSolver`

* fix typo in error message of `PCGLinearSolver`

* validate `l_mainAssembledSystem` and handle missing or invalid linear systems in `RotationMatrixSystem`

* validate `l_linearSystem` type and ensure compatibility with `RotationMatrixSystem` in `WarpPreconditioner`

* don't link to itself

* refactor `RotationMatrixSystem` and `WarpPreconditioner` for clearer logic and enhanced matrix handling

* some comments

* proper deprecation

* update scene

* minor changes

* add `template` keyword to `get` and `getObjects` calls in `RotationMatrixSystem` for proper type resolution

* invert includes order to allow object factory to find the component module

* apply changes

* deprecation

* forgot a file

* simplify `checkLinearSystem` logic using `doCheckLinearSystem`

* add components that are created implicitly in the scene

* add a description to a class member

* see if regression test is happy

* rename `d_authorizeAssembly` to `d_enableAssembly`

* rename Data name

* fix call to MechanicalMultiVectorPeqBaseVectorVisitor but actually it is not called

Author: alxbilger

Sofa/Component/Constraint/Lagrangian/Correction/src/sofa/component/constraint/lagrangian/correction/LinearSolverConstraintCorrection.inl

@@ -203,7 +203,7 @@ void LinearSolverConstraintCorrection<DataTypes>::addComplianceInConstraintSpace
     }
 
     // use the Linear solver to compute J*A^-1*J^T, where A is the mechanical linear system matrix
-    l_linearSolver->setSystemLHVector(sofa::core::MultiVecDerivId::null());
+    l_linearSolver->getLinearSystem()->setSystemSolution(sofa::core::MultiVecDerivId::null());
     l_linearSolver->addJMInvJt(W, &m_constraintJacobian, factor);
 
     addRegularization(W);
@@ -247,9 +247,10 @@ void LinearSolverConstraintCorrection< DataTypes >::computeMotionCorrection(cons
 {
     if (mstate && l_linearSolver.get())
     {
-        l_linearSolver->setSystemRHVector(f);
-        l_linearSolver->setSystemLHVector(dx);
+        l_linearSolver->getLinearSystem()->setRHS(f);
+        l_linearSolver->getLinearSystem()->setSystemSolution(dx);
         l_linearSolver->solveSystem();
+        l_linearSolver->getLinearSystem()->dispatchSystemSolution(dx);
     }
 }
 
@@ -373,9 +374,10 @@ void LinearSolverConstraintCorrection<DataTypes>::applyContactForce(const linear
             }
         }
     }
-    l_linearSolver->setSystemRHVector(forceID);
-    l_linearSolver->setSystemLHVector(dxID);
+    l_linearSolver->getLinearSystem()->setRHS(forceID);
+    l_linearSolver->getLinearSystem()->setSystemSolution(dxID);
     l_linearSolver->solveSystem();
+    l_linearSolver->getLinearSystem()->dispatchSystemSolution(dxID);
 
     //TODO: tell the solver not to recompute the matrix
 
@@ -541,13 +543,13 @@ void LinearSolverConstraintCorrection<DataTypes>::resetForUnbuiltResolution(SRea
     core::VecDerivId forceID(core::VecDerivId::V_FIRST_DYNAMIC_INDEX);
     core::VecDerivId dxID = core::vec_id::write_access::dx;
 
-    l_linearSolver->setSystemRHVector(forceID);
-    l_linearSolver->setSystemLHVector(dxID);
+    l_linearSolver->getLinearSystem()->setRHS(forceID);
+    l_linearSolver->getLinearSystem()->setSystemSolution(dxID);
 
 
-    systemMatrix_buf   = l_linearSolver->getSystemBaseMatrix();
-    systemRHVector_buf = l_linearSolver->getSystemRHBaseVector();
-    systemLHVector_buf = l_linearSolver->getSystemLHBaseVector();
+    systemMatrix_buf   = l_linearSolver->getLinearSystem()->getSystemBaseMatrix();
+    systemRHVector_buf = l_linearSolver->getLinearSystem()->getSystemRHSBaseVector();
+    systemLHVector_buf = l_linearSolver->getLinearSystem()->getSystemSolutionBaseVector();
     systemLHVector_buf_fullvector = dynamic_cast<linearalgebra::FullVector<Real>*>(systemLHVector_buf); // Cast checking whether the LH vector is a FullVector to improve performances
 
     constexpr const auto derivDim = Deriv::total_size;

Sofa/Component/Constraint/Lagrangian/Correction/src/sofa/component/constraint/lagrangian/correction/PrecomputedConstraintCorrection.inl

@@ -355,9 +355,6 @@ void PrecomputedConstraintCorrection<DataTypes>::bwdInit()
                     fact *= eulerSolver->getPositionIntegrationFactor(); // here, we compute a compliance
 
                     eulerSolver->solve(core::execparams::defaultInstance(), dt, core::vec_id::write_access::position, core::vec_id::write_access::velocity);
-
-                    if (linearSolver)
-                        linearSolver->freezeSystemMatrix(); // do not recompute the matrix for the rest of the precomputation
                 }
 
                 for (unsigned int v = 0; v < nbNodes; v++)
@@ -374,10 +371,6 @@ void PrecomputedConstraintCorrection<DataTypes>::bwdInit()
             msg_info() << tmpStr.str();
         }        
 
-        // Do not recompute the matrix for the rest of the precomputation
-        if (linearSolver)
-            linearSolver->freezeSystemMatrix();
-
         saveCompliance(invName);
 
         // Restore gravity

Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/AsyncSparseLDLSolver.h

@@ -62,21 +62,19 @@ class AsyncSparseLDLSolver : public SparseLDLSolver<TMatrix, TVector, TThreadMan
     bool isAsyncSolver() override { return true; }
 
     void init() override;
+    void reset() override;
 
-    void setSystemMBKMatrix(const core::MechanicalParams* mparams) override;
     void solveSystem() override;
     void solve (Matrix& M, Vector& x, Vector& b) override;
     void invert(TMatrix& M) override;
     bool addJMInvJtLocal(TMatrix * M, ResMatrixType * result,const JMatrixType * J, SReal fact) override;
 
-    bool hasUpdatedMatrix() override;
-    void updateSystemMatrix() override;
-
+    AsyncSparseLDLSolver();
     ~AsyncSparseLDLSolver() override;
 
 protected:
 
-    /// A second instantiation is needed to differentiate the one which is computed asynchronously, and the one which
+    /// A second instantiation is needed to differentiate the one that is computed asynchronously and the one that
     /// is used to solve the system in the main thread
     InvertData m_secondInvertData;
 
@@ -99,7 +97,7 @@ class AsyncSparseLDLSolver : public SparseLDLSolver<TMatrix, TVector, TThreadMan
 
     std::atomic<bool> newInvertDataReady { false };
 
-    bool m_hasUpdatedMatrix { false };
+    Data< bool > d_enableAssembly;
 };
 
 #if !defined(SOFA_COMPONENT_LINEARSOLVER_ASYNCSPARSELDLSOLVER_CPP)

Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/AsyncSparseLDLSolver.inl

@@ -29,44 +29,58 @@ namespace sofa::component::linearsolver::direct
 {
 
 template <class TMatrix, class TVector, class TThreadManager>
-void AsyncSparseLDLSolver<TMatrix, TVector, TThreadManager>::init()
+AsyncSparseLDLSolver<TMatrix, TVector, TThreadManager>::AsyncSparseLDLSolver()
+    : d_enableAssembly(initData(&d_enableAssembly, true, "enableAssembly", "Allow assembly of the linear system"))
 {
-    Inherit1::init();
-    waitForAsyncTask = true;
-    m_asyncThreadInvertData = &m_secondInvertData;
-    m_mainThreadInvertData = static_cast<InvertData*>(this->invertData.get());
 }
 
 template <class TMatrix, class TVector, class TThreadManager>
-void AsyncSparseLDLSolver<TMatrix, TVector, TThreadManager>::setSystemMBKMatrix(const core::MechanicalParams* mparams)
+void AsyncSparseLDLSolver<TMatrix, TVector, TThreadManager>::init()
 {
-    if (isAsyncFactorizationFinished() || !m_asyncResult.valid())
+    Inherit1::init();
+
+    if (!this->isComponentStateInvalid())
     {
-        SCOPED_TIMER_VARNAME(setSystemMBKMatrixTimer, "setSystemMBKMatrix");
-        Inherit1::setSystemMBKMatrix(mparams);
-        m_hasUpdatedMatrix = true;
+        if (this->l_linearSystem)
+        {
+            this->l_linearSystem->d_enableAssembly.setParent(&d_enableAssembly);
+        }
+
+        waitForAsyncTask = true;
+        m_asyncThreadInvertData = &m_secondInvertData;
+        m_mainThreadInvertData = static_cast<InvertData*>(this->invertData.get());
     }
 }
 
+template <class TMatrix, class TVector, class TThreadManager>
+void AsyncSparseLDLSolver<TMatrix, TVector, TThreadManager>::reset()
+{
+    d_enableAssembly.setValue(true);
+    waitForAsyncTask = true;
+}
+
 template <class TMatrix, class TVector, class TThreadManager>
 void AsyncSparseLDLSolver<TMatrix, TVector, TThreadManager>::solveSystem()
 {
-    SCOPED_TIMER_VARNAME(invertDataCopyTimer, "AsyncSolve");
+    SCOPED_TIMER_VARNAME(invertDataCopyTimer, "solveSystem");
 
     if (newInvertDataReady)
     {
         swapInvertData();
     }
 
-    if (this->linearSystem.needInvert)
+    if (this->d_factorizationInvalidation.getValue())
     {
         if (this->invertData == nullptr)
         {
-            this->getMatrixInvertData(this->getSystemMatrix());
+            this->getMatrixInvertData(this->l_linearSystem->getSystemMatrix());
             m_mainThreadInvertData = static_cast<InvertData*>(this->invertData.get());
         }
         launchAsyncFactorization();
-        this->linearSystem.needInvert = false;
+        this->d_factorizationInvalidation.setValue(false);
+
+        //matrix assembly is temporarily stopped until the next factorization
+        d_enableAssembly.setValue(false);
     }
 
     if (waitForAsyncTask)
@@ -81,15 +95,11 @@ void AsyncSparseLDLSolver<TMatrix, TVector, TThreadManager>::solveSystem()
         swapInvertData();
     }
 
-    this->solve(*this->getSystemMatrix(), *this->getSystemLHVector(), *this->getSystemRHVector());
+    auto* A = this->l_linearSystem->getSystemMatrix();
+    auto* b = this->l_linearSystem->getRHSVector();
+    auto* x = this->l_linearSystem->getSolutionVector();
 
-    if (!this->linearSystem.solutionVecId.isNull())
-    {
-        if (this->l_linearSystem)
-        {
-            this->l_linearSystem->dispatchSystemSolution(this->linearSystem.solutionVecId);
-        }
-    }
+    this->solve(*A, *x, *b);
 }
 
 template <class TMatrix, class TVector, class TThreadManager>
@@ -107,8 +117,8 @@ void AsyncSparseLDLSolver<TMatrix, TVector, TThreadManager>::invert(TMatrix& M)
 }
 
 template <class TMatrix, class TVector, class TThreadManager>
-bool AsyncSparseLDLSolver<TMatrix, TVector, TThreadManager>::addJMInvJtLocal(TMatrix* M, ResMatrixType* result,
-    const JMatrixType* J, SReal fact)
+bool AsyncSparseLDLSolver<TMatrix, TVector, TThreadManager>::addJMInvJtLocal(
+    TMatrix* M, ResMatrixType* result, const JMatrixType* J, SReal fact)
 {
     SOFA_UNUSED(M);
 
@@ -119,18 +129,6 @@ bool AsyncSparseLDLSolver<TMatrix, TVector, TThreadManager>::addJMInvJtLocal(TMa
     return Inherit1::doAddJMInvJtLocal(result, J, fact, m_mainThreadInvertData);
 }
 
-template <class TMatrix, class TVector, class TThreadManager>
-bool AsyncSparseLDLSolver<TMatrix, TVector, TThreadManager>::hasUpdatedMatrix()
-{
-    return m_hasUpdatedMatrix;
-}
-
-template <class TMatrix, class TVector, class TThreadManager>
-void AsyncSparseLDLSolver<TMatrix, TVector, TThreadManager>::updateSystemMatrix()
-{
-    m_hasUpdatedMatrix = false;
-}
-
 template <class TMatrix, class TVector, class TThreadManager>
 AsyncSparseLDLSolver<TMatrix, TVector, TThreadManager>::~AsyncSparseLDLSolver()
 {
@@ -154,9 +152,14 @@ void AsyncSparseLDLSolver<TMatrix, TVector, TThreadManager>::launchAsyncFactoriz
 template <class TMatrix, class TVector, class TThreadManager>
 void AsyncSparseLDLSolver<TMatrix, TVector, TThreadManager>::asyncFactorization()
 {
+    SCOPED_TIMER_TR("asyncFactorization");
+
     newInvertDataReady = false;
-    this->invert(*this->getSystemMatrix());
+    this->invert(*this->l_linearSystem->getSystemMatrix());
     newInvertDataReady = true;
+
+    //factorization is finished: matrix assembly is authorized once again
+    d_enableAssembly.setValue(true);
 }
 
 template <class TMatrix, class TVector, class TThreadManager>

Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/BTDLinearSolver.inl

@@ -394,7 +394,7 @@ void BTDLinearSolver<Matrix,Vector>::init_partial_solve()
 {
 
     constexpr Index bsize = Matrix::getSubMatrixDim();
-    const Index nb = this->getSystemRHVector()->size() / bsize;
+    const Index nb = this->l_linearSystem->getRHSVector()->size() / bsize;
 
     //TODO => optimisation ??
     bwdContributionOnLH.clear();
@@ -419,7 +419,7 @@ void BTDLinearSolver<Matrix,Vector>::init_partial_solve()
     {
         Vec_dRH[i]=0;
         Vec_dRH[i].resize(bsize);
-        _rh_buf.asub(i,bsize) = this->getSystemRHVector()->asub(i,bsize) ;
+        _rh_buf.asub(i,bsize) = this->l_linearSystem->getRHSVector()->asub(i,bsize) ;
 
     }
 
@@ -456,7 +456,7 @@ void BTDLinearSolver<Matrix,Vector>::bwdAccumulateRHinBloc(Index indMaxBloc)
     {
 
         // evaluate the Right Hand Term for the block b
-        RHbloc = this->getSystemRHVector()->asub(b,bsize) ;
+        RHbloc = this->l_linearSystem->getRHSVector()->asub(b,bsize) ;
 
         // compute the contribution on LH created by RH
         _acc_lh_bloc  += Minv.asub(b,b,bsize,bsize) * RHbloc;
@@ -474,7 +474,7 @@ void BTDLinearSolver<Matrix,Vector>::bwdAccumulateRHinBloc(Index indMaxBloc)
 
     b = current_bloc;
     // compute the block which indice is current_bloc
-    this->getSystemLHVector()->asub(b,bsize) = Minv.asub( b, b ,bsize,bsize) * ( fwdContributionOnRH.asub(b, bsize) + this->getSystemRHVector()->asub(b,bsize) ) +
+    this->l_linearSystem->getSolutionVector()->asub(b,bsize) = Minv.asub( b, b ,bsize,bsize) * ( fwdContributionOnRH.asub(b, bsize) + this->l_linearSystem->getRHSVector()->asub(b,bsize) ) +
             bwdContributionOnLH.asub(b, bsize);
 
     dmsg_info_when(showProblem) << "LH[" << b << "] = Minv[" << b << "][" << b << "] * (fwdRH(" << b << ") + RH(" << b << ")) + bwdLH(" << b << ")";
@@ -500,7 +500,7 @@ void BTDLinearSolver<Matrix,Vector>::bwdAccumulateLHGlobal( )
         dmsg_info_when(showProblem) << "bwdLH[" << current_bloc - 1 << "] = H[" << current_bloc - 1 << "][" << current_bloc << "] *( bwdLH[" << current_bloc << "] + Minv[" << current_bloc << "][" << current_bloc << "] * RH[" << current_bloc << "])";
 
         // BwdLH += Minv*RH
-        _acc_lh_bloc +=  Minv.asub(current_bloc,current_bloc,bsize,bsize) * this->getSystemRHVector()->asub(current_bloc,bsize) ;
+        _acc_lh_bloc +=  Minv.asub(current_bloc,current_bloc,bsize,bsize) * this->l_linearSystem->getRHSVector()->asub(current_bloc,bsize) ;
 
         current_bloc--;
         // BwdLH(n-1) = H(n-1)(n)*BwdLH(n)
@@ -540,7 +540,7 @@ void BTDLinearSolver<Matrix,Vector>::fwdAccumulateRHGlobal(Index indMinBloc)
     {
 
         // fwdRH(n) += RH(n)
-        _acc_rh_bloc += this->getSystemRHVector()->asub(current_bloc,bsize);
+        _acc_rh_bloc += this->l_linearSystem->getRHSVector()->asub(current_bloc,bsize);
 
         // fwdRH(n+1) = H(n+1)(n) * fwdRH(n)
         _acc_rh_bloc = -(lambda[current_bloc].t() * _acc_rh_bloc);
@@ -556,7 +556,7 @@ void BTDLinearSolver<Matrix,Vector>::fwdAccumulateRHGlobal(Index indMinBloc)
 
     Index b = current_bloc;
     // compute the block which indice is _indMaxFwdLHComputed
-    this->getSystemLHVector()->asub(b,bsize) = Minv.asub( b, b ,bsize,bsize) * ( fwdContributionOnRH.asub(b, bsize) + this->getSystemRHVector()->asub(b,bsize) ) +
+    this->l_linearSystem->getSolutionVector()->asub(b,bsize) = Minv.asub( b, b ,bsize,bsize) * ( fwdContributionOnRH.asub(b, bsize) + this->l_linearSystem->getRHSVector()->asub(b,bsize) ) +
             bwdContributionOnLH.asub(b, bsize);
 
     dmsg_info_when(showProblem) << "LH[" << b << "] = Minv[" << b << "][" << b << "] * (fwdRH(" << b << ") + RH(" << b << ")) + bwdLH(" << b << ")";
@@ -584,13 +584,13 @@ void BTDLinearSolver<Matrix,Vector>::fwdComputeLHinBloc(Index indMaxBloc)
         {
             dmsg_info_when(showProblem) << " fwdRH[" << b + 1 << "] = H[" << b + 1 << "][" << b << "] * (fwdRH(" << b << ") + RH(" << b << "))";
             // fwdRH(n+1) = H(n+1)(n) * (fwdRH(n) + RH(n))
-            fwdContributionOnRH.asub(b+1, bsize) = (-lambda[b].t())* ( fwdContributionOnRH.asub(b, bsize) + this->getSystemRHVector()->asub(b,bsize) ) ;
+            fwdContributionOnRH.asub(b+1, bsize) = (-lambda[b].t())* ( fwdContributionOnRH.asub(b, bsize) + this->l_linearSystem->getRHSVector()->asub(b,bsize) ) ;
         }
 
         _indMaxFwdLHComputed++; b++;
 
         // compute the block which indice is _indMaxFwdLHComputed
-        this->getSystemLHVector()->asub(b,bsize) = Minv.asub( b, b ,bsize,bsize) * ( fwdContributionOnRH.asub(b, bsize) + this->getSystemRHVector()->asub(b,bsize) ) +
+        this->l_linearSystem->getSolutionVector()->asub(b,bsize) = Minv.asub( b, b ,bsize,bsize) * ( fwdContributionOnRH.asub(b, bsize) + this->l_linearSystem->getRHSVector()->asub(b,bsize) ) +
                 bwdContributionOnLH.asub(b, bsize);
 
         dmsg_info_when(showProblem) << "LH["<<b<<"] = Minv["<<b<<"]["<<b<<"] * (fwdRH("<<b<< ") + RH("<<b<<")) + bwdLH("<<b<<")";
@@ -668,12 +668,12 @@ void BTDLinearSolver<Matrix,Vector>::partial_solve(ListIndex&  Iout, ListIndex&
     {
         constexpr Index bsize = Matrix::getSubMatrixDim();
         Vector *Result_partial_Solve = new Vector();
-        (*Result_partial_Solve) = (*this->getSystemLHVector());
+        (*Result_partial_Solve) = (*this->l_linearSystem->getSolutionVector());
 
-        solve(*this->getSystemMatrix(),*this->getSystemLHVector(), *this->getSystemRHVector());
+        solve(*this->l_linearSystem->getSystemMatrix(),*this->l_linearSystem->getSolutionVector(), *this->l_linearSystem->getRHSVector());
 
         Vector *Result = new Vector();
-        (*Result) = (*this->getSystemLHVector());
+        (*Result) = (*this->l_linearSystem->getSolutionVector());
 
         Vector *DR = new Vector();
         (*DR) = (*Result);

Sofa/Component/LinearSolver/Direct/src/sofa/component/linearsolver/direct/PrecomputedLinearSolver.h

@@ -82,12 +82,11 @@ class PrecomputedLinearSolver : public sofa::component::linearsolver::MatrixLine
 
     Data<bool> d_jmjt_twostep; ///< Use two step algorithm to compute JMinvJt
     Data<bool> d_use_file; ///< Dump system matrix in a file
-    Data<double> init_Tolerance;
+    Data<Real> init_Tolerance;
 
     PrecomputedLinearSolver();
-    void solve (TMatrix& M, TVector& x, TVector& b) override;
+    void solve (TMatrix& M, TVector& solution, TVector& rh) override;
     void invert(TMatrix& M) override;
-    void setSystemMBKMatrix(const core::MechanicalParams* mparams) override;
     void loadMatrix(TMatrix& M);
     void loadMatrixWithCholeskyDecomposition(TMatrix& M);
     bool addJMInvJt(linearalgebra::BaseMatrix* result, linearalgebra::BaseMatrix* J, SReal fact) override;
@@ -120,8 +119,8 @@ protected :
 private :
     bool first;
     unsigned systemSize;
-    double dt;
-    double factInt;
+    Real dt;
+    Real factInt;
     std::vector<bool> isActiveDofs;
 };