Finalize NRSE

power_grid_model_c/power_grid_model/include/power_grid_model/grouped_index_vector.hpp

@@ -260,7 +260,7 @@ class DenseGroupedIdxVector {
         : DenseGroupedIdxVector{std::move(dense_group_elements), num_groups} {}
 
   private:
-    Idx num_groups_;
+    Idx num_groups_{};
     IdxVector dense_vector_;
 };
 

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/block_matrix.hpp

@@ -15,22 +15,18 @@
 namespace power_grid_model::math_solver {
 
 template <scalar_value T, bool sym, bool is_tensor, int n_sub_block> struct block_trait {
-    static constexpr int n_row = sym ? n_sub_block : n_sub_block * 3;
-    static constexpr int n_col = [] {
-        if constexpr (is_tensor) {
-            return sym ? n_sub_block : n_sub_block * 3;
-        } else {
-            return 1;
-        }
-    }();
+    static constexpr int sub_block_size = sym ? 1 : 3;
+    static constexpr int n_row = n_sub_block * sub_block_size;
+    static constexpr int n_col = is_tensor ? n_sub_block * sub_block_size : 1;
 
     using ArrayType = Eigen::Array<T, n_row, n_col, Eigen::ColMajor>;
 };
 
 template <class T, bool sym, bool is_tensor, int n_sub_block>
 class Block : public block_trait<T, sym, is_tensor, n_sub_block>::ArrayType {
   public:
-    using ArrayType = typename block_trait<T, sym, is_tensor, n_sub_block>::ArrayType;
+    using block_traits = block_trait<T, sym, is_tensor, n_sub_block>;
+    using ArrayType = typename block_traits::ArrayType;
     using ArrayType::operator();
 
     // default zero
@@ -42,24 +38,52 @@ class Block : public block_trait<T, sym, is_tensor, n_sub_block>::ArrayType {
         return *this;
     }
 
-    template <int r> static auto get_asym_row_idx() { return Eigen::seqN(Eigen::fix<r * 3>, Eigen::fix<3>); }
-    template <int c> static auto get_asym_col_idx() {
+    template <int start, int size, int n_total> static auto get_sequence() {
+        static_assert(start >= 0);
+        static_assert(size >= 0);
+        static_assert(start + size <= n_total);
+        return Eigen::seqN(Eigen::fix<start>, Eigen::fix<size>);
+    }
+    template <int block_start, int block_size, int n_total> static auto get_block_sequence() {
+        constexpr auto size = block_size * block_traits::sub_block_size;
+        constexpr auto start = block_start * size;
+        return Block::get_sequence<start, size, n_total>();
+    }
+    template <int r, int r_size> static auto get_block_row_idx() {
+        return Block::get_block_sequence<r, r_size, block_traits::n_row>();
+    }
+    template <int c, int c_size> static auto get_block_col_idx() {
         if constexpr (is_tensor) {
-            return Eigen::seqN(Eigen::fix<c * 3>, Eigen::fix<3>);
+            return Block::get_block_sequence<c, c_size, block_traits::n_col>();
         } else {
-            return Eigen::seqN(Eigen::fix<0>, Eigen::fix<1>);
+            return Block::get_sequence<c, c_size, block_traits::n_col>();
         }
     }
+    template <int r> static auto get_row_idx() { return Block::get_block_row_idx<r, 1>(); }
+    template <int c> static auto get_col_idx() { return Block::get_block_col_idx<c, 1>(); }
 
     template <int r, int c>
     using GetterType =
-        std::conditional_t<sym, T&, decltype(std::declval<ArrayType>()(get_asym_row_idx<r>(), get_asym_col_idx<c>()))>;
+        std::conditional_t<sym, T&, decltype(std::declval<ArrayType>()(get_row_idx<r>(), get_col_idx<c>()))>;
 
     template <int r, int c> GetterType<r, c> get_val() {
         if constexpr (sym) {
             return (*this)(r, c);
         } else {
-            return (*this)(get_asym_row_idx<r>(), get_asym_col_idx<c>());
+            return (*this)(get_row_idx<r>(), get_col_idx<c>());
+        }
+    }
+
+    template <int r, int c, int r_size, int c_size>
+    using BlockGetterType = std::conditional_t<r_size == 1 && c_size == 1, GetterType<r, c>,
+                                               decltype(std::declval<ArrayType>()(get_block_row_idx<r, r_size>(),
+                                                                                  get_block_col_idx<c, c_size>()))>;
+
+    template <int r, int c, int r_size, int c_size> BlockGetterType<r, c, r_size, c_size> get_block_val() {
+        if constexpr (r_size == 1 && c_size == 1) {
+            return get_val<r, c>();
+        } else {
+            return (*this)(get_block_row_idx<r, r_size>(), get_block_col_idx<c, c_size>());
         }
     }
 

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp

@@ -51,6 +51,8 @@ template <bool sym> using NRSERhs = NRSEUnknown<sym>;
 template <bool sym> class NRSEGainBlock : public Block<double, sym, true, 4> {
   public:
     template <int r, int c> using GetterType = typename Block<double, sym, true, 4>::template GetterType<r, c>;
+    template <int r, int c, int r_size, int c_size>
+    using BlockGetterType = typename Block<double, sym, true, 4>::template BlockGetterType<r, c, r_size, c_size>;
 
     // eigen expression
     using Block<double, sym, true, 4>::Block;
@@ -61,6 +63,7 @@ template <bool sym> class NRSEGainBlock : public Block<double, sym, true, 4> {
     GetterType<1, 0> g_Q_theta() { return this->template get_val<1, 0>(); }
     GetterType<1, 1> g_Q_v() { return this->template get_val<1, 1>(); }
 
+    // (Q^T)(X, Y)
     GetterType<0, 2> qt_P_theta() { return this->template get_val<0, 2>(); }
     GetterType<0, 3> qt_P_v() { return this->template get_val<0, 3>(); }
     GetterType<1, 2> qt_Q_theta() { return this->template get_val<1, 2>(); }
@@ -75,6 +78,11 @@ template <bool sym> class NRSEGainBlock : public Block<double, sym, true, 4> {
     GetterType<2, 3> r_P_v() { return this->template get_val<2, 3>(); }
     GetterType<3, 2> r_Q_theta() { return this->template get_val<3, 2>(); }
     GetterType<3, 3> r_Q_v() { return this->template get_val<3, 3>(); }
+
+    BlockGetterType<0, 0, 2, 2> g() { return this->template get_block_val<0, 0, 2, 2>(); }
+    BlockGetterType<0, 1, 2, 2> qt() { return this->template get_block_val<0, 1, 2, 2>(); }
+    BlockGetterType<1, 0, 2, 2> q() { return this->template get_block_val<1, 0, 2, 2>(); }
+    BlockGetterType<1, 1, 2, 2> r() { return this->template get_block_val<1, 1, 2, 2>(); }
 };
 
 // solver
@@ -159,8 +167,8 @@ template <bool sym> class NewtonRaphsonSESolver {
             sub_timer = Timer(calculation_info, 2224, "Prepare LHS rhs");
             prepare_matrix_and_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_with_prefactorized_matrix(data_gain_, perm_, delta_x_rhs_, delta_x_rhs_);
+            sub_timer = Timer(calculation_info, 2225, "Solve sparse linear equation");
+            sparse_solver_.prefactorize_and_solve(data_gain_, perm_, delta_x_rhs_, delta_x_rhs_);
             sub_timer = Timer(calculation_info, 2226, "Iterate unknown");
             max_dev = iterate_unknown(output.u, measured_values);
         };
@@ -245,20 +253,9 @@ template <bool sym> class NewtonRaphsonSESolver {
 
             for (Idx data_idx_lu = row_indptr[row]; data_idx_lu != row_indptr[row + 1]; ++data_idx_lu) {
                 // get data idx of y bus,
-                // skip for a fill-in
                 Idx const data_idx = y_bus.map_lu_y_bus()[data_idx_lu];
-                if (data_idx == -1) {
-                    continue;
-                }
-
                 Idx const col = col_indices[data_idx_lu];
 
-                u_state.uj = current_u[col];
-                u_state.abs_uj_inv = diagonal_inverse(x_[col].v());
-                u_state.uj_uj_conj = vector_outer_product(u_state.uj, conj(u_state.uj));
-                u_state.ui_uj_conj = vector_outer_product(u_state.ui, conj(u_state.uj));
-                u_state.uj_ui_conj = vector_outer_product(u_state.uj, conj(u_state.ui));
-
                 // get a reference and reset block to zero
                 NRSEGainBlock<sym>& block = data_gain_[data_idx_lu];
                 if (row == col) {
@@ -269,6 +266,17 @@ template <bool sym> class NewtonRaphsonSESolver {
                     block.clear();
                 }
 
+                // skip anything else for a fill-in
+                if (data_idx == -1) {
+                    continue;
+                }
+
+                u_state.uj = current_u[col];
+                u_state.abs_uj_inv = diagonal_inverse(x_[col].v());
+                u_state.uj_uj_conj = vector_outer_product(u_state.uj, conj(u_state.uj));
+                u_state.ui_uj_conj = vector_outer_product(u_state.ui, conj(u_state.uj));
+                u_state.uj_ui_conj = vector_outer_product(u_state.uj, conj(u_state.ui));
+
                 // fill block with branch, shunt measurement
                 for (Idx element_idx = y_bus.y_bus_entry_indptr()[data_idx];
                      element_idx != y_bus.y_bus_entry_indptr()[data_idx + 1]; ++element_idx) {
@@ -323,9 +331,6 @@ template <bool sym> class NewtonRaphsonSESolver {
 
         fill_qt(y_bus);
         process_lagrange_multiplier(y_bus);
-
-        // prefactorize
-        sparse_solver_.prefactorize(data_gain_, perm_);
     }
 
     /// Q_ij = 0
@@ -462,12 +467,7 @@ template <bool sym> class NewtonRaphsonSESolver {
     void fill_qt(YBus<sym> const& y_bus) {
         iterate_matrix_skip_fills(
             [this](Idx /* row */, Idx /* col */, Idx data_idx, Idx data_idx_transpose) {
-                auto& block = data_gain_[data_idx];
-
-                block.qt_P_theta() = data_gain_[data_idx_transpose].q_P_theta();
-                block.qt_P_v() = data_gain_[data_idx_transpose].q_Q_theta();
-                block.qt_Q_theta() = data_gain_[data_idx_transpose].q_P_v();
-                block.qt_Q_v() = data_gain_[data_idx_transpose].q_Q_v();
+                data_gain_[data_idx].qt() = hermitian_transpose(data_gain_[data_idx_transpose].q());
             },
             y_bus);
     }

tests/cpp_unit_tests/test_grouped_index_vector.cpp

@@ -55,7 +55,13 @@ TEST_CASE_TEMPLATE("Grouped idx data structure", IdxVectorConstructor, TypePair<
 
     auto const idx_vector = construct_from<IdxVectorType, ConstructFromTag>(groups, num_groups);
 
-    SUBCASE("Empty grouped idx vector") {
+    SUBCASE("Empty grouped idx vector - no explicit initialization") {
+        IdxVectorType const indices;
+        CHECK(indices.element_size() == 0);
+        CHECK(indices.size() == 0);
+    }
+
+    SUBCASE("Empty grouped idx vector - explicit initialization") {
         IdxVectorType const indices{};
         CHECK(indices.element_size() == 0);
         CHECK(indices.size() == 0);
@@ -124,6 +130,13 @@ TEST_CASE_TEMPLATE("Enumerated zip iterator for grouped index data structures",
     auto const idx_vector_b = construct_from<B, from_natural_t>(groups_b, num_groups);
     auto const idx_vector_c = construct_from<C, from_natural_t>(groups_c, num_groups);
 
+    SUBCASE("empty input") {
+        auto const empty_idx_vector = A{};
+        for ([[maybe_unused]] auto [index, element_range] : enumerated_zip_sequence(empty_idx_vector)) {
+            FAIL("this code should not be reached");
+        }
+    }
+
     SUBCASE("1 input") {
         // Test single zipped iteration
         IdxRanges actual_ranges_a{};