From 6d68c6f48c40eebddd27145bd148924ac9c4e965 Mon Sep 17 00:00:00 2001 From: Alex Bilger Date: Thu, 2 Apr 2026 17:02:20 +0200 Subject: [PATCH] add some unit tests --- .../Direct/tests/SparseLDLSolver_test.cpp | 300 ++++++++++++++++++ 1 file changed, 300 insertions(+) diff --git a/Sofa/Component/LinearSolver/Direct/tests/SparseLDLSolver_test.cpp b/Sofa/Component/LinearSolver/Direct/tests/SparseLDLSolver_test.cpp index 3d72df3da06..d00efa006f4 100644 --- a/Sofa/Component/LinearSolver/Direct/tests/SparseLDLSolver_test.cpp +++ b/Sofa/Component/LinearSolver/Direct/tests/SparseLDLSolver_test.cpp @@ -21,6 +21,7 @@ ******************************************************************************/ #include #include +#include #include #include #include @@ -28,6 +29,8 @@ #include #include +#include +#include TEST(SparseLDLSolver, EmptySystem) @@ -140,3 +143,300 @@ TEST(SparseLDLSolver, AssociatedLinearSystem) EXPECT_EQ(MatrixSystem::GetCustomTemplateName(), MatrixType::Name()); } + +TEST(SparseLDLSolver, Scalar1x1) +{ + using MatrixType = sofa::linearalgebra::CompressedRowSparseMatrix; + using VectorType = sofa::linearalgebra::FullVector; + using Solver = sofa::component::linearsolver::direct::SparseLDLSolver; + const Solver::SPtr solver = sofa::core::objectmodel::New(); + + // Explicitly set NaturalOrderingMethod to avoid uninitialized permutation + using NaturalOrdering = sofa::component::linearsolver::ordering::NaturalOrderingMethod; + const NaturalOrdering::SPtr ordering = sofa::core::objectmodel::New(); + solver->l_orderingMethod.set(ordering.get()); + + MatrixType A(1, 1); + A.add(0, 0, 2.0); + A.compress(); + + VectorType b(1), x(1); + b[0] = 4.0; + + solver->init(); + solver->invert(A); + solver->solve(A, x, b); + + msg_info("Test") << "x[0] = " << x[0]; + std::cout << "[DEBUG_LOG] Scalar1x1: x[0] = " << x[0] << std::endl; + EXPECT_NEAR(x[0], 2.0, 1e-10) << "Expected 2.0, got " << x[0]; +} + +TEST(SparseLDLSolver, IdentityMatrix) +{ + using MatrixType = sofa::linearalgebra::CompressedRowSparseMatrix; + using VectorType = sofa::linearalgebra::FullVector; + using Solver = sofa::component::linearsolver::direct::SparseLDLSolver; + const Solver::SPtr solver = sofa::core::objectmodel::New(); + + // Explicitly set NaturalOrderingMethod + using NaturalOrdering = sofa::component::linearsolver::ordering::NaturalOrderingMethod; + const NaturalOrdering::SPtr ordering = sofa::core::objectmodel::New(); + solver->l_orderingMethod.set(ordering.get()); + + constexpr int n = 5; + MatrixType A(n, n); + for (int i = 0; i < n; ++i) + { + A.add(i, i, 1.0); + } + A.compress(); + + VectorType b(n), x(n); + for (int i = 0; i < n; ++i) + { + b[i] = (SReal)(i + 1); + } + + solver->init(); + solver->invert(A); + solver->solve(A, x, b); + + for (int i = 0; i < n; ++i) + { + EXPECT_NEAR(x[i], b[i], 1e-10) << "At index " << i << ", expected " << b[i] << ", got " << x[i]; + } +} + +TEST(SparseLDLSolver, BlockMatrix3x3) +{ + using Block = sofa::type::Mat<3, 3, SReal>; + using MatrixType = sofa::linearalgebra::CompressedRowSparseMatrix; + using VectorType = sofa::linearalgebra::FullVector; + using Solver = sofa::component::linearsolver::direct::SparseLDLSolver; + const Solver::SPtr solver = sofa::core::objectmodel::New(); + + // Explicitly set NaturalOrderingMethod + using NaturalOrdering = sofa::component::linearsolver::ordering::NaturalOrderingMethod; + const NaturalOrdering::SPtr ordering = sofa::core::objectmodel::New(); + solver->l_orderingMethod.set(ordering.get()); + + constexpr int nBlocks = 2; + MatrixType A(nBlocks * 3, nBlocks * 3); + + // [B0 0] + // [0 B1] + // where B0 = [1 0 0] + // [0 1 0] + // [0 0 1] + // and B1 = [2 0 0] + // [0 4 0] + // [0 0 8] + + // Block 0: Identity + Block B0; B0.identity(); + A.add(0, 0, B0); + + // Block 1: Diagonal + Block B1; B1.clear(); + B1[0][0] = 2.0; B1[1][1] = 4.0; B1[2][2] = 8.0; + A.add(3, 3, B1); + + A.compress(); + + // The system size is nBlocks * 3 = 6 + VectorType b(6), x(6); + b[0] = 1.0; b[1] = 2.0; b[2] = 3.0; + b[3] = 4.0; b[4] = 8.0; b[5] = 16.0; + + solver->init(); + solver->invert(A); + solver->solve(A, x, b); + + // Expected x = [1, 2, 3, 2, 2, 2] + + EXPECT_NEAR(x[0], 1.0, 1e-10); + EXPECT_NEAR(x[1], 2.0, 1e-10); + EXPECT_NEAR(x[2], 3.0, 1e-10); + + EXPECT_NEAR(x[3], 2.0, 1e-10); + EXPECT_NEAR(x[4], 2.0, 1e-10); + EXPECT_NEAR(x[5], 2.0, 1e-10); +} + +TEST(SparseLDLSolver, BlockMatrix3x3_NonDiagonal) +{ + using Block = sofa::type::Mat<3, 3, SReal>; + using MatrixType = sofa::linearalgebra::CompressedRowSparseMatrix; + using VectorType = sofa::linearalgebra::FullVector; + using Solver = sofa::component::linearsolver::direct::SparseLDLSolver; + const Solver::SPtr solver = sofa::core::objectmodel::New(); + + // Explicitly set NaturalOrderingMethod + using NaturalOrdering = sofa::component::linearsolver::ordering::NaturalOrderingMethod; + const NaturalOrdering::SPtr ordering = sofa::core::objectmodel::New(); + solver->l_orderingMethod.set(ordering.get()); + + constexpr int nBlocks = 1; + MatrixType A(nBlocks * 3, nBlocks * 3); + + // Symmetric 3x3 block + // [ 2 1 0 ] + // [ 1 2 1 ] + // [ 0 1 2 ] + Block B; + B.clear(); + B[0][0] = 2.0; B[0][1] = 1.0; B[0][2] = 0.0; + B[1][0] = 1.0; B[1][1] = 2.0; B[1][2] = 1.0; + B[2][0] = 0.0; B[2][1] = 1.0; B[2][2] = 2.0; + A.addBlock(0, 0, B); + + A.compress(); + + VectorType b(3), x(3); + // b = A * [1, 1, 1] = [3, 4, 3] + b[0] = 3.0; b[1] = 4.0; b[2] = 3.0; + + solver->init(); + solver->invert(A); + solver->solve(A, x, b); + + EXPECT_NEAR(x[0], 1.0, 1e-10); + EXPECT_NEAR(x[1], 1.0, 1e-10); + EXPECT_NEAR(x[2], 1.0, 1e-10); +} + +TEST(SparseLDLSolver, DiagonalMatrix) +{ + using MatrixType = sofa::linearalgebra::CompressedRowSparseMatrix; + using VectorType = sofa::linearalgebra::FullVector; + using Solver = sofa::component::linearsolver::direct::SparseLDLSolver; + const Solver::SPtr solver = sofa::core::objectmodel::New(); + + const int n = 5; + MatrixType A(n, n); + for (int i = 0; i < n; ++i) + A.add(i, i, (SReal)(i + 1)); + A.compress(); + + VectorType b(n), x(n); + for (int i = 0; i < n; ++i) + b[i] = 1.0; + + solver->init(); + solver->invert(A); + solver->solve(A, x, b); + + for (int i = 0; i < n; ++i) + EXPECT_NEAR(x[i], 1.0 / (i + 1), 1e-10); +} + +TEST(SparseLDLSolver, SimpleSPD2x2) +{ + using MatrixType = sofa::linearalgebra::CompressedRowSparseMatrix; + using VectorType = sofa::linearalgebra::FullVector; + using Solver = sofa::component::linearsolver::direct::SparseLDLSolver; + const Solver::SPtr solver = sofa::core::objectmodel::New(); + + // A = [ 2 -1 ] + // [ -1 2 ] + MatrixType A(2, 2); + A.add(0, 0, 2.0); + A.add(0, 1, -1.0); + A.add(1, 0, -1.0); + A.add(1, 1, 2.0); + A.compress(); + + VectorType b(2), x(2); + b[0] = 1.0; + b[1] = 0.0; + // Expected x = [2/3, 1/3] + + solver->init(); + solver->invert(A); + solver->solve(A, x, b); + + EXPECT_NEAR(x[0], 2.0/3.0, 1e-10); + EXPECT_NEAR(x[1], 1.0/3.0, 1e-10); +} + +TEST(SparseLDLSolver, SimpleSPD3x3) +{ + using MatrixType = sofa::linearalgebra::CompressedRowSparseMatrix; + using VectorType = sofa::linearalgebra::FullVector; + using Solver = sofa::component::linearsolver::direct::SparseLDLSolver; + const Solver::SPtr solver = sofa::core::objectmodel::New(); + + // A = [ 4 12 -16 ] + // [ 12 37 -43 ] + // [-16 -43 98 ] + MatrixType A(3, 3); + A.add(0, 0, 4.0); A.add(0, 1, 12.0); A.add(0, 2, -16.0); + A.add(1, 0, 12.0); A.add(1, 1, 37.0); A.add(1, 2, -43.0); + A.add(2, 0, -16.0); A.add(2, 1, -43.0); A.add(2, 2, 98.0); + A.compress(); + + VectorType b(3), x(3); + b[0] = 1.0; b[1] = 2.0; b[2] = 3.0; + + solver->init(); + solver->invert(A); + solver->solve(A, x, b); + + // Verify A*x = b + VectorType Ax(3); + for(int i=0; i<3; ++i) { + Ax[i] = 0; + for(int j=0; j<3; ++j) Ax[i] += A(i,j) * x[j]; + } + + EXPECT_NEAR(Ax[0], b[0], 1e-10); + EXPECT_NEAR(Ax[1], b[1], 1e-10); + EXPECT_NEAR(Ax[2], b[2], 1e-10); +} + +TEST(SparseLDLSolver, MultiStepFactorization) +{ + using MatrixType = sofa::linearalgebra::CompressedRowSparseMatrix; + using VectorType = sofa::linearalgebra::FullVector; + using Solver = sofa::component::linearsolver::direct::SparseLDLSolver; + const Solver::SPtr solver = sofa::core::objectmodel::New(); + + const int n = 3; + MatrixType A(n, n); + A.add(0, 0, 2.0); A.add(0, 1, -1.0); + A.add(1, 0, -1.0); A.add(1, 1, 2.0); A.add(1, 2, -1.0); + A.add(2, 1, -1.0); A.add(2, 2, 2.0); + A.compress(); + + VectorType b(n), x(n); + for(int i=0; iinit(); + + // First solve + solver->invert(A); + solver->solve(A, x, b); + + // Change values but NOT shape + { + auto& values = const_cast(A.getColsValue()); + values[0] = 4.0; // A(0,0) + } + solver->invert(A); + solver->solve(A, x, b); + + VectorType Ax(n); + for(int i=0; i