Fix incomplete cholesky.

This commit is contained in:
Antonio Sánchez
2024-03-08 19:18:10 +00:00
committed by Rasmus Munk Larsen
parent f1adb0ccc2
commit 352ede96e4
3 changed files with 67 additions and 19 deletions

View File

@@ -217,15 +217,18 @@ class SparseMatrix : public SparseCompressedBase<SparseMatrix<Scalar_, Options_,
return m_data.atInRange(m_outerIndex[outer], end, inner);
}
/** \returns a non-const reference to the value of the matrix at position \a i, \a j
/** \returns a non-const reference to the value of the matrix at position \a i, \a j.
*
* If the element does not exist then it is inserted via the insert(Index,Index) function
* which itself turns the matrix into a non compressed form if that was not the case.
* The output parameter `inserted` is set to true.
*
* Otherwise, if the element does exist, `inserted` will be set to false.
*
* This is a O(log(nnz_j)) operation (binary search) plus the cost of insert(Index,Index)
* function if the element does not already exist.
*/
inline Scalar& coeffRef(Index row, Index col) {
inline Scalar& findOrInsertCoeff(Index row, Index col, bool* inserted) {
eigen_assert(row >= 0 && row < rows() && col >= 0 && col < cols());
const Index outer = IsRowMajor ? row : col;
const Index inner = IsRowMajor ? col : row;
@@ -240,17 +243,37 @@ class SparseMatrix : public SparseCompressedBase<SparseMatrix<Scalar_, Options_,
m_innerNonZeros[outer]++;
m_data.index(end) = StorageIndex(inner);
m_data.value(end) = Scalar(0);
if (inserted != nullptr) {
*inserted = true;
}
return m_data.value(end);
}
}
if ((dst < end) && (m_data.index(dst) == inner))
if ((dst < end) && (m_data.index(dst) == inner)) {
// this coefficient exists, return a refernece to it
if (inserted != nullptr) {
*inserted = false;
}
return m_data.value(dst);
else
} else {
if (inserted != nullptr) {
*inserted = true;
}
// insertion will require reconfiguring the buffer
return insertAtByOuterInner(outer, inner, dst);
}
}
/** \returns a non-const reference to the value of the matrix at position \a i, \a j
*
* If the element does not exist then it is inserted via the insert(Index,Index) function
* which itself turns the matrix into a non compressed form if that was not the case.
*
* This is a O(log(nnz_j)) operation (binary search) plus the cost of insert(Index,Index)
* function if the element does not already exist.
*/
inline Scalar& coeffRef(Index row, Index col) { return findOrInsertCoeff(row, col, nullptr); }
/** \returns a reference to a novel non zero coefficient with coordinates \a row x \a col.
* The non zero coefficient must \b not already exist.
*