From a4cb33e96b3745f79dd813be15508335f188a0ff Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 28 Aug 2026 11:30:10 +0900 Subject: [PATCH 01/12] add unit tests on problematic cases --- ...mpressedRowSparseMatrixConstraint_test.cpp | 73 ++++++ .../CompressedRowSparseMatrixGeneric_test.cpp | 172 +++++++++++++ ...mpressedRowSparseMatrixMechanical_test.cpp | 230 ++++++++++++++++++ 3 files changed, 475 insertions(+) diff --git a/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixConstraint_test.cpp b/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixConstraint_test.cpp index 20f425d537b..ce5d4220500 100644 --- a/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixConstraint_test.cpp +++ b/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixConstraint_test.cpp @@ -26,6 +26,8 @@ #include +#include + #include #include #include @@ -2422,4 +2424,75 @@ TEST(CompressedRowSparseMatrixConstraint, multTransposeBaseVector) EXPECT_NEAR(res[2][2], 0.0, tol); } +// ==================== Regression tests ==================== +// +// The tests below currently FAIL. Each one documents a defect found by review; +// the comment above it names the offending line. They are expected to pass once +// the corresponding fix lands. + +// RowConstIterator's defaulted default constructor leaves m_internal and +// m_matrix without initialisers, so a default-constructed iterator is not in the +// invalid state its own isInvalid() is meant to report. Adding NSDMIs +// (= s_invalidIndex / = nullptr) fixes it. +// CompressedRowSparseMatrixConstraint.h:255,431-433 +TEST(CompressedRowSparseMatrixConstraint, DefaultConstructedRowIteratorIsInvalid) +{ + using Vec3 = sofa::type::Vec3; + using Matrix = sofa::linearalgebra::CompressedRowSparseMatrixConstraint; + + Matrix::RowConstIterator it{}; + EXPECT_TRUE(it.isInvalid()) << "a default-constructed row iterator must not alias row 0"; +} + +// clearRowBlock() (inherited from CompressedRowSparseMatrixGeneric) calls +// rowIndex.back() before checking that rowIndex is non-empty, which segfaults +// instead of failing an assertion. CompressedRowSparseMatrixGeneric.h:984-985 +// +// EXPECT_EXIT runs the body in a forked child, so the crash is contained and the +// test binary survives. The child exits 0 only when the call both returns and +// leaves the matrix untouched, so this reports FAILED today ("Terminated by +// signal 11") and PASSED once the guard is added. The suite is named *DeathTest +// per the googletest convention: suites whose name ends in DeathTest are run +// before all others, because forking is only safe before any test starts threads. +#if GTEST_HAS_DEATH_TEST + +TEST(CompressedRowSparseMatrixConstraintDeathTest, ClearRowBlockOnEmptyMatrix) +{ + using Vec3 = sofa::type::Vec3; + using Matrix = sofa::linearalgebra::CompressedRowSparseMatrixConstraint; + + EXPECT_EXIT( + { + Matrix m; + m.clearRowBlock(0); + std::exit(m.empty() ? 0 : 2); + }, + ::testing::ExitedWithCode(0), ""); +} + +#endif // GTEST_HAS_DEATH_TEST + +// Guard (passes today): setLine() on an absent row must not trip over the same +// unguarded back(). +TEST(CompressedRowSparseMatrixConstraint, SetLineOnEmptyMatrix) +{ + using Vec3 = sofa::type::Vec3; + using Matrix = sofa::linearalgebra::CompressedRowSparseMatrixConstraint; + + Matrix src; + src.writeLine(0).addCol(2, Vec3(1, 0, 0)); + src.compress(); + + Matrix dst; + EXPECT_NO_THROW(dst.setLine(0, src.readLine(0).row())); + dst.compress(); + + auto row = dst.readLine(0); + ASSERT_NE(row, dst.end()); + auto col = row.begin(); + ASSERT_NE(col, row.end()); + EXPECT_EQ(col.index(), 2); + EXPECT_EQ(col.val(), Vec3(1, 0, 0)); +} + } // namespace sofa diff --git a/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixGeneric_test.cpp b/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixGeneric_test.cpp index f11671c8805..207069960c0 100644 --- a/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixGeneric_test.cpp +++ b/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixGeneric_test.cpp @@ -23,6 +23,8 @@ #include #include +#include + using CRS = sofa::linearalgebra::CompressedRowSparseMatrixGeneric; using Mat3 = sofa::type::Mat<3, 3, double>; using CRSMat3 = sofa::linearalgebra::CompressedRowSparseMatrixGeneric; @@ -646,3 +648,173 @@ TEST(CompressedRowSparseMatrixGeneric, Mat3x3dMul) EXPECT_NEAR(result(i, j), expected, kTol); } } + +// ==================== Regression tests: known defects ==================== +// +// The tests in this section FAIL today. Each comment names the offending line. +// They are expected to pass once the corresponding fix lands. +// +// The clearRowBlock() cases segfault rather than fail cleanly, so they live in +// the *DeathTest suite further down and run the call in a forked child. + +// block() short-circuits on "is j the first / last column of this row?" without +// first checking that the row range is non-empty. fullRows() registers rows with +// an empty range, after which colsIndex[rowRange.first] reads the *next* row's +// first entry and colsIndex[rowRange.second - 1] reads the *previous* row's last +// entry (index -1 for row 0). CompressedRowSparseMatrixGeneric.h:768-769 +TEST(CompressedRowSparseMatrixGeneric, BlockOnEmptyRowAfterFullRows) +{ + CRS m(6, 6); + m.setBlock(4, 4, 7.0); + m.compress(); + m.fullRows(); + + ASSERT_EQ(m.getRowIndex().size(), 6u); + + for (CRS::Index i = 0; i < 6; ++i) + { + for (CRS::Index j = 0; j < 6; ++j) + { + const double expected = (i == 4 && j == 4) ? 7.0 : 0.0; + EXPECT_NEAR(m.block(i, j), expected, kTol) << "block(" << i << "," << j << ")"; + } + } +} + +// Narrower reproduction of the same defect: reading column 4 of the empty row 0 +// returns the value stored at (4,4). +TEST(CompressedRowSparseMatrixGeneric, BlockOnEmptyRowDoesNotBorrowNeighbourValue) +{ + CRS m(6, 6); + m.setBlock(4, 4, 7.0); + m.compress(); + m.fullRows(); + + EXPECT_NEAR(m.block(0, 4), 0.0, kTol) << "row 0 is empty, must not return row 4's value"; + EXPECT_NEAR(m.block(0, 0), 0.0, kTol) << "row 0 is empty, must not read colsIndex[-1]"; +} + +// clearRowColBlock() computes foundRowId but never uses it: rowRange is built +// from the rowId left behind by a failed sortedFind, and deleteRow(rowId) runs +// whenever *either* index was found. On a matrix where row i is absent but +// column i exists, this deletes an unrelated row. +// CompressedRowSparseMatrixGeneric.h:1074-1108 +TEST(CompressedRowSparseMatrixGeneric, ClearRowColBlockWithAbsentRow) +{ + CRS m(4, 4); + // row 0 is absent, but column 0 exists (in row 1) + m.setBlock(1, 0, 7.0); + m.setBlock(1, 1, 3.0); + m.setBlock(2, 2, 5.0); + m.compress(); + ASSERT_EQ(m.getRowIndex().size(), 2u); + + m.clearRowColBlock(0); + + EXPECT_NEAR(m.block(1, 0), 0.0, kTol) << "column 0 must be cleared"; + EXPECT_NEAR(m.block(1, 1), 3.0, kTol) << "row 1 must not be deleted"; + EXPECT_NEAR(m.block(2, 2), 5.0, kTol) << "row 2 must be untouched"; +} + +// clearRowBlock() calls rowIndex.back() / rowIndex.front() before checking that +// rowIndex is non-empty, which segfaults instead of failing an assertion. +// CompressedRowSparseMatrixGeneric.h:984-985 +// +// EXPECT_EXIT runs the body in a forked child, so the crash is contained and the +// test binary survives to run everything after it. The child exits 0 only when +// the call both returns and leaves the matrix untouched, so this reports FAILED +// today ("Terminated by signal 11") and PASSED once the guard is added. +// The suite is named *DeathTest per the googletest convention: suites whose name +// ends in DeathTest are run before all others, because forking is only safe +// before any test has started threads. +#if GTEST_HAS_DEATH_TEST + +TEST(CompressedRowSparseMatrixGenericDeathTest, ClearRowBlockOnEmptyMatrix) +{ + EXPECT_EXIT( + { + CRS m(4, 4); + m.compress(); + m.clearRowBlock(1); + std::exit(m.getRowIndex().empty() ? 0 : 2); + }, + ::testing::ExitedWithCode(0), ""); +} + +TEST(CompressedRowSparseMatrixGenericDeathTest, ClearRowBlockOnDefaultConstructedMatrix) +{ + EXPECT_EXIT( + { + CRS m; + m.clearRowBlock(0); + std::exit(m.getRowIndex().empty() ? 0 : 2); + }, + ::testing::ExitedWithCode(0), ""); +} + +#endif // GTEST_HAS_DEATH_TEST + +// ============ Regression guards: undefined behaviour, host-dependent ============ +// +// These pass on this host but exercise genuine UB, so they are kept as guards +// rather than as demonstrations of a wrong result: +// - the out-of-range colsIndex / rowBegin reads throw std::logic_error in a +// Debug build, where sofa::type::vector bounds-checks (NDEBUG undefined); +// - the divide-by-zero ones are only caught by UBSan, or by the hardware on +// x86_64 where integer division by zero raises SIGFPE. ARM's UDIV silently +// returns 0, which is why they pass here. + +// clearColBlock() repeats block()'s unguarded first/last-column fast path over +// empty row ranges. CompressedRowSparseMatrixGeneric.h:1017-1018 +TEST(CompressedRowSparseMatrixGeneric, ClearColBlockAfterFullRows) +{ + CRS m(4, 4); + m.setBlock(1, 1, 3.0); + m.compress(); + m.fullRows(); + + m.clearColBlock(1); + + for (CRS::Index i = 0; i < 4; ++i) + EXPECT_NEAR(m.block(i, 1), 0.0, kTol) << "block(" << i << ",1)"; +} + +// Neither the row nor the column exists: clearRowColBlock() still builds a range +// from a stale rowId before reporting the error. +TEST(CompressedRowSparseMatrixGeneric, ClearRowColBlockWithAbsentRowAndCol) +{ + CRS m(4, 4); + m.setBlock(1, 1, 3.0); + m.setBlock(2, 1, 4.0); + m.compress(); + + m.clearRowColBlock(3); + + EXPECT_NEAR(m.block(1, 1), 3.0, kTol); + EXPECT_NEAR(m.block(2, 1), 4.0, kTol); +} + +// The hinted wblock() overload divides by nBlockRow / nBlockCol without the zero +// guard its unhinted sibling has. CompressedRowSparseMatrixGeneric.h:895,904 +TEST(CompressedRowSparseMatrixGeneric, HintedWblockOnEmptyMatrix) +{ + CRS m; + CRS::Index rowId = 0; + CRS::Index colId = 0; + EXPECT_NO_THROW(m.setBlock(0, 0, rowId, colId, 1.0)); +} + +// getMaxColIndex() reads colsIndex[rowBegin[rowId + 1] - 1] for every registered +// row, which points into the previous row when a row is empty. +// CompressedRowSparseMatrixGeneric.h:420-427 +TEST(CompressedRowSparseMatrixGeneric, MaxColIndexWithEmptyRows) +{ + CRS m(4, 4); + m.setBlock(1, 2, 3.0); + m.compress(); + m.fullRows(); + + // exercised through clearColBlock, which calls getMaxColIndex under AutoSize + EXPECT_NO_THROW(m.clearColBlock(2)); + EXPECT_NEAR(m.block(1, 2), 0.0, kTol); +} diff --git a/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixMechanical_test.cpp b/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixMechanical_test.cpp index bc578192f23..d7c90315088 100644 --- a/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixMechanical_test.cpp +++ b/Sofa/framework/LinearAlgebra/test/CompressedRowSparseMatrixMechanical_test.cpp @@ -526,3 +526,233 @@ TEST(CompressedRowSparseMatrixMechanical, Mat3x3dMulVector) EXPECT_NEAR(res[1], 2.0, kTol); EXPECT_NEAR(res[2], 3.0, kTol); } + +// ==================== Regression tests: known defects ==================== +// +// The tests in this section FAIL today. Each comment names the offending line. +// They are expected to pass once the corresponding fix lands. +// +// WARNING: CompressedRowSparseMatrixGeneric.ClearRowBlockOnEmptyMatrix and its +// siblings segfault rather than fail cleanly. Until the fixes land, run with +// --gtest_filter=-*OnEmptyMatrix* +// if you need the rest of the suite to complete. + +// clearRowCol() indexes rowBegin with a block-column number, but rowBegin is +// indexed by the *position* of a row inside rowIndex. The two only coincide when +// every row is present. Here rows 0,1,6,7 are absent, so the search for the +// symmetric block (3,2) is run over row 5's range and finds nothing -- entry +// (3,2) is silently left in column 2. All rowBegin accesses stay in bounds, so +// this reproduces deterministically rather than depending on adjacent memory. +// CompressedRowSparseMatrixMechanical.h:411 +TEST(CompressedRowSparseMatrixMechanical, ClearRowColMissesSymmetricBlockWhenRowsAreSparse) +{ + CRSMech m(8, 8); + m.set(2, 2, 22.0); + m.set(2, 3, 23.0); + m.set(3, 2, 32.0); + m.set(3, 3, 33.0); + m.set(4, 4, 44.0); + m.set(5, 5, 55.0); + m.compress(); + ASSERT_EQ(m.getRowIndex().size(), 4u) << "test needs rowIndex != identity"; + ASSERT_EQ(m.getRowIndex()[0], 2); + + m.clearRowCol(2); + + EXPECT_NEAR(m.element(2, 2), 0.0, kTol) << "row 2"; + EXPECT_NEAR(m.element(2, 3), 0.0, kTol) << "row 2"; + EXPECT_NEAR(m.element(3, 2), 0.0, kTol) << "column 2 must be cleared too"; + + // everything outside row 2 / column 2 must survive + EXPECT_NEAR(m.element(3, 3), 33.0, kTol); + EXPECT_NEAR(m.element(4, 4), 44.0, kTol); + EXPECT_NEAR(m.element(5, 5), 55.0, kTol); +} + +// In clearRowCol(), when the symmetric block (j,i) does not exist the local +// pointer `b` is never reset, so it still points at block (i,j) and the second +// loop zeroes a *column* of that block. Only observable for NL > 1. +// CompressedRowSparseMatrixMechanical.h:417-422 +TEST(CompressedRowSparseMatrixMechanical, ClearRowColAsymmetricPatternMat3) +{ + // 2x2 grid of 3x3 blocks; block (1,0) is deliberately absent so that the + // search for the symmetric block fails. + CRSMechMat3 m(6, 6); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + m.set(i, j, 1.0); // block (0,0) + for (int i = 0; i < 3; ++i) + for (int j = 3; j < 6; ++j) + m.set(i, j, 2.0); // block (0,1) + for (int i = 3; i < 6; ++i) + for (int j = 3; j < 6; ++j) + m.set(i, j, 3.0); // block (1,1) + m.compress(); + + m.clearRowCol(0); + + for (int j = 0; j < 6; ++j) + EXPECT_NEAR(m.element(0, j), 0.0, kTol) << "element(0," << j << ")"; + for (int i = 0; i < 6; ++i) + EXPECT_NEAR(m.element(i, 0), 0.0, kTol) << "element(" << i << ",0)"; + + // block (0,1) rows 1 and 2 are neither in row 0 nor in column 0 and must survive + EXPECT_NEAR(m.element(1, 3), 2.0, kTol); + EXPECT_NEAR(m.element(2, 3), 2.0, kTol); + EXPECT_NEAR(m.element(1, 4), 2.0, kTol); + EXPECT_NEAR(m.element(2, 5), 2.0, kTol); + + // block (1,1) is untouched + EXPECT_NEAR(m.element(3, 3), 3.0, kTol); + EXPECT_NEAR(m.element(5, 5), 3.0, kTol); +} + +// element() delegates to block(), whose "is j the first / last column of this +// row?" fast paths dereference colsIndex without first checking that the row +// range is non-empty. fullRows() registers exactly such empty ranges, after +// which element() returns a neighbouring row's value or reads colsIndex[-1]. +// CompressedRowSparseMatrixGeneric.h:768-769 +TEST(CompressedRowSparseMatrixMechanical, ElementAfterFullRows) +{ + CRSMech m(6, 6); + m.set(4, 4, 7.0); + m.compress(); + m.fullRows(); + + for (int i = 0; i < 6; ++i) + { + for (int j = 0; j < 6; ++j) + { + const double expected = (i == 4 && j == 4) ? 7.0 : 0.0; + EXPECT_NEAR(m.element(i, j), expected, kTol) << "element(" << i << "," << j << ")"; + } + } +} + +// filterValues() only emits a row when it produced at least one value, so the +// keepEmptyRows flag threaded through every copy*() wrapper has no effect. The +// compensating pop_back() at line 591 is dead: rowBegin.back() == vid cannot +// hold there. CompressedRowSparseMatrixMechanical.h:585,591 +// +// The flag concerns destination rows emptied *by the filter*, not source rows +// that were never registered: filterValues only visits srcMatrix.rowIndex. +// CRSMechanicalPolicy has CompressZeros == false, so the explicitly stored zero +// at (2,2) survives compression on the source side and gives the nonzeros filter +// a row to empty out. +TEST(CompressedRowSparseMatrixMechanical, CopyNonZerosKeepEmptyRows) +{ + CRSMech src(4, 4); + src.set(1, 1, 2.0); + src.set(2, 2, 0.0); + src.compress(); + ASSERT_EQ(src.getRowIndex().size(), 2u) << "the stored zero must keep row 2 on the source"; + + CRSMech dst; + dst.copyNonZeros(src, /*keepEmptyRows*/ true); + + ASSERT_EQ(dst.getRowIndex().size(), 2u) << "row 2 was emptied by the filter but must be kept"; + EXPECT_EQ(dst.getRowIndex()[0], 1); + EXPECT_EQ(dst.getRowIndex()[1], 2); + + // the kept row must be genuinely empty, and the surviving value intact + const auto range = dst.getRowRange(1); + EXPECT_TRUE(range.empty()) << "row 2 should hold no block"; + EXPECT_NEAR(dst.element(1, 1), 2.0, kTol); +} + +// Same source, default flag: the row emptied by the filter is dropped. +TEST(CompressedRowSparseMatrixMechanical, CopyNonZerosDropsRowEmptiedByFilter) +{ + CRSMech src(4, 4); + src.set(1, 1, 2.0); + src.set(2, 2, 0.0); + src.compress(); + + CRSMech dst; + dst.copyNonZeros(src); + + ASSERT_EQ(dst.getRowIndex().size(), 1u); + EXPECT_EQ(dst.getRowIndex()[0], 1); + EXPECT_NEAR(dst.element(1, 1), 2.0, kTol); +} + +// ============ Regression guards: undefined behaviour, host-dependent ============ +// +// These pass on this host but exercise genuine UB, so they are kept as guards +// rather than as demonstrations of a wrong result: +// - the out-of-range colsIndex / rowBegin reads throw std::logic_error in a +// Debug build, where sofa::type::vector bounds-checks (NDEBUG undefined); +// - the divide-by-zero ones are only caught by UBSan, or by the hardware on +// x86_64 where integer division by zero raises SIGFPE. ARM's UDIV silently +// returns 0, which is why they pass here. + +// The `i * rowIndex.size() / nBlockRow` search hint is guarded against +// nBlockRow == 0 in block() and wblock(), but the guard was dropped in every +// copy below. Confirmed with UBSan at Mechanical.h:343, :748 and :841. +TEST(CompressedRowSparseMatrixMechanical, ClearRowOnEmptyMatrix) +{ + CRSMech m; + EXPECT_NO_THROW(m.clearRow(3)); + EXPECT_EQ(m.getRowIndex().size(), 0u); +} + +TEST(CompressedRowSparseMatrixMechanical, ClearColOnEmptyMatrix) +{ + CRSMech m; + EXPECT_NO_THROW(m.clearCol(3)); + EXPECT_EQ(m.getRowIndex().size(), 0u); +} + +TEST(CompressedRowSparseMatrixMechanical, ClearRowColOnEmptyMatrix) +{ + CRSMech m; + EXPECT_NO_THROW(m.clearRowCol(3)); + EXPECT_EQ(m.getRowIndex().size(), 0u); +} + +TEST(CompressedRowSparseMatrixMechanical, BlockAccessorsOnEmptyMatrix) +{ + CRSMech m; + EXPECT_NO_THROW((void) m.blockGet(3, 3)); + EXPECT_NO_THROW((void) m.blockGetW(3, 3)); + EXPECT_NO_THROW((void) m.blockCreate(3, 3)); +} + +TEST(CompressedRowSparseMatrixMechanical, BRowIteratorsOnEmptyMatrix) +{ + CRSMech m; + EXPECT_NO_THROW((void) m.bRowBegin(3)); + EXPECT_NO_THROW((void) m.bRowEnd(3)); + EXPECT_NO_THROW((void) m.bRowRange(3)); +} + +// clearCol() -> clearColBlock() repeats block()'s unguarded fast path. Under the +// ClearByZeros policy the stray index still lands on an entry of column j, so +// the result happens to be correct; only the out-of-range read is wrong. +// CompressedRowSparseMatrixGeneric.h:1017-1018 +TEST(CompressedRowSparseMatrixMechanical, ClearColAfterFullRows) +{ + CRSMech m(4, 4); + m.set(1, 1, 3.0); + m.compress(); + m.fullRows(); + + m.clearCol(1); + + for (int i = 0; i < 4; ++i) + EXPECT_NEAR(m.element(i, 1), 0.0, kTol) << "element(" << i << ",1)"; +} + +// Control case: the default must still drop empty rows. Passes today. +TEST(CompressedRowSparseMatrixMechanical, CopyNonZerosDropsEmptyRowsByDefault) +{ + CRSMech src(4, 4); + src.set(1, 1, 2.0); + src.compress(); + + CRSMech dst; + dst.copyNonZeros(src); + + EXPECT_EQ(dst.getRowIndex().size(), 1u); + EXPECT_NEAR(dst.element(1, 1), 2.0, kTol); +} From 11e0466d9b912ef50fbc9b2413341d374cd4a131 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 28 Aug 2026 12:20:26 +0900 Subject: [PATCH 02/12] [LinearAlgebra] Fix block() returning a neighbouring row's value block() short-circuits on "is j the first or the last column of this row?" without first checking that the row range is non-empty. A registered row may hold no block at all -- fullRows() and fullDiagonal() both create such rows -- and then colsIndex[rowRange.first] reads the next row's first entry while colsIndex[rowRange.second - 1] reads the previous row's last one, which is index -1 for the first row. Reading any empty row therefore returned another row's value, or read out of bounds. On a 6x6 matrix holding only (4,4) = 7, block(0,4) returned 7 and block(0,0) returned uninitialised memory. The same row lookup and column lookup are open-coded at a dozen sites, so extract them into findRow() and findColInRange() and use them here. Both helpers fold in the guards the copies were missing: the empty-range check, and the nBlockRow == 0 / nBlockCol == 0 division guard. Fixes CompressedRowSparseMatrixGeneric.BlockOnEmptyRowAfterFullRows, CompressedRowSparseMatrixGeneric.BlockOnEmptyRowDoesNotBorrowNeighbourValue and CompressedRowSparseMatrixMechanical.ElementAfterFullRows. Co-Authored-By: Claude Opus 5 (cherry picked from commit 2b234b6ec19dc962d0a24590c59e390869c2c27e) --- .../CompressedRowSparseMatrixGeneric.h | 73 +++++++++++++++---- 1 file changed, 58 insertions(+), 15 deletions(-) diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h index 066e72b8604..8f288f022be 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h @@ -428,6 +428,62 @@ public : return maxColIndex; } + /** + * \brief Look for column j among the blocks of a single row. + * + * The first and last registered columns of the row are checked directly, as + * they are by far the most common queries; anything else falls back to a + * binary search. The range must be checked for emptiness first: a registered + * row may hold no block at all (fullRows() and fullDiagonal() both create + * such rows), in which case rowRange.begin() addresses the next row's first + * block and rowRange.end() - 1 the previous row's last one. + * + * @param rowRange : range of this row inside colsIndex / colsValue + * @param j : column index to look for + * @param colId : on success, position of the block in colsIndex / colsValue + * @return true if the column holds a block in this row + **/ + bool findColInRange(const Range& rowRange, Index j, Index& colId) const + { + if (rowRange.empty()) return false; + if (j == colsIndex[rowRange.begin()]) + { + colId = rowRange.begin(); + return true; + } + if (j == colsIndex[rowRange.end() - 1]) + { + colId = rowRange.end() - 1; + return true; + } + colId = (nBlockCol == 0) ? rowRange.begin() + : rowRange.begin() + j * rowRange.size() / nBlockCol; + return sortedFind(colsIndex, rowRange, j, colId); + } + + /** + * \brief Look for row i and return its position inside rowIndex. + * @param i : row index + * @param rowId : on success, position of the row inside rowIndex + * @return true if the row is registered + **/ + bool findRow(Index i, Index& rowId) const + { + if (rowIndex.empty()) return false; + if (i == rowIndex.back()) + { + rowId = Index(rowIndex.size() - 1); + return true; + } + if (i == rowIndex.front()) + { + rowId = 0; + return true; + } + rowId = (nBlockRow == 0) ? 0 : Index(i * rowIndex.size() / nBlockRow); + return sortedFind(rowIndex, i, rowId); + } + /** * \brief Method to easy delete row given position in rowIndex. * @param RowId position on line in rowIndex @@ -755,23 +811,10 @@ public : if constexpr (Policy::AutoSize) if (j > nBlockCol) return empty; /// Matrix is auto sized so requested column could not exist Index rowId = 0; - if (i == rowIndex.back()) rowId = Index(rowIndex.size() - 1); /// Optimization to avoid do a find when looking for the last line registred - else if (i == rowIndex.front()) rowId = 0; /// Optimization to avoid do a find when looking for the first line registred - else - { - rowId = (nBlockRow == 0) ? 0 : Index(i * rowIndex.size() / nBlockRow); - if (!sortedFind(rowIndex, i, rowId)) return empty; - } + if (!findRow(i, rowId)) return empty; - Range rowRange(rowBegin[rowId], rowBegin[rowId+1]); Index colId = 0; - if (j == colsIndex[rowRange.first]) colId = rowRange.first; /// Optimization to avoid do a find when looking for the first column registred for specific column - else if (j == colsIndex[rowRange.second - 1]) colId = rowRange.second - 1; /// Optimization to avoid do a find when looking for the last column registred for specific column - else - { - colId = (nBlockCol == 0) ? 0 : rowRange.begin() + j * rowRange.size() / nBlockCol; - if (!sortedFind(colsIndex, rowRange, j, colId)) return empty; - } + if (!findColInRange(Range(rowBegin[rowId], rowBegin[rowId+1]), j, colId)) return empty; return colsValue[colId]; } From ccdbb2607b76dab108ad326d768cc11704b0cbe8 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 28 Aug 2026 12:21:36 +0900 Subject: [PATCH 03/12] [LinearAlgebra] Fix clearRowBlock() segfaulting on an empty matrix clearRowBlock() called rowIndex.back() and rowIndex.front() as fast paths before checking that rowIndex was non-empty, so clearing a row of a matrix with no registered row dereferenced past a null pointer and segfaulted. findRow() already reports an empty matrix as "row not found", so use it here instead of the open-coded lookup. Fixes CompressedRowSparseMatrixGenericDeathTest.ClearRowBlockOnEmptyMatrix, CompressedRowSparseMatrixGenericDeathTest.ClearRowBlockOnDefaultConstructedMatrix and CompressedRowSparseMatrixConstraintDeathTest.ClearRowBlockOnEmptyMatrix. Co-Authored-By: Claude Opus 5 (cherry picked from commit 14d101ebcb2fa13c352e2ec298e54c8137c0b40b) --- .../sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h index 8f288f022be..376901a220c 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h @@ -1024,13 +1024,7 @@ public : } Index rowId = 0; - if (i == rowIndex.back()) rowId = Index(rowIndex.size() - 1); /// Optimization to avoid do a find when looking for the last line registred - else if (i == rowIndex.front()) rowId = 0; /// Optimization to avoid do a find when looking for the first line registred - else - { - rowId = (nBlockRow == 0) ? 0 : Index(i * rowIndex.size() / nBlockRow); - if (!sortedFind(rowIndex, i, rowId)) return; - } + if (!findRow(i, rowId)) return; /// Nothing to clear: the matrix is empty or the row holds no block deleteRow(rowId); From 6d10c501abf1129c9de8f906364c2420dc2aa048 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 28 Aug 2026 12:22:58 +0900 Subject: [PATCH 04/12] [LinearAlgebra] Fix clearRowColBlock() deleting an unrelated row clearRowColBlock() computed foundRowId but never acted on it. When the row lookup failed, rowId kept whatever value the failed binary search left behind, and that value was then used to build rowRange -- reading out of bounds -- and passed straight to deleteRow(). On a matrix where row i is absent but column i exists, this dropped a completely unrelated row: clearing row/column 0 of a matrix holding (1,0), (1,1) and (2,2) deleted row 1, losing (1,1). An absent row is normal in a sparse matrix and column i may still hold blocks in other rows, so it is not an error either. Reserve the diagnostic for an index that is genuinely outside the matrix, delete the row only when findRow() locates it, and always clear the column. clearColBlock() carried its own copy of the unguarded column lookup, so route it through findColInRange() as well since clearRowColBlock() delegates to it. Fixes CompressedRowSparseMatrixGeneric.ClearRowColBlockWithAbsentRow. Co-Authored-By: Claude Opus 5 (cherry picked from commit 3275851fb0e1c693bb1f4828221867d26aa44ece) --- .../CompressedRowSparseMatrixGeneric.h | 43 ++++++------------- 1 file changed, 12 insertions(+), 31 deletions(-) diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h index 376901a220c..898e0551593 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h @@ -1050,15 +1050,8 @@ public : { Range rowRange(rowBegin[rowId], rowBegin[rowId+1]); - Index colId = -1; - if (j == colsIndex[rowRange.first]) colId = rowRange.first; /// Optimization to avoid do a find when looking for the first column registred for specific column - else if (j == colsIndex[rowRange.second - 1]) colId = rowRange.second - 1; /// Optimization to avoid do a find when looking for the last column registred for specific column - else - { - colId = (nBlockCol == 0) ? 0 : rowRange.begin() + j * rowRange.size() / nBlockCol; - if (!sortedFind(colsIndex, rowRange, j, colId)) colId = -1; - } - if (colId != -1) /// Means col exist in this line + Index colId = 0; + if (findColInRange(rowRange, j, colId)) /// Means col exist in this line { if constexpr (Policy::ClearByZeros) { @@ -1113,34 +1106,22 @@ public : /// If AutoCompress policy is activated, we neeed to be sure not missing btemp registered value. if constexpr (Policy::AutoCompress) compress(); - bool foundRowId = true; - Index rowId = 0; - if (i == rowIndex.back()) rowId = rowIndex.size() - 1; /// Optimization to avoid do a find when looking for the last line registred - else if (i == rowIndex.front()) rowId = 0; /// Optimization to avoid do a find when looking for the first line registred - else - { - rowId = (nBlockRow == 0) ? 0 : i * rowIndex.size() / nBlockRow; - if (!sortedFind(rowIndex, i, rowId)) foundRowId = false; - } - - bool foundColId = true; - Range rowRange(rowBegin[rowId], rowBegin[rowId+1]); - Index colId = 0; - if (i == colsIndex[rowRange.first]) colId = rowRange.first; /// Optimization to avoid do a find when looking for the first column registred for specific column - else if (i == colsIndex[rowRange.second - 1]) colId = rowRange.second - 1; /// Optimization to avoid do a find when looking for the last column registred for specific column - else + if (i < 0 || i >= nBlockRow || i >= nBlockCol) { - colId = (nBlockCol == 0) ? 0 : rowRange.begin() + i * rowRange.size() / nBlockCol; - if (!sortedFind(colsIndex, rowRange, i, colId)) foundColId = false;; + msg_error("CompressedRowSparseMatrixGeneric") << "invalid write access to row and column "<Name() << " of size ("<Name() << " of size ("< Date: Fri, 28 Aug 2026 12:23:42 +0900 Subject: [PATCH 05/12] [LinearAlgebra] Fix clearRowCol() indexing rowBegin with a column number Two defects in the symmetric half of clearRowCol(). First, having read the block column j out of colsIndex, it built the symmetric block's range as rowBegin[j] .. rowBegin[j + 1]. rowBegin is indexed by the position of a row inside rowIndex, not by the row number; the two coincide only when every row is present. On a sparse matrix this searched the wrong row -- reading out of bounds when j exceeded the number of registered rows -- and so failed to clear entries that are in column i. Clearing row/column 2 of an 8x8 matrix holding rows 2..5 searched row 5's range for the symmetric block and left (3,2) untouched. Second, when the symmetric block (j,i) was not found, the local pointer was never reset, so the second loop zeroed a column of block (i,j) instead -- entries that lie in neither row i nor column i. Only observable for NL > 1; with 3x3 blocks, clearing row/column 0 wiped (1,3) and (2,3). Look row j up with findRow() before taking its range, locate the symmetric block with findColInRange(), and skip the second clear entirely when it does not exist. Handle the diagonal block explicitly, since it is its own symmetric counterpart. This also drops an unguarded division by nBlockRow. Fixes CompressedRowSparseMatrixMechanical.ClearRowColMissesSymmetricBlockWhenRowsAreSparse and CompressedRowSparseMatrixMechanical.ClearRowColAsymmetricPatternMat3. Co-Authored-By: Claude Opus 5 (cherry picked from commit 230d96cf49e3d0d895468d481e4f728cd20ab517) --- .../CompressedRowSparseMatrixMechanical.h | 58 +++++++++++-------- 1 file changed, 33 insertions(+), 25 deletions(-) diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixMechanical.h b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixMechanical.h index bfa7345be82..503abf8a2cf 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixMechanical.h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixMechanical.h @@ -392,36 +392,44 @@ class CompressedRowSparseMatrixMechanical final // final is used to allow the co if constexpr (Policy::AutoCompress) this->compress(); Index bi=0; split_row_index(i, bi); - Index rowId = Index(i * this->rowIndex.size() / this->nBlockRow); - if (this->sortedFind(this->rowIndex, i, rowId)) - { - Range rowRange(this->rowBegin[rowId], this->rowBegin[rowId+1]); - for (Index xj = rowRange.begin(); xj < rowRange.end(); ++xj) - { - Block* b = &this->colsValue[xj]; - // first clear (i,j) - for (Index bj = 0; bj < (Index)NC; ++bj) - traits::vset(*b, bi, bj, 0); - // then clear (j,i) - Index j = this->colsIndex[xj]; - - if (j != i) - { - Range jrowRange(this->rowBegin[j], this->rowBegin[j + 1]); - Index colId = 0; + Index rowId = 0; + if (!this->findRow(i, rowId)) return; - // look for column i - if (this->sortedFind(this->colsIndex, jrowRange, i, colId)) - { - b = &this->colsValue[colId]; - } - } + const Range rowRange(this->rowBegin[rowId], this->rowBegin[rowId+1]); + for (Index xj = rowRange.begin(); xj < rowRange.end(); ++xj) + { + // first clear (i,j) + Block& b = this->colsValue[xj]; + for (Index bj = 0; bj < (Index)NC; ++bj) + traits::vset(b, bi, bj, 0); + const Index j = this->colsIndex[xj]; + + // the diagonal block is its own symmetric counterpart + if (j == i) + { for (Index bj = 0; bj < (Index)NL; ++bj) - traits::vset(*b, bj, bi, 0); - + traits::vset(b, bj, bi, 0); + continue; } + + // then clear (j,i), when that block exists. rowBegin is indexed + // by the position of a row inside rowIndex, not by the row + // number, so row j has to be looked up rather than used as an + // index -- the two only coincide when every row is present. + Index jRowId = 0; + if (!this->findRow(j, jRowId)) continue; + + Index colId = 0; + if (!this->findColInRange(Range(this->rowBegin[jRowId], this->rowBegin[jRowId+1]), i, colId)) continue; + + // never fall back to block (i,j) when (j,i) is missing: clearing + // a column of it would zero entries that lie in neither row i + // nor column i + Block& bSym = this->colsValue[colId]; + for (Index bj = 0; bj < (Index)NL; ++bj) + traits::vset(bSym, bj, bi, 0); } } } From 8bf0533aa8f7a1df24ce3e295db1cea37433aa46 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 28 Aug 2026 12:24:31 +0900 Subject: [PATCH 06/12] [LinearAlgebra] Make filterValues() honour keepEmptyRows filterValues() only pushed a destination row when that row had produced at least one value, so the keepEmptyRows flag threaded through all eight copy*() wrappers did nothing. The compensating pop_back() meant to undo the push could never fire either: rows are only pushed once oldVid != vid, so rowBegin.back() is always strictly less than vid at that point. Push the row when the caller asked to keep empty ones, and drop the dead pop_back(). The flag concerns destination rows emptied by the filter itself, not source rows that were never registered: filterValues() only iterates srcMatrix.rowIndex. No caller in the tree passes keepEmptyRows; every call uses the default false. Fixes CompressedRowSparseMatrixMechanical.CopyNonZerosKeepEmptyRows. Co-Authored-By: Claude Opus 5 --- .../CompressedRowSparseMatrixMechanical.h | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixMechanical.h b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixMechanical.h index 503abf8a2cf..c95c3acfdc1 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixMechanical.h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixMechanical.h @@ -589,18 +589,14 @@ class CompressedRowSparseMatrixMechanical final // final is used to allow the co } } - if (oldVid != vid) //check in case all sub-blocks have been filtered out + // a destination row left empty by the filter is dropped unless + // the caller asked to keep it + if (oldVid != vid || keepEmptyRows) { this->rowIndex.push_back(scalarRowId / DstBlockRows + subRow); this->rowBegin.push_back(oldVid); } } - - if (!keepEmptyRows && !this->rowBegin.empty() && this->rowBegin.back() == vid) // row was empty - { - this->rowIndex.pop_back(); - this->rowBegin.pop_back(); - } } this->rowBegin.push_back(vid); // end of last row } From a20dc8d9c3edf9796cf70fbefe49ba5fc30c3f01 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 28 Aug 2026 12:25:41 +0900 Subject: [PATCH 07/12] [LinearAlgebra] Initialise RowConstIterator's members RowConstIterator's defaulted default constructor left m_internal and m_matrix without initialisers. A default-constructed iterator therefore held an indeterminate row position and a dangling matrix pointer, and its own isInvalid() reported false -- so it read as a valid iterator onto row 0 of an unspecified matrix. Give both members a default member initialiser so the default-constructed state is the invalid one the class already knows how to describe. Fixes CompressedRowSparseMatrixConstraint.DefaultConstructedRowIteratorIsInvalid. Co-Authored-By: Claude Opus 5 --- .../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 2e8ae05181a..0187e3fe7ad 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixConstraint.h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixConstraint.h @@ -430,8 +430,10 @@ class CompressedRowSparseMatrixConstraint : public sofa::linearalgebra::Compress private: - Index m_internal; - const CompressedRowSparseMatrixConstraint* m_matrix; + /// a default-constructed iterator must report itself as invalid rather + /// than alias row 0 of an unspecified matrix + Index m_internal = s_invalidIndex; + const CompressedRowSparseMatrixConstraint* m_matrix = nullptr; }; /// Get the iterator corresponding to the beginning of the rows of blocks From 93ba262013cbf9a2f26eaebb26b2fe8f6f4e3b5c Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 28 Aug 2026 12:27:46 +0900 Subject: [PATCH 08/12] [LinearAlgebra] Route the remaining row/column lookups through the guarded helpers The "find row i" and "find column j inside a row" idioms were open-coded at a dozen more sites. Each copy carried the same two hazards the previous commits fixed at individual call sites: an unguarded division by nBlockRow / nBlockCol, undefined on a default-constructed matrix, and the first/last-column fast paths dereferencing colsIndex without checking the row range for emptiness. The division silently yields 0 on ARM, which is why the accompanying guard tests pass there, but integer division by zero raises SIGFPE on x86_64. Confirmed with UBSan at Mechanical.h:343 (clearRow), :748 (blockGet) and :841 (bRowBegin) before this change, clean after. Converted: clearRow, blockGet, blockGetW, blockCreate, bRowBegin, bRowEnd and bRowRange in the mechanical matrix; both wblock() overloads -- including the hinted one, whose divisions had no guard at all -- and getMaxColIndex in the generic one. The ordered-insertion append path now treats an empty range on the last row as "append here" rather than reading colsIndex[-1]. No behaviour change on well-formed input; this is the undefined behaviour left over once the six defects with failing tests were fixed. Co-Authored-By: Claude Opus 5 --- .../CompressedRowSparseMatrixGeneric.h | 48 +++++++------------ .../CompressedRowSparseMatrixMechanical.h | 40 ++++++++-------- 2 files changed, 37 insertions(+), 51 deletions(-) diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h index 898e0551593..21234772623 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h @@ -422,6 +422,9 @@ public : Index maxColIndex = 0; for (Index rowId = 0; rowId < static_cast(rowIndex.size()); rowId++) { + /// a registered row may hold no block, in which case rowBegin[rowId+1] - 1 + /// would address the previous row's last column + if (rowBegin[rowId] == rowBegin[rowId+1]) continue; Index lastColIndex = colsIndex[rowBegin[rowId+1] - 1]; if (lastColIndex > maxColIndex) maxColIndex = lastColIndex; } @@ -850,11 +853,10 @@ public : { Index rowId = Index(rowIndex.size() - 1); Range rowRange(rowBegin[rowId], rowBegin[rowId+1]); - if (j == colsIndex[rowRange.second - 1]) /// In this case, we are trying to write on last registered column, directly return ref on it - { - return &colsValue[rowRange.second - 1]; - } - else if (j > colsIndex[rowRange.second - 1]) /// Optimization we are trying to write on last line et upper of last column, directly create it. + /// A registered row may hold no block at all, in which case + /// rowRange.end() - 1 would address the previous row's last one. + /// Appending is then the right move, as this is the last row. + if (rowRange.empty() || j > colsIndex[rowRange.end() - 1]) /// Optimization we are trying to write on last line et upper of last column, directly create it. { if (!create) return nullptr; colsIndex.push_back(j); @@ -868,8 +870,8 @@ public : } else { - Index colId = (nBlockCol == 0) ? 0 : rowRange.begin() + j * rowRange.size() / nBlockCol; - if (!sortedFind(colsIndex, rowRange, j, colId)) return create ? insertBtemp(i,j) : nullptr; + Index colId = 0; + if (!findColInRange(rowRange, j, colId)) return create ? insertBtemp(i,j) : nullptr; return &colsValue[colId]; } } @@ -877,34 +879,20 @@ public : if constexpr (Policy::AutoSize) if (j > nBlockCol) return create ? insertBtemp(i,j) : nullptr; /// Matrix is auto sized so requested column could not exist Index rowId = 0; - if (i == rowIndex.back()) rowId = Index(rowIndex.size() - 1); /// Optimization to avoid do a find when looking for the last line registred - else if (i == rowIndex.front()) rowId = 0; /// Optimization to avoid do a find when looking for the first line registred - else - { - rowId = (nBlockRow == 0) ? 0 : Index(i * rowIndex.size() / nBlockRow); - if (!sortedFind(rowIndex, i, rowId)) return create ? insertBtemp(i,j) : nullptr; - } + if (!findRow(i, rowId)) return create ? insertBtemp(i,j) : nullptr; - Range rowRange(rowBegin[rowId], rowBegin[rowId+1]); Index colId = 0; - if (j == colsIndex[rowRange.first]) colId = rowRange.first; /// Optimization to avoid do a find when looking for the first column registred for specific column - else if (j == colsIndex[rowRange.second - 1]) colId = rowRange.second - 1; /// Optimization to avoid do a find when looking for the last column registred for specific column - else - { - colId = (nBlockCol == 0) ? 0 : rowRange.begin() + j * rowRange.size() / nBlockCol; - if (!sortedFind(colsIndex, rowRange, j, colId)) return create ? insertBtemp(i,j) : nullptr; - } + if (!findColInRange(Range(rowBegin[rowId], rowBegin[rowId+1]), j, colId)) return create ? insertBtemp(i,j) : nullptr; return &colsValue[colId]; } else { - Index rowId = (nBlockRow == 0) ? 0 : Index(i * rowIndex.size() / nBlockRow); - if (sortedFind(rowIndex, i, rowId)) + Index rowId = 0; + if (findRow(i, rowId)) { - Range rowRange(rowBegin[rowId], rowBegin[rowId+1]); - Index colId = (nBlockCol == 0) ? 0 : rowRange.begin() + j * rowRange.size() / nBlockCol; - if (sortedFind(colsIndex, rowRange, j, colId)) + Index colId = 0; + if (findColInRange(Range(rowBegin[rowId], rowBegin[rowId+1]), j, colId)) { return &colsValue[colId]; } @@ -935,8 +923,7 @@ public : bool rowFound = true; if (rowId < 0 || rowId >= static_cast(rowIndex.size()) || rowIndex[rowId] != i) { - rowId = Index(i * rowIndex.size() / nBlockRow); - rowFound = sortedFind(rowIndex, i, rowId); + rowFound = findRow(i, rowId); } if (rowFound) { @@ -944,8 +931,7 @@ public : Range rowRange(rowBegin[rowId], rowBegin[rowId+1]); if (colId < rowRange.begin() || colId >= rowRange.end() || colsIndex[colId] != j) { - colId = rowRange.begin() + j * rowRange.size() / nBlockCol; - colFound = sortedFind(colsIndex, rowRange, j, colId); + colFound = findColInRange(rowRange, j, colId); } if (colFound) { diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixMechanical.h b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixMechanical.h index c95c3acfdc1..2b07219eb95 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixMechanical.h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixMechanical.h @@ -340,8 +340,8 @@ class CompressedRowSparseMatrixMechanical final // final is used to allow the co if constexpr (Policy::AutoCompress) this->compress(); /// If AutoCompress policy is activated, we neeed to be sure not missing btemp registered value. Index bi=0; split_row_index(i, bi); - Index rowId = Index(i * this->rowIndex.size() / this->nBlockRow); - if (this->sortedFind(this->rowIndex, i, rowId)) + Index rowId = 0; + if (this->findRow(i, rowId)) { Range rowRange(this->rowBegin[rowId], this->rowBegin[rowId+1]); for (Index xj = rowRange.begin(); xj < rowRange.end(); ++xj) @@ -749,12 +749,12 @@ 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 ! - Index rowId = Index(i * this->rowIndex.size() / this->nBlockRow); - if (this->sortedFind(this->rowIndex, i, rowId)) + Index rowId = 0; + if (this->findRow(i, rowId)) { Range rowRange(this->rowBegin[rowId], this->rowBegin[rowId+1]); - Index colId = rowRange.begin() + j * rowRange.size() / this->nBlockCol; - if (this->sortedFind(this->colsIndex, rowRange, j, colId)) + Index colId = 0; + if (this->findColInRange(rowRange, j, colId)) { return createBlockConstAccessor(i, j, colId); } @@ -767,12 +767,12 @@ class CompressedRowSparseMatrixMechanical final // final is used to allow the co { if constexpr (Policy::AutoCompress) compress(); - Index rowId = Index(i * this->rowIndex.size() / this->nBlockRow); - if (this->sortedFind(this->rowIndex, i, rowId)) + Index rowId = 0; + if (this->findRow(i, rowId)) { Range rowRange(this->rowBegin[rowId], this->rowBegin[rowId+1]); - Index colId = rowRange.begin() + j * rowRange.size() / this->nBlockCol; - if (this->sortedFind(this->colsIndex, rowRange, j, colId)) + Index colId = 0; + if (this->findColInRange(rowRange, j, colId)) { return createBlockAccessor(i, j, colId); } @@ -783,12 +783,12 @@ class CompressedRowSparseMatrixMechanical final // final is used to allow the co /// Get write access to a block, possibly creating it BlockAccessor blockCreate(Index i, Index j) { - Index rowId = Index(i * this->rowIndex.size() / this->nBlockRow); - if (this->sortedFind(this->rowIndex, i, rowId)) + Index rowId = 0; + if (this->findRow(i, rowId)) { Range rowRange(this->rowBegin[rowId], this->rowBegin[rowId+1]); - Index colId = rowRange.begin() + j * rowRange.size() / this->nBlockCol; - if (this->sortedFind(this->colsIndex, rowRange, j, colId)) + Index colId = 0; + if (this->findColInRange(rowRange, j, colId)) { return createBlockAccessor(i, j, colId); } @@ -842,9 +842,9 @@ class CompressedRowSparseMatrixMechanical final // final is used to allow the co ColBlockConstIterator bRowBegin(Index ib) const override { if constexpr (Policy::AutoCompress) const_cast(this)->compress(); /// \warning this violates the const-ness of the method ! - Index rowId = Index(ib * this->rowIndex.size() / this->nBlockRow); + Index rowId = 0; Index index = 0; - if (this->sortedFind(this->rowIndex, ib, rowId)) + if (this->findRow(ib, rowId)) { index = this->rowBegin[rowId]; } @@ -855,9 +855,9 @@ class CompressedRowSparseMatrixMechanical final // final is used to allow the co ColBlockConstIterator bRowEnd(Index ib) const override { if constexpr (Policy::AutoCompress) const_cast(this)->compress(); /// \warning this violates the const-ness of the method ! - Index rowId = Index(ib * this->rowIndex.size() / this->nBlockRow); + Index rowId = 0; Index index2 = 0; - if (this->sortedFind(this->rowIndex, ib, rowId)) + if (this->findRow(ib, rowId)) { index2 = this->rowBegin[rowId+1]; } @@ -868,9 +868,9 @@ class CompressedRowSparseMatrixMechanical final // final is used to allow the co std::pair bRowRange(Index ib) const override { if constexpr (Policy::AutoCompress) const_cast(this)->compress(); /// \warning this violates the const-ness of the method ! - Index rowId = Index(ib * this->rowIndex.size() / this->nBlockRow); + Index rowId = 0; Index index = 0, index2 = 0; - if (this->sortedFind(this->rowIndex, ib, rowId)) + if (this->findRow(ib, rowId)) { index = this->rowBegin[rowId]; index2 = this->rowBegin[rowId+1]; From d3aac9086ede8df4745672afd0e8dfd9f928a0aa Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 28 Aug 2026 13:28:28 +0900 Subject: [PATCH 09/12] [LinearAlgebra] Remove the unused and non-functional check_matrix() check_matrix() has no callers anywhere in the tree, and would not work if it had any: - the non-static overload takes &rowBegin[0], &colsIndex[0] and &colsValue[0] without checking the vectors are non-empty; - it walks a_p from 1 to m, where m is rowBSize(), but rowBegin holds rowIndex.size() + 1 entries -- the two only match after fullRows(); - it requires a_p to increase strictly, so any empty row is reported as an error; - the column loop increments i both in the for-header and in the body, so most entries are never checked; - one diagnostic prints a_p[i] where it means a_i[i]; - and success is reported through msg_error. Delete it rather than repair a validator nothing calls. Co-Authored-By: Claude Opus 5 --- .../CompressedRowSparseMatrixGeneric.h | 86 ------------------- 1 file changed, 86 deletions(-) diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h index 21234772623..c937fd5bd6f 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h @@ -1475,92 +1475,6 @@ public : return name.c_str(); } - bool check_matrix() - { - return check_matrix( - Index(this->getColsValue().size()), - this->rowBSize(), - this->colBSize(), - static_cast (&(rowBegin[0])), - static_cast (&(colsIndex[0])), - static_cast (&(colsValue[0])) - ); - } - - static bool check_matrix( - Index nzmax, // nb values - Index m, // number of row - Index n, // number of columns - Index * a_p, // column pointers (size n+1) or col indices (size nzmax) - Index * a_i, // row indices, size nzmax - Block * a_x // numerical values, size nzmax - ) - { - // check ap, size m beecause ther is at least the diagonal value wich is different of 0 - if (a_p[0]!=0) - { - msg_error("CompressedRowSparseMatrixGeneric") << "First value of row indices (a_p) should be 0"; - return false; - } - - for (Index i=1; i<=m; i++) - { - if (a_p[i]<=a_p[i-1]) - { - msg_error("CompressedRowSparseMatrixGeneric") << "Row (a_p) indices are not sorted index " << i-1 << " : " << a_p[i-1] << " , " << i << " : " << a_p[i]; - return false; - } - } - if (nzmax == -1) - { - nzmax = a_p[m]; - } - else if (a_p[m]!=nzmax) - { - msg_error("CompressedRowSparseMatrixGeneric") << "Last value of row indices (a_p) should be " << nzmax << " and is " << a_p[m]; - return false; - } - - - Index k=1; - for (Index i=0; i=n) - { - msg_error("CompressedRowSparseMatrixGeneric") << "Column (a_i) indices are not correct " << i << " : " << a_i[i]; - return false; - } - } - k++; - } - - for (Index i=0; i Date: Fri, 28 Aug 2026 13:29:59 +0900 Subject: [PATCH 10/12] [LinearAlgebra] Remove the write-only touchedBlock member touchedBlock was declared alongside the CSR arrays and cleared once in compressBtemp(), but never written to and never read -- by this class or by anything else in the tree. Its documented purpose, tracking which blocks were touched since the last compression, was never implemented. It was also the one data member swap() did not exchange, so swapping two matrices left it behind; removing it makes that discrepancy moot. The VecFlag alias it used is kept: it is now unused in-tree, but it is public API on both CRSBlockTraits and the matrix classes, so out-of-tree code may still name it. Co-Authored-By: Claude Opus 5 --- .../src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h | 2 -- 1 file changed, 2 deletions(-) diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h index c937fd5bd6f..194856b01e0 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h @@ -221,7 +221,6 @@ public : VecIndex rowBegin; ///< column indices of non-empty blocks in each row. The column indices of the non-empty block within the i-th non-empty row are all the colsIndex[j], j in [rowBegin[i],rowBegin[i+1]) VecIndex colsIndex; ///< column indices of all the non-empty blocks, sorted by increasing row index and column index VecBlock colsValue; ///< values of the non-empty blocks, in the same order as in colsIndex - VecFlag touchedBlock; ///< boolean vector, i-th value is true if block has been touched since last compression. /// Additional storage to make block insertion more efficient VecIndexedBlock btemp; ///< unsorted blocks and their indices @@ -581,7 +580,6 @@ public : rowBegin.clear(); colsIndex.clear(); colsValue.clear(); - touchedBlock.clear(); rowIndex.reserve(oldRowIndex.size()); rowBegin.reserve(oldRowIndex.size() + 1); From 073c5f8fc6c4f7173bf7804974dff723c4f7c8ed Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 28 Aug 2026 13:33:11 +0900 Subject: [PATCH 11/12] [LinearAlgebra] Remove the unused readVector()/writeVector() helpers Neither helper has ever been called: write()/read() serialise the CSR arrays through the stream operators instead, and nothing outside the class can reach these two since they are protected. readVector() would also have been wrong for most instantiations. It parses every entry through safeStrToInt into an int before pushing it back, so it cannot round-trip a vector of blocks or of any real type -- only the index vectors. Dropping them leaves sofa/type/hardening.h unused in this header, so remove that include too. Verified that no transitive consumer depended on it by rebuilding Sofa.Core, Sofa.Component.Constraint.Projective, Sofa.Component.LinearSolver.Direct, Sofa.Component.LinearSystem, Sofa.Component.StateContainer, Sofa.Component.SolidMechanics.FEM.Elastic and Sofa.Component.Mapping.Linear. Co-Authored-By: Claude Opus 5 --- .../CompressedRowSparseMatrixGeneric.h | 27 ------------------- 1 file changed, 27 deletions(-) diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h index 194856b01e0..16677756060 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h @@ -29,7 +29,6 @@ #include #include -#include #include #include #include @@ -1509,32 +1508,6 @@ public : return is; } -protected: - - template - void writeVector(const TVec& vec, std::ostream& os) - { - for (auto& v : vec) - os < - void readVector(TVec& vec, std::istream& in) - { - std::string temp; - while (std::getline(in, temp, ';')) - { - int val{}; - if(sofa::type::hardening::safeStrToInt(temp, val)) - { - vec.push_back(val); - } - else - { - msg_warning("CompressedRowSparseMatrixGeneric") << "could not parse " << temp << " ; skipping entry."; - } - } - } }; #if !defined(SOFA_COMPONENT_LINEARSOLVER_COMPRESSEDROWSPARSEMATRIXGENERIC_CPP) From d90377fed3f772a1993de51ede9e5cde335c8257 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 28 Aug 2026 13:35:16 +0900 Subject: [PATCH 12/12] [LinearAlgebra] Remove the uncallable RowConstIterator::operator* RowConstIterator declared template Real operator*(const VecDeriv& v) const but Real appears only in the return type, so it can be neither deduced from the argument nor defaulted, and operator syntax offers no way to supply it explicitly. No call to it can compile: overload resolution does not even consider the member, and a *rowIt * v expression fails with "indirection requires pointer operand". The functionality it duplicated is available, and actually used, through the free CompressedRowSparseMatrixVecDerivMult(row, vec), which defaults Real to VecDeriv::value_type::Real. Removed rather than repaired: reinstating it would mean giving Real that same default, which turns a member nothing has ever been able to call into new API. That is a deliberate addition, not dead-code cleanup, so it is left out of this series. Co-Authored-By: Claude Opus 5 --- .../linearalgebra/CompressedRowSparseMatrixConstraint.h | 6 ------ 1 file changed, 6 deletions(-) diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixConstraint.h b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixConstraint.h index 0187e3fe7ad..4c649a0de9e 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixConstraint.h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixConstraint.h @@ -422,12 +422,6 @@ class CompressedRowSparseMatrixConstraint : public sofa::linearalgebra::Compress return !(*this < other); } - template - Real operator*(const VecDeriv& v) const - { - return CompressedRowSparseMatrixVecDerivMult(row(), v); - } - private: /// a default-constructed iterator must report itself as invalid rather