From a807c110eb9939e351306fbf0024f727f00c16e6 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Mon, 18 Oct 2021 21:58:20 +0200
Subject: [PATCH] CholmodSolver: Implement rowAdd method
This method allows to add new matrix rows to the Cholesky factor.
See the CHOLMOD manual for more details.
---
dune/solvers/solvers/cholmodsolver.hh | 57 ++++++++++++++++++++++++++
dune/solvers/test/cholmodsolvertest.cc | 13 ++++++
2 files changed, 70 insertions(+)
diff --git a/dune/solvers/solvers/cholmodsolver.hh b/dune/solvers/solvers/cholmodsolver.hh
index 99d9bd72..87b65225 100644
--- a/dune/solvers/solvers/cholmodsolver.hh
+++ b/dune/solvers/solvers/cholmodsolver.hh
@@ -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_;
diff --git a/dune/solvers/test/cholmodsolvertest.cc b/dune/solvers/test/cholmodsolvertest.cc
index faf1cb94..a0d7733c 100644
--- a/dune/solvers/test/cholmodsolvertest.cc
+++ b/dune/solvers/test/cholmodsolvertest.cc
@@ -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;
}
--
GitLab