Skip to content

[LinearSolver.Direct] Add some unit tests #6068

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
alxbilger wants to merge 1 commit into
base: master
Choose a base branch
from
+300 −0
Open

[LinearSolver.Direct] Add some unit tests #6068

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
300 changes: 300 additions & 0 deletions Sofa/Component/LinearSolver/Direct/tests/SparseLDLSolver_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,13 +21,16 @@
******************************************************************************/
#include <sofa/testing/BaseTest.h>
#include <sofa/component/linearsolver/direct/SparseLDLSolver.h>
#include <sofa/component/linearsolver/ordering/NaturalOrderingMethod.h>
#include <sofa/component/linearsolver/direct/SparseCommon.h>
#include <sofa/component/linearsystem/MatrixLinearSystem.h>
#include <sofa/simulation/Node.h>
#include <sofa/simulation/graph/DAGSimulation.h>
#include <sofa/simpleapi/SimpleApi.h>

#include <sofa/testing/NumericTest.h>
#include <sofa/type/Mat.h>
#include <sofa/type/Vec.h>


TEST(SparseLDLSolver, EmptySystem)
Expand Down Expand Up @@ -140,3 +143,300 @@ TEST(SparseLDLSolver, AssociatedLinearSystem)

EXPECT_EQ(MatrixSystem::GetCustomTemplateName(), MatrixType::Name());
}

TEST(SparseLDLSolver, Scalar1x1)
{
using MatrixType = sofa::linearalgebra::CompressedRowSparseMatrix<SReal>;
using VectorType = sofa::linearalgebra::FullVector<SReal>;
using Solver = sofa::component::linearsolver::direct::SparseLDLSolver<MatrixType, VectorType>;
const Solver::SPtr solver = sofa::core::objectmodel::New<Solver>();

// Explicitly set NaturalOrderingMethod to avoid uninitialized permutation
using NaturalOrdering = sofa::component::linearsolver::ordering::NaturalOrderingMethod;
const NaturalOrdering::SPtr ordering = sofa::core::objectmodel::New<NaturalOrdering>();
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<SReal>;
using VectorType = sofa::linearalgebra::FullVector<SReal>;
using Solver = sofa::component::linearsolver::direct::SparseLDLSolver<MatrixType, VectorType>;
const Solver::SPtr solver = sofa::core::objectmodel::New<Solver>();

// Explicitly set NaturalOrderingMethod
using NaturalOrdering = sofa::component::linearsolver::ordering::NaturalOrderingMethod;
const NaturalOrdering::SPtr ordering = sofa::core::objectmodel::New<NaturalOrdering>();
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<Block>;
using VectorType = sofa::linearalgebra::FullVector<SReal>;
using Solver = sofa::component::linearsolver::direct::SparseLDLSolver<MatrixType, VectorType>;
const Solver::SPtr solver = sofa::core::objectmodel::New<Solver>();

// Explicitly set NaturalOrderingMethod
using NaturalOrdering = sofa::component::linearsolver::ordering::NaturalOrderingMethod;
const NaturalOrdering::SPtr ordering = sofa::core::objectmodel::New<NaturalOrdering>();
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<Block>;
using VectorType = sofa::linearalgebra::FullVector<SReal>;
using Solver = sofa::component::linearsolver::direct::SparseLDLSolver<MatrixType, VectorType>;
const Solver::SPtr solver = sofa::core::objectmodel::New<Solver>();

// Explicitly set NaturalOrderingMethod
using NaturalOrdering = sofa::component::linearsolver::ordering::NaturalOrderingMethod;
const NaturalOrdering::SPtr ordering = sofa::core::objectmodel::New<NaturalOrdering>();
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<SReal>;
using VectorType = sofa::linearalgebra::FullVector<SReal>;
using Solver = sofa::component::linearsolver::direct::SparseLDLSolver<MatrixType, VectorType>;
const Solver::SPtr solver = sofa::core::objectmodel::New<Solver>();

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<SReal>;
using VectorType = sofa::linearalgebra::FullVector<SReal>;
using Solver = sofa::component::linearsolver::direct::SparseLDLSolver<MatrixType, VectorType>;
const Solver::SPtr solver = sofa::core::objectmodel::New<Solver>();

// 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<SReal>;
using VectorType = sofa::linearalgebra::FullVector<SReal>;
using Solver = sofa::component::linearsolver::direct::SparseLDLSolver<MatrixType, VectorType>;
const Solver::SPtr solver = sofa::core::objectmodel::New<Solver>();

// 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<SReal>;
using VectorType = sofa::linearalgebra::FullVector<SReal>;
using Solver = sofa::component::linearsolver::direct::SparseLDLSolver<MatrixType, VectorType>;
const Solver::SPtr solver = sofa::core::objectmodel::New<Solver>();

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; i<n; ++i) b[i] = 1.0;

solver->init();

// First solve
solver->invert(A);
solver->solve(A, x, b);

// Change values but NOT shape
{
auto& values = const_cast<MatrixType::VecBlock&>(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<n; ++i)
{
Ax[i] = 0;
for(int j=0; j<n; ++j)
{
Ax[i] += A(i,j) * x[j];
}
}
for(int i=0; i<n; ++i)
{
EXPECT_NEAR(Ax[i], b[i], 1e-10);
}
}
Loading