revise based on comments

Signed-off-by: Tony Xiang <[email protected]>

Author: TonyXiang8787

include/power_grid_model/math_solver/iterative_current_pf_solver.hpp

@@ -171,7 +171,7 @@ class IterativeCurrentPFSolver : public IterativePFSolver<sym, IterativeCurrentP
     // Solve the linear equations I_inj = YU
     // inplace
     void solve_matrix() {
-        sparse_solver_.solve(*mat_data_, *perm_, rhs_u_, rhs_u_);
+        sparse_solver_.solve_with_prefactorization(*mat_data_, *perm_, rhs_u_, rhs_u_);
     }
 
     // Find maximum deviation in voltage among all buses

include/power_grid_model/math_solver/iterative_linear_se_solver.hpp

@@ -621,7 +621,7 @@ class IterativeLinearSESolver {
             prepare_rhs(y_bus, measured_values, output.u);
             // solve with prefactorization
             sub_timer = Timer(calculation_info, 2225, "Solve sparse linear equation (pre-factorized)");
-            sparse_solver_.solve((std::vector<SEGainBlock<sym>> const&)data_gain_, perm_, x_rhs_, x_rhs_);
+            sparse_solver_.solve_with_prefactorization(data_gain_, perm_, x_rhs_, x_rhs_);
             sub_timer = Timer(calculation_info, 2226, "Iterate unknown");
             max_dev = iterate_unknown(output.u, measured_values.has_angle_measurement());
         }

include/power_grid_model/math_solver/linear_pf_solver.hpp

@@ -105,7 +105,7 @@ class LinearPFSolver {
         // solve
         // u vector will have I_injection for slack bus for now
         sub_timer = Timer(calculation_info, 2222, "Solve sparse linear equation");
-        sparse_solver_.solve(mat_data_, perm_, output.u, output.u);
+        sparse_solver_.prefactorize_and_solve(mat_data_, perm_, output.u, output.u);
 
         // calculate math result
         sub_timer = Timer(calculation_info, 2223, "Calculate Math Result");

include/power_grid_model/math_solver/newton_raphson_pf_solver.hpp

@@ -363,7 +363,7 @@ class NewtonRaphsonPFSolver : public IterativePFSolver<sym, NewtonRaphsonPFSolve
 
     // Solve the linear Equations
     void solve_matrix() {
-        sparse_solver_.solve(data_jac_, perm_, del_x_pq_, del_x_pq_);
+        sparse_solver_.prefactorize_and_solve(data_jac_, perm_, del_x_pq_, del_x_pq_);
     }
 
     // Get maximum deviation among all bus voltages

include/power_grid_model/math_solver/sparse_lu_solver.hpp

@@ -64,7 +64,7 @@ struct sparse_lu_entry_trait<Tensor, RHSVector, XVector, enable_tensor_lu_t<Tens
     struct BlockPerm {
         typename LUFactor::PermutationPType p;
         typename LUFactor::PermutationQType q;
-    };  //Extract permutation matrices p and q from LUFactor
+    };  // Extract permutation matrices p and q from LUFactor
     using BlockPermArray = std::vector<BlockPerm>;
 };
 
@@ -90,18 +90,19 @@ class SparseLUSolver {
     }
 
     // solve with new matrix data, need to factorize first
-    void solve(std::vector<Tensor>& data,         // matrix data, factorize in-place
-               BlockPermArray& block_perm_array,  // pre-allocated permutation array, will be overwritten
-               std::vector<RHSVector> const& rhs, std::vector<XVector>& x) {
+    void prefactorize_and_solve(
+        std::vector<Tensor>& data,         // matrix data, factorize in-place
+        BlockPermArray& block_perm_array,  // pre-allocated permutation array, will be overwritten
+        std::vector<RHSVector> const& rhs, std::vector<XVector>& x) {
         prefactorize(data, block_perm_array);
         // call solve with const method
-        solve((std::vector<Tensor> const&)data, block_perm_array, rhs, x);
+        solve_with_prefactorization((std::vector<Tensor> const&)data, block_perm_array, rhs, x);
     }
 
     // solve with existing pre-factorization
-    void solve(std::vector<Tensor> const& data,         // pre-factoirzed data, const ref
-               BlockPermArray const& block_perm_array,  // pre-calculated permutation, const ref
-               std::vector<RHSVector> const& rhs, std::vector<XVector>& x) {
+    void solve_with_prefactorization(std::vector<Tensor> const& data,         // pre-factoirzed data, const ref
+                                     BlockPermArray const& block_perm_array,  // pre-calculated permutation, const ref
+                                     std::vector<RHSVector> const& rhs, std::vector<XVector>& x) {
         // local reference
         auto const& row_indptr = *row_indptr_;
         auto const& col_indices = *col_indices_;
@@ -333,7 +334,7 @@ class SparseLUSolver {
 
    private:
     Idx size_;
-    Idx nnz_;
+    Idx nnz_;  // number of non zeroes (in block)
     std::shared_ptr<IdxVector const> row_indptr_;
     std::shared_ptr<IdxVector const> col_indices_;
     std::shared_ptr<IdxVector const> diag_lu_;

tests/cpp_unit_tests/test_sparse_lu_solver.cpp

@@ -46,7 +46,7 @@ TEST_CASE("Test Sparse LU solver") {
     auto col_indices = std::make_shared<IdxVector const>(IdxVector{0, 1, 2, 0, 1, 2, 0, 1, 2});
     auto diag_lu = std::make_shared<IdxVector const>(IdxVector{0, 4, 8});
 
-    SECTION("Test scalar(double) calculation") {
+    SECTION("Scalar(double) calculation") {
         // [4 1 5        3          21
         //  3 7 f     * [-1]   =  [ 2 ]
         //  2 f 6]       2          18
@@ -62,23 +62,31 @@ TEST_CASE("Test Sparse LU solver") {
         SparseLUSolver<double, double, double>::BlockPermArray block_perm{};
 
         SECTION("Test calculation") {
-            solver.solve(data, block_perm, rhs, x);
+            solver.prefactorize_and_solve(data, block_perm, rhs, x);
             check_result(x, x_ref);
         }
 
         SECTION("Test (pseudo) singular") {
             data[0] = 0.0;
-            CHECK_THROWS_AS(solver.solve(data, block_perm, rhs, x), SparseMatrixError);
+            CHECK_THROWS_AS(solver.prefactorize_and_solve(data, block_perm, rhs, x), SparseMatrixError);
         }
 
         SECTION("Test prefactorize") {
             solver.prefactorize(data, block_perm);
-            solver.solve((std::vector<double> const&)data, block_perm, rhs, x);
+            solver.solve_with_prefactorization((std::vector<double> const&)data, block_perm, rhs, x);
             check_result(x, x_ref);
         }
+
+        SECTION("Data is prefactorized by solve") {
+            auto prefactorized_data = data;
+            auto prefactorized_block_perm = block_perm;
+            solver.prefactorize(prefactorized_data, prefactorized_block_perm);
+            solver.prefactorize_and_solve(data, block_perm, rhs, x);
+            CHECK(prefactorized_data == data);
+        }
     }
 
-    SECTION("Test block(double 2*2) calculation") {
+    SECTION("Block(double 2*2) calculation") {
         // [  0 1   1   2   3   4           3             38
         //  100 0   7  -1   5   6           4            356
         //    1 2   0 200   f   f       * [ -1 ]   =  [ -389 ]
@@ -105,17 +113,17 @@ TEST_CASE("Test Sparse LU solver") {
         SparseLUSolver<Tensor, Array, Array>::BlockPermArray block_perm(3);
 
         SECTION("Test calculation") {
-            solver.solve(data, block_perm, rhs, x);
+            solver.prefactorize_and_solve(data, block_perm, rhs, x);
             check_result(x, x_ref);
         }
         SECTION("Test (pseudo) singular") {
             data[0](0, 1) = 0.0;
-            CHECK_THROWS_AS(solver.solve(data, block_perm, rhs, x), SparseMatrixError);
+            CHECK_THROWS_AS(solver.prefactorize_and_solve(data, block_perm, rhs, x), SparseMatrixError);
         }
 
         SECTION("Test prefactorize") {
             solver.prefactorize(data, block_perm);
-            solver.solve((std::vector<Tensor> const&)data, block_perm, rhs, x);
+            solver.solve_with_prefactorization((std::vector<Tensor> const&)data, block_perm, rhs, x);
             check_result(x, x_ref);
         }
     }

tests/cpp_unit_tests/test_y_bus.cpp

@@ -211,7 +211,7 @@ TEST_CASE("Test one bus system") {
 TEST_CASE("Test fill-in y bus") {
     /*
      * struct
-     * [1] --0--> [0] --[1]--> [2]
+     * [1] --0--> [0] --1--> [2]
      * extra fill-in: (1, 2) by removing node 0
      *
      * [