From acadc261a9637c7e45ac77fbae650fe15a0efce3 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 27 Feb 2026 11:29:48 +0900 Subject: [PATCH 01/11] fix substraction operation (copy paste?) --- .../LinearAlgebra/src/sofa/linearalgebra/SparseMatrix.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrix.h b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrix.h index 7c0e4e169d2..71fa262a410 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrix.h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/SparseMatrix.h @@ -357,9 +357,9 @@ class SparseMatrix : public linearalgebra::BaseMatrix } template - MatrixExpr< MatrixAddition< SparseMatrix, SparseMatrix > > operator-(const SparseMatrix& m) const + MatrixExpr< MatrixSubtraction< SparseMatrix, SparseMatrix > > operator-(const SparseMatrix& m) const { - return MatrixExpr { MatrixAddition< SparseMatrix, SparseMatrix >(*this, m) }; + return MatrixExpr { MatrixSubtraction< SparseMatrix, SparseMatrix >(*this, m) }; } void swap(SparseMatrix& m) From 6558a85df0c465c1e60621d45341f37d2cbedb5f Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 27 Feb 2026 11:31:12 +0900 Subject: [PATCH 02/11] correct test in CRS::resizeBlock --- .../src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h index 35f38db1dcf..119b451abb8 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h @@ -279,7 +279,7 @@ public : virtual void resizeBlock(Index nbBRow, Index nbBCol) { - if (nBlockRow == nbBRow && nBlockRow == nbBCol) + if (nBlockRow == nbBRow && nBlockCol == nbBCol) { /// Just clear the matrix for (Index i = 0; i < static_cast(colsValue.size()); ++i) From 0022eccfba495705dd6221300efdb35bd40b37eb Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 27 Feb 2026 12:14:23 +0900 Subject: [PATCH 03/11] fix access indexing in blockfullmatrix --- .../LinearAlgebra/src/sofa/linearalgebra/BlockFullMatrix.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/BlockFullMatrix.h b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/BlockFullMatrix.h index b174d29f7b5..4fb1ed3d493 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/BlockFullMatrix.h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/BlockFullMatrix.h @@ -180,7 +180,7 @@ class BlockFullMatrix : public linearalgebra::BaseMatrix Real r = 0; for (Index j=0; j Date: Tue, 3 Mar 2026 13:10:08 +0900 Subject: [PATCH 04/11] add unit tests (with expected fails from clear()) --- .../test/BlockFullMatrix_test.cpp | 502 ++++++++++++++++++ .../LinearAlgebra/test/CMakeLists.txt | 1 + 2 files changed, 503 insertions(+) create mode 100644 Sofa/framework/LinearAlgebra/test/BlockFullMatrix_test.cpp diff --git a/Sofa/framework/LinearAlgebra/test/BlockFullMatrix_test.cpp b/Sofa/framework/LinearAlgebra/test/BlockFullMatrix_test.cpp new file mode 100644 index 00000000000..64446c9a9cf --- /dev/null +++ b/Sofa/framework/LinearAlgebra/test/BlockFullMatrix_test.cpp @@ -0,0 +1,502 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#include +#include + +#include + +#include + +namespace sofa +{ + +using BFM2 = linearalgebra::BlockFullMatrix<2, double>; +using Block2 = BFM2::Block; +using Vec2 = type::Vec<2, double>; + +constexpr double tol = 1e-10; + +// ============================================================ +// Construction & Dimensions +// ============================================================ + +TEST(BlockFullMatrix, DefaultConstruction) +{ + BFM2 m; + EXPECT_EQ(m.rowSize(), 0); + EXPECT_EQ(m.colSize(), 0); + EXPECT_EQ(m.ptr(), nullptr); +} + +TEST(BlockFullMatrix, ConstructWithDimensions) +{ + BFM2 m(4, 4); + EXPECT_EQ(m.rowSize(), 4); + EXPECT_EQ(m.colSize(), 4); + EXPECT_NE(m.ptr(), nullptr); +} + +TEST(BlockFullMatrix, Resize) +{ + // Use nBCol >= 4 so that clear() (which iterates 3*nBRow) doesn't go OOB. + // With BSIZE=2, 8 cols → nBCol=4. + BFM2 m; + m.resize(4, 8); + EXPECT_EQ(m.rowSize(), 4); + EXPECT_EQ(m.colSize(), 8); + EXPECT_NE(m.ptr(), nullptr); + + m.resize(6, 8); + EXPECT_EQ(m.rowSize(), 6); + EXPECT_EQ(m.colSize(), 8); +} + +TEST(BlockFullMatrix, ResizeClearsData) +{ + // Use 4x8 (nBRow=2, nBCol=4) so clear() under-clears rather than OOB. + // clear() iterates 3*nBRow=6 blocks out of nBRow*nBCol=8. + // Blocks at indices 6,7 (block row 1, columns 2-3) won't be cleared. + BFM2 m(4, 8); + m.set(0, 0, 7.0); + m.set(3, 7, 9.0); // block (1,3) = data[7], won't be cleared by buggy clear() + + m.resize(4, 8); + + for (SignedIndex i = 0; i < 4; ++i) + for (SignedIndex j = 0; j < 8; ++j) + EXPECT_NEAR(m.element(i, j), 0.0, tol) + << "element(" << i << "," << j << ") not zero after resize"; +} + +// ============================================================ +// Element Access (set / element / add / clear) +// ============================================================ + +TEST(BlockFullMatrix, SetAndElement) +{ + BFM2 m(4, 4); + + // element in block (0,0) + m.set(0, 1, 1.5); + EXPECT_NEAR(m.element(0, 1), 1.5, tol); + + // element in block (1,0) + m.set(2, 0, -3.0); + EXPECT_NEAR(m.element(2, 0), -3.0, tol); + + // element in block (0,1) + m.set(1, 3, 4.25); + EXPECT_NEAR(m.element(1, 3), 4.25, tol); + + // element in block (1,1) + m.set(3, 3, 100.0); + EXPECT_NEAR(m.element(3, 3), 100.0, tol); +} + +TEST(BlockFullMatrix, Add) +{ + BFM2 m(4, 4); + m.set(0, 0, 3.0); + m.add(0, 0, 2.0); + EXPECT_NEAR(m.element(0, 0), 5.0, tol); +} + +TEST(BlockFullMatrix, ClearElement) +{ + BFM2 m(4, 4); + m.set(1, 1, 5.0); + m.set(1, 0, 2.0); + m.clear(1, 1); + EXPECT_NEAR(m.element(1, 1), 0.0, tol); + EXPECT_NEAR(m.element(1, 0), 2.0, tol); +} + +// ============================================================ +// Block Access +// ============================================================ + +TEST(BlockFullMatrix, BlocAccess) +{ + BFM2 m(4, 4); + + // Write through block reference + Block2& blk = m.bloc(0, 1); + blk(0, 0) = 10.0; + blk(1, 1) = 20.0; + + // Read via element() + EXPECT_NEAR(m.element(0, 2), 10.0, tol); + EXPECT_NEAR(m.element(1, 3), 20.0, tol); +} + +TEST(BlockFullMatrix, SubAndAsub) +{ + BFM2 m(4, 4); + m.set(2, 2, 42.0); + + // sub() takes element indices, returns block reference + const Block2& s = m.sub(2, 2, 2, 2); + EXPECT_NEAR(s.element(0, 0), 42.0, tol); + + // asub() takes block indices + const Block2& a = m.asub(1, 1, 2, 2); + EXPECT_NEAR(a.element(0, 0), 42.0, tol); +} + +TEST(BlockFullMatrix, GetSetSubMatrix) +{ + BFM2 m(4, 4); + Block2 input; + input.set(0, 0, 1.0); + input.set(0, 1, 2.0); + input.set(1, 0, 3.0); + input.set(1, 1, 4.0); + + m.setSubMatrix(2, 2, 2, 2, input); + + Block2 output; + m.getSubMatrix(2, 2, 2, 2, output); + + for (SignedIndex i = 0; i < 2; ++i) + for (SignedIndex j = 0; j < 2; ++j) + EXPECT_NEAR(output.element(i, j), input.element(i, j), tol); +} + +// ============================================================ +// Clearing Operations +// ============================================================ + +TEST(BlockFullMatrix, ClearAll) +{ + // 4x8 matrix with BSIZE=2 => nBRow=2, nBCol=4. + // clear() iterates 3*nBRow=6 instead of nBRow*nBCol=8, so the last 2 blocks + // (block row 1, columns 2-3) are NOT cleared. This exposes the bug without + // causing an out-of-bounds write (which would happen with nBCol < 3). + BFM2 m(4, 8); + m.set(0, 0, 1.0); // block (0,0) = data[0] + m.set(0, 2, 2.0); // block (0,1) = data[1] + m.set(0, 4, 3.0); // block (0,2) = data[2] + m.set(0, 6, 4.0); // block (0,3) = data[3] + m.set(2, 0, 5.0); // block (1,0) = data[4] + m.set(2, 2, 6.0); // block (1,1) = data[5] + m.set(2, 4, 7.0); // block (1,2) = data[6] — NOT cleared by buggy code + m.set(2, 6, 8.0); // block (1,3) = data[7] — NOT cleared by buggy code + + m.clear(); + + for (SignedIndex i = 0; i < 4; ++i) + for (SignedIndex j = 0; j < 8; ++j) + EXPECT_NEAR(m.element(i, j), 0.0, tol) + << "element(" << i << "," << j << ") not zero after clear()"; +} + +TEST(BlockFullMatrix, ClearRow) +{ + BFM2 m(4, 4); + for (SignedIndex i = 0; i < 4; ++i) + for (SignedIndex j = 0; j < 4; ++j) + m.set(i, j, static_cast(i * 4 + j + 1)); + + m.clearRow(1); + + for (SignedIndex j = 0; j < 4; ++j) + EXPECT_NEAR(m.element(1, j), 0.0, tol) + << "row 1, col " << j << " should be zero"; + + // Other rows unchanged + EXPECT_NEAR(m.element(0, 0), 1.0, tol); + EXPECT_NEAR(m.element(2, 2), 11.0, tol); + EXPECT_NEAR(m.element(3, 3), 16.0, tol); +} + +TEST(BlockFullMatrix, ClearCol) +{ + BFM2 m(4, 4); + for (SignedIndex i = 0; i < 4; ++i) + for (SignedIndex j = 0; j < 4; ++j) + m.set(i, j, static_cast(i * 4 + j + 1)); + + m.clearCol(2); + + for (SignedIndex i = 0; i < 4; ++i) + EXPECT_NEAR(m.element(i, 2), 0.0, tol) + << "row " << i << ", col 2 should be zero"; + + // Other columns unchanged + EXPECT_NEAR(m.element(0, 0), 1.0, tol); + EXPECT_NEAR(m.element(1, 1), 6.0, tol); + EXPECT_NEAR(m.element(3, 3), 16.0, tol); +} + +TEST(BlockFullMatrix, ClearRowCol) +{ + BFM2 m(4, 4); + for (SignedIndex i = 0; i < 4; ++i) + for (SignedIndex j = 0; j < 4; ++j) + m.set(i, j, static_cast(i * 4 + j + 1)); + + m.clearRowCol(1); + + // Row 1 all zeros + for (SignedIndex j = 0; j < 4; ++j) + EXPECT_NEAR(m.element(1, j), 0.0, tol); + + // Col 1 all zeros + for (SignedIndex i = 0; i < 4; ++i) + EXPECT_NEAR(m.element(i, 1), 0.0, tol); + + // Other elements unchanged (except those on row/col 1) + EXPECT_NEAR(m.element(0, 0), 1.0, tol); + EXPECT_NEAR(m.element(2, 2), 11.0, tol); + EXPECT_NEAR(m.element(3, 3), 16.0, tol); +} + +// ============================================================ +// Matrix-Vector Multiplication +// ============================================================ + +TEST(BlockFullMatrix, IdentityTimesVector) +{ + BFM2 m(4, 4); + // Set to identity + for (SignedIndex i = 0; i < 4; ++i) + m.set(i, i, 1.0); + + linearalgebra::FullVector v(4); + v[0] = 1.0; v[1] = 2.0; v[2] = 3.0; v[3] = 4.0; + + linearalgebra::FullVector res = m * v; + + for (SignedIndex i = 0; i < 4; ++i) + EXPECT_NEAR(res[i], v[i], tol); +} + +TEST(BlockFullMatrix, GeneralProduct) +{ + // 4x4 matrix, BSIZE=2 + // M = [1 2 0 0] v = [1] + // [3 4 0 0] [2] + // [0 0 5 6] [3] + // [0 0 7 8] [4] + // + // M*v = [1*1+2*2, 3*1+4*2, 5*3+6*4, 7*3+8*4] = [5, 11, 39, 53] + BFM2 m(4, 4); + m.set(0, 0, 1.0); m.set(0, 1, 2.0); + m.set(1, 0, 3.0); m.set(1, 1, 4.0); + m.set(2, 2, 5.0); m.set(2, 3, 6.0); + m.set(3, 2, 7.0); m.set(3, 3, 8.0); + + linearalgebra::FullVector v(4); + v[0] = 1.0; v[1] = 2.0; v[2] = 3.0; v[3] = 4.0; + + linearalgebra::FullVector res = m * v; + + EXPECT_NEAR(res[0], 5.0, tol); + EXPECT_NEAR(res[1], 11.0, tol); + EXPECT_NEAR(res[2], 39.0, tol); + EXPECT_NEAR(res[3], 53.0, tol); +} + +TEST(BlockFullMatrix, SingleBlockColumn) +{ + // Matrix with nBCol=1: 4 rows, 2 cols (BSIZE=2 => nBRow=2, nBCol=1) + // This tests the split-loop boundary in operator* + // M = [1 2] v = [3] + // [3 4] [5] + // [5 6] + // [7 8] + // + // M*v = [1*3+2*5, 3*3+4*5, 5*3+6*5, 7*3+8*5] = [13, 29, 45, 61] + BFM2 m(4, 2); + m.set(0, 0, 1.0); m.set(0, 1, 2.0); + m.set(1, 0, 3.0); m.set(1, 1, 4.0); + m.set(2, 0, 5.0); m.set(2, 1, 6.0); + m.set(3, 0, 7.0); m.set(3, 1, 8.0); + + linearalgebra::FullVector v(2); + v[0] = 3.0; v[1] = 5.0; + + linearalgebra::FullVector res = m * v; + + EXPECT_NEAR(res[0], 13.0, tol); + EXPECT_NEAR(res[1], 29.0, tol); + EXPECT_NEAR(res[2], 45.0, tol); + EXPECT_NEAR(res[3], 61.0, tol); +} + +// ============================================================ +// Block Nested Class +// ============================================================ + +TEST(BlockFullMatrix, BlockNrowsNcols) +{ + Block2 b; + EXPECT_EQ(b.Nrows(), 2); + EXPECT_EQ(b.Ncols(), 2); +} + +TEST(BlockFullMatrix, BlockSetAddElement) +{ + Block2 b; + b.clear(); + b.set(0, 0, 3.0); + EXPECT_NEAR(b.element(0, 0), 3.0, tol); + + b.add(0, 0, 2.0); + EXPECT_NEAR(b.element(0, 0), 5.0, tol); + + b.set(1, 0, -1.0); + EXPECT_NEAR(b.element(1, 0), -1.0, tol); +} + +TEST(BlockFullMatrix, BlockTranspose) +{ + Block2 b; + b.clear(); + b.set(0, 0, 1.0); b.set(0, 1, 2.0); + b.set(1, 0, 3.0); b.set(1, 1, 4.0); + + auto tb = b.t(); + + // tb * v should compute b^T * v + // b^T = [1 3] v = [1] => [1*1+3*2, 2*1+4*2] = [7, 10] + // [2 4] [2] + Vec2 v(1.0, 2.0); + Vec2 r = tb * v; + + EXPECT_NEAR(r[0], 7.0, tol); + EXPECT_NEAR(r[1], 10.0, tol); +} + +TEST(BlockFullMatrix, BlockInverse) +{ + Block2 b; + b.clear(); + b.set(0, 0, 4.0); b.set(0, 1, 7.0); + b.set(1, 0, 2.0); b.set(1, 1, 6.0); + + Block2 inv = b.i(); + + // b * inv should be identity + auto product = b * inv; + EXPECT_NEAR(product(0, 0), 1.0, tol); + EXPECT_NEAR(product(0, 1), 0.0, tol); + EXPECT_NEAR(product(1, 0), 0.0, tol); + EXPECT_NEAR(product(1, 1), 1.0, tol); +} + +TEST(BlockFullMatrix, BlockNegation) +{ + Block2 b; + b.clear(); + b.set(0, 0, 1.0); b.set(0, 1, -2.0); + b.set(1, 0, 3.0); b.set(1, 1, -4.0); + + auto neg = -b; + + EXPECT_NEAR(neg(0, 0), -1.0, tol); + EXPECT_NEAR(neg(0, 1), 2.0, tol); + EXPECT_NEAR(neg(1, 0), -3.0, tol); + EXPECT_NEAR(neg(1, 1), 4.0, tol); +} + +TEST(BlockFullMatrix, BlockSubtraction) +{ + Block2 a; + a.clear(); + a.set(0, 0, 5.0); a.set(0, 1, 6.0); + a.set(1, 0, 7.0); a.set(1, 1, 8.0); + + Block2 b; + b.clear(); + b.set(0, 0, 1.0); b.set(0, 1, 2.0); + b.set(1, 0, 3.0); b.set(1, 1, 4.0); + + auto diff = a - b; + + EXPECT_NEAR(diff(0, 0), 4.0, tol); + EXPECT_NEAR(diff(0, 1), 4.0, tol); + EXPECT_NEAR(diff(1, 0), 4.0, tol); + EXPECT_NEAR(diff(1, 1), 4.0, tol); +} + +TEST(BlockFullMatrix, BlockMatVec) +{ + Block2 b; + b.clear(); + b.set(0, 0, 1.0); b.set(0, 1, 2.0); + b.set(1, 0, 3.0); b.set(1, 1, 4.0); + + Vec2 v(5.0, 6.0); + Vec2 r = b * v; + + // [1 2] * [5] = [1*5+2*6, 3*5+4*6] = [17, 39] + // [3 4] [6] + EXPECT_NEAR(r[0], 17.0, tol); + EXPECT_NEAR(r[1], 39.0, tol); +} + +// ============================================================ +// Name +// ============================================================ + +TEST(BlockFullMatrix, Name) +{ + std::string name = BFM2::Name(); + EXPECT_EQ(name, "BlockFullMatrix2d"); +} + +// ============================================================ +// Spot-Check with Production Type (BlockFullMatrix<6, SReal>) +// ============================================================ + +TEST(BlockFullMatrix, ProductionType) +{ + using BFM6 = linearalgebra::BlockFullMatrix<6, SReal>; + + BFM6 m(12, 12); + EXPECT_EQ(m.rowSize(), 12); + EXPECT_EQ(m.colSize(), 12); + + m.set(0, 0, 1.0); + m.set(5, 5, 2.0); + m.set(6, 6, 3.0); + m.set(11, 11, 4.0); + EXPECT_NEAR(m.element(0, 0), 1.0, tol); + EXPECT_NEAR(m.element(5, 5), 2.0, tol); + EXPECT_NEAR(m.element(6, 6), 3.0, tol); + EXPECT_NEAR(m.element(11, 11), 4.0, tol); + + // Identity matrix-vector product + linearalgebra::FullVector v(12); + for (SignedIndex i = 0; i < 12; ++i) + { + m.set(i, i, 1.0); + v[i] = static_cast(i + 1); + } + + linearalgebra::FullVector res = m * v; + for (SignedIndex i = 0; i < 12; ++i) + EXPECT_NEAR(res[i], v[i], tol); +} + +} // namespace sofa diff --git a/Sofa/framework/LinearAlgebra/test/CMakeLists.txt b/Sofa/framework/LinearAlgebra/test/CMakeLists.txt index b786e3df8a9..96c6ad817ec 100644 --- a/Sofa/framework/LinearAlgebra/test/CMakeLists.txt +++ b/Sofa/framework/LinearAlgebra/test/CMakeLists.txt @@ -3,6 +3,7 @@ cmake_minimum_required(VERSION 3.22) project(Sofa.LinearAlgebra_test) set(SOURCE_FILES + BlockFullMatrix_test.cpp BTDMatrix_test.cpp BaseMatrix_test.cpp CompressedRowSparseMatrix_test.cpp From d3020df14a025fed109f8e0ab4adb3ad9dce2deb Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Tue, 3 Mar 2026 13:11:26 +0900 Subject: [PATCH 05/11] fix for clear() replacing 3*nBRow with nBRow*nBCol so that clear() iterates over the actual number of allocated blocks instead of a hardcoded factor of 3. --- .../LinearAlgebra/src/sofa/linearalgebra/BlockFullMatrix.inl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/BlockFullMatrix.inl b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/BlockFullMatrix.inl index 7974d3dc590..5c99bccdfdf 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/BlockFullMatrix.inl +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/BlockFullMatrix.inl @@ -280,7 +280,7 @@ void BlockFullMatrix::clearRowCol(Index i) template void BlockFullMatrix::clear() { - for (Index i=0; i<3*nBRow; ++i) + for (Index i=0; i Date: Tue, 3 Mar 2026 13:31:49 +0900 Subject: [PATCH 06/11] add unit tests for CompressedRowSparseMatrixGeneric The ResizeBlockNonSquare test specifically covers the regression for commit 6558a85df0 (the nBlockRow == nbBRow && nBlockCol == nbBCol fix), verifying that resizing with same rows but different cols (or vice versa) triggers a full reset rather than just zeroing values. --- .../LinearAlgebra/test/CMakeLists.txt | 1 + .../CompressedRowSparseMatrixGeneric_test.cpp | 648 ++++++++++++++++++ 2 files changed, 649 insertions(+) create mode 100644 Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixGeneric_test.cpp diff --git a/Sofa/framework/LinearAlgebra/test/CMakeLists.txt b/Sofa/framework/LinearAlgebra/test/CMakeLists.txt index 96c6ad817ec..bdcc27aac9f 100644 --- a/Sofa/framework/LinearAlgebra/test/CMakeLists.txt +++ b/Sofa/framework/LinearAlgebra/test/CMakeLists.txt @@ -7,6 +7,7 @@ set(SOURCE_FILES BTDMatrix_test.cpp BaseMatrix_test.cpp CompressedRowSparseMatrix_test.cpp + CompressedRowSparseMatrixGeneric_test.cpp CompressedRowSparseMatrixConstraint_test.cpp Matrix_test.cpp RotationMatrix_test.cpp diff --git a/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixGeneric_test.cpp b/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixGeneric_test.cpp new file mode 100644 index 00000000000..f11671c8805 --- /dev/null +++ b/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixGeneric_test.cpp @@ -0,0 +1,648 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#include +#include +#include + +using CRS = sofa::linearalgebra::CompressedRowSparseMatrixGeneric; +using Mat3 = sofa::type::Mat<3, 3, double>; +using CRSMat3 = sofa::linearalgebra::CompressedRowSparseMatrixGeneric; + +constexpr double kTol = 1e-10; + +// ==================== Construction & Sizing ==================== + +TEST(CompressedRowSparseMatrixGeneric, DefaultConstruction) +{ + CRS m; + EXPECT_EQ(m.rowBSize(), 0); + EXPECT_EQ(m.colBSize(), 0); + EXPECT_TRUE(m.getRowIndex().empty()); + EXPECT_TRUE(m.getRowBegin().empty()); + EXPECT_TRUE(m.getColsIndex().empty()); + EXPECT_TRUE(m.getColsValue().empty()); +} + +TEST(CompressedRowSparseMatrixGeneric, ConstructWithDimensions) +{ + CRS m(3, 4); + EXPECT_EQ(m.rowBSize(), 3); + EXPECT_EQ(m.colBSize(), 4); +} + +TEST(CompressedRowSparseMatrixGeneric, ResizeBlock) +{ + CRS m; + m.resizeBlock(3, 4); + EXPECT_EQ(m.rowBSize(), 3); + EXPECT_EQ(m.colBSize(), 4); + + m.resizeBlock(5, 6); + EXPECT_EQ(m.rowBSize(), 5); + EXPECT_EQ(m.colBSize(), 6); +} + +TEST(CompressedRowSparseMatrixGeneric, ResizeBlockSameSize) +{ + CRS m(3, 4); + m.setBlock(0, 0, 1.0); + m.setBlock(1, 2, 5.0); + m.compress(); + + const auto numValues = m.getColsValue().size(); + m.resizeBlock(3, 4); + + EXPECT_EQ(m.rowBSize(), 3); + EXPECT_EQ(m.colBSize(), 4); + EXPECT_EQ(m.getColsValue().size(), numValues); + for (const auto& v : m.getColsValue()) + { + EXPECT_NEAR(v, 0.0, kTol); + } +} + +TEST(CompressedRowSparseMatrixGeneric, ResizeBlockNonSquare) +{ + // Regression test for the nBlockRow==nbBRow && nBlockCol==nbBCol fix + CRS m(3, 4); + m.setBlock(0, 0, 1.0); + m.setBlock(1, 3, 7.0); + m.compress(); + + // Same rows, different cols: must fully reset, not just zero + m.resizeBlock(3, 6); + EXPECT_EQ(m.rowBSize(), 3); + EXPECT_EQ(m.colBSize(), 6); + EXPECT_TRUE(m.getRowIndex().empty()); + EXPECT_TRUE(m.getColsValue().empty()); + + // Same cols, different rows: must fully reset + CRS m2(3, 4); + m2.setBlock(0, 0, 2.0); + m2.compress(); + m2.resizeBlock(5, 4); + EXPECT_EQ(m2.rowBSize(), 5); + EXPECT_EQ(m2.colBSize(), 4); + EXPECT_TRUE(m2.getRowIndex().empty()); + EXPECT_TRUE(m2.getColsValue().empty()); +} + +// ==================== Block Insertion & Retrieval ==================== + +TEST(CompressedRowSparseMatrixGeneric, SetBlockAndBlock) +{ + CRS m(4, 4); + m.setBlock(1, 2, 3.14); + EXPECT_NEAR(m.block(1, 2), 3.14, kTol); +} + +TEST(CompressedRowSparseMatrixGeneric, AddBlock) +{ + CRS m(4, 4); + m.setBlock(0, 0, 3.0); + m.compress(); + m.addBlock(0, 0, 2.0); + EXPECT_NEAR(m.block(0, 0), 5.0, kTol); +} + +TEST(CompressedRowSparseMatrixGeneric, AddBlockMultipleToSamePosition) +{ + CRS m(4, 4); + m.addBlock(1, 0, 2.0); + m.addBlock(1, 1, 3.0); + m.addBlock(1, 0, 1.0); // separate btemp entry since last was (1,1) + m.compress(); + EXPECT_NEAR(m.block(1, 0), 3.0, kTol); // 2.0 + 1.0 + EXPECT_NEAR(m.block(1, 1), 3.0, kTol); +} + +TEST(CompressedRowSparseMatrixGeneric, WblockCreate) +{ + CRS m(4, 4); + double* ptr = m.wblock(2, 3, true); + ASSERT_NE(ptr, nullptr); + *ptr = 42.0; + EXPECT_NEAR(m.block(2, 3), 42.0, kTol); +} + +TEST(CompressedRowSparseMatrixGeneric, WblockNoCreate) +{ + CRS m(4, 4); + double* ptr = m.wblock(2, 3, false); + EXPECT_EQ(ptr, nullptr); +} + +TEST(CompressedRowSparseMatrixGeneric, GetBlock) +{ + CRS m(4, 4); + m.setBlock(1, 2, 7.5); + const double& ref = m.getBlock(1, 2); + EXPECT_NEAR(ref, 7.5, kTol); +} + +TEST(CompressedRowSparseMatrixGeneric, InsertionOrder) +{ + CRS m(6, 6); + m.setBlock(4, 1, 10.0); + m.setBlock(0, 3, 20.0); + m.setBlock(2, 0, 30.0); + m.setBlock(4, 5, 40.0); + m.setBlock(0, 0, 50.0); + + // AutoCompress on read verifies all blocks + EXPECT_NEAR(m.block(0, 0), 50.0, kTol); + EXPECT_NEAR(m.block(0, 3), 20.0, kTol); + EXPECT_NEAR(m.block(2, 0), 30.0, kTol); + EXPECT_NEAR(m.block(4, 1), 10.0, kTol); + EXPECT_NEAR(m.block(4, 5), 40.0, kTol); + EXPECT_NEAR(m.block(1, 1), 0.0, kTol); +} + +// ==================== Compression ==================== + +TEST(CompressedRowSparseMatrixGeneric, CompressAfterInsert) +{ + CRS m(4, 5); + m.setBlock(0, 1, 1.0); + m.setBlock(0, 3, 2.0); + m.setBlock(2, 0, 3.0); + m.setBlock(2, 4, 4.0); + m.compress(); + + const auto& ri = m.getRowIndex(); + const auto& rb = m.getRowBegin(); + const auto& ci = m.getColsIndex(); + const auto& cv = m.getColsValue(); + + ASSERT_EQ(ri.size(), 2u); + EXPECT_EQ(ri[0], 0); + EXPECT_EQ(ri[1], 2); + + ASSERT_EQ(rb.size(), 3u); + EXPECT_EQ(rb[0], 0); + EXPECT_EQ(rb[1], 2); + EXPECT_EQ(rb[2], 4); + + ASSERT_EQ(ci.size(), 4u); + EXPECT_EQ(ci[0], 1); + EXPECT_EQ(ci[1], 3); + EXPECT_EQ(ci[2], 0); + EXPECT_EQ(ci[3], 4); + + ASSERT_EQ(cv.size(), 4u); + EXPECT_NEAR(cv[0], 1.0, kTol); + EXPECT_NEAR(cv[1], 2.0, kTol); + EXPECT_NEAR(cv[2], 3.0, kTol); + EXPECT_NEAR(cv[3], 4.0, kTol); +} + +TEST(CompressedRowSparseMatrixGeneric, CompressRemovesZeros) +{ + CRS m(4, 4); + m.setBlock(0, 0, 5.0); + m.setBlock(0, 1, 3.0); + m.setBlock(1, 0, 1.0); + m.addBlock(0, 0, -5.0); // separate btemp entry, net (0,0) = 0 + + m.compress(); + + // The zero at (0,0) should not be in colsValue; only (0,1)=3 and (1,0)=1 remain + EXPECT_EQ(m.getColsValue().size(), 2u); + EXPECT_NEAR(m.block(0, 0), 0.0, kTol); + EXPECT_NEAR(m.block(0, 1), 3.0, kTol); + EXPECT_NEAR(m.block(1, 0), 1.0, kTol); +} + +TEST(CompressedRowSparseMatrixGeneric, CompressMergesBtemp) +{ + CRS m(4, 5); + m.setBlock(0, 0, 1.0); + m.setBlock(0, 1, 2.0); + m.compress(); + + // These go to btemp since the positions don't exist in CSR + m.setBlock(1, 0, 5.0); + m.setBlock(0, 3, 3.0); + m.compress(); + + EXPECT_NEAR(m.block(0, 0), 1.0, kTol); + EXPECT_NEAR(m.block(0, 1), 2.0, kTol); + EXPECT_NEAR(m.block(0, 3), 3.0, kTol); + EXPECT_NEAR(m.block(1, 0), 5.0, kTol); +} + +TEST(CompressedRowSparseMatrixGeneric, AutoCompressOnRead) +{ + CRS m(4, 4); + m.setBlock(1, 2, 9.0); + // No explicit compress(); block() auto-compresses + EXPECT_NEAR(m.block(1, 2), 9.0, kTol); + EXPECT_FALSE(m.getRowIndex().empty()); +} + +TEST(CompressedRowSparseMatrixGeneric, CountEmptyBlocks) +{ + CRS m(4, 4); + m.setBlock(0, 0, 1.0); + m.setBlock(0, 1, 2.0); + m.setBlock(1, 0, 3.0); + m.compress(); + + EXPECT_EQ(m.countEmptyBlocks(), 0u); + + // Zero out one entry in-place (existing CSR entry) + m.setBlock(0, 0, 0.0); + EXPECT_EQ(m.countEmptyBlocks(), 1u); +} + +// ==================== CSR Structure Access ==================== + +TEST(CompressedRowSparseMatrixGeneric, GetRowIndex) +{ + CRS m(6, 6); + m.setBlock(0, 0, 1.0); + m.setBlock(2, 1, 2.0); + m.setBlock(5, 3, 3.0); + m.compress(); + + const auto& ri = m.getRowIndex(); + ASSERT_EQ(ri.size(), 3u); + EXPECT_EQ(ri[0], 0); + EXPECT_EQ(ri[1], 2); + EXPECT_EQ(ri[2], 5); +} + +TEST(CompressedRowSparseMatrixGeneric, GetRowBegin) +{ + CRS m(6, 6); + m.setBlock(0, 0, 1.0); + m.setBlock(0, 2, 2.0); + m.setBlock(2, 1, 3.0); + m.compress(); + + const auto& rb = m.getRowBegin(); + ASSERT_EQ(rb.size(), 3u); + EXPECT_EQ(rb[0], 0); + EXPECT_EQ(rb[1], 2); + EXPECT_EQ(rb[2], 3); +} + +TEST(CompressedRowSparseMatrixGeneric, GetRowRange) +{ + CRS m(6, 6); + m.setBlock(0, 0, 1.0); + m.setBlock(0, 2, 2.0); + m.setBlock(2, 1, 3.0); + m.compress(); + + auto r0 = m.getRowRange(0); + EXPECT_EQ(r0.begin(), 0); + EXPECT_EQ(r0.end(), 2); + EXPECT_EQ(r0.size(), 2); + + auto r1 = m.getRowRange(1); + EXPECT_EQ(r1.begin(), 2); + EXPECT_EQ(r1.end(), 3); + EXPECT_EQ(r1.size(), 1); + + auto r2 = m.getRowRange(2); + EXPECT_TRUE(r2.isInvalid()); +} + +TEST(CompressedRowSparseMatrixGeneric, GetRowRangeEmpty) +{ + CRS m; + auto r = m.getRowRange(0); + EXPECT_TRUE(r.isInvalid()); +} + +TEST(CompressedRowSparseMatrixGeneric, GetColsIndexAndValue) +{ + CRS m(4, 5); + m.setBlock(0, 1, 1.0); + m.setBlock(0, 3, 2.0); + m.setBlock(2, 0, 3.0); + m.compress(); + + const auto& ci = m.getColsIndex(); + const auto& cv = m.getColsValue(); + + ASSERT_EQ(ci.size(), 3u); + ASSERT_EQ(cv.size(), 3u); + + EXPECT_EQ(ci[0], 1); + EXPECT_NEAR(cv[0], 1.0, kTol); + EXPECT_EQ(ci[1], 3); + EXPECT_NEAR(cv[1], 2.0, kTol); + EXPECT_EQ(ci[2], 0); + EXPECT_NEAR(cv[2], 3.0, kTol); +} + +TEST(CompressedRowSparseMatrixGeneric, SortedFind) +{ + CRS::VecIndex v; + v.push_back(1); + v.push_back(3); + v.push_back(5); + v.push_back(8); + v.push_back(12); + + CRS::Index result = 0; + + EXPECT_TRUE(CRS::sortedFind(v, 5, result)); + EXPECT_EQ(result, 2); + + EXPECT_TRUE(CRS::sortedFind(v, 1, result)); + EXPECT_EQ(result, 0); + + EXPECT_TRUE(CRS::sortedFind(v, 12, result)); + EXPECT_EQ(result, 4); + + EXPECT_FALSE(CRS::sortedFind(v, 7, result)); + EXPECT_FALSE(CRS::sortedFind(v, 0, result)); + + // Sub-range search + result = 0; + CRS::Range subRange(1, 4); // elements at indices 1..3: {3, 5, 8} + EXPECT_TRUE(CRS::sortedFind(v, subRange, 5, result)); + EXPECT_EQ(result, 2); + EXPECT_FALSE(CRS::sortedFind(v, subRange, 1, result)); + EXPECT_FALSE(CRS::sortedFind(v, subRange, 12, result)); +} + +// ==================== Clearing Operations ==================== + +TEST(CompressedRowSparseMatrixGeneric, Clear) +{ + CRS m(4, 4); + m.setBlock(0, 0, 1.0); + m.setBlock(1, 2, 5.0); + m.compress(); + EXPECT_FALSE(m.getColsValue().empty()); + + m.clear(); + // ClearByZeros + CompressZeros: zeros all values then compresses them away + EXPECT_TRUE(m.getColsValue().empty()); + EXPECT_TRUE(m.getRowIndex().empty()); +} + +TEST(CompressedRowSparseMatrixGeneric, ClearRowBlock) +{ + CRS m(4, 4); + m.setBlock(0, 0, 1.0); + m.setBlock(0, 1, 2.0); + m.setBlock(1, 0, 3.0); + m.setBlock(2, 2, 4.0); + m.compress(); + + m.clearRowBlock(0); + + EXPECT_NEAR(m.block(0, 0), 0.0, kTol); + EXPECT_NEAR(m.block(0, 1), 0.0, kTol); + EXPECT_NEAR(m.block(1, 0), 3.0, kTol); + EXPECT_NEAR(m.block(2, 2), 4.0, kTol); +} + +TEST(CompressedRowSparseMatrixGeneric, ClearColBlock) +{ + CRS m(4, 4); + m.setBlock(0, 1, 1.0); + m.setBlock(1, 1, 2.0); + m.setBlock(2, 0, 3.0); + m.setBlock(2, 1, 4.0); + m.compress(); + + m.clearColBlock(1); + + EXPECT_NEAR(m.block(0, 1), 0.0, kTol); + EXPECT_NEAR(m.block(1, 1), 0.0, kTol); + EXPECT_NEAR(m.block(2, 1), 0.0, kTol); + EXPECT_NEAR(m.block(2, 0), 3.0, kTol); +} + +// ==================== Row Operations ==================== + +TEST(CompressedRowSparseMatrixGeneric, FullRows) +{ + CRS m(6, 6); + m.setBlock(0, 0, 1.0); + m.setBlock(5, 1, 2.0); + m.compress(); + + EXPECT_EQ(m.getRowIndex().size(), 2u); + + m.fullRows(); + + const auto& ri = m.getRowIndex(); + EXPECT_EQ(ri.size(), 6u); + for (CRS::Index i = 0; i < 6; ++i) + EXPECT_EQ(ri[i], i); +} + +TEST(CompressedRowSparseMatrixGeneric, FullRowsEmpty) +{ + CRS m(4, 4); + m.fullRows(); + EXPECT_EQ(m.getRowIndex().size(), 4u); + EXPECT_TRUE(m.getColsValue().empty()); +} + +TEST(CompressedRowSparseMatrixGeneric, ShiftIndices) +{ + CRS m(4, 4); + m.setBlock(0, 1, 1.0); + m.setBlock(2, 3, 2.0); + m.compress(); + + const auto riBefore = m.getRowIndex(); + const auto rbBefore = m.getRowBegin(); + const auto ciBefore = m.getColsIndex(); + + m.shiftIndices(1); + + const auto& ri = m.getRowIndex(); + const auto& rb = m.getRowBegin(); + const auto& ci = m.getColsIndex(); + + for (std::size_t i = 0; i < ri.size(); ++i) + EXPECT_EQ(ri[i], riBefore[i] + 1); + for (std::size_t i = 0; i < rb.size(); ++i) + EXPECT_EQ(rb[i], rbBefore[i] + 1); + for (std::size_t i = 0; i < ci.size(); ++i) + EXPECT_EQ(ci[i], ciBefore[i] + 1); +} + +// ==================== Matrix Operations ==================== + +TEST(CompressedRowSparseMatrixGeneric, Mul) +{ + // A(3x2) * B(2x2) = C(3x2) + CRS A(3, 2), B(2, 2), C; + + A.setBlock(0, 0, 1.0); A.setBlock(0, 1, 2.0); + A.setBlock(1, 0, 3.0); A.setBlock(1, 1, 4.0); + A.setBlock(2, 0, 5.0); A.setBlock(2, 1, 6.0); + A.compress(); + + B.setBlock(0, 0, 7.0); B.setBlock(0, 1, 8.0); + B.setBlock(1, 0, 9.0); B.setBlock(1, 1, 10.0); + B.compress(); + + A.mul(C, B); + + EXPECT_EQ(C.rowBSize(), 3); + EXPECT_EQ(C.colBSize(), 2); + + EXPECT_NEAR(C.block(0, 0), 25.0, kTol); // 1*7 + 2*9 + EXPECT_NEAR(C.block(0, 1), 28.0, kTol); // 1*8 + 2*10 + EXPECT_NEAR(C.block(1, 0), 57.0, kTol); // 3*7 + 4*9 + EXPECT_NEAR(C.block(1, 1), 64.0, kTol); // 3*8 + 4*10 + EXPECT_NEAR(C.block(2, 0), 89.0, kTol); // 5*7 + 6*9 + EXPECT_NEAR(C.block(2, 1), 100.0, kTol); // 5*8 + 6*10 +} + +TEST(CompressedRowSparseMatrixGeneric, MulTranspose) +{ + // C = A^T(2x3) * B(3x2) = (2x2) + CRS A(3, 2), B(3, 2), C; + + A.setBlock(0, 0, 1.0); A.setBlock(0, 1, 2.0); + A.setBlock(1, 0, 3.0); A.setBlock(1, 1, 4.0); + A.setBlock(2, 0, 5.0); A.setBlock(2, 1, 6.0); + A.compress(); + + B.setBlock(0, 0, 1.0); B.setBlock(0, 1, 2.0); + B.setBlock(1, 0, 3.0); B.setBlock(1, 1, 4.0); + B.setBlock(2, 0, 5.0); B.setBlock(2, 1, 6.0); + B.compress(); + + A.mulTranspose(C, B); + + EXPECT_EQ(C.rowBSize(), 2); + EXPECT_EQ(C.colBSize(), 2); + + EXPECT_NEAR(C.block(0, 0), 35.0, kTol); // 1*1 + 3*3 + 5*5 + EXPECT_NEAR(C.block(0, 1), 44.0, kTol); // 1*2 + 3*4 + 5*6 + EXPECT_NEAR(C.block(1, 0), 44.0, kTol); // 2*1 + 4*3 + 6*5 + EXPECT_NEAR(C.block(1, 1), 56.0, kTol); // 2*2 + 4*4 + 6*6 +} + +TEST(CompressedRowSparseMatrixGeneric, TransposeFullRows) +{ + CRS A(3, 2), AT; + + A.setBlock(0, 0, 1.0); A.setBlock(0, 1, 2.0); + A.setBlock(1, 0, 3.0); A.setBlock(1, 1, 4.0); + A.setBlock(2, 0, 5.0); + A.compress(); + A.fullRows(); + + A.transposeFullRows(AT); + + EXPECT_EQ(AT.rowBSize(), 2); + EXPECT_EQ(AT.colBSize(), 3); + + EXPECT_NEAR(AT.block(0, 0), 1.0, kTol); + EXPECT_NEAR(AT.block(0, 1), 3.0, kTol); + EXPECT_NEAR(AT.block(0, 2), 5.0, kTol); + EXPECT_NEAR(AT.block(1, 0), 2.0, kTol); + EXPECT_NEAR(AT.block(1, 1), 4.0, kTol); + EXPECT_NEAR(AT.block(1, 2), 0.0, kTol); +} + +// ==================== Swap ==================== + +TEST(CompressedRowSparseMatrixGeneric, Swap) +{ + CRS a(3, 3), b(4, 4); + a.setBlock(0, 0, 1.0); + a.compress(); + b.setBlock(2, 3, 7.0); + b.compress(); + + a.swap(b); + + EXPECT_EQ(a.rowBSize(), 4); + EXPECT_EQ(a.colBSize(), 4); + EXPECT_NEAR(a.block(2, 3), 7.0, kTol); + + EXPECT_EQ(b.rowBSize(), 3); + EXPECT_EQ(b.colBSize(), 3); + EXPECT_NEAR(b.block(0, 0), 1.0, kTol); +} + +// ==================== Name ==================== + +TEST(CompressedRowSparseMatrixGeneric, Name) +{ + EXPECT_STREQ(CRS::Name(), "CompressedRowSparseMatrixd"); +} + +// ==================== Block Type Spot-Check (Mat3x3d) ==================== + +TEST(CompressedRowSparseMatrixGeneric, Mat3x3dInsertAndRetrieve) +{ + CRSMat3 m(2, 2); + Mat3 block; + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + block(i, j) = static_cast(i * 3 + j + 1); + + m.setBlock(0, 1, block); + const auto& retrieved = m.block(0, 1); + + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + EXPECT_NEAR(retrieved(i, j), block(i, j), kTol); +} + +TEST(CompressedRowSparseMatrixGeneric, Mat3x3dMul) +{ + // A(1x2 blocks) * B(2x1 blocks) = C(1x1 block) + CRSMat3 A(1, 2), B(2, 1), C; + + Mat3 I; + I.identity(); + + Mat3 D; + D.clear(); + for (int i = 0; i < 3; ++i) D(i, i) = 2.0; + + A.setBlock(0, 0, I); + A.setBlock(0, 1, D); + A.compress(); + + B.setBlock(0, 0, I); + B.setBlock(1, 0, I); + B.compress(); + + A.mul(C, B); + + // C(0,0) = I*I + D*I = I + D = diag(3,3,3) + const auto& result = C.block(0, 0); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + { + const double expected = (i == j) ? 3.0 : 0.0; + EXPECT_NEAR(result(i, j), expected, kTol); + } +} From ff72d63acb994a038a07625a6737eadb624aede6 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Tue, 3 Mar 2026 14:07:04 +0900 Subject: [PATCH 07/11] add unit tests for CompressedRowSparseMatrixMechanical --- .../LinearAlgebra/test/CMakeLists.txt | 1 + ...mpressedRowSparseMatrixMechanical_test.cpp | 503 ++++++++++++++++++ 2 files changed, 504 insertions(+) create mode 100644 Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixMechanical_test.cpp diff --git a/Sofa/framework/LinearAlgebra/test/CMakeLists.txt b/Sofa/framework/LinearAlgebra/test/CMakeLists.txt index bdcc27aac9f..4c1b3d8c07e 100644 --- a/Sofa/framework/LinearAlgebra/test/CMakeLists.txt +++ b/Sofa/framework/LinearAlgebra/test/CMakeLists.txt @@ -9,6 +9,7 @@ set(SOURCE_FILES CompressedRowSparseMatrix_test.cpp CompressedRowSparseMatrixGeneric_test.cpp CompressedRowSparseMatrixConstraint_test.cpp + CompressedRowSparseMatrixMechanical_test.cpp Matrix_test.cpp RotationMatrix_test.cpp SparseMatrixProduct_test.cpp diff --git a/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixMechanical_test.cpp b/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixMechanical_test.cpp new file mode 100644 index 00000000000..c0aff0f3fe6 --- /dev/null +++ b/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixMechanical_test.cpp @@ -0,0 +1,503 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#include +#include +#include + +using CRSMech = sofa::linearalgebra::CompressedRowSparseMatrixMechanical; +using CRSMechMat3 = sofa::linearalgebra::CompressedRowSparseMatrixMechanical>; +using FVec = sofa::linearalgebra::FullVector; + +constexpr double kTol = 1e-10; + +// ==================== Construction & Sizing ==================== + +TEST(CompressedRowSparseMatrixMechanical, DefaultConstruction) +{ + CRSMech m; + EXPECT_EQ(m.rowSize(), 0); + EXPECT_EQ(m.colSize(), 0); + EXPECT_EQ(m.rowBSize(), 0); + EXPECT_EQ(m.colBSize(), 0); + EXPECT_EQ(m.nRow, 0); + EXPECT_EQ(m.nCol, 0); +} + +TEST(CompressedRowSparseMatrixMechanical, ConstructWithScalarDimensions) +{ + CRSMech m(6, 8); + EXPECT_EQ(m.rowSize(), 6); + EXPECT_EQ(m.colSize(), 8); + EXPECT_EQ(m.rowBSize(), 6); + EXPECT_EQ(m.colBSize(), 8); +} + +TEST(CompressedRowSparseMatrixMechanical, Resize) +{ + CRSMech m(10, 10); + m.set(0, 0, 5.0); + m.compress(); + + m.resize(5, 7); + EXPECT_EQ(m.rowSize(), 5); + EXPECT_EQ(m.colSize(), 7); + EXPECT_EQ(m.rowBSize(), 5); + EXPECT_EQ(m.colBSize(), 7); +} + +TEST(CompressedRowSparseMatrixMechanical, ResizeBlock) +{ + CRSMech m; + m.resizeBlock(3, 4); + EXPECT_EQ(m.rowBSize(), 3); + EXPECT_EQ(m.colBSize(), 4); + EXPECT_EQ(m.nRow, 3); + EXPECT_EQ(m.nCol, 4); +} + +TEST(CompressedRowSparseMatrixMechanical, Extend) +{ + CRSMech m(4, 4); + m.set(0, 0, 1.0); + m.set(1, 1, 2.0); + m.compress(); + + m.extend(10, 10); + EXPECT_EQ(m.rowSize(), 10); + EXPECT_EQ(m.colSize(), 10); + EXPECT_EQ(m.rowBSize(), 10); + EXPECT_EQ(m.colBSize(), 10); + + // data is preserved + EXPECT_NEAR(m.element(0, 0), 1.0, kTol); + EXPECT_NEAR(m.element(1, 1), 2.0, kTol); +} + +// ==================== Scalar Element Access ==================== + +TEST(CompressedRowSparseMatrixMechanical, SetAndElement) +{ + CRSMech m(5, 5); + m.set(2, 3, 7.5); + EXPECT_NEAR(m.element(2, 3), 7.5, kTol); +} + +TEST(CompressedRowSparseMatrixMechanical, AddScalar) +{ + CRSMech m(5, 5); + m.set(0, 0, 3.0); + m.compress(); + m.add(0, 0, 2.0); + EXPECT_NEAR(m.element(0, 0), 5.0, kTol); +} + +TEST(CompressedRowSparseMatrixMechanical, ClearScalar) +{ + CRSMech m(5, 5); + m.set(1, 2, 5.0); + m.compress(); + m.clear(1, 2); + EXPECT_NEAR(m.element(1, 2), 0.0, kTol); +} + +TEST(CompressedRowSparseMatrixMechanical, ElementNonExistent) +{ + CRSMech m(5, 5); + m.set(0, 0, 1.0); + m.compress(); + EXPECT_NEAR(m.element(3, 4), 0.0, kTol); +} + +// ==================== Row/Column Clear ==================== + +TEST(CompressedRowSparseMatrixMechanical, ClearRow) +{ + CRSMech m(4, 4); + m.set(1, 0, 1.0); + m.set(1, 1, 2.0); + m.set(1, 2, 3.0); + m.set(1, 3, 4.0); + m.set(0, 0, 10.0); + m.set(2, 2, 20.0); + m.compress(); + + m.clearRow(1); + + EXPECT_NEAR(m.element(1, 0), 0.0, kTol); + EXPECT_NEAR(m.element(1, 1), 0.0, kTol); + EXPECT_NEAR(m.element(1, 2), 0.0, kTol); + EXPECT_NEAR(m.element(1, 3), 0.0, kTol); + // other rows unaffected + EXPECT_NEAR(m.element(0, 0), 10.0, kTol); + EXPECT_NEAR(m.element(2, 2), 20.0, kTol); +} + +TEST(CompressedRowSparseMatrixMechanical, ClearCol) +{ + CRSMech m(4, 4); + m.set(0, 2, 1.0); + m.set(1, 2, 2.0); + m.set(2, 2, 3.0); + m.set(3, 2, 4.0); + m.set(0, 0, 10.0); + m.set(3, 3, 20.0); + m.compress(); + + m.clearCol(2); + + EXPECT_NEAR(m.element(0, 2), 0.0, kTol); + EXPECT_NEAR(m.element(1, 2), 0.0, kTol); + EXPECT_NEAR(m.element(2, 2), 0.0, kTol); + EXPECT_NEAR(m.element(3, 2), 0.0, kTol); + // other columns unaffected + EXPECT_NEAR(m.element(0, 0), 10.0, kTol); + EXPECT_NEAR(m.element(3, 3), 20.0, kTol); +} + +TEST(CompressedRowSparseMatrixMechanical, ClearRowCol) +{ + CRSMech m(4, 4); + // set up a symmetric-looking matrix with fullDiagonal + m.set(0, 0, 1.0); + m.set(0, 1, 2.0); + m.set(1, 0, 2.0); + m.set(1, 1, 3.0); + m.set(1, 2, 4.0); + m.set(2, 1, 4.0); + m.set(2, 2, 5.0); + m.set(3, 3, 6.0); + m.compress(); + m.fullDiagonal(); + + m.clearRowCol(1); + + // row 1 and column 1 should be zero + for (int j = 0; j < 4; ++j) + EXPECT_NEAR(m.element(1, j), 0.0, kTol) << "element(1," << j << ")"; + for (int i = 0; i < 4; ++i) + EXPECT_NEAR(m.element(i, 1), 0.0, kTol) << "element(" << i << ",1)"; + + // other entries unaffected + EXPECT_NEAR(m.element(0, 0), 1.0, kTol); + EXPECT_NEAR(m.element(2, 2), 5.0, kTol); + EXPECT_NEAR(m.element(3, 3), 6.0, kTol); +} + +TEST(CompressedRowSparseMatrixMechanical, ClearOverride) +{ + CRSMech m(4, 4); + m.set(0, 0, 1.0); + m.set(1, 2, 5.0); + m.set(2, 3, 7.0); + m.compress(); + + const auto numValuesBefore = m.getColsValue().size(); + EXPECT_GT(numValuesBefore, 0u); + + m.clear(); + + // ClearByZeros=true zeros values, CompressZeros=false preserves structure + EXPECT_EQ(m.getColsValue().size(), numValuesBefore); + for (const auto& v : m.getColsValue()) + { + EXPECT_NEAR(v, 0.0, kTol); + } +} + +// ==================== Diagonal ==================== + +TEST(CompressedRowSparseMatrixMechanical, FullDiagonal) +{ + CRSMech m(4, 4); + m.set(0, 1, 2.0); + m.set(2, 3, 5.0); + m.compress(); + + m.fullDiagonal(); + + const auto& ri = m.getRowIndex(); + ASSERT_EQ(ri.size(), 4u); + for (int i = 0; i < 4; ++i) + EXPECT_EQ(ri[i], i); + + // diagonal entries exist (value is 0 where not explicitly set) + EXPECT_NEAR(m.element(0, 0), 0.0, kTol); + EXPECT_NEAR(m.element(1, 1), 0.0, kTol); + EXPECT_NEAR(m.element(2, 2), 0.0, kTol); + EXPECT_NEAR(m.element(3, 3), 0.0, kTol); + // off-diagonal preserved + EXPECT_NEAR(m.element(0, 1), 2.0, kTol); + EXPECT_NEAR(m.element(2, 3), 5.0, kTol); +} + +// ==================== Matrix-Vector Products ==================== + +TEST(CompressedRowSparseMatrixMechanical, MulVector) +{ + // M = [1 2 0] + // [0 3 4] + // [5 0 6] + CRSMech m(3, 3); + m.set(0, 0, 1.0); m.set(0, 1, 2.0); + m.set(1, 1, 3.0); m.set(1, 2, 4.0); + m.set(2, 0, 5.0); m.set(2, 2, 6.0); + m.compress(); + + FVec v(3); + v[0] = 1.0; v[1] = 2.0; v[2] = 3.0; + + FVec res(3); + m.mul(res, v); + + // M*v = [1+4, 6+12, 5+18] = [5, 18, 23] + EXPECT_NEAR(res[0], 5.0, kTol); + EXPECT_NEAR(res[1], 18.0, kTol); + EXPECT_NEAR(res[2], 23.0, kTol); +} + +TEST(CompressedRowSparseMatrixMechanical, AddMulVector) +{ + // M = [2 1] + // [3 0] + CRSMech m(2, 2); + m.set(0, 0, 2.0); m.set(0, 1, 1.0); + m.set(1, 0, 3.0); + m.compress(); + + FVec v(2); + v[0] = 1.0; v[1] = 2.0; + + FVec res(2); + res[0] = 0.0; res[1] = 0.0; + m.addMul(res, v); + + // res = M * v = [2+2, 3+0] = [4, 3] + EXPECT_NEAR(res[0], 4.0, kTol); + EXPECT_NEAR(res[1], 3.0, kTol); +} + +// ==================== Arithmetic Operators ==================== + +TEST(CompressedRowSparseMatrixMechanical, OperatorPlus) +{ + CRSMech A(3, 3), B(3, 3); + A.set(0, 0, 1.0); A.set(0, 1, 2.0); A.set(1, 1, 3.0); + B.set(0, 0, 4.0); B.set(1, 0, 1.0); B.set(1, 1, 2.0); + A.compress(); + B.compress(); + + CRSMech C = A + B; + EXPECT_NEAR(C.element(0, 0), 5.0, kTol); + EXPECT_NEAR(C.element(0, 1), 2.0, kTol); + EXPECT_NEAR(C.element(1, 0), 1.0, kTol); + EXPECT_NEAR(C.element(1, 1), 5.0, kTol); +} + +TEST(CompressedRowSparseMatrixMechanical, OperatorPlusEquals) +{ + CRSMech A(3, 3), B(3, 3); + A.set(0, 0, 1.0); A.set(1, 1, 2.0); + B.set(0, 0, 3.0); B.set(1, 1, 4.0); B.set(2, 2, 5.0); + A.compress(); + B.compress(); + + A += B; + EXPECT_NEAR(A.element(0, 0), 4.0, kTol); + EXPECT_NEAR(A.element(1, 1), 6.0, kTol); + EXPECT_NEAR(A.element(2, 2), 5.0, kTol); +} + +TEST(CompressedRowSparseMatrixMechanical, OperatorMinusEquals) +{ + CRSMech A(3, 3), B(3, 3); + A.set(0, 0, 5.0); A.set(1, 1, 8.0); + B.set(0, 0, 2.0); B.set(1, 1, 3.0); + A.compress(); + B.compress(); + + A -= B; + EXPECT_NEAR(A.element(0, 0), 3.0, kTol); + EXPECT_NEAR(A.element(1, 1), 5.0, kTol); +} + +// ==================== Info Methods ==================== + +TEST(CompressedRowSparseMatrixMechanical, GetCategory) +{ + CRSMech m; + EXPECT_EQ(m.getCategory(), sofa::linearalgebra::BaseMatrix::MATRIX_SPARSE); +} + +TEST(CompressedRowSparseMatrixMechanical, GetBlockRowsCols) +{ + CRSMech mScalar; + EXPECT_EQ(mScalar.getBlockRows(), 1); + EXPECT_EQ(mScalar.getBlockCols(), 1); + + CRSMechMat3 mBlock; + EXPECT_EQ(mBlock.getBlockRows(), 3); + EXPECT_EQ(mBlock.getBlockCols(), 3); +} + +TEST(CompressedRowSparseMatrixMechanical, GetBandWidth) +{ + CRSMech mScalar; + EXPECT_EQ(mScalar.getBandWidth(), 0); // NC - 1 = 0 + + CRSMechMat3 mBlock; + EXPECT_EQ(mBlock.getBandWidth(), 2); // NC - 1 = 2 +} + +TEST(CompressedRowSparseMatrixMechanical, RowColSize) +{ + CRSMech m(6, 8); + EXPECT_EQ(m.rowSize(), 6); + EXPECT_EQ(m.colSize(), 8); + EXPECT_EQ(m.bRowSize(), 6); + EXPECT_EQ(m.bColSize(), 8); +} + +// ==================== Swap ==================== + +TEST(CompressedRowSparseMatrixMechanical, Swap) +{ + CRSMech a(4, 4), b(6, 6); + a.set(0, 0, 1.0); + a.set(1, 1, 2.0); + a.compress(); + + b.set(3, 3, 7.0); + b.set(5, 5, 9.0); + b.compress(); + + a.swap(b); + + EXPECT_EQ(a.rowSize(), 6); + EXPECT_EQ(a.colSize(), 6); + EXPECT_NEAR(a.element(3, 3), 7.0, kTol); + EXPECT_NEAR(a.element(5, 5), 9.0, kTol); + + EXPECT_EQ(b.rowSize(), 4); + EXPECT_EQ(b.colSize(), 4); + EXPECT_NEAR(b.element(0, 0), 1.0, kTol); + EXPECT_NEAR(b.element(1, 1), 2.0, kTol); +} + +// ==================== Filtering ==================== + +TEST(CompressedRowSparseMatrixMechanical, CopyUpper) +{ + CRSMech src(3, 3); + src.set(0, 0, 1.0); src.set(0, 1, 2.0); src.set(0, 2, 3.0); + src.set(1, 0, 4.0); src.set(1, 1, 5.0); src.set(1, 2, 6.0); + src.set(2, 0, 7.0); src.set(2, 1, 8.0); src.set(2, 2, 9.0); + src.compress(); + + CRSMech dst; + dst.copyUpper(src); + + // upper triangle: (0,0), (0,1), (0,2), (1,1), (1,2), (2,2) + EXPECT_NEAR(dst.element(0, 0), 1.0, kTol); + EXPECT_NEAR(dst.element(0, 1), 2.0, kTol); + EXPECT_NEAR(dst.element(0, 2), 3.0, kTol); + EXPECT_NEAR(dst.element(1, 1), 5.0, kTol); + EXPECT_NEAR(dst.element(1, 2), 6.0, kTol); + EXPECT_NEAR(dst.element(2, 2), 9.0, kTol); + + // lower triangle should not be present + EXPECT_NEAR(dst.element(1, 0), 0.0, kTol); + EXPECT_NEAR(dst.element(2, 0), 0.0, kTol); + EXPECT_NEAR(dst.element(2, 1), 0.0, kTol); +} + +TEST(CompressedRowSparseMatrixMechanical, CopyLower) +{ + CRSMech src(3, 3); + src.set(0, 0, 1.0); src.set(0, 1, 2.0); src.set(0, 2, 3.0); + src.set(1, 0, 4.0); src.set(1, 1, 5.0); src.set(1, 2, 6.0); + src.set(2, 0, 7.0); src.set(2, 1, 8.0); src.set(2, 2, 9.0); + src.compress(); + + CRSMech dst; + dst.copyLower(src); + + // lower triangle: (0,0), (1,0), (1,1), (2,0), (2,1), (2,2) + EXPECT_NEAR(dst.element(0, 0), 1.0, kTol); + EXPECT_NEAR(dst.element(1, 0), 4.0, kTol); + EXPECT_NEAR(dst.element(1, 1), 5.0, kTol); + EXPECT_NEAR(dst.element(2, 0), 7.0, kTol); + EXPECT_NEAR(dst.element(2, 1), 8.0, kTol); + EXPECT_NEAR(dst.element(2, 2), 9.0, kTol); + + // upper triangle should not be present + EXPECT_NEAR(dst.element(0, 1), 0.0, kTol); + EXPECT_NEAR(dst.element(0, 2), 0.0, kTol); + EXPECT_NEAR(dst.element(1, 2), 0.0, kTol); +} + +// ==================== Mat3x3d Block Spot-Check ==================== + +TEST(CompressedRowSparseMatrixMechanical, Mat3x3dScalarAccess) +{ + CRSMechMat3 m(9, 9); + + // set scalar entries within the block at block position (0, 0) + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + m.set(i, j, static_cast(i * 3 + j + 1)); + + m.compress(); + + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + EXPECT_NEAR(m.element(i, j), static_cast(i * 3 + j + 1), kTol) + << "element(" << i << "," << j << ")"; +} + +TEST(CompressedRowSparseMatrixMechanical, Mat3x3dResize) +{ + CRSMechMat3 m(9, 9); + EXPECT_EQ(m.rowSize(), 9); + EXPECT_EQ(m.colSize(), 9); + EXPECT_EQ(m.rowBSize(), 3); + EXPECT_EQ(m.colBSize(), 3); +} + +TEST(CompressedRowSparseMatrixMechanical, Mat3x3dMulVector) +{ + // Single 3x3 block = diagonal [1, 2, 3] + CRSMechMat3 m(3, 3); + m.set(0, 0, 1.0); + m.set(1, 1, 2.0); + m.set(2, 2, 3.0); + m.compress(); + + FVec v(3); + v[0] = 1.0; v[1] = 1.0; v[2] = 1.0; + + FVec res(3); + m.mul(res, v); + + EXPECT_NEAR(res[0], 1.0, kTol); + EXPECT_NEAR(res[1], 2.0, kTol); + EXPECT_NEAR(res[2], 3.0, kTol); +} From f905d32fa274ac19d957639cbb41feba07d44280 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Tue, 3 Mar 2026 14:18:50 +0900 Subject: [PATCH 08/11] CRSM, alwayymmetric : fix typo in code and subsequent compilation error (missing template) --- .../sofa/linearalgebra/CompressedRowSparseMatrixMechanical.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixMechanical.h b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixMechanical.h index 003f2920571..f5e003aeff4 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixMechanical.h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixMechanical.h @@ -1141,9 +1141,9 @@ class CompressedRowSparseMatrixMechanical final // final is used to allow the co if constexpr (Policy::AutoCompress) const_cast(this)->compress(); /// \warning this violates the const-ness of the method ! vresize( res, this->colBSize(), colSize() ); - if constexpr (Policy::IsAlwaysSymetric) /// In symetric case this^T = this + if constexpr (Policy::IsAlwaysSymmetric) /// In symmetric case this^T = this { - taddMul(res, vec); + taddMul(res, vec); return; } From db4c1d702542b1dcf6d1c86c5d7689fd4c6b4e60 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Tue, 3 Mar 2026 14:19:25 +0900 Subject: [PATCH 09/11] add AddMultTranspose with symmetric matrix (not possible before) --- ...mpressedRowSparseMatrixMechanical_test.cpp | 25 +++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixMechanical_test.cpp b/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixMechanical_test.cpp index c0aff0f3fe6..bc578192f23 100644 --- a/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixMechanical_test.cpp +++ b/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixMechanical_test.cpp @@ -295,6 +295,31 @@ TEST(CompressedRowSparseMatrixMechanical, AddMulVector) EXPECT_NEAR(res[1], 3.0, kTol); } +TEST(CompressedRowSparseMatrixMechanical, AddMultTranspose) +{ + // For CRSMechanicalPolicy, IsAlwaysSymmetric=true so addMultTranspose + // delegates to addMul (this^T == this for a symmetric matrix). + // Use a symmetric matrix to verify. + // + // M = [1 2] + // [2 3] + // v = [1, 2], expected M^T * v = M * v = [1+4, 2+6] = [5, 8] + CRSMech m(2, 2); + m.set(0, 0, 1.0); m.set(0, 1, 2.0); + m.set(1, 0, 2.0); m.set(1, 1, 3.0); + m.compress(); + + FVec v(2); + v[0] = 1.0; v[1] = 2.0; + + FVec res(2); + res[0] = 0.0; res[1] = 0.0; + m.addMultTranspose(res, v); + + EXPECT_NEAR(res[0], 5.0, kTol); + EXPECT_NEAR(res[1], 8.0, kTol); +} + // ==================== Arithmetic Operators ==================== TEST(CompressedRowSparseMatrixMechanical, OperatorPlus) From bec9987071fdb9338420d82107ee83e0d00e6edb Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Tue, 3 Mar 2026 14:34:24 +0900 Subject: [PATCH 10/11] CompressedRowSparseMatrixConstraint: fix bug when calling cbegin/cend on empty matrix --- .../linearalgebra/CompressedRowSparseMatrixConstraint.h | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixConstraint.h b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixConstraint.h index e3a6d3db154..51f583d017f 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixConstraint.h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixConstraint.h @@ -454,14 +454,16 @@ class CompressedRowSparseMatrixConstraint : public sofa::linearalgebra::Compress RowConstIterator cbegin() const { if constexpr(Policy::AutoCompress) const_cast(this)->compress(); /// \warning this violates the const-ness of the method ! - return RowConstIterator(this, 0); + return RowConstIterator(this, + this->rowIndex.empty() ? s_invalidIndex : 0); } /// Get the iterator corresponding to the end of the rows of blocks RowConstIterator cend() const { if constexpr(Policy::AutoCompress) const_cast(this)->compress(); /// \warning this violates the const-ness of the method ! - return RowConstIterator(this, Index(this->rowIndex.size())); + return RowConstIterator(this, + this->rowIndex.empty() ? s_invalidIndex : Index(this->rowIndex.size())); } class RowWriteAccessor From 893c1b315eb39fd97875e207f337587452c7c7ac Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Tue, 3 Mar 2026 14:34:41 +0900 Subject: [PATCH 11/11] add more unit tests for CompressedRowSparseMatrixConstraint --- ...mpressedRowSparseMatrixConstraint_test.cpp | 231 ++++++++++++++++++ 1 file changed, 231 insertions(+) diff --git a/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixConstraint_test.cpp b/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixConstraint_test.cpp index 9b164ce6bd1..20f425d537b 100644 --- a/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixConstraint_test.cpp +++ b/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixConstraint_test.cpp @@ -28,6 +28,7 @@ #include #include +#include #include #include #include @@ -2191,4 +2192,234 @@ Constraint ID : 4 dof ID : 5 value : 0 0 0.474108 dof ID : 7 value : 0 0.983 EXPECT_EQ(oss.str(), expectedOutput); } +// ==================== cbegin / cend ==================== + +TEST(CompressedRowSparseMatrixConstraint, cbeginCendEmptyMatrix) +{ + const sofa::linearalgebra::CompressedRowSparseMatrixConstraint m; + + // cbegin/cend must be consistent with begin/end on an empty matrix + EXPECT_EQ(m.cbegin(), m.cend()); + EXPECT_EQ(m.cbegin(), m.begin()); + EXPECT_EQ(m.cend(), m.end()); +} + +TEST(CompressedRowSparseMatrixConstraint, cbeginCendNonEmpty) +{ + sofa::linearalgebra::CompressedRowSparseMatrixConstraint m; + m.writeLine(0).addCol(0, sofa::type::Vec3(1, 2, 3)); + m.writeLine(5).addCol(1, sofa::type::Vec3(4, 5, 6)); + m.compress(); + + EXPECT_EQ(m.cbegin(), m.begin()); + EXPECT_EQ(m.cend(), m.end()); + + int count = 0; + for (auto it = m.cbegin(); it != m.cend(); ++it) + ++count; + EXPECT_EQ(count, 2); +} + +// ==================== Constructor with dimensions ==================== + +TEST(CompressedRowSparseMatrixConstraint, ConstructWithDimensions) +{ + sofa::linearalgebra::CompressedRowSparseMatrixConstraint m(5, 10); + EXPECT_EQ(m.rowBSize(), 5); + EXPECT_EQ(m.colBSize(), 10); + EXPECT_TRUE(m.empty()); +} + +// ==================== Name ==================== + +TEST(CompressedRowSparseMatrixConstraint, Name) +{ + using Matrix = sofa::linearalgebra::CompressedRowSparseMatrixConstraint; + const std::string name = Matrix::Name(); + EXPECT_FALSE(name.empty()); + EXPECT_NE(name.find("CompressedRowSparseMatrixConstraint"), std::string::npos); +} + +// ==================== RowType::find ==================== + +TEST(CompressedRowSparseMatrixConstraint, RowTypeFind) +{ + using Vec3 = sofa::type::Vec3; + using Matrix = sofa::linearalgebra::CompressedRowSparseMatrixConstraint; + + Matrix m; + auto row = m.writeLine(0); + row.addCol(5, Vec3(1, 2, 3)); + row.addCol(10, Vec3(4, 5, 6)); + row.addCol(15, Vec3(7, 8, 9)); + m.compress(); + + auto rowType = m.readLine(0).row(); + + // find existing column + auto found = rowType.find(10); + ASSERT_NE(found, rowType.end()); + EXPECT_EQ(found.index(), 10u); + EXPECT_EQ(found.val(), Vec3(4, 5, 6)); + + // find non-existing column returns end + auto notFound = rowType.find(7); + EXPECT_EQ(notFound, rowType.end()); +} + +// ==================== Row iterator arithmetic ==================== + +TEST(CompressedRowSparseMatrixConstraint, RowIteratorArithmetic) +{ + using Vec3 = sofa::type::Vec3; + using Matrix = sofa::linearalgebra::CompressedRowSparseMatrixConstraint; + + Matrix m; + m.writeLine(0).addCol(0, Vec3(1, 0, 0)); + m.writeLine(1).addCol(0, Vec3(0, 1, 0)); + m.writeLine(2).addCol(0, Vec3(0, 0, 1)); + m.compress(); + + auto it = m.begin(); + + // operator+ + auto it2 = it + 2; + EXPECT_EQ(it2.index(), 2u); + + // operator-(difference) + auto it0 = it2 - 2; + EXPECT_EQ(it0, it); + + // iterator difference + EXPECT_EQ(it2 - it, 2); + + // operator+=, operator-= + auto it3 = it; + it3 += 1; + EXPECT_EQ(it3.index(), 1u); + it3 -= 1; + EXPECT_EQ(it3, it); + + // comparison operators + EXPECT_TRUE(it < it2); + EXPECT_TRUE(it <= it2); + EXPECT_TRUE(it2 > it); + EXPECT_TRUE(it2 >= it); + EXPECT_TRUE(it <= it); + EXPECT_TRUE(it >= it); + EXPECT_FALSE(it > it2); + EXPECT_FALSE(it2 < it); +} + +// ==================== Column iterator arithmetic ==================== + +TEST(CompressedRowSparseMatrixConstraint, ColIteratorArithmetic) +{ + using Vec3 = sofa::type::Vec3; + using Matrix = sofa::linearalgebra::CompressedRowSparseMatrixConstraint; + + Matrix m; + auto row = m.writeLine(0); + row.addCol(5, Vec3(1, 0, 0)); + row.addCol(10, Vec3(0, 1, 0)); + row.addCol(15, Vec3(0, 0, 1)); + m.compress(); + + auto rowIt = m.readLine(0); + auto begin = rowIt.begin(); + auto end = rowIt.end(); + + // operator+= + auto it = begin; + it += 2; + EXPECT_EQ(it.index(), 15u); + + // operator-= + it -= 2; + EXPECT_EQ(it, begin); + + // prefix decrement + it += 2; + --it; + EXPECT_EQ(it.index(), 10u); + + // postfix decrement + auto prev = it--; + EXPECT_EQ(prev.index(), 10u); + EXPECT_EQ(it.index(), 5u); + + // comparison operators + EXPECT_TRUE(begin < end); + EXPECT_TRUE(end > begin); + EXPECT_TRUE(begin <= begin); + EXPECT_TRUE(begin >= begin); + EXPECT_TRUE(begin <= end); + EXPECT_TRUE(end >= begin); + EXPECT_FALSE(begin > end); + EXPECT_FALSE(end < begin); +} + +// ==================== RowConstIterator::empty ==================== + +TEST(CompressedRowSparseMatrixConstraint, RowConstIteratorEmpty) +{ + using Vec3 = sofa::type::Vec3; + using Matrix = sofa::linearalgebra::CompressedRowSparseMatrixConstraint; + + Matrix m; + // Row 0 has a column entry + m.writeLine(0).addCol(5, Vec3(1, 2, 3)); + // Row 1 is created in btemp but has no columns: since writeLine + // returns a RowWriteAccessor without adding data, row 1 won't exist + // after compress. Instead, add and then clear to get an empty row. + m.writeLine(1).addCol(0, Vec3(0, 0, 0)); + m.compress(); + + auto it0 = m.readLine(0); + EXPECT_FALSE(it0.empty()); +} + +// ==================== multTransposeBaseVector ==================== + +TEST(CompressedRowSparseMatrixConstraint, multTransposeBaseVector) +{ + using Vec3 = sofa::type::Vec3; + using Matrix = sofa::linearalgebra::CompressedRowSparseMatrixConstraint; + + Matrix m; + // Row 0: col 1 = (1, 0, 0), col 2 = (0, 1, 0) + auto row0 = m.writeLine(0); + row0.addCol(1, Vec3(1, 0, 0)); + row0.addCol(2, Vec3(0, 1, 0)); + // Row 1: col 0 = (0, 0, 1) + auto row1 = m.writeLine(1); + row1.addCol(0, Vec3(0, 0, 1)); + m.compress(); + + // lambda = [2.0, 3.0] + sofa::linearalgebra::FullVector lambda(2); + lambda[0] = 2.0; + lambda[1] = 3.0; + + // res starts as zero + sofa::type::vector res; + res.resize(3, Vec3(0, 0, 0)); + + m.multTransposeBaseVector(res, &lambda); + + // res[1] += (1,0,0) * 2 = (2,0,0) + // res[2] += (0,1,0) * 2 = (0,2,0) + // res[0] += (0,0,1) * 3 = (0,0,3) + constexpr double tol = 1e-10; + EXPECT_NEAR(res[0][0], 0.0, tol); + EXPECT_NEAR(res[0][1], 0.0, tol); + EXPECT_NEAR(res[0][2], 3.0, tol); + EXPECT_NEAR(res[1][0], 2.0, tol); + EXPECT_NEAR(res[1][1], 0.0, tol); + EXPECT_NEAR(res[1][2], 0.0, tol); + EXPECT_NEAR(res[2][0], 0.0, tol); + EXPECT_NEAR(res[2][1], 2.0, tol); + EXPECT_NEAR(res[2][2], 0.0, tol); +} + } // namespace sofa