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 diff --git a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h index 066e72b8604..2f76c20ece0 100644 --- a/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h +++ b/Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h @@ -422,12 +422,97 @@ 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; } 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 +840,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]; } @@ -807,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); @@ -825,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]; } } @@ -834,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]; } @@ -892,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) { @@ -901,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) { @@ -981,13 +1036,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); @@ -1013,15 +1062,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) { @@ -1076,34 +1118,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 ("<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) @@ -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); } } } @@ -581,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 } @@ -745,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); } @@ -763,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); } @@ -779,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); } @@ -838,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]; } @@ -851,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]; } @@ -864,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]; 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); +}