Skip to content

[LinearAlgebra] CompressedRowSparseMatrix family: Fix problematic cases #6270

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
fredroy wants to merge 8 commits into
base: master
Choose a base branch
from
Open

[LinearAlgebra] CompressedRowSparseMatrix family: Fix problematic cases #6270

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
5cbf6de
add unit tests on problematic cases
fredroy Aug 28, 2026
a8573a5
[LinearAlgebra] Fix block() returning a neighbouring row's value
fredroy Aug 28, 2026
9cd39c4
[LinearAlgebra] Fix clearRowBlock() segfaulting on an empty matrix
fredroy Aug 28, 2026
613ca8f
[LinearAlgebra] Fix clearRowColBlock() deleting an unrelated row
fredroy Aug 28, 2026
ddd7314
[LinearAlgebra] Fix clearRowCol() indexing rowBegin with a column number
fredroy Aug 28, 2026
31497dd
[LinearAlgebra] Make filterValues() honour keepEmptyRows
fredroy Aug 28, 2026
7763f50
[LinearAlgebra] Initialise RowConstIterator's members
fredroy Aug 28, 2026
96ca07a
[LinearAlgebra] Route the remaining row/column lookups through the gu…
fredroy Aug 28, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 4 additions & 2 deletions Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixConstraint.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
198 changes: 114 additions & 84 deletions Sofa/framework/LinearAlgebra/src/sofa/linearalgebra/CompressedRowSparseMatrixGeneric.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -422,12 +422,97 @@ public :
Index maxColIndex = 0;
for (Index rowId = 0; rowId < static_cast<Index>(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
Expand Down Expand Up @@ -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];
}
Expand Down Expand Up @@ -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);
Expand All @@ -825,43 +896,29 @@ 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];
}
}

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];
}
Expand Down Expand Up @@ -892,17 +949,15 @@ public :
bool rowFound = true;
if (rowId < 0 || rowId >= static_cast<Index>(rowIndex.size()) || rowIndex[rowId] != i)
{
rowId = Index(i * rowIndex.size() / nBlockRow);
rowFound = sortedFind(rowIndex, i, rowId);
rowFound = searchRow(i, rowId);
}
if (rowFound)
{
bool colFound = true;
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)
{
Expand Down Expand Up @@ -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);

Expand All @@ -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)
{
Expand Down Expand Up @@ -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 "<<i<<" in "<< this->Name() << " of size ("<<rowBSize()<<","<<colBSize()<<")";
return;
}

if (!foundRowId && !foundColId)
/// An absent row is not an error: the matrix is sparse, so row i may hold
/// no block while column i still does in other rows. Only delete the row
/// when it is actually registered, otherwise rowId is meaningless and
/// deleteRow() would drop an unrelated row.
Index rowId = 0;
if (findRow(i, rowId))
{
msg_error("CompressedRowSparseMatrixGeneric") << "invalid write access to row and column "<<i<<" in "<< this->Name() << " of size ("<<rowBSize()<<","<<colBSize()<<")";
return;
deleteRow(rowId); /// Do not call clearRowBlock to only compress zero if activated once.
}

deleteRow(rowId); /// Do not call clearRow to only compress zero if activated once.
clearColBlock(i);
}

Expand Down
Loading
Loading