From 5cbf6de4a4e08f13d78419afb5d6c056250b5139 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 28 Aug 2026 11:30:10 +0900 Subject: [PATCH 1/8] 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 a8573a563186c3d7a49360c054149e8e408f6d61 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 28 Aug 2026 12:20:26 +0900 Subject: [PATCH 2/8] [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 here, in two variants: - searchRow() / searchColInRange() are the interpolated binary search on its own, with the guards the copies were missing -- the nBlockRow == 0 and nBlockCol == 0 division guards, and sortedFind()'s own "empty range means not found" behaviour; - findRow() / findColInRange() add the first/last-row and first/last-column checks on top, and are the ones block() uses. Keeping the two apart matters for performance: those end checks pay for themselves when the queried entry is usually an end one, as it is here, but on the insertion paths they almost never hit and cost two extra loads and branches per call. Later commits pick whichever variant matches the call site. Fixes CompressedRowSparseMatrixGeneric.BlockOnEmptyRowAfterFullRows, CompressedRowSparseMatrixGeneric.BlockOnEmptyRowDoesNotBorrowNeighbourValue and CompressedRowSparseMatrixMechanical.ElementAfterFullRows. --- .../CompressedRowSparseMatrixGeneric.h | 99 ++++++++++++++++--- 1 file changed, 84 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..3c4ba674475 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h @@ -428,6 +428,88 @@ 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; + } + return searchColInRange(rowRange, j, colId); + } + + /** + * \brief Binary search for column j inside a row, from an interpolated guess. + * + * findColInRange() without the first/last-column checks. Those two loads pay + * for themselves when the queried column is usually an end one, and cost more + * than they save on the insertion paths, where they almost never hit. + * sortedFind() already reports an empty range as "not found", so this is safe + * on a row that holds no block. + **/ + bool searchColInRange(const Range& rowRange, Index j, Index& colId) const + { + 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; + } + return searchRow(i, rowId); + } + + /** + * \brief Binary search for row i, from an interpolated guess. + * + * findRow() without the first/last-row checks, for the same reason as + * searchColInRange(). sortedFind() reports an empty rowIndex as "not found", + * and the division is guarded, so this is safe on an empty matrix. + **/ + bool searchRow(Index i, Index& rowId) const + { + 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 +837,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 9cd39c4040d4b83bda8cfe5454da41f81858481b Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 28 Aug 2026 12:21:36 +0900 Subject: [PATCH 3/8] [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. --- .../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 3c4ba674475..0222cf33fcc 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h @@ -1050,13 +1050,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 613ca8f459594fbd159ec76076f13b087918e669 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 28 Aug 2026 12:22:58 +0900 Subject: [PATCH 4/8] [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. --- .../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 0222cf33fcc..81f36b9500c 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h @@ -1076,15 +1076,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) { @@ -1139,34 +1132,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 5/8] [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 searchRow() before taking its range, locate the symmetric block with searchColInRange(), 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. The lean lookup variants are the right ones here: this code walks a row's entries and the queried row is not usually the first or last registered one. Correctness still costs something -- finding the symmetric block is now a real row lookup per entry, where the buggy version used a wrong index and found nothing. Measured at about 75 ns per clearRowCol() call on a 5000-block-row matrix, and it is called once per constrained scalar DOF per step. Fixes CompressedRowSparseMatrixMechanical.ClearRowColMissesSymmetricBlockWhenRowsAreSparse and CompressedRowSparseMatrixMechanical.ClearRowColAsymmetricPatternMat3. --- .../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..75722da7b5e 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->searchRow(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->searchRow(j, jRowId)) continue; + + Index colId = 0; + if (!this->searchColInRange(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 31497dd28f883a77a5c67530718210ca305c50fb Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 28 Aug 2026 12:24:31 +0900 Subject: [PATCH 6/8] [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. --- .../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 75722da7b5e..04a5ed3abca 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 7763f507c97fb2c2a3efe12a28486aa539056e12 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 28 Aug 2026 12:25:41 +0900 Subject: [PATCH 7/8] [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. --- .../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 96ca07ab31bfb04c4fbc42aa0bbd80edd2a69fc8 Mon Sep 17 00:00:00 2001 From: Frederick Roy Date: Fri, 28 Aug 2026 12:27:46 +0900 Subject: [PATCH 8/8] [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]. All of these use the lean searchRow() / searchColInRange() variants, matching what they did before the series: none of them had the first/last-entry checks, and adding them here would slow assembly down measurably. wblock() in particular was already safe -- it guarded both divisions, and sortedFind() reports an empty range as not found -- so it gains nothing from those checks and would only pay for them. Measured on a 5000-block-row FEM-shaped matrix, giving wblock() the end checks cost 28% on add(); with the lean variants it is within noise of 503ce3a9, while element() misses improve 19% and block() 9%. No behaviour change on well-formed input; this is the undefined behaviour left over once the six defects with failing tests were fixed. --- .../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 81f36b9500c..2f76c20ece0 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; } @@ -876,11 +879,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); @@ -894,8 +896,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]; } } @@ -903,34 +905,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 (searchRow(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 (searchColInRange(Range(rowBegin[rowId], rowBegin[rowId+1]), j, colId)) { return &colsValue[colId]; } @@ -961,8 +949,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 = searchRow(i, rowId); } if (rowFound) { @@ -970,8 +957,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 = searchColInRange(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 04a5ed3abca..d907177f4e8 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->searchRow(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->searchRow(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->searchColInRange(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->searchRow(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->searchColInRange(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->searchRow(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->searchColInRange(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->searchRow(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->searchRow(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->searchRow(ib, rowId)) { index = this->rowBegin[rowId]; index2 = this->rowBegin[rowId+1];