[LinearSolver] Untangle linear solver and linear system

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;
 };