Repository: PowerGridModel/power-grid-model (language: cpp)

Issue:
[FEATURE] Add prefactorization for EIGEN superLU solver



Discussion / hints:

Fix patch:
diff --git a/include/power_grid_model/math_solver/eigen_superlu_solver.hpp b/include/power_grid_model/math_solver/eigen_superlu_solver.hpp
index 3f43c5804..25ede6191 100644
--- a/include/power_grid_model/math_solver/eigen_superlu_solver.hpp
+++ b/include/power_grid_model/math_solver/eigen_superlu_solver.hpp
@@ -30,96 +30,131 @@ class EigenSuperLUSolver final {
     using BufferVector = Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, 1>, Eigen::Aligned8>;
     using SparseIdxEntry = Eigen::Triplet<Idx, Idx>;
 
+    struct BSRHandle {
+        Idx matrix_size_in_block;
+        Idx block_size;
+        Idx nnz_block;
+        Idx block_capacity;
+        Idx nnz;
+        Idx matrix_size;
+        std::shared_ptr<IdxVector const> data_mapping;
+    };
+
    public:
     EigenSuperLUSolver(Idx matrix_size_in_block,                    // size of matrix in block form
                        Idx block_size,                              // size of block
                        std::shared_ptr<IdxVector const> const& ia,  // indptr
                        std::shared_ptr<IdxVector const> const& ja   // column indices
-                       )
-        : matrix_size_in_block_{matrix_size_in_block},
-          block_size_{block_size},
-          nnz_block_{ia->back()},
-          block_capacity_{block_size_ * block_size_},
-          nnz_{nnz_block_ * block_capacity_},
-          matrix_size_{matrix_size_in_block_ * block_size_} {
-        initialize_solver(*ia, *ja);
+    ) {
+        bsr_handle_.matrix_size_in_block = matrix_size_in_block;
+        bsr_handle_.block_size = block_size;
+        bsr_handle_.nnz_block = ia->back();
+        bsr_handle_.block_capacity = block_size * block_size;
+        bsr_handle_.nnz = bsr_handle_.nnz_block * bsr_handle_.block_capacity;
+        bsr_handle_.matrix_size = matrix_size_in_block * block_size;
+        initialize_matrix(*ia, *ja);
+        initialize_solver();
+    }
+
+    // copy construction
+    EigenSuperLUSolver(EigenSuperLUSolver const& other)
+        : bsr_handle_{other.bsr_handle_}, sparse_matrix_{other.sparse_matrix_} {
+        initialize_solver();
+    }
+
+    // copy assignment
+    EigenSuperLUSolver& operator=(EigenSuperLUSolver const& other) {
+        // copy new handle
+        bsr_handle_ = other.bsr_handle_;
+        // copy new matrix
+        sparse_matrix_ = other.sparse_matrix_;
+        // re initialize
+        initialize_solver();
+        prefactorized_ = false;
+        return *this;
     }
 
     void solve(void const* data, void* b, void* x, bool use_prefactorization = false) {
-        (void)use_prefactorization;  // suppress unused variable warning
-        // copy data
-        copy_matrix_data(data);
-        BufferVector bm{std::launder(reinterpret_cast<T*>(b)), matrix_size_};
-        BufferVector xm{std::launder(reinterpret_cast<T*>(x)), matrix_size_};
-        SparseSolver sparse_solver;
-        sparse_solver.analyzePattern(sparse_matrix_);
-        sparse_solver.factorize(sparse_matrix_);
-        if (sparse_solver.info() != Eigen::Success) {
-            throw SparseMatrixError{sparse_solver.info(), sparse_solver.lastErrorMessage()};
+        // reset possible pre-factorization if we are not using prefactorization
+        prefactorized_ = prefactorized_ && use_prefactorization;
+        // run factorization
+        if (!prefactorized_) {
+            prefactorize(data);
         }
-        xm = sparse_solver.solve(bm);
-        if (sparse_solver.info() != Eigen::Success) {
-            throw SparseMatrixError{sparse_solver.info(), sparse_solver.lastErrorMessage()};
+        // run solve
+        BufferVector bm{std::launder(reinterpret_cast<T*>(b)), bsr_handle_.matrix_size};
+        BufferVector xm{std::launder(reinterpret_cast<T*>(x)), bsr_handle_.matrix_size};
+        xm = sparse_solver_.solve(bm);
+        if (sparse_solver_.info() != Eigen::Success) {
+            throw SparseMatrixError{sparse_solver_.info(), sparse_solver_.lastErrorMessage()};
         }
     }
 
-    // empty prefactorize function
     void prefactorize(void const* data) {
-        (void)data;
+        // copy data
+        copy_matrix_data(data);
+        // factorize
+        sparse_solver_.factorize(sparse_matrix_);
+        if (sparse_solver_.info() != Eigen::Success) {
+            throw SparseMatrixError{sparse_solver_.info(), sparse_solver_.lastErrorMessage()};
+        }
+        prefactorized_ = true;
     }
 
     void invalidate_prefactorization() {
+        prefactorized_ = false;
     }
 
    private:
-    Idx matrix_size_in_block_;
-    Idx block_size_;
-    Idx nnz_block_;
-    Idx block_capacity_;
-    Idx nnz_;
-    Idx matrix_size_;
+    BSRHandle bsr_handle_;
     SparseMatrix sparse_matrix_;
-    // SparseSolver sparse_solver_;  // not copyable or movable
-    std::shared_ptr<IdxVector const> data_mapping_;
+    SparseSolver sparse_solver_;  // not copyable or movable
+    bool prefactorized_{false};
 
-    void initialize_solver(IdxVector const& ia,  // indptr
+    void initialize_matrix(IdxVector const& ia,  // indptr
                            IdxVector const& ja   // column indices
     ) {
         // loop to add triplet
         std::vector<SparseIdxEntry> idx_list;
-        idx_list.reserve(nnz_);
+        idx_list.reserve(bsr_handle_.nnz);
         // two loops for blocks in csr
-        for (Idx bi = 0; bi != matrix_size_in_block_; ++bi) {
+        for (Idx bi = 0; bi != bsr_handle_.matrix_size_in_block; ++bi) {
             for (Idx block_ind = ia[bi]; block_ind != ia[bi + 1]; ++block_ind) {
                 Idx const bj = ja[block_ind];
                 // two loops for rows/columns inside block
-                for (Idx ci = 0; ci != block_size_; ++ci) {
-                    for (Idx cj = 0; cj != block_size_; ++cj) {
-                        Idx const i = bi * block_size_ + ci;
-                        Idx const j = bj * block_size_ + cj;
-                        Idx const data_ind = block_ind * block_capacity_ + ci * block_size_ + cj;
+                for (Idx ci = 0; ci != bsr_handle_.block_size; ++ci) {
+                    for (Idx cj = 0; cj != bsr_handle_.block_size; ++cj) {
+                        Idx const i = bi * bsr_handle_.block_size + ci;
+                        Idx const j = bj * bsr_handle_.block_size + cj;
+                        Idx const data_ind = block_ind * bsr_handle_.block_capacity + ci * bsr_handle_.block_size + cj;
                         idx_list.emplace_back(i, j, data_ind);
                     }
                 }
             }
         }
         // idx matrix
-        SparseIndexMatrix idx_matrix{matrix_size_, matrix_size_};
+        SparseIndexMatrix idx_matrix{bsr_handle_.matrix_size, bsr_handle_.matrix_size};
         idx_matrix.setFromTriplets(idx_list.begin(), idx_list.end());
         assert(idx_matrix.isCompressed());
         // copy idx data
-        data_mapping_ = std::make_shared<IdxVector const>(idx_matrix.valuePtr(), idx_matrix.valuePtr() + nnz_);
+        bsr_handle_.data_mapping =
+            std::make_shared<IdxVector const>(idx_matrix.valuePtr(), idx_matrix.valuePtr() + bsr_handle_.nnz);
         // initialize sparse matrix
         sparse_matrix_ = idx_matrix.cast<T>();
         assert(sparse_matrix_.isCompressed());
     }
 
+    void initialize_solver() {
+        sparse_solver_.analyzePattern(sparse_matrix_);
+    }
+
     void copy_matrix_data(void const* data) {
         T const* const input_ptr = std::launder(reinterpret_cast<T const*>(data));
         T* const data_ptr = sparse_matrix_.valuePtr();
-        std::transform(data_mapping_->begin(), data_mapping_->end(), data_ptr, [&input_ptr](Idx j) {
-            return input_ptr[j];
-        });
+        std::transform(bsr_handle_.data_mapping->begin(), bsr_handle_.data_mapping->end(), data_ptr,
+                       [&input_ptr](Idx j) {
+                           return input_ptr[j];
+                       });
     }
 };
 


Test patch:
diff --git a/tests/cpp_unit_tests/test_bsr_solver.cpp b/tests/cpp_unit_tests/test_bsr_solver.cpp
index f11e6d766..285d9af80 100644
--- a/tests/cpp_unit_tests/test_bsr_solver.cpp
+++ b/tests/cpp_unit_tests/test_bsr_solver.cpp
@@ -104,11 +104,6 @@ TEST_CASE("Test BSR solver") {
         solver.solve(data.data(), rhs.data(), x_solver.data(), true);
         check_result(x, x_solver);
 
-#if defined(POWER_GRID_MODEL_USE_MKL_AT_RUNTIME)
-        if (!get_pardiso_handle().has_pardiso) {
-            return;
-        }
-#endif  // POWER_GRID_MODEL_USE_MKL_AT_RUNTIME
         // basically data * 2
         std::vector<double> other_data = {
             2.0, 0.0, 0.0, 4.0,  // 0, 0
@@ -122,12 +117,10 @@ TEST_CASE("Test BSR solver") {
         // basically x / 2
         std::vector<double> other_x = {0.5, 0.5, 1.0, 1.0, 0.5, 0.5, 1.0, 1.0};
 
-#if defined(POWER_GRID_MODEL_USE_MKL)
         // because all data should be prefactorized, changing the data should not
         // change the result when use_prefactorization = true
         solver.solve(other_data.data(), rhs.data(), x_solver.data(), true);
         check_result(x, x_solver);
-#endif  // POWER_GRID_MODEL_USE_MKL
 
         // prefactorize other_data, then solve and compare with other_x
         solver.prefactorize(other_data.data());
@@ -137,6 +130,13 @@ TEST_CASE("Test BSR solver") {
         // solve and compare with other_x without using prefactorization
         solver.solve(other_data.data(), rhs.data(), x_solver.data(), false);
         check_result(other_x, x_solver);
+
+        // invalidate pre-factorization
+        // and re-run with original data with pre-factorization enabled
+        // it should still re-do the factorization
+        solver.invalidate_prefactorization();
+        solver.solve(data.data(), rhs.data(), x_solver.data(), true);
+        check_result(x, x_solver);
     }
 }