Skip to content
Snippets Groups Projects
Commit a807c110 authored by oliver.sander_at_tu-dresden.de's avatar oliver.sander_at_tu-dresden.de
Browse files

CholmodSolver: Implement rowAdd method

This method allows to add new matrix rows to the Cholesky factor.
See the CHOLMOD manual for more details.
parent 89f27925
No related branches found
No related tags found
No related merge requests found
Pipeline #41706 failed
......@@ -268,6 +268,63 @@ public:
errorCode_ = cholmodCommonObject().status;
}
/** \brief Adds row and column to an LDL^T factorization
*
* \param k The scalar(!) row and column to add
*
* The pattern and values are taken from the original matrix given to CholmodSolver
*
* \todo Currently only works if the matrix entries are scalars or 1x1 FieldMatrix
*/
void rowAdd(size_t k)
{
// Length of a row
auto n = istlCholmodSolver_.cholmodFactor()->n;
size_t nonZeros = 0;
for ([[maybe_unused]] auto&& entry : sparseRange((*this->matrix_)[k]))
nonZeros++;
// Create a Cholmod data structure with the row to be added
const auto deleter = [c = &cholmodCommonObject()](auto* p) {
cholmod_free_sparse(&p, c);
};
auto R = std::unique_ptr<cholmod_sparse, decltype(deleter)>(
cholmod_allocate_sparse(n, // # rows
1, // # cols
nonZeros, // # of nonzeroes
1, // indices are sorted ( 1 = true)
1, // matrix is "packed" ( 1 = true)
0, // stype of matrix ( 0 = consider the entire matrix )
CHOLMOD_REAL, // xtype of matrix ( CHOLMOD_REAL = single array, no complex numbers)
&cholmodCommonObject()
), deleter);
// copy the data of BCRS matrix to Cholmod Sparse matrix
int* Ap = static_cast<int*>(R->p);
int* Ai = static_cast<int*>(R->i);
double* Ax = static_cast<double*>(R->x);
Ap[0] = 0; // We have a single column, which starts at 0
Ap[1] = nonZeros;
for (auto [entry, col] : sparseRange((*this->matrix_)[k]))
{
*Ax++ = entry;
*Ai++ = invPermutation_[col];
}
cholmod_factor *L = istlCholmodSolver_.cholmodFactor();
cholmod_rowadd(invPermutation_[k], // row/column index to add
R.get(), // row/column of matrix to factorize (n x 1)
L, // Factor to modify
&cholmodCommonObject());
// get the error code from Cholmod in case something happened
errorCode_ = cholmodCommonObject().status;
}
//! dune-istl Cholmod solver object
Dune::Cholmod<VectorType> istlCholmodSolver_;
......
......@@ -144,6 +144,19 @@ bool checkMultipleSolves()
passed = false;
}
// Put the deleted row/column back in
cholmodSolver.rowAdd(1);
reference = {1,2,1,2};
cholmodSolver.solve();
diff = reference;
diff -= x;
if (std::isnan(diff.infinity_norm()) || diff.infinity_norm() > 1e-8)
{
std::cerr << "CholmodSolver didn't produce the correct result!" << std::endl;
passed = false;
}
return passed;
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment