diff --git a/Sofa/Component/LinearSolver/Direct/tests/SparseLDLSolver_test.cpp b/Sofa/Component/LinearSolver/Direct/tests/SparseLDLSolver_test.cpp
index 3d72df3da06..882c0cb34af 100644
--- a/Sofa/Component/LinearSolver/Direct/tests/SparseLDLSolver_test.cpp
+++ b/Sofa/Component/LinearSolver/Direct/tests/SparseLDLSolver_test.cpp
@@ -19,14 +19,14 @@
 *                                                                             *
 * Contact information: contact@sofa-framework.org                             *
 ******************************************************************************/
-#include <sofa/testing/BaseTest.h>
-#include <sofa/component/linearsolver/direct/SparseLDLSolver.h>
 #include <sofa/component/linearsolver/direct/SparseCommon.h>
+#include <sofa/component/linearsolver/direct/SparseLDLSolver.h>
 #include <sofa/component/linearsystem/MatrixLinearSystem.h>
+#include <sofa/helper/RandomGenerator.h>
+#include <sofa/simpleapi/SimpleApi.h>
 #include <sofa/simulation/Node.h>
 #include <sofa/simulation/graph/DAGSimulation.h>
-#include <sofa/simpleapi/SimpleApi.h>
-
+#include <sofa/testing/BaseTest.h>
 #include <sofa/testing/NumericTest.h>
 
 
@@ -140,3 +140,76 @@ TEST(SparseLDLSolver, AssociatedLinearSystem)
 
     EXPECT_EQ(MatrixSystem::GetCustomTemplateName(), MatrixType::Name());
 }
+
+
+TEST(SparseLDLSolver, TestInvertingRandomMatrix)
+{
+    using MatrixType = sofa::linearalgebra::CompressedRowSparseMatrix<SReal>;
+    using Solver = sofa::component::linearsolver::direct::SparseLDLSolver<MatrixType, sofa::linearalgebra::FullVector<SReal> >;
+    const Solver::SPtr solver = sofa::core::objectmodel::New<Solver>();
+
+    solver->init();
+
+    unsigned nbRows = 300;
+    unsigned nbCols = 300;
+    SReal reg = 5;
+    float sparsity = 0.05;
+    const auto nbNonZero = static_cast<sofa::SignedIndex>(sparsity * static_cast<SReal>(nbRows*nbCols));
+
+
+    sofa::linearalgebra::FullMatrix<SReal> tempMatrix(nbRows,nbCols), finalTempMatrix;
+    tempMatrix.clear();
+
+    sofa::helper::RandomGenerator randomGenerator;
+    randomGenerator.initSeed(2807);
+
+
+    for (sofa::SignedIndex i = 0; i < nbNonZero; ++i)
+    {
+        const auto value = static_cast<SReal>(fabs(sofa::helper::drand(2)));
+        const auto row = randomGenerator.random<sofa::Index>(0, nbRows);
+        const auto col = randomGenerator.random<sofa::Index>(0, nbCols);
+        tempMatrix.set(row,col,value);
+    }
+
+    for (sofa::SignedIndex i = 0; i < nbRows; ++i)
+    {
+        tempMatrix.set(i,i,tempMatrix(i,i) + reg);
+    }
+    tempMatrix.mulT(finalTempMatrix, tempMatrix);
+
+    sofa::linearalgebra::CompressedRowSparseMatrix<SReal> matrix;
+    matrix.resize(nbRows, nbCols);
+
+    unsigned nbNZ = 0;
+    for (unsigned i=0; i<nbRows; ++i)
+    {
+        for (unsigned j=0; j<nbCols; ++j)
+        {
+            if (finalTempMatrix(i,j) > 1e-8)
+                matrix.set(i,j,finalTempMatrix(i,j));
+            else
+                ++nbNZ;
+        }
+
+    }
+    msg_info("TestInvertingRandomMatrix") << "REAL SPARSITY (#zeros/#elements) : " <<(nbNZ)/static_cast<double>(nbRows*nbCols) ;
+    matrix.compress();
+
+    sofa::linearalgebra::FullVector<SReal> known(nbCols), unknown(nbCols), rhs(nbRows);
+    for (unsigned i=0; i<nbRows; ++i)
+    {
+        const auto value = static_cast<SReal>(sofa::helper::drand(1));
+        known.set(i,value);
+    }
+    matrix.mul(rhs, known);
+
+    solver->invert(matrix);
+    solver->solve(matrix,  unknown, rhs);
+
+    for (unsigned i=0; i<nbCols; ++i)
+    {
+        EXPECT_NEAR(unknown[i], known[i], 1e-12);
+    }
+
+}