Merge links on nodes: included set solution system

power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp

@@ -224,8 +224,7 @@ inline void backward_substitution(EliminationResult& elimination_result) {
 }
 
 // reduced echelon form based on custom forward elimination and backward substitution procedures
-[[nodiscard]] inline EliminationResult reduced_echelon_form(std::vector<BranchIdx> edges,
-                                                            std::vector<DoubleComplex> node_loads) {
+inline EliminationResult reduced_echelon_form(std::vector<BranchIdx> edges, std::vector<DoubleComplex> node_loads) {
     EliminationResult result{};
     auto const edge_number{edges.size()};
     auto const node_number{node_loads.size()};
@@ -241,4 +240,50 @@ inline void backward_substitution(EliminationResult& elimination_result) {
     backward_substitution(result);
     return result;
 }
+
+// The degrees of freedom matrix (dfs_matrix) is associated with the degrees of freedom vector according
+// internal_loads = extended_rhs - dfs_matrix * lambda
+struct SolutionSet {
+    CooSparseMatrix dfs_matrix{};
+    std::vector<DoubleComplex> extended_rhs{};
+};
+
+// Constructs the dfs_matrix and the extended_rhs for the set of solutions.
+// The dfs_matrix is constructed from the matrix rows associated with the free column indices
+// and from negative unit elements at locations where the variables are associated with the free parameters.
+inline SolutionSet set_solution_system(EliminationResult& result) {
+    SolutionSet solution_set{};
+
+    auto& [dfs_matrix, extended_rhs] = solution_set;
+    auto const pivot_indices_size = narrow_cast<Idx>(result.pivot_edge_indices.size());
+    auto const free_indices_size = narrow_cast<Idx>(result.free_edge_indices.size());
+    Idx const total_indices_size = pivot_indices_size + free_indices_size;
+    dfs_matrix.prepare(total_indices_size, free_indices_size);
+    extended_rhs.resize(total_indices_size);
+    constexpr auto const free_matrix_element = IntS{-1};
+    Idx free_edge_idx;
+
+    // The part constructed from result.matrix and result.rhs.
+    for (auto matrix_row : std::views::iota(Idx{}, pivot_indices_size)) {
+        auto const pivot_edge_idx = result.pivot_edge_indices[matrix_row];
+        for (auto dfs_matrix_col : std::views::iota(Idx{}, free_indices_size)) {
+            free_edge_idx = result.free_edge_indices[dfs_matrix_col];
+            result.matrix.get_value(matrix_row, free_edge_idx)
+                .transform([&dfs_matrix, pivot_edge_idx, dfs_matrix_col](IntS matrix_element) {
+                    dfs_matrix.set_value(matrix_element, pivot_edge_idx, dfs_matrix_col);
+                    return matrix_element;
+                });
+        }
+        extended_rhs[pivot_edge_idx] = result.rhs[matrix_row];
+    }
+
+    // The part constructed from negative unit elements.
+    for (auto dfs_matrix_col : std::views::iota(Idx{}, free_indices_size)) {
+        free_edge_idx = result.free_edge_indices[dfs_matrix_col];
+        dfs_matrix.set_value(free_matrix_element, free_edge_idx, dfs_matrix_col);
+        extended_rhs[result.free_edge_indices[dfs_matrix_col]] = DoubleComplex{};
+    }
+    return solution_set;
+};
+
 } // namespace power_grid_model::link_solver::detail

tests/cpp_unit_tests/test_link_solver.cpp

@@ -8,7 +8,9 @@
 
 #include <doctest/doctest.h>
 
+#include <cstddef>
 #include <ranges>
+#include <span>
 #include <unordered_map>
 #include <unordered_set>
 #include <vector>
@@ -346,7 +348,7 @@ TEST_CASE("Test the link solver algorithm") {
             result.matrix.set_value(1, 3, 5);
             result.matrix.set_value(1, 3, 6);
 
-            result.rhs = std::vector<DoubleComplex>{{-1.0, -1.0}, {-1.0, -1.0}, {-2.0, -2.0}, {0.0, 0.0}, {0.0, 0.0}};
+            result.rhs = std::vector<DoubleComplex>{{-1.0, -1.0}, {-1.0, -1.0}, {-2.0, -2.0}, {0.0, 0.0}};
 
             result.free_edge_indices = std::vector<Idx>{3, 4, 6};
             result.pivot_edge_indices = std::vector<Idx>{0, 1, 2, 5};
@@ -380,9 +382,86 @@ TEST_CASE("Test the link solver algorithm") {
             CHECK(-1 == result.matrix.get_value(2, 4));
             CHECK(1 == result.matrix.get_value(3, 5));
             CHECK(1 == result.matrix.get_value(3, 6));
-            REQUIRE(result.rhs.size() == 5);
-            CHECK(result.rhs ==
-                  std::vector<DoubleComplex>{{0.0, 0.0}, {1.0, 1.0}, {-2.0, -2.0}, {0.0, 0.0}, {0.0, 0.0}});
+            REQUIRE(result.rhs.size() == 4);
+
+            CHECK(result.rhs == std::vector<DoubleComplex>{{0.0, 0.0}, {1.0, 1.0}, {-2.0, -2.0}, {0.0, 0.0}});
+        }
+    }
+    SUBCASE("Testing the set_solution_system routine") {
+        auto generate_input_result = [](std::span<const IntS> data, std::span<const Idx> row, std::span<const Idx> col,
+                                        Idx row_number, Idx col_number) {
+            EliminationResult result{};
+            result.matrix.prepare(row_number, col_number);
+            for (auto idx = size_t{0}; idx < data.size(); ++idx) {
+                result.matrix.set_value(data[idx], row[idx], col[idx]);
+            }
+            return result;
+        };
+
+        SUBCASE("Complex case with complex loads") {
+            std::vector<IntS> data = {1, 0, 0, 1, 0, 0, -1, 1, -1, -1, 0, 0, 1, 1, 1, 1, 0, 1};
+            std::vector<Idx> row = {0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 1, 0, 0, 0};
+            std::vector<Idx> col = {0, 1, 2, 1, 2, 3, 6, 2, 3, 4, 5, 6, 5, 6, 4, 3, 4, 6};
+
+            auto result = generate_input_result(data, row, col, Idx{5}, Idx{7});
+
+            result.rhs = {{0, 0}, {1, 1}, {-2, -2}, {-0, -0}};
+            result.free_edge_indices = {3, 4, 6};
+            result.pivot_edge_indices = {0, 1, 2, 5};
+
+            SolutionSet const solution_set = set_solution_system(result);
+
+            CHECK(1 == solution_set.dfs_matrix.get_value(0, 0));
+            CHECK(1 == solution_set.dfs_matrix.get_value(0, 2));
+            CHECK(1 == solution_set.dfs_matrix.get_value(1, 1));
+            CHECK(-1 == solution_set.dfs_matrix.get_value(1, 2));
+            CHECK(-1 == solution_set.dfs_matrix.get_value(2, 0));
+            CHECK(-1 == solution_set.dfs_matrix.get_value(2, 1));
+            CHECK(-1 == solution_set.dfs_matrix.get_value(3, 0));
+            CHECK(-1 == solution_set.dfs_matrix.get_value(4, 1));
+            CHECK(1 == solution_set.dfs_matrix.get_value(5, 2));
+            CHECK(-1 == solution_set.dfs_matrix.get_value(6, 2));
+            CHECK(solution_set.extended_rhs ==
+                  std::vector<DoubleComplex>({{0, 0}, {1, 1}, {-2, -2}, {0, 0}, {0, 0}, {0, 0}, {0, 0}}));
+        }
+
+        SUBCASE("Two edges, two nodes, two real loads") {
+
+            std::vector<IntS> data = {1, 0, 1, -1, -1};
+            std::vector<Idx> row = {0, 0, 1, 1, 0};
+            std::vector<Idx> col = {0, 1, 1, 2, 2};
+
+            auto result = generate_input_result(data, row, col, Idx{1}, Idx{3});
+
+            result.rhs = {{1, 0}, {0, 0}};
+            result.free_edge_indices = {2};
+            result.pivot_edge_indices = {0, 1};
+
+            SolutionSet const solution_set = set_solution_system(result);
+
+            CHECK(-1 == solution_set.dfs_matrix.get_value(0, 0));
+            CHECK(-1 == solution_set.dfs_matrix.get_value(1, 0));
+            CHECK(-1 == solution_set.dfs_matrix.get_value(2, 0));
+            CHECK(solution_set.extended_rhs == std::vector<DoubleComplex>({{1, 0}, {0, 0}, {0, 0}}));
+        }
+
+        SUBCASE("Four edges, four nodes, two real loads") {
+            std::vector<IntS> data = {1, 0, 0, 1, 1, 1, 1};
+            std::vector<Idx> row = {0, 0, 0, 1, 1, 2, 2};
+            std::vector<Idx> col = {0, 1, 3, 1, 3, 2, 3};
+
+            auto result = generate_input_result(data, row, col, Idx{1}, Idx{4});
+
+            result.rhs = {{1, 0}, {-1, 0}, {-1, 0}};
+            result.free_edge_indices = {3};
+            result.pivot_edge_indices = {0, 1, 2};
+
+            const SolutionSet solution_set = set_solution_system(result);
+
+            CHECK(1 == solution_set.dfs_matrix.get_value(1, 0));
+            CHECK(1 == solution_set.dfs_matrix.get_value(2, 0));
+            CHECK(-1 == solution_set.dfs_matrix.get_value(3, 0));
+            CHECK(solution_set.extended_rhs == std::vector<DoubleComplex>({{1, 0}, {-1, 0}, {-1, 0}, {0, 0}}));
         }
     }
 }