From e3e5be814ac55f4e28a4305f869b5127aaec913d Mon Sep 17 00:00:00 2001 From: Santiago Figueroa Manrique Date: Tue, 10 Mar 2026 14:49:57 +0100 Subject: [PATCH 01/12] forward elimination implementation Signed-off-by: Santiago Figueroa Manrique --- .../include/power_grid_model/link_solver.hpp | 145 ++++++++++++++++++ tests/cpp_unit_tests/CMakeLists.txt | 1 + 2 files changed, 146 insertions(+) create mode 100644 power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp new file mode 100644 index 0000000000..5e5e70edd4 --- /dev/null +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp @@ -0,0 +1,145 @@ +// SPDX-FileCopyrightText: Contributors to the Power Grid Model project +// +// SPDX-License-Identifier: MPL-2.0 + +#pragma once + +#include "calculation_parameters.hpp" +#include "common/common.hpp" +#include "common/counting_iterator.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace power_grid_model::link_solver { + +namespace detail { +constexpr uint8_t starting_row{0}; + +enum class EdgeEvent : std::uint8_t { + Deleted = 0, // pivot edge - used as pivot row + Replaced = 1, // reattached to a different node + ContractedToPoint = 2 // from == to, becomes a degree of freedom +}; + +enum class EdgeDirection : IntS { Away = -1, Towards = 1 }; + +struct COOSparseMatrix { + std::vector data; + std::vector row; + std::vector col; + + void append(IntS value, uint64_t row_idx, uint64_t col_idx) { + data.push_back(value); + row.push_back(row_idx); + col.push_back(col_idx); + } +}; + +struct EdgeHistory { + std::vector rows{}; + std::vector events{}; +}; + +struct EliminationResult { + COOSparseMatrix matrix; + std::vector rhs; // RHS value at each pivot row + std::vector free_edge_indices; // index of degrees of freedom (self loop edges) + std::vector edge_history; // edges elimination history +}; + +// convention: from node at position 0, to node at position 1 +[[nodiscard]] constexpr Idx from_node(BranchIdx const& edge) { return edge[0]; } +[[nodiscard]] constexpr Idx& from_node(BranchIdx& edge) { return edge[0]; } +[[nodiscard]] constexpr Idx to_node(BranchIdx const& edge) { return edge[1]; } +[[nodiscard]] constexpr Idx& to_node(BranchIdx& edge) { return edge[1]; } + +// map from node index to the set of adjacent edge indices +// unordered_map because node IDs may be sparse/non-contiguous +// unordered_set for O(1) insert/erase during reattachment +[[nodiscard]] auto build_adjacency_list(std::span edges) { + std::unordered_map> adjacency_list{}; + for (auto const& [raw_index, edge] : enumerate(edges)) { + auto const index = static_cast(raw_index); + auto const [from_node, to_node] = edge; + adjacency_list[from_node].insert(index); + adjacency_list[to_node].insert(index); + } + return adjacency_list; +} + +void write_edge_history(EdgeHistory& edge_history, EdgeEvent event, uint64_t row) { + edge_history.events.push_back(event); + edge_history.rows.push_back(row); +} + +// forward elimination is performed via a modified Gaussian elimination procedure +// we name this procedure elimination-game +[[nodiscard]] EliminationResult forward_elimination(std::vector edges, + std::vector node_loads) { + EliminationResult result{.edge_history = std::vector(edges.size())}; + using enum EdgeEvent; + using enum EdgeDirection; + + auto& [matrix, rhs, free_edge_indices, edge_history] = result; + + std::unordered_set edges_contracted_to_point{}; + uint64_t matrix_row{starting_row}; + auto adjacency_list = build_adjacency_list(edges); + + for (auto const& [raw_index, edge] : enumerate(edges)) { + auto const index = static_cast(raw_index); + if (edges_contracted_to_point.contains(index)) { + free_edge_indices.push_back(index); + } else { + write_edge_history(edge_history[index], Deleted, matrix_row); // Delete edge -> pivot there + matrix.append(std::to_underlying(Towards), matrix_row, index); + + Idx const from_node_idx = from_node(edge); + Idx const to_node_idx = to_node(edge); + + // update adjacency list for deleted edge + adjacency_list[from_node_idx].erase(index); + adjacency_list[to_node_idx].erase(index); + + // Gaussian elimination like steps + node_loads[from_node_idx] += node_loads[to_node_idx]; + rhs.push_back(node_loads[to_node_idx]); + + auto const adjacent_edges_snapshot = + std::vector(adjacency_list[to_node_idx].begin(), adjacency_list[to_node_idx].end()); + + for (uint64_t const adjacent_edge_idx : adjacent_edges_snapshot) { + if (from_node(edges[adjacent_edge_idx]) == to_node_idx) { + from_node(edges[adjacent_edge_idx]) = from_node_idx; // re-attachment step + matrix.append(std::to_underlying(Away), matrix_row, adjacent_edge_idx); + } else { + to_node(edges[adjacent_edge_idx]) = from_node_idx; // re-attachment step + matrix.append(std::to_underlying(Towards), matrix_row, adjacent_edge_idx); + } + + // update adjacency list after re-attachment + adjacency_list[to_node_idx].erase(adjacent_edge_idx); + adjacency_list[from_node_idx].insert(adjacent_edge_idx); + + if (from_node(edges[adjacent_edge_idx]) == to_node(edges[adjacent_edge_idx])) { + write_edge_history(edge_history[adjacent_edge_idx], ContractedToPoint, matrix_row); + edges_contracted_to_point.insert(adjacent_edge_idx); + } else { + write_edge_history(edge_history[adjacent_edge_idx], Replaced, matrix_row); + } + } + ++matrix_row; + } + } + + return result; +} +} // namespace detail +} // namespace power_grid_model::link_solver diff --git a/tests/cpp_unit_tests/CMakeLists.txt b/tests/cpp_unit_tests/CMakeLists.txt index d3a4fee8ad..7257ca9f08 100644 --- a/tests/cpp_unit_tests/CMakeLists.txt +++ b/tests/cpp_unit_tests/CMakeLists.txt @@ -130,6 +130,7 @@ add_executable( "test_index_mapping.cpp" "test_job_dispatch.cpp" "test_calculation_preparation.cpp" + "test_link_solver.cpp" ) target_link_libraries( From 12e9f5acb2702b479fe3fdf74ebd8d1950373dae Mon Sep 17 00:00:00 2001 From: Santiago Figueroa Manrique Date: Tue, 10 Mar 2026 14:50:13 +0100 Subject: [PATCH 02/12] test forward elimination Signed-off-by: Santiago Figueroa Manrique --- tests/cpp_unit_tests/test_link_solver.cpp | 138 ++++++++++++++++++++++ 1 file changed, 138 insertions(+) create mode 100644 tests/cpp_unit_tests/test_link_solver.cpp diff --git a/tests/cpp_unit_tests/test_link_solver.cpp b/tests/cpp_unit_tests/test_link_solver.cpp new file mode 100644 index 0000000000..013a7801b8 --- /dev/null +++ b/tests/cpp_unit_tests/test_link_solver.cpp @@ -0,0 +1,138 @@ +// SPDX-FileCopyrightText: Contributors to the Power Grid Model project +// +// SPDX-License-Identifier: MPL-2.0 + +#include + +#include + +namespace power_grid_model::link_solver { +TEST_CASE("Test the link solver algorithm") { + using namespace detail; + + SUBCASE("Test forward elimination - elimination game") { + using enum EdgeEvent; + + SUBCASE("One edge, two nodes, two real loads") { + auto const edges = std::vector{{0, 1}}; + auto const node_loads = std::vector{{-1.0, 0.0}, {1.0, 0.0}}; + auto const result = forward_elimination(edges, node_loads); + + REQUIRE(result.matrix.data.size() == 1); + CHECK(result.matrix.data == std::vector{1}); + REQUIRE(result.matrix.row.size() == 1); + CHECK(result.matrix.row == std::vector{0}); + REQUIRE(result.matrix.col.size() == 1); + CHECK(result.matrix.col == std::vector{0}); + REQUIRE(result.rhs.size() == 1); + CHECK(result.rhs == std::vector{{1.0, 0.0}}); + CHECK(result.free_edge_indices.empty()); + REQUIRE(result.edge_history.size() == 1); + CHECK(result.edge_history[0].events == + std::vector{Deleted}); // Should this be deleted or contracted to point? + CHECK(result.edge_history[0].rows == std::vector{0}); + } + + SUBCASE("Two edges, three nodes, two real loads") { + auto const edges = std::vector{{1, 0}, {1, 2}}; + auto const node_loads = std::vector{{-1.0, 0.0}, {1.0, 0.0}, {0.0, 0.0}}; + auto const result = forward_elimination(edges, node_loads); + + REQUIRE(result.matrix.data.size() == 2); + CHECK(result.matrix.data == std::vector{1, 1}); + REQUIRE(result.matrix.row.size() == 2); + CHECK(result.matrix.row == std::vector{0, 1}); + REQUIRE(result.matrix.col.size() == 2); + CHECK(result.matrix.col == std::vector{0, 1}); + REQUIRE(result.rhs.size() == 2); + CHECK(result.rhs == std::vector{{-1.0, 0.0}, {0.0, 0.0}}); + CHECK(result.free_edge_indices.empty()); + REQUIRE(result.edge_history.size() == 2); + CHECK(result.edge_history[0].events == std::vector{Deleted}); + CHECK(result.edge_history[0].rows == std::vector{0}); + CHECK(result.edge_history[1].events == std::vector{Deleted}); + CHECK(result.edge_history[1].rows == std::vector{1}); + } + + SUBCASE("Three edges, three nodes, two real loads") { + auto const edges = std::vector{{0, 1}, {1, 2}, {2, 0}}; + auto const node_loads = std::vector{{-1.0, 0.0}, {1.0, 0.0}, {0.0, 0.0}}; + auto const result = forward_elimination(edges, node_loads); + + REQUIRE(result.matrix.data.size() == 4); + CHECK(result.matrix.data == std::vector{1, -1, 1, -1}); + REQUIRE(result.matrix.row.size() == 4); + CHECK(result.matrix.row == std::vector{0, 0, 1, 1}); + REQUIRE(result.matrix.col.size() == 4); + CHECK(result.matrix.col == std::vector{0, 1, 1, 2}); + REQUIRE(result.rhs.size() == 2); + CHECK(result.rhs == std::vector{{1.0, 0.0}, {0.0, 0.0}}); + REQUIRE(result.free_edge_indices.size() == 1); + CHECK(result.free_edge_indices == std::vector{2}); + REQUIRE(result.edge_history.size() == 3); + CHECK(result.edge_history[0].events == std::vector{Deleted}); + CHECK(result.edge_history[0].rows == std::vector{0}); + CHECK(result.edge_history[1].events == std::vector{Replaced, Deleted}); + CHECK(result.edge_history[1].rows == std::vector{0, 1}); + CHECK(result.edge_history[2].events == std::vector{ContractedToPoint}); + CHECK(result.edge_history[2].rows == std::vector{1}); + } + + SUBCASE("Two edges, two nodes, two real loads") { + auto const edges = std::vector{{0, 1}, {0, 1}}; + auto const node_loads = std::vector{{-1.0, 0.0}, {1.0, 0.0}}; + auto const result = forward_elimination(edges, node_loads); + + REQUIRE(result.matrix.data.size() == 2); + CHECK(result.matrix.data == std::vector{1, 1}); + REQUIRE(result.matrix.row.size() == 2); + CHECK(result.matrix.row == std::vector{0, 0}); + REQUIRE(result.matrix.col.size() == 2); + CHECK(result.matrix.col == std::vector{0, 1}); + REQUIRE(result.rhs.size() == 1); + CHECK(result.rhs == std::vector{{1.0, 0.0}}); + REQUIRE(result.free_edge_indices.size() == 1); + CHECK(result.free_edge_indices == std::vector{1}); + REQUIRE(result.edge_history.size() == 2); + CHECK(result.edge_history[0].events == std::vector{Deleted}); + CHECK(result.edge_history[0].rows == std::vector{0}); + CHECK(result.edge_history[1].events == std::vector{ContractedToPoint}); + CHECK(result.edge_history[1].rows == std::vector{0}); + } + + SUBCASE("Complex case with complex loads") { + auto const edges = std::vector{{3, 0}, {1, 0}, {2, 0}, {3, 2}, {1, 2}, {1, 4}, {3, 4}}; + auto const node_loads = + std::vector{{-1.0, -1.0}, {-1.0, -1.0}, {2.0, 2.0}, {0.0, 0.0}, {0.0, 0.0}}; + auto const result = forward_elimination(edges, node_loads); + + REQUIRE(result.matrix.data.size() == 14); + CHECK(result.matrix.data == std::vector{1, 1, 1, 1, 1, -1, -1, 1, -1, -1, -1, -1, 1, 1}); + REQUIRE(result.matrix.row.size() == 14); + CHECK(result.matrix.row == std::vector{0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3}); + REQUIRE(result.matrix.col.size() == 14); + CHECK(result.matrix.col == std::vector{0, 2, 1, 1, 2, 6, 3, 2, 3, 6, 5, 4, 5, 6}); + REQUIRE(result.rhs.size() == 4); + CHECK(result.rhs == std::vector{{-1.0, -1.0}, {-1.0, -1.0}, {-2.0, -2.0}, {0.0, 0.0}}); + REQUIRE(result.free_edge_indices.size() == 3); + CHECK(result.free_edge_indices == std::vector{3, 4, 6}); + REQUIRE(result.edge_history.size() == 7); + CHECK(result.edge_history[0].events == std::vector{Deleted}); + CHECK(result.edge_history[0].rows == std::vector{0}); + CHECK(result.edge_history[1].events == std::vector{Replaced, Deleted}); + CHECK(result.edge_history[1].rows == std::vector{0, 1}); + CHECK(result.edge_history[2].events == std::vector{Replaced, Replaced, Deleted}); + CHECK(result.edge_history[2].rows == std::vector{0, 1, 2}); + CHECK(result.edge_history[3].events == std::vector{Replaced, ContractedToPoint}); + CHECK(result.edge_history[3].rows == std::vector{1, 2}); + CHECK(result.edge_history[4].events == std::vector{ContractedToPoint}); + CHECK(result.edge_history[4].rows == std::vector{2}); + CHECK(result.edge_history[5].events == std::vector{Replaced, Deleted}); + CHECK(result.edge_history[5].rows == std::vector{2, 3}); + CHECK(result.edge_history[6].events == std::vector{Replaced, Replaced, ContractedToPoint}); + CHECK(result.edge_history[6].rows == std::vector{1, 2, 3}); + } + } +} + +} // namespace power_grid_model::link_solver From c185dbe0144e090e6145c8f162fa0ec99a64faa5 Mon Sep 17 00:00:00 2001 From: Santiago Figueroa Manrique Date: Tue, 10 Mar 2026 15:14:14 +0100 Subject: [PATCH 03/12] minor Signed-off-by: Santiago Figueroa Manrique --- .../include/power_grid_model/link_solver.hpp | 17 +++++++++-------- tests/cpp_unit_tests/test_link_solver.cpp | 8 ++++++-- 2 files changed, 15 insertions(+), 10 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp index 5e5e70edd4..a8b58ba100 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp @@ -31,9 +31,9 @@ enum class EdgeEvent : std::uint8_t { enum class EdgeDirection : IntS { Away = -1, Towards = 1 }; struct COOSparseMatrix { - std::vector data; - std::vector row; - std::vector col; + std::vector data{}; + std::vector row{}; + std::vector col{}; void append(IntS value, uint64_t row_idx, uint64_t col_idx) { data.push_back(value); @@ -48,10 +48,10 @@ struct EdgeHistory { }; struct EliminationResult { - COOSparseMatrix matrix; - std::vector rhs; // RHS value at each pivot row - std::vector free_edge_indices; // index of degrees of freedom (self loop edges) - std::vector edge_history; // edges elimination history + COOSparseMatrix matrix{}; + std::vector rhs{}; // RHS value at each pivot row + std::vector free_edge_indices{}; // index of degrees of freedom (self loop edges) + std::vector edge_history{}; // edges elimination history }; // convention: from node at position 0, to node at position 1 @@ -83,7 +83,8 @@ void write_edge_history(EdgeHistory& edge_history, EdgeEvent event, uint64_t row // we name this procedure elimination-game [[nodiscard]] EliminationResult forward_elimination(std::vector edges, std::vector node_loads) { - EliminationResult result{.edge_history = std::vector(edges.size())}; + EliminationResult result{}; + result.edge_history.resize(edges.size()); using enum EdgeEvent; using enum EdgeDirection; diff --git a/tests/cpp_unit_tests/test_link_solver.cpp b/tests/cpp_unit_tests/test_link_solver.cpp index 013a7801b8..07d6fe0df1 100644 --- a/tests/cpp_unit_tests/test_link_solver.cpp +++ b/tests/cpp_unit_tests/test_link_solver.cpp @@ -107,11 +107,15 @@ TEST_CASE("Test the link solver algorithm") { auto const result = forward_elimination(edges, node_loads); REQUIRE(result.matrix.data.size() == 14); - CHECK(result.matrix.data == std::vector{1, 1, 1, 1, 1, -1, -1, 1, -1, -1, -1, -1, 1, 1}); + // TODO(figueroa1395): Some checks are commented out because the use of unordered_set in adjacency list + // makes the order of re-attachments non-deterministic (per-platform - locally it's fine), which affects the + // order of entries in the COO matrix should we use an ordered set instead to test this fine grained + // details, or is it enough to test overall structure and result correctness? CHECK(result.matrix.data == + // std::vector{1, 1, 1, 1, 1, -1, -1, 1, -1, -1, -1, -1, 1, 1}); REQUIRE(result.matrix.row.size() == 14); CHECK(result.matrix.row == std::vector{0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3}); REQUIRE(result.matrix.col.size() == 14); - CHECK(result.matrix.col == std::vector{0, 2, 1, 1, 2, 6, 3, 2, 3, 6, 5, 4, 5, 6}); + // CHECK(result.matrix.col == std::vector{0, 2, 1, 1, 2, 6, 3, 2, 3, 6, 5, 4, 5, 6}); REQUIRE(result.rhs.size() == 4); CHECK(result.rhs == std::vector{{-1.0, -1.0}, {-1.0, -1.0}, {-2.0, -2.0}, {0.0, 0.0}}); REQUIRE(result.free_edge_indices.size() == 3); From 5333b44fb2a55b438b39d44528054780c52e71c5 Mon Sep 17 00:00:00 2001 From: Santiago Figueroa Manrique Date: Wed, 11 Mar 2026 09:28:39 +0100 Subject: [PATCH 04/12] clang tidy and sonar Signed-off-by: Santiago Figueroa Manrique --- .../include/power_grid_model/link_solver.hpp | 85 ++++++++++--------- tests/cpp_unit_tests/test_link_solver.cpp | 4 + 2 files changed, 51 insertions(+), 38 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp index a8b58ba100..7047dbd780 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp @@ -9,17 +9,14 @@ #include "common/counting_iterator.hpp" #include -#include -#include -#include +#include #include #include #include +#include #include -namespace power_grid_model::link_solver { - -namespace detail { +namespace power_grid_model::link_solver::detail { constexpr uint8_t starting_row{0}; enum class EdgeEvent : std::uint8_t { @@ -54,6 +51,8 @@ struct EliminationResult { std::vector edge_history{}; // edges elimination history }; +using AdjacencyList = std::unordered_map>; + // convention: from node at position 0, to node at position 1 [[nodiscard]] constexpr Idx from_node(BranchIdx const& edge) { return edge[0]; } [[nodiscard]] constexpr Idx& from_node(BranchIdx& edge) { return edge[0]; } @@ -63,9 +62,9 @@ struct EliminationResult { // map from node index to the set of adjacent edge indices // unordered_map because node IDs may be sparse/non-contiguous // unordered_set for O(1) insert/erase during reattachment -[[nodiscard]] auto build_adjacency_list(std::span edges) { - std::unordered_map> adjacency_list{}; - for (auto const& [raw_index, edge] : enumerate(edges)) { +[[nodiscard]] inline auto build_adjacency_list(std::span edges) -> AdjacencyList { + AdjacencyList adjacency_list{}; + for (auto const& [raw_index, edge] : enumerate(std::as_const(edges))) { auto const index = static_cast(raw_index); auto const [from_node, to_node] = edge; adjacency_list[from_node].insert(index); @@ -74,15 +73,47 @@ struct EliminationResult { return adjacency_list; } -void write_edge_history(EdgeHistory& edge_history, EdgeEvent event, uint64_t row) { +inline void write_edge_history(EdgeHistory& edge_history, EdgeEvent event, uint64_t row) { edge_history.events.push_back(event); edge_history.rows.push_back(row); } +inline void re_attachment_step(Idx from_node_idx, Idx to_node_idx, uint64_t matrix_row, std::vector& edges, + AdjacencyList& adjacency_list, COOSparseMatrix& matrix, + std::vector& edge_history, + std::unordered_set& edges_contracted_to_point) { + using enum EdgeDirection; + using enum EdgeEvent; + + auto const adjacent_edges_snapshot = + std::vector(adjacency_list[to_node_idx].begin(), adjacency_list[to_node_idx].end()); + + for (uint64_t const adjacent_edge_idx : adjacent_edges_snapshot) { + if (from_node(edges[adjacent_edge_idx]) == to_node_idx) { + from_node(edges[adjacent_edge_idx]) = from_node_idx; // re-attachment step + matrix.append(std::to_underlying(Away), matrix_row, adjacent_edge_idx); + } else { + to_node(edges[adjacent_edge_idx]) = from_node_idx; // re-attachment step + matrix.append(std::to_underlying(Towards), matrix_row, adjacent_edge_idx); + } + + // update adjacency list after re-attachment + adjacency_list[to_node_idx].erase(adjacent_edge_idx); + adjacency_list[from_node_idx].insert(adjacent_edge_idx); + + if (from_node(edges[adjacent_edge_idx]) == to_node(edges[adjacent_edge_idx])) { + write_edge_history(edge_history[adjacent_edge_idx], ContractedToPoint, matrix_row); + edges_contracted_to_point.insert(adjacent_edge_idx); + } else { + write_edge_history(edge_history[adjacent_edge_idx], Replaced, matrix_row); + } + } +} + // forward elimination is performed via a modified Gaussian elimination procedure // we name this procedure elimination-game -[[nodiscard]] EliminationResult forward_elimination(std::vector edges, - std::vector node_loads) { +[[nodiscard]] inline EliminationResult forward_elimination(std::vector edges, + std::vector node_loads) { EliminationResult result{}; result.edge_history.resize(edges.size()); using enum EdgeEvent; @@ -94,7 +125,7 @@ void write_edge_history(EdgeHistory& edge_history, EdgeEvent event, uint64_t row uint64_t matrix_row{starting_row}; auto adjacency_list = build_adjacency_list(edges); - for (auto const& [raw_index, edge] : enumerate(edges)) { + for (auto const& [raw_index, edge] : enumerate(std::as_const(edges))) { auto const index = static_cast(raw_index); if (edges_contracted_to_point.contains(index)) { free_edge_indices.push_back(index); @@ -113,34 +144,12 @@ void write_edge_history(EdgeHistory& edge_history, EdgeEvent event, uint64_t row node_loads[from_node_idx] += node_loads[to_node_idx]; rhs.push_back(node_loads[to_node_idx]); - auto const adjacent_edges_snapshot = - std::vector(adjacency_list[to_node_idx].begin(), adjacency_list[to_node_idx].end()); - - for (uint64_t const adjacent_edge_idx : adjacent_edges_snapshot) { - if (from_node(edges[adjacent_edge_idx]) == to_node_idx) { - from_node(edges[adjacent_edge_idx]) = from_node_idx; // re-attachment step - matrix.append(std::to_underlying(Away), matrix_row, adjacent_edge_idx); - } else { - to_node(edges[adjacent_edge_idx]) = from_node_idx; // re-attachment step - matrix.append(std::to_underlying(Towards), matrix_row, adjacent_edge_idx); - } - - // update adjacency list after re-attachment - adjacency_list[to_node_idx].erase(adjacent_edge_idx); - adjacency_list[from_node_idx].insert(adjacent_edge_idx); - - if (from_node(edges[adjacent_edge_idx]) == to_node(edges[adjacent_edge_idx])) { - write_edge_history(edge_history[adjacent_edge_idx], ContractedToPoint, matrix_row); - edges_contracted_to_point.insert(adjacent_edge_idx); - } else { - write_edge_history(edge_history[adjacent_edge_idx], Replaced, matrix_row); - } - } + re_attachment_step(from_node_idx, to_node_idx, matrix_row, edges, adjacency_list, matrix, edge_history, + edges_contracted_to_point); ++matrix_row; } } return result; } -} // namespace detail -} // namespace power_grid_model::link_solver +} // namespace power_grid_model::link_solver::detail diff --git a/tests/cpp_unit_tests/test_link_solver.cpp b/tests/cpp_unit_tests/test_link_solver.cpp index 07d6fe0df1..19b5d405a9 100644 --- a/tests/cpp_unit_tests/test_link_solver.cpp +++ b/tests/cpp_unit_tests/test_link_solver.cpp @@ -2,10 +2,14 @@ // // SPDX-License-Identifier: MPL-2.0 +#include #include #include +#include +#include + namespace power_grid_model::link_solver { TEST_CASE("Test the link solver algorithm") { using namespace detail; From 49dae7c002e7112f89db18f81c4a25f7fd8e1d3c Mon Sep 17 00:00:00 2001 From: Santiago Figueroa Manrique Date: Wed, 11 Mar 2026 14:11:43 +0100 Subject: [PATCH 05/12] test adjacency list building Signed-off-by: Santiago Figueroa Manrique --- tests/cpp_unit_tests/test_link_solver.cpp | 57 ++++++++++++++++++++++- 1 file changed, 55 insertions(+), 2 deletions(-) diff --git a/tests/cpp_unit_tests/test_link_solver.cpp b/tests/cpp_unit_tests/test_link_solver.cpp index 19b5d405a9..7133ff56c1 100644 --- a/tests/cpp_unit_tests/test_link_solver.cpp +++ b/tests/cpp_unit_tests/test_link_solver.cpp @@ -2,6 +2,7 @@ // // SPDX-License-Identifier: MPL-2.0 +#include #include #include @@ -14,6 +15,59 @@ namespace power_grid_model::link_solver { TEST_CASE("Test the link solver algorithm") { using namespace detail; + SUBCASE("Test build adjacency list") { + + SUBCASE("One edge, two nodes") { + auto edges = std::vector{{0, 1}}; + AdjacencyList const adjacency_list = build_adjacency_list(edges); + + REQUIRE(adjacency_list.size() == 2); + CHECK(adjacency_list.at(0) == std::unordered_set{0}); + CHECK(adjacency_list.at(1) == std::unordered_set{0}); + } + + SUBCASE("Two edges, three nodes") { + auto edges = std::vector{{1, 0}, {1, 2}}; + AdjacencyList const adjacency_list = build_adjacency_list(edges); + + REQUIRE(adjacency_list.size() == 3); + CHECK(adjacency_list.at(0) == std::unordered_set{0}); + CHECK(adjacency_list.at(1) == std::unordered_set{0, 1}); + CHECK(adjacency_list.at(2) == std::unordered_set{1}); + } + + SUBCASE("Three edges, three nodes") { + auto edges = std::vector{{0, 1}, {1, 2}, {2, 0}}; + AdjacencyList const adjacency_list = build_adjacency_list(edges); + + REQUIRE(adjacency_list.size() == 3); + CHECK(adjacency_list.at(0) == std::unordered_set{0, 2}); + CHECK(adjacency_list.at(1) == std::unordered_set{0, 1}); + CHECK(adjacency_list.at(2) == std::unordered_set{1, 2}); + } + + SUBCASE("Two edges, two nodes") { + auto edges = std::vector{{0, 1}, {0, 1}}; + AdjacencyList const adjacency_list = build_adjacency_list(edges); + + REQUIRE(adjacency_list.size() == 2); + CHECK(adjacency_list.at(0) == std::unordered_set{0, 1}); + CHECK(adjacency_list.at(1) == std::unordered_set{0, 1}); + } + + SUBCASE("Seven edges, five nodes") { + auto edges = std::vector{{3, 0}, {1, 0}, {2, 0}, {3, 2}, {1, 2}, {1, 4}, {3, 4}}; + AdjacencyList const adjacency_list = build_adjacency_list(edges); + + REQUIRE(adjacency_list.size() == 5); + CHECK(adjacency_list.at(0) == std::unordered_set{0, 1, 2}); + CHECK(adjacency_list.at(1) == std::unordered_set{1, 4, 5}); + CHECK(adjacency_list.at(2) == std::unordered_set{2, 3, 4}); + CHECK(adjacency_list.at(3) == std::unordered_set{0, 3, 6}); + CHECK(adjacency_list.at(4) == std::unordered_set{5, 6}); + } + } + SUBCASE("Test forward elimination - elimination game") { using enum EdgeEvent; @@ -32,8 +86,7 @@ TEST_CASE("Test the link solver algorithm") { CHECK(result.rhs == std::vector{{1.0, 0.0}}); CHECK(result.free_edge_indices.empty()); REQUIRE(result.edge_history.size() == 1); - CHECK(result.edge_history[0].events == - std::vector{Deleted}); // Should this be deleted or contracted to point? + CHECK(result.edge_history[0].events == std::vector{Deleted}); CHECK(result.edge_history[0].rows == std::vector{0}); } From e0846b1e9cf2588703af7b4a020de19a62740770 Mon Sep 17 00:00:00 2001 From: Santiago Figueroa Manrique Date: Wed, 11 Mar 2026 15:42:43 +0100 Subject: [PATCH 06/12] address comments Signed-off-by: Santiago Figueroa Manrique --- .../include/power_grid_model/link_solver.hpp | 82 ++++++++++--------- tests/cpp_unit_tests/test_link_solver.cpp | 80 +++++++++--------- 2 files changed, 84 insertions(+), 78 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp index 7047dbd780..ac35681a11 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp @@ -19,6 +19,9 @@ namespace power_grid_model::link_solver::detail { constexpr uint8_t starting_row{0}; +using AdjacencyMap = std::unordered_map>; +using ContractedEdgesSet = std::unordered_set; + enum class EdgeEvent : std::uint8_t { Deleted = 0, // pivot edge - used as pivot row Replaced = 1, // reattached to a different node @@ -48,11 +51,9 @@ struct EliminationResult { COOSparseMatrix matrix{}; std::vector rhs{}; // RHS value at each pivot row std::vector free_edge_indices{}; // index of degrees of freedom (self loop edges) - std::vector edge_history{}; // edges elimination history + std::vector edges_history{}; // edges elimination history }; -using AdjacencyList = std::unordered_map>; - // convention: from node at position 0, to node at position 1 [[nodiscard]] constexpr Idx from_node(BranchIdx const& edge) { return edge[0]; } [[nodiscard]] constexpr Idx& from_node(BranchIdx& edge) { return edge[0]; } @@ -62,8 +63,8 @@ using AdjacencyList = std::unordered_map>; // map from node index to the set of adjacent edge indices // unordered_map because node IDs may be sparse/non-contiguous // unordered_set for O(1) insert/erase during reattachment -[[nodiscard]] inline auto build_adjacency_list(std::span edges) -> AdjacencyList { - AdjacencyList adjacency_list{}; +[[nodiscard]] inline auto build_adjacency_list(std::span edges) -> AdjacencyMap { + AdjacencyMap adjacency_list{}; for (auto const& [raw_index, edge] : enumerate(std::as_const(edges))) { auto const index = static_cast(raw_index); auto const [from_node, to_node] = edge; @@ -73,40 +74,33 @@ using AdjacencyList = std::unordered_map>; return adjacency_list; } -inline void write_edge_history(EdgeHistory& edge_history, EdgeEvent event, uint64_t row) { - edge_history.events.push_back(event); - edge_history.rows.push_back(row); +inline void write_edge_history(EdgeHistory& edges_history, EdgeEvent event, uint64_t row) { + edges_history.events.push_back(event); + edges_history.rows.push_back(row); } -inline void re_attachment_step(Idx from_node_idx, Idx to_node_idx, uint64_t matrix_row, std::vector& edges, - AdjacencyList& adjacency_list, COOSparseMatrix& matrix, - std::vector& edge_history, - std::unordered_set& edges_contracted_to_point) { +inline void replace_and_write(uint64_t edge_idx, Idx from_node_idx, Idx to_node_idx, uint64_t matrix_row, + std::vector& edges, COOSparseMatrix& matrix) { using enum EdgeDirection; - using enum EdgeEvent; - - auto const adjacent_edges_snapshot = - std::vector(adjacency_list[to_node_idx].begin(), adjacency_list[to_node_idx].end()); - for (uint64_t const adjacent_edge_idx : adjacent_edges_snapshot) { - if (from_node(edges[adjacent_edge_idx]) == to_node_idx) { - from_node(edges[adjacent_edge_idx]) = from_node_idx; // re-attachment step - matrix.append(std::to_underlying(Away), matrix_row, adjacent_edge_idx); - } else { - to_node(edges[adjacent_edge_idx]) = from_node_idx; // re-attachment step - matrix.append(std::to_underlying(Towards), matrix_row, adjacent_edge_idx); - } + if (from_node(edges[edge_idx]) == to_node_idx) { + from_node(edges[edge_idx]) = from_node_idx; // re-attachment step + matrix.append(std::to_underlying(Away), matrix_row, edge_idx); + } else { + to_node(edges[edge_idx]) = from_node_idx; // re-attachment step + matrix.append(std::to_underlying(Towards), matrix_row, edge_idx); + } +} - // update adjacency list after re-attachment - adjacency_list[to_node_idx].erase(adjacent_edge_idx); - adjacency_list[from_node_idx].insert(adjacent_edge_idx); +inline void update_edge_info(uint64_t edge_idx, uint64_t matrix_row, std::vector& edges, + std::vector& edges_history, ContractedEdgesSet& edges_contracted_to_point) { + using enum EdgeEvent; - if (from_node(edges[adjacent_edge_idx]) == to_node(edges[adjacent_edge_idx])) { - write_edge_history(edge_history[adjacent_edge_idx], ContractedToPoint, matrix_row); - edges_contracted_to_point.insert(adjacent_edge_idx); - } else { - write_edge_history(edge_history[adjacent_edge_idx], Replaced, matrix_row); - } + if (from_node(edges[edge_idx]) == to_node(edges[edge_idx])) { + write_edge_history(edges_history[edge_idx], ContractedToPoint, matrix_row); + edges_contracted_to_point.insert(edge_idx); + } else { + write_edge_history(edges_history[edge_idx], Replaced, matrix_row); } } @@ -115,11 +109,11 @@ inline void re_attachment_step(Idx from_node_idx, Idx to_node_idx, uint64_t matr [[nodiscard]] inline EliminationResult forward_elimination(std::vector edges, std::vector node_loads) { EliminationResult result{}; - result.edge_history.resize(edges.size()); + result.edges_history.resize(edges.size()); using enum EdgeEvent; using enum EdgeDirection; - auto& [matrix, rhs, free_edge_indices, edge_history] = result; + auto& [matrix, rhs, free_edge_indices, edges_history] = result; std::unordered_set edges_contracted_to_point{}; uint64_t matrix_row{starting_row}; @@ -130,7 +124,7 @@ inline void re_attachment_step(Idx from_node_idx, Idx to_node_idx, uint64_t matr if (edges_contracted_to_point.contains(index)) { free_edge_indices.push_back(index); } else { - write_edge_history(edge_history[index], Deleted, matrix_row); // Delete edge -> pivot there + write_edge_history(edges_history[index], Deleted, matrix_row); // Delete edge -> pivot there matrix.append(std::to_underlying(Towards), matrix_row, index); Idx const from_node_idx = from_node(edge); @@ -144,8 +138,20 @@ inline void re_attachment_step(Idx from_node_idx, Idx to_node_idx, uint64_t matr node_loads[from_node_idx] += node_loads[to_node_idx]; rhs.push_back(node_loads[to_node_idx]); - re_attachment_step(from_node_idx, to_node_idx, matrix_row, edges, adjacency_list, matrix, edge_history, - edges_contracted_to_point); + auto const adjacent_edges_snapshot = + std::vector(adjacency_list[to_node_idx].begin(), adjacency_list[to_node_idx].end()); + + for (uint64_t const adjacent_edge_idx : adjacent_edges_snapshot) { + // re-attach edge and write to matrix + replace_and_write(adjacent_edge_idx, from_node_idx, to_node_idx, matrix_row, edges, matrix); + + // update adjacency list after re-attachment + adjacency_list[to_node_idx].erase(adjacent_edge_idx); + adjacency_list[from_node_idx].insert(adjacent_edge_idx); + + // update edges_history and edges_contracted_to_point (if needed) after re-attachment + update_edge_info(adjacent_edge_idx, matrix_row, edges, edges_history, edges_contracted_to_point); + } ++matrix_row; } } diff --git a/tests/cpp_unit_tests/test_link_solver.cpp b/tests/cpp_unit_tests/test_link_solver.cpp index 7133ff56c1..0e0dad7e71 100644 --- a/tests/cpp_unit_tests/test_link_solver.cpp +++ b/tests/cpp_unit_tests/test_link_solver.cpp @@ -19,7 +19,7 @@ TEST_CASE("Test the link solver algorithm") { SUBCASE("One edge, two nodes") { auto edges = std::vector{{0, 1}}; - AdjacencyList const adjacency_list = build_adjacency_list(edges); + AdjacencyMap const adjacency_list = build_adjacency_list(edges); REQUIRE(adjacency_list.size() == 2); CHECK(adjacency_list.at(0) == std::unordered_set{0}); @@ -28,7 +28,7 @@ TEST_CASE("Test the link solver algorithm") { SUBCASE("Two edges, three nodes") { auto edges = std::vector{{1, 0}, {1, 2}}; - AdjacencyList const adjacency_list = build_adjacency_list(edges); + AdjacencyMap const adjacency_list = build_adjacency_list(edges); REQUIRE(adjacency_list.size() == 3); CHECK(adjacency_list.at(0) == std::unordered_set{0}); @@ -38,7 +38,7 @@ TEST_CASE("Test the link solver algorithm") { SUBCASE("Three edges, three nodes") { auto edges = std::vector{{0, 1}, {1, 2}, {2, 0}}; - AdjacencyList const adjacency_list = build_adjacency_list(edges); + AdjacencyMap const adjacency_list = build_adjacency_list(edges); REQUIRE(adjacency_list.size() == 3); CHECK(adjacency_list.at(0) == std::unordered_set{0, 2}); @@ -48,7 +48,7 @@ TEST_CASE("Test the link solver algorithm") { SUBCASE("Two edges, two nodes") { auto edges = std::vector{{0, 1}, {0, 1}}; - AdjacencyList const adjacency_list = build_adjacency_list(edges); + AdjacencyMap const adjacency_list = build_adjacency_list(edges); REQUIRE(adjacency_list.size() == 2); CHECK(adjacency_list.at(0) == std::unordered_set{0, 1}); @@ -57,7 +57,7 @@ TEST_CASE("Test the link solver algorithm") { SUBCASE("Seven edges, five nodes") { auto edges = std::vector{{3, 0}, {1, 0}, {2, 0}, {3, 2}, {1, 2}, {1, 4}, {3, 4}}; - AdjacencyList const adjacency_list = build_adjacency_list(edges); + AdjacencyMap const adjacency_list = build_adjacency_list(edges); REQUIRE(adjacency_list.size() == 5); CHECK(adjacency_list.at(0) == std::unordered_set{0, 1, 2}); @@ -85,9 +85,9 @@ TEST_CASE("Test the link solver algorithm") { REQUIRE(result.rhs.size() == 1); CHECK(result.rhs == std::vector{{1.0, 0.0}}); CHECK(result.free_edge_indices.empty()); - REQUIRE(result.edge_history.size() == 1); - CHECK(result.edge_history[0].events == std::vector{Deleted}); - CHECK(result.edge_history[0].rows == std::vector{0}); + REQUIRE(result.edges_history.size() == 1); + CHECK(result.edges_history[0].events == std::vector{Deleted}); + CHECK(result.edges_history[0].rows == std::vector{0}); } SUBCASE("Two edges, three nodes, two real loads") { @@ -104,11 +104,11 @@ TEST_CASE("Test the link solver algorithm") { REQUIRE(result.rhs.size() == 2); CHECK(result.rhs == std::vector{{-1.0, 0.0}, {0.0, 0.0}}); CHECK(result.free_edge_indices.empty()); - REQUIRE(result.edge_history.size() == 2); - CHECK(result.edge_history[0].events == std::vector{Deleted}); - CHECK(result.edge_history[0].rows == std::vector{0}); - CHECK(result.edge_history[1].events == std::vector{Deleted}); - CHECK(result.edge_history[1].rows == std::vector{1}); + REQUIRE(result.edges_history.size() == 2); + CHECK(result.edges_history[0].events == std::vector{Deleted}); + CHECK(result.edges_history[0].rows == std::vector{0}); + CHECK(result.edges_history[1].events == std::vector{Deleted}); + CHECK(result.edges_history[1].rows == std::vector{1}); } SUBCASE("Three edges, three nodes, two real loads") { @@ -126,13 +126,13 @@ TEST_CASE("Test the link solver algorithm") { CHECK(result.rhs == std::vector{{1.0, 0.0}, {0.0, 0.0}}); REQUIRE(result.free_edge_indices.size() == 1); CHECK(result.free_edge_indices == std::vector{2}); - REQUIRE(result.edge_history.size() == 3); - CHECK(result.edge_history[0].events == std::vector{Deleted}); - CHECK(result.edge_history[0].rows == std::vector{0}); - CHECK(result.edge_history[1].events == std::vector{Replaced, Deleted}); - CHECK(result.edge_history[1].rows == std::vector{0, 1}); - CHECK(result.edge_history[2].events == std::vector{ContractedToPoint}); - CHECK(result.edge_history[2].rows == std::vector{1}); + REQUIRE(result.edges_history.size() == 3); + CHECK(result.edges_history[0].events == std::vector{Deleted}); + CHECK(result.edges_history[0].rows == std::vector{0}); + CHECK(result.edges_history[1].events == std::vector{Replaced, Deleted}); + CHECK(result.edges_history[1].rows == std::vector{0, 1}); + CHECK(result.edges_history[2].events == std::vector{ContractedToPoint}); + CHECK(result.edges_history[2].rows == std::vector{1}); } SUBCASE("Two edges, two nodes, two real loads") { @@ -150,11 +150,11 @@ TEST_CASE("Test the link solver algorithm") { CHECK(result.rhs == std::vector{{1.0, 0.0}}); REQUIRE(result.free_edge_indices.size() == 1); CHECK(result.free_edge_indices == std::vector{1}); - REQUIRE(result.edge_history.size() == 2); - CHECK(result.edge_history[0].events == std::vector{Deleted}); - CHECK(result.edge_history[0].rows == std::vector{0}); - CHECK(result.edge_history[1].events == std::vector{ContractedToPoint}); - CHECK(result.edge_history[1].rows == std::vector{0}); + REQUIRE(result.edges_history.size() == 2); + CHECK(result.edges_history[0].events == std::vector{Deleted}); + CHECK(result.edges_history[0].rows == std::vector{0}); + CHECK(result.edges_history[1].events == std::vector{ContractedToPoint}); + CHECK(result.edges_history[1].rows == std::vector{0}); } SUBCASE("Complex case with complex loads") { @@ -177,21 +177,21 @@ TEST_CASE("Test the link solver algorithm") { CHECK(result.rhs == std::vector{{-1.0, -1.0}, {-1.0, -1.0}, {-2.0, -2.0}, {0.0, 0.0}}); REQUIRE(result.free_edge_indices.size() == 3); CHECK(result.free_edge_indices == std::vector{3, 4, 6}); - REQUIRE(result.edge_history.size() == 7); - CHECK(result.edge_history[0].events == std::vector{Deleted}); - CHECK(result.edge_history[0].rows == std::vector{0}); - CHECK(result.edge_history[1].events == std::vector{Replaced, Deleted}); - CHECK(result.edge_history[1].rows == std::vector{0, 1}); - CHECK(result.edge_history[2].events == std::vector{Replaced, Replaced, Deleted}); - CHECK(result.edge_history[2].rows == std::vector{0, 1, 2}); - CHECK(result.edge_history[3].events == std::vector{Replaced, ContractedToPoint}); - CHECK(result.edge_history[3].rows == std::vector{1, 2}); - CHECK(result.edge_history[4].events == std::vector{ContractedToPoint}); - CHECK(result.edge_history[4].rows == std::vector{2}); - CHECK(result.edge_history[5].events == std::vector{Replaced, Deleted}); - CHECK(result.edge_history[5].rows == std::vector{2, 3}); - CHECK(result.edge_history[6].events == std::vector{Replaced, Replaced, ContractedToPoint}); - CHECK(result.edge_history[6].rows == std::vector{1, 2, 3}); + REQUIRE(result.edges_history.size() == 7); + CHECK(result.edges_history[0].events == std::vector{Deleted}); + CHECK(result.edges_history[0].rows == std::vector{0}); + CHECK(result.edges_history[1].events == std::vector{Replaced, Deleted}); + CHECK(result.edges_history[1].rows == std::vector{0, 1}); + CHECK(result.edges_history[2].events == std::vector{Replaced, Replaced, Deleted}); + CHECK(result.edges_history[2].rows == std::vector{0, 1, 2}); + CHECK(result.edges_history[3].events == std::vector{Replaced, ContractedToPoint}); + CHECK(result.edges_history[3].rows == std::vector{1, 2}); + CHECK(result.edges_history[4].events == std::vector{ContractedToPoint}); + CHECK(result.edges_history[4].rows == std::vector{2}); + CHECK(result.edges_history[5].events == std::vector{Replaced, Deleted}); + CHECK(result.edges_history[5].rows == std::vector{2, 3}); + CHECK(result.edges_history[6].events == std::vector{Replaced, Replaced, ContractedToPoint}); + CHECK(result.edges_history[6].rows == std::vector{1, 2, 3}); } } } From 9992a9e0a2f735dd1c01b1bdcccc77a3bed43b8a Mon Sep 17 00:00:00 2001 From: Santiago Figueroa Manrique Date: Mon, 16 Mar 2026 16:10:22 +0100 Subject: [PATCH 07/12] changed COO sparse matrix struct Signed-off-by: Santiago Figueroa Manrique --- .../include/power_grid_model/link_solver.hpp | 70 ++++++---- tests/cpp_unit_tests/test_link_solver.cpp | 129 +++++++++--------- 2 files changed, 111 insertions(+), 88 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp index ac35681a11..15087f7aef 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp @@ -10,6 +10,7 @@ #include #include +#include #include #include #include @@ -17,7 +18,6 @@ #include namespace power_grid_model::link_solver::detail { -constexpr uint8_t starting_row{0}; using AdjacencyMap = std::unordered_map>; using ContractedEdgesSet = std::unordered_set; @@ -31,14 +31,27 @@ enum class EdgeEvent : std::uint8_t { enum class EdgeDirection : IntS { Away = -1, Towards = 1 }; struct COOSparseMatrix { - std::vector data{}; - std::vector row{}; - std::vector col{}; - - void append(IntS value, uint64_t row_idx, uint64_t col_idx) { - data.push_back(value); - row.push_back(row_idx); - col.push_back(col_idx); + uint64_t row_number{}; + std::unordered_map data_map{}; + + void prepare(uint64_t row_size, uint64_t col_size) { + row_number = col_size; + data_map.reserve(row_size * + col_size); // worst case to guarantee O(1) insertion/retrieval, but typically much sparser + } + void set_value(IntS value, uint64_t row_idx, uint64_t col_idx) { + assert(row_number != 0 && "row_number must be set before setting values in data_map"); + data_map[row_idx * row_number + col_idx] = value; + } + // TODO(figueroa1395): maybe this is overkill - revisit after implementation is further ahead + [[nodiscard]] bool get_value(IntS& value, uint64_t row_idx, uint64_t col_idx) const { + assert(row_number != 0 && "row_number must be set before getting values from data_map"); + auto const it = data_map.find(row_idx * row_number + col_idx); + if (it != data_map.end()) { + value = it->second; + return true; + } + return false; } }; @@ -63,15 +76,15 @@ struct EliminationResult { // map from node index to the set of adjacent edge indices // unordered_map because node IDs may be sparse/non-contiguous // unordered_set for O(1) insert/erase during reattachment -[[nodiscard]] inline auto build_adjacency_list(std::span edges) -> AdjacencyMap { - AdjacencyMap adjacency_list{}; +[[nodiscard]] inline auto build_adjacency_map(std::span edges) -> AdjacencyMap { + AdjacencyMap adjacency_map{}; for (auto const& [raw_index, edge] : enumerate(std::as_const(edges))) { auto const index = static_cast(raw_index); auto const [from_node, to_node] = edge; - adjacency_list[from_node].insert(index); - adjacency_list[to_node].insert(index); + adjacency_map[from_node].insert(index); + adjacency_map[to_node].insert(index); } - return adjacency_list; + return adjacency_map; } inline void write_edge_history(EdgeHistory& edges_history, EdgeEvent event, uint64_t row) { @@ -85,10 +98,10 @@ inline void replace_and_write(uint64_t edge_idx, Idx from_node_idx, Idx to_node_ if (from_node(edges[edge_idx]) == to_node_idx) { from_node(edges[edge_idx]) = from_node_idx; // re-attachment step - matrix.append(std::to_underlying(Away), matrix_row, edge_idx); + matrix.set_value(std::to_underlying(Away), matrix_row, edge_idx); } else { to_node(edges[edge_idx]) = from_node_idx; // re-attachment step - matrix.append(std::to_underlying(Towards), matrix_row, edge_idx); + matrix.set_value(std::to_underlying(Towards), matrix_row, edge_idx); } } @@ -106,18 +119,24 @@ inline void update_edge_info(uint64_t edge_idx, uint64_t matrix_row, std::vector // forward elimination is performed via a modified Gaussian elimination procedure // we name this procedure elimination-game +// node_loads convention: caller passes negated external loads (as in RHS of power flow equations)) [[nodiscard]] inline EliminationResult forward_elimination(std::vector edges, std::vector node_loads) { EliminationResult result{}; - result.edges_history.resize(edges.size()); + uint64_t edge_number{edges.size()}; + uint64_t node_number{node_loads.size()}; + + result.edges_history.resize(edge_number); using enum EdgeEvent; using enum EdgeDirection; auto& [matrix, rhs, free_edge_indices, edges_history] = result; + matrix.prepare(edge_number, node_number); - std::unordered_set edges_contracted_to_point{}; + constexpr uint8_t starting_row{}; uint64_t matrix_row{starting_row}; - auto adjacency_list = build_adjacency_list(edges); + auto adjacency_map = build_adjacency_map(edges); + ContractedEdgesSet edges_contracted_to_point{}; // TODO(figueroa1395): is this really needed? for (auto const& [raw_index, edge] : enumerate(std::as_const(edges))) { auto const index = static_cast(raw_index); @@ -125,29 +144,29 @@ inline void update_edge_info(uint64_t edge_idx, uint64_t matrix_row, std::vector free_edge_indices.push_back(index); } else { write_edge_history(edges_history[index], Deleted, matrix_row); // Delete edge -> pivot there - matrix.append(std::to_underlying(Towards), matrix_row, index); + matrix.set_value(std::to_underlying(Towards), matrix_row, index); Idx const from_node_idx = from_node(edge); Idx const to_node_idx = to_node(edge); // update adjacency list for deleted edge - adjacency_list[from_node_idx].erase(index); - adjacency_list[to_node_idx].erase(index); + adjacency_map[from_node_idx].erase(index); + adjacency_map[to_node_idx].erase(index); // Gaussian elimination like steps node_loads[from_node_idx] += node_loads[to_node_idx]; rhs.push_back(node_loads[to_node_idx]); auto const adjacent_edges_snapshot = - std::vector(adjacency_list[to_node_idx].begin(), adjacency_list[to_node_idx].end()); + std::vector(adjacency_map[to_node_idx].begin(), adjacency_map[to_node_idx].end()); for (uint64_t const adjacent_edge_idx : adjacent_edges_snapshot) { // re-attach edge and write to matrix replace_and_write(adjacent_edge_idx, from_node_idx, to_node_idx, matrix_row, edges, matrix); // update adjacency list after re-attachment - adjacency_list[to_node_idx].erase(adjacent_edge_idx); - adjacency_list[from_node_idx].insert(adjacent_edge_idx); + adjacency_map[to_node_idx].erase(adjacent_edge_idx); + adjacency_map[from_node_idx].insert(adjacent_edge_idx); // update edges_history and edges_contracted_to_point (if needed) after re-attachment update_edge_info(adjacent_edge_idx, matrix_row, edges, edges_history, edges_contracted_to_point); @@ -155,7 +174,6 @@ inline void update_edge_info(uint64_t edge_idx, uint64_t matrix_row, std::vector ++matrix_row; } } - return result; } } // namespace power_grid_model::link_solver::detail diff --git a/tests/cpp_unit_tests/test_link_solver.cpp b/tests/cpp_unit_tests/test_link_solver.cpp index 0e0dad7e71..a1ca651f71 100644 --- a/tests/cpp_unit_tests/test_link_solver.cpp +++ b/tests/cpp_unit_tests/test_link_solver.cpp @@ -9,9 +9,16 @@ #include #include +#include +#include #include namespace power_grid_model::link_solver { +namespace { +[[nodiscard]] IntS data_map_key_gen(uint64_t row_idx, uint64_t col_idx, uint64_t row_number) { + return row_idx * row_number + col_idx; +} +} // namespace TEST_CASE("Test the link solver algorithm") { using namespace detail; @@ -19,52 +26,52 @@ TEST_CASE("Test the link solver algorithm") { SUBCASE("One edge, two nodes") { auto edges = std::vector{{0, 1}}; - AdjacencyMap const adjacency_list = build_adjacency_list(edges); + AdjacencyMap const adjacency_map = build_adjacency_map(edges); - REQUIRE(adjacency_list.size() == 2); - CHECK(adjacency_list.at(0) == std::unordered_set{0}); - CHECK(adjacency_list.at(1) == std::unordered_set{0}); + REQUIRE(adjacency_map.size() == 2); + CHECK(adjacency_map.at(0) == std::unordered_set{0}); + CHECK(adjacency_map.at(1) == std::unordered_set{0}); } SUBCASE("Two edges, three nodes") { auto edges = std::vector{{1, 0}, {1, 2}}; - AdjacencyMap const adjacency_list = build_adjacency_list(edges); + AdjacencyMap const adjacency_map = build_adjacency_map(edges); - REQUIRE(adjacency_list.size() == 3); - CHECK(adjacency_list.at(0) == std::unordered_set{0}); - CHECK(adjacency_list.at(1) == std::unordered_set{0, 1}); - CHECK(adjacency_list.at(2) == std::unordered_set{1}); + REQUIRE(adjacency_map.size() == 3); + CHECK(adjacency_map.at(0) == std::unordered_set{0}); + CHECK(adjacency_map.at(1) == std::unordered_set{0, 1}); + CHECK(adjacency_map.at(2) == std::unordered_set{1}); } SUBCASE("Three edges, three nodes") { auto edges = std::vector{{0, 1}, {1, 2}, {2, 0}}; - AdjacencyMap const adjacency_list = build_adjacency_list(edges); + AdjacencyMap const adjacency_map = build_adjacency_map(edges); - REQUIRE(adjacency_list.size() == 3); - CHECK(adjacency_list.at(0) == std::unordered_set{0, 2}); - CHECK(adjacency_list.at(1) == std::unordered_set{0, 1}); - CHECK(adjacency_list.at(2) == std::unordered_set{1, 2}); + REQUIRE(adjacency_map.size() == 3); + CHECK(adjacency_map.at(0) == std::unordered_set{0, 2}); + CHECK(adjacency_map.at(1) == std::unordered_set{0, 1}); + CHECK(adjacency_map.at(2) == std::unordered_set{1, 2}); } SUBCASE("Two edges, two nodes") { auto edges = std::vector{{0, 1}, {0, 1}}; - AdjacencyMap const adjacency_list = build_adjacency_list(edges); + AdjacencyMap const adjacency_map = build_adjacency_map(edges); - REQUIRE(adjacency_list.size() == 2); - CHECK(adjacency_list.at(0) == std::unordered_set{0, 1}); - CHECK(adjacency_list.at(1) == std::unordered_set{0, 1}); + REQUIRE(adjacency_map.size() == 2); + CHECK(adjacency_map.at(0) == std::unordered_set{0, 1}); + CHECK(adjacency_map.at(1) == std::unordered_set{0, 1}); } SUBCASE("Seven edges, five nodes") { auto edges = std::vector{{3, 0}, {1, 0}, {2, 0}, {3, 2}, {1, 2}, {1, 4}, {3, 4}}; - AdjacencyMap const adjacency_list = build_adjacency_list(edges); - - REQUIRE(adjacency_list.size() == 5); - CHECK(adjacency_list.at(0) == std::unordered_set{0, 1, 2}); - CHECK(adjacency_list.at(1) == std::unordered_set{1, 4, 5}); - CHECK(adjacency_list.at(2) == std::unordered_set{2, 3, 4}); - CHECK(adjacency_list.at(3) == std::unordered_set{0, 3, 6}); - CHECK(adjacency_list.at(4) == std::unordered_set{5, 6}); + AdjacencyMap const adjacency_map = build_adjacency_map(edges); + + REQUIRE(adjacency_map.size() == 5); + CHECK(adjacency_map.at(0) == std::unordered_set{0, 1, 2}); + CHECK(adjacency_map.at(1) == std::unordered_set{1, 4, 5}); + CHECK(adjacency_map.at(2) == std::unordered_set{2, 3, 4}); + CHECK(adjacency_map.at(3) == std::unordered_set{0, 3, 6}); + CHECK(adjacency_map.at(4) == std::unordered_set{5, 6}); } } @@ -75,14 +82,10 @@ TEST_CASE("Test the link solver algorithm") { auto const edges = std::vector{{0, 1}}; auto const node_loads = std::vector{{-1.0, 0.0}, {1.0, 0.0}}; auto const result = forward_elimination(edges, node_loads); + uint64_t const node_number = node_loads.size(); - REQUIRE(result.matrix.data.size() == 1); - CHECK(result.matrix.data == std::vector{1}); - REQUIRE(result.matrix.row.size() == 1); - CHECK(result.matrix.row == std::vector{0}); - REQUIRE(result.matrix.col.size() == 1); - CHECK(result.matrix.col == std::vector{0}); - REQUIRE(result.rhs.size() == 1); + REQUIRE(result.matrix.data_map.size() == 1); + CHECK(result.matrix.data_map.at(data_map_key_gen(0, 0, node_number)) == 1); CHECK(result.rhs == std::vector{{1.0, 0.0}}); CHECK(result.free_edge_indices.empty()); REQUIRE(result.edges_history.size() == 1); @@ -94,13 +97,11 @@ TEST_CASE("Test the link solver algorithm") { auto const edges = std::vector{{1, 0}, {1, 2}}; auto const node_loads = std::vector{{-1.0, 0.0}, {1.0, 0.0}, {0.0, 0.0}}; auto const result = forward_elimination(edges, node_loads); + uint64_t const node_number = node_loads.size(); - REQUIRE(result.matrix.data.size() == 2); - CHECK(result.matrix.data == std::vector{1, 1}); - REQUIRE(result.matrix.row.size() == 2); - CHECK(result.matrix.row == std::vector{0, 1}); - REQUIRE(result.matrix.col.size() == 2); - CHECK(result.matrix.col == std::vector{0, 1}); + REQUIRE(result.matrix.data_map.size() == 2); + CHECK(result.matrix.data_map.at(data_map_key_gen(0, 0, node_number)) == 1); + CHECK(result.matrix.data_map.at(data_map_key_gen(1, 1, node_number)) == 1); REQUIRE(result.rhs.size() == 2); CHECK(result.rhs == std::vector{{-1.0, 0.0}, {0.0, 0.0}}); CHECK(result.free_edge_indices.empty()); @@ -115,13 +116,13 @@ TEST_CASE("Test the link solver algorithm") { auto const edges = std::vector{{0, 1}, {1, 2}, {2, 0}}; auto const node_loads = std::vector{{-1.0, 0.0}, {1.0, 0.0}, {0.0, 0.0}}; auto const result = forward_elimination(edges, node_loads); + uint64_t const node_number = node_loads.size(); - REQUIRE(result.matrix.data.size() == 4); - CHECK(result.matrix.data == std::vector{1, -1, 1, -1}); - REQUIRE(result.matrix.row.size() == 4); - CHECK(result.matrix.row == std::vector{0, 0, 1, 1}); - REQUIRE(result.matrix.col.size() == 4); - CHECK(result.matrix.col == std::vector{0, 1, 1, 2}); + REQUIRE(result.matrix.data_map.size() == 4); + CHECK(result.matrix.data_map.at(data_map_key_gen(0, 0, node_number)) == 1); + CHECK(result.matrix.data_map.at(data_map_key_gen(0, 1, node_number)) == -1); + CHECK(result.matrix.data_map.at(data_map_key_gen(1, 1, node_number)) == 1); + CHECK(result.matrix.data_map.at(data_map_key_gen(1, 2, node_number)) == -1); REQUIRE(result.rhs.size() == 2); CHECK(result.rhs == std::vector{{1.0, 0.0}, {0.0, 0.0}}); REQUIRE(result.free_edge_indices.size() == 1); @@ -139,13 +140,11 @@ TEST_CASE("Test the link solver algorithm") { auto const edges = std::vector{{0, 1}, {0, 1}}; auto const node_loads = std::vector{{-1.0, 0.0}, {1.0, 0.0}}; auto const result = forward_elimination(edges, node_loads); + uint64_t const node_number = node_loads.size(); - REQUIRE(result.matrix.data.size() == 2); - CHECK(result.matrix.data == std::vector{1, 1}); - REQUIRE(result.matrix.row.size() == 2); - CHECK(result.matrix.row == std::vector{0, 0}); - REQUIRE(result.matrix.col.size() == 2); - CHECK(result.matrix.col == std::vector{0, 1}); + REQUIRE(result.matrix.data_map.size() == 2); + CHECK(result.matrix.data_map.at(data_map_key_gen(0, 0, node_number)) == 1); + CHECK(result.matrix.data_map.at(data_map_key_gen(0, 1, node_number)) == 1); REQUIRE(result.rhs.size() == 1); CHECK(result.rhs == std::vector{{1.0, 0.0}}); REQUIRE(result.free_edge_indices.size() == 1); @@ -162,17 +161,23 @@ TEST_CASE("Test the link solver algorithm") { auto const node_loads = std::vector{{-1.0, -1.0}, {-1.0, -1.0}, {2.0, 2.0}, {0.0, 0.0}, {0.0, 0.0}}; auto const result = forward_elimination(edges, node_loads); - - REQUIRE(result.matrix.data.size() == 14); - // TODO(figueroa1395): Some checks are commented out because the use of unordered_set in adjacency list - // makes the order of re-attachments non-deterministic (per-platform - locally it's fine), which affects the - // order of entries in the COO matrix should we use an ordered set instead to test this fine grained - // details, or is it enough to test overall structure and result correctness? CHECK(result.matrix.data == - // std::vector{1, 1, 1, 1, 1, -1, -1, 1, -1, -1, -1, -1, 1, 1}); - REQUIRE(result.matrix.row.size() == 14); - CHECK(result.matrix.row == std::vector{0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3}); - REQUIRE(result.matrix.col.size() == 14); - // CHECK(result.matrix.col == std::vector{0, 2, 1, 1, 2, 6, 3, 2, 3, 6, 5, 4, 5, 6}); + uint64_t const node_number = node_loads.size(); + + REQUIRE(result.matrix.data_map.size() == 14); + CHECK(result.matrix.data_map.at(data_map_key_gen(0, 0, node_number)) == 1); + CHECK(result.matrix.data_map.at(data_map_key_gen(0, 2, node_number)) == 1); + CHECK(result.matrix.data_map.at(data_map_key_gen(0, 1, node_number)) == 1); + CHECK(result.matrix.data_map.at(data_map_key_gen(1, 1, node_number)) == 1); + CHECK(result.matrix.data_map.at(data_map_key_gen(1, 2, node_number)) == 1); + CHECK(result.matrix.data_map.at(data_map_key_gen(1, 6, node_number)) == -1); + CHECK(result.matrix.data_map.at(data_map_key_gen(1, 3, node_number)) == -1); + CHECK(result.matrix.data_map.at(data_map_key_gen(2, 2, node_number)) == 1); + CHECK(result.matrix.data_map.at(data_map_key_gen(2, 3, node_number)) == -1); + CHECK(result.matrix.data_map.at(data_map_key_gen(2, 6, node_number)) == -1); + CHECK(result.matrix.data_map.at(data_map_key_gen(2, 5, node_number)) == -1); + CHECK(result.matrix.data_map.at(data_map_key_gen(2, 4, node_number)) == -1); + CHECK(result.matrix.data_map.at(data_map_key_gen(3, 5, node_number)) == 1); + CHECK(result.matrix.data_map.at(data_map_key_gen(3, 6, node_number)) == 1); REQUIRE(result.rhs.size() == 4); CHECK(result.rhs == std::vector{{-1.0, -1.0}, {-1.0, -1.0}, {-2.0, -2.0}, {0.0, 0.0}}); REQUIRE(result.free_edge_indices.size() == 3); From 95d6a4847f400b54b18f6dcb7c303af614326cc4 Mon Sep 17 00:00:00 2001 From: Santiago Figueroa Manrique Date: Tue, 17 Mar 2026 11:26:33 +0100 Subject: [PATCH 08/12] sonar cloud Signed-off-by: Santiago Figueroa Manrique --- .../power_grid_model/include/power_grid_model/link_solver.hpp | 3 +-- tests/cpp_unit_tests/test_link_solver.cpp | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp index 15087f7aef..9f21f14e1b 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp @@ -46,8 +46,7 @@ struct COOSparseMatrix { // TODO(figueroa1395): maybe this is overkill - revisit after implementation is further ahead [[nodiscard]] bool get_value(IntS& value, uint64_t row_idx, uint64_t col_idx) const { assert(row_number != 0 && "row_number must be set before getting values from data_map"); - auto const it = data_map.find(row_idx * row_number + col_idx); - if (it != data_map.end()) { + if (auto const it = data_map.find(row_idx * row_number + col_idx); it != data_map.end()) { value = it->second; return true; } diff --git a/tests/cpp_unit_tests/test_link_solver.cpp b/tests/cpp_unit_tests/test_link_solver.cpp index a1ca651f71..4c9fd5d07f 100644 --- a/tests/cpp_unit_tests/test_link_solver.cpp +++ b/tests/cpp_unit_tests/test_link_solver.cpp @@ -15,7 +15,7 @@ namespace power_grid_model::link_solver { namespace { -[[nodiscard]] IntS data_map_key_gen(uint64_t row_idx, uint64_t col_idx, uint64_t row_number) { +[[nodiscard]] uint64_t data_map_key_gen(uint64_t row_idx, uint64_t col_idx, uint64_t row_number) { return row_idx * row_number + col_idx; } } // namespace From 7c99d73648e7a8dfeae030907290590467b13f1d Mon Sep 17 00:00:00 2001 From: Santiago Figueroa Manrique Date: Tue, 17 Mar 2026 19:37:16 +0100 Subject: [PATCH 09/12] backward substitution implementation Signed-off-by: Santiago Figueroa Manrique --- .../include/power_grid_model/link_solver.hpp | 60 ++++++++++++++++--- 1 file changed, 52 insertions(+), 8 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp index 9f21f14e1b..b93979ca7d 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp @@ -43,8 +43,7 @@ struct COOSparseMatrix { assert(row_number != 0 && "row_number must be set before setting values in data_map"); data_map[row_idx * row_number + col_idx] = value; } - // TODO(figueroa1395): maybe this is overkill - revisit after implementation is further ahead - [[nodiscard]] bool get_value(IntS& value, uint64_t row_idx, uint64_t col_idx) const { + bool get_value(IntS& value, uint64_t row_idx, uint64_t col_idx) const { assert(row_number != 0 && "row_number must be set before getting values from data_map"); if (auto const it = data_map.find(row_idx * row_number + col_idx); it != data_map.end()) { value = it->second; @@ -52,6 +51,18 @@ struct COOSparseMatrix { } return false; } + void add_to_value(IntS value, uint64_t row_idx, uint64_t col_idx) { + assert(row_number != 0 && "row_number must be set before adding values in data_map"); + auto const key = row_idx * row_number + col_idx; + if (auto const it = data_map.find(key); it != data_map.end()) { + it->second += value; + if (it->second == 0) { + data_map.erase(it); // maintain sparsity by erasing zero entries + } + } else if (value != 0) { + data_map[key] = value; + } + } }; struct EdgeHistory { @@ -61,9 +72,10 @@ struct EdgeHistory { struct EliminationResult { COOSparseMatrix matrix{}; - std::vector rhs{}; // RHS value at each pivot row - std::vector free_edge_indices{}; // index of degrees of freedom (self loop edges) - std::vector edges_history{}; // edges elimination history + std::vector rhs{}; // RHS value at each pivot row + std::vector free_edge_indices{}; // index of degrees of freedom (self loop edges) + std::vector pivot_edge_indices{}; // index of pivot edges + std::vector edges_history{}; // edges elimination history }; // convention: from node at position 0, to node at position 1 @@ -122,14 +134,14 @@ inline void update_edge_info(uint64_t edge_idx, uint64_t matrix_row, std::vector [[nodiscard]] inline EliminationResult forward_elimination(std::vector edges, std::vector node_loads) { EliminationResult result{}; - uint64_t edge_number{edges.size()}; - uint64_t node_number{node_loads.size()}; + uint64_t const edge_number{edges.size()}; + uint64_t const node_number{node_loads.size()}; result.edges_history.resize(edge_number); using enum EdgeEvent; using enum EdgeDirection; - auto& [matrix, rhs, free_edge_indices, edges_history] = result; + auto& [matrix, rhs, free_edge_indices, pivot_edge_indices, edges_history] = result; matrix.prepare(edge_number, node_number); constexpr uint8_t starting_row{}; @@ -143,6 +155,7 @@ inline void update_edge_info(uint64_t edge_idx, uint64_t matrix_row, std::vector free_edge_indices.push_back(index); } else { write_edge_history(edges_history[index], Deleted, matrix_row); // Delete edge -> pivot there + pivot_edge_indices.push_back(index); matrix.set_value(std::to_underlying(Towards), matrix_row, index); Idx const from_node_idx = from_node(edge); @@ -175,4 +188,35 @@ inline void update_edge_info(uint64_t edge_idx, uint64_t matrix_row, std::vector } return result; } + +// backward substitution is performed in a sparse way +// using the result from the elimination game +inline void backward_substitution(EliminationResult& elimination_result) { + auto pivot_col_indices = elimination_result.pivot_edge_indices | std::views::drop(1) | std::ranges::views::reverse; + auto free_col_indices = std::span(elimination_result.free_edge_indices); + + for (auto const pivot_col_idx : pivot_col_indices) { + auto const& edge_history = elimination_result.edges_history[pivot_col_idx]; + uint64_t const pivot_row_idx = edge_history.rows.back(); + + for (auto const row_idx : edge_history.rows | std::views::reverse | std::views::drop(1)) { + IntS multiplier_value{}; + elimination_result.matrix.get_value(multiplier_value, row_idx, pivot_col_idx); + elimination_result.matrix.add_to_value( + -multiplier_value, row_idx, + pivot_col_idx); // the initial pivot is always one because of the forward sweep + + for (auto const backward_col_idx : std::ranges::subrange( + std::ranges::upper_bound(free_col_indices, pivot_col_idx), free_col_indices.end())) { + IntS pivot_value{}; + if (elimination_result.matrix.get_value(pivot_value, pivot_row_idx, backward_col_idx)) { + elimination_result.matrix.add_to_value(-multiplier_value * pivot_value, row_idx, backward_col_idx); + } + } + + elimination_result.rhs[row_idx] += + -static_cast(multiplier_value) * elimination_result.rhs[pivot_row_idx]; + } + } +} } // namespace power_grid_model::link_solver::detail From 0f8249c36c33470600e45726f3a315a771c2cbc2 Mon Sep 17 00:00:00 2001 From: Santiago Figueroa Manrique Date: Wed, 18 Mar 2026 08:35:48 +0100 Subject: [PATCH 10/12] test backward substitution Signed-off-by: Santiago Figueroa Manrique --- .../include/power_grid_model/link_solver.hpp | 7 +- tests/cpp_unit_tests/test_link_solver.cpp | 118 +++++++++++++----- 2 files changed, 92 insertions(+), 33 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp index b93979ca7d..14f4b441e7 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp @@ -55,7 +55,7 @@ struct COOSparseMatrix { assert(row_number != 0 && "row_number must be set before adding values in data_map"); auto const key = row_idx * row_number + col_idx; if (auto const it = data_map.find(key); it != data_map.end()) { - it->second += value; + it->second = static_cast(it->second + value); if (it->second == 0) { data_map.erase(it); // maintain sparsity by erasing zero entries } @@ -203,14 +203,15 @@ inline void backward_substitution(EliminationResult& elimination_result) { IntS multiplier_value{}; elimination_result.matrix.get_value(multiplier_value, row_idx, pivot_col_idx); elimination_result.matrix.add_to_value( - -multiplier_value, row_idx, + static_cast(-multiplier_value), row_idx, pivot_col_idx); // the initial pivot is always one because of the forward sweep for (auto const backward_col_idx : std::ranges::subrange( std::ranges::upper_bound(free_col_indices, pivot_col_idx), free_col_indices.end())) { IntS pivot_value{}; if (elimination_result.matrix.get_value(pivot_value, pivot_row_idx, backward_col_idx)) { - elimination_result.matrix.add_to_value(-multiplier_value * pivot_value, row_idx, backward_col_idx); + elimination_result.matrix.add_to_value(static_cast(-multiplier_value * pivot_value), row_idx, + backward_col_idx); } } diff --git a/tests/cpp_unit_tests/test_link_solver.cpp b/tests/cpp_unit_tests/test_link_solver.cpp index 4c9fd5d07f..0462eace76 100644 --- a/tests/cpp_unit_tests/test_link_solver.cpp +++ b/tests/cpp_unit_tests/test_link_solver.cpp @@ -15,8 +15,10 @@ namespace power_grid_model::link_solver { namespace { -[[nodiscard]] uint64_t data_map_key_gen(uint64_t row_idx, uint64_t col_idx, uint64_t row_number) { - return row_idx * row_number + col_idx; +void check_value(IntS expected, uint64_t row_idx, uint64_t col_idx, detail::COOSparseMatrix const& matrix) { + IntS value{}; + CHECK(matrix.get_value(value, row_idx, col_idx)); + CHECK(value == expected); } } // namespace TEST_CASE("Test the link solver algorithm") { @@ -82,10 +84,9 @@ TEST_CASE("Test the link solver algorithm") { auto const edges = std::vector{{0, 1}}; auto const node_loads = std::vector{{-1.0, 0.0}, {1.0, 0.0}}; auto const result = forward_elimination(edges, node_loads); - uint64_t const node_number = node_loads.size(); REQUIRE(result.matrix.data_map.size() == 1); - CHECK(result.matrix.data_map.at(data_map_key_gen(0, 0, node_number)) == 1); + check_value(1, 0, 0, result.matrix); CHECK(result.rhs == std::vector{{1.0, 0.0}}); CHECK(result.free_edge_indices.empty()); REQUIRE(result.edges_history.size() == 1); @@ -97,11 +98,10 @@ TEST_CASE("Test the link solver algorithm") { auto const edges = std::vector{{1, 0}, {1, 2}}; auto const node_loads = std::vector{{-1.0, 0.0}, {1.0, 0.0}, {0.0, 0.0}}; auto const result = forward_elimination(edges, node_loads); - uint64_t const node_number = node_loads.size(); REQUIRE(result.matrix.data_map.size() == 2); - CHECK(result.matrix.data_map.at(data_map_key_gen(0, 0, node_number)) == 1); - CHECK(result.matrix.data_map.at(data_map_key_gen(1, 1, node_number)) == 1); + check_value(1, 0, 0, result.matrix); + check_value(1, 1, 1, result.matrix); REQUIRE(result.rhs.size() == 2); CHECK(result.rhs == std::vector{{-1.0, 0.0}, {0.0, 0.0}}); CHECK(result.free_edge_indices.empty()); @@ -116,13 +116,12 @@ TEST_CASE("Test the link solver algorithm") { auto const edges = std::vector{{0, 1}, {1, 2}, {2, 0}}; auto const node_loads = std::vector{{-1.0, 0.0}, {1.0, 0.0}, {0.0, 0.0}}; auto const result = forward_elimination(edges, node_loads); - uint64_t const node_number = node_loads.size(); REQUIRE(result.matrix.data_map.size() == 4); - CHECK(result.matrix.data_map.at(data_map_key_gen(0, 0, node_number)) == 1); - CHECK(result.matrix.data_map.at(data_map_key_gen(0, 1, node_number)) == -1); - CHECK(result.matrix.data_map.at(data_map_key_gen(1, 1, node_number)) == 1); - CHECK(result.matrix.data_map.at(data_map_key_gen(1, 2, node_number)) == -1); + check_value(1, 0, 0, result.matrix); + check_value(-1, 0, 1, result.matrix); + check_value(1, 1, 1, result.matrix); + check_value(-1, 1, 2, result.matrix); REQUIRE(result.rhs.size() == 2); CHECK(result.rhs == std::vector{{1.0, 0.0}, {0.0, 0.0}}); REQUIRE(result.free_edge_indices.size() == 1); @@ -140,11 +139,10 @@ TEST_CASE("Test the link solver algorithm") { auto const edges = std::vector{{0, 1}, {0, 1}}; auto const node_loads = std::vector{{-1.0, 0.0}, {1.0, 0.0}}; auto const result = forward_elimination(edges, node_loads); - uint64_t const node_number = node_loads.size(); REQUIRE(result.matrix.data_map.size() == 2); - CHECK(result.matrix.data_map.at(data_map_key_gen(0, 0, node_number)) == 1); - CHECK(result.matrix.data_map.at(data_map_key_gen(0, 1, node_number)) == 1); + check_value(1, 0, 0, result.matrix); + check_value(1, 0, 1, result.matrix); REQUIRE(result.rhs.size() == 1); CHECK(result.rhs == std::vector{{1.0, 0.0}}); REQUIRE(result.free_edge_indices.size() == 1); @@ -161,23 +159,22 @@ TEST_CASE("Test the link solver algorithm") { auto const node_loads = std::vector{{-1.0, -1.0}, {-1.0, -1.0}, {2.0, 2.0}, {0.0, 0.0}, {0.0, 0.0}}; auto const result = forward_elimination(edges, node_loads); - uint64_t const node_number = node_loads.size(); REQUIRE(result.matrix.data_map.size() == 14); - CHECK(result.matrix.data_map.at(data_map_key_gen(0, 0, node_number)) == 1); - CHECK(result.matrix.data_map.at(data_map_key_gen(0, 2, node_number)) == 1); - CHECK(result.matrix.data_map.at(data_map_key_gen(0, 1, node_number)) == 1); - CHECK(result.matrix.data_map.at(data_map_key_gen(1, 1, node_number)) == 1); - CHECK(result.matrix.data_map.at(data_map_key_gen(1, 2, node_number)) == 1); - CHECK(result.matrix.data_map.at(data_map_key_gen(1, 6, node_number)) == -1); - CHECK(result.matrix.data_map.at(data_map_key_gen(1, 3, node_number)) == -1); - CHECK(result.matrix.data_map.at(data_map_key_gen(2, 2, node_number)) == 1); - CHECK(result.matrix.data_map.at(data_map_key_gen(2, 3, node_number)) == -1); - CHECK(result.matrix.data_map.at(data_map_key_gen(2, 6, node_number)) == -1); - CHECK(result.matrix.data_map.at(data_map_key_gen(2, 5, node_number)) == -1); - CHECK(result.matrix.data_map.at(data_map_key_gen(2, 4, node_number)) == -1); - CHECK(result.matrix.data_map.at(data_map_key_gen(3, 5, node_number)) == 1); - CHECK(result.matrix.data_map.at(data_map_key_gen(3, 6, node_number)) == 1); + check_value(1, 0, 0, result.matrix); + check_value(1, 0, 2, result.matrix); + check_value(1, 0, 1, result.matrix); + check_value(1, 1, 1, result.matrix); + check_value(1, 1, 2, result.matrix); + check_value(-1, 1, 6, result.matrix); + check_value(-1, 1, 3, result.matrix); + check_value(1, 2, 2, result.matrix); + check_value(-1, 2, 3, result.matrix); + check_value(-1, 2, 6, result.matrix); + check_value(-1, 2, 5, result.matrix); + check_value(-1, 2, 4, result.matrix); + check_value(1, 3, 5, result.matrix); + check_value(1, 3, 6, result.matrix); REQUIRE(result.rhs.size() == 4); CHECK(result.rhs == std::vector{{-1.0, -1.0}, {-1.0, -1.0}, {-2.0, -2.0}, {0.0, 0.0}}); REQUIRE(result.free_edge_indices.size() == 3); @@ -199,6 +196,67 @@ TEST_CASE("Test the link solver algorithm") { CHECK(result.edges_history[6].rows == std::vector{1, 2, 3}); } } + + SUBCASE("Test backward substitution") { + using enum EdgeEvent; + + SUBCASE("Complex case with complex loads") { + EliminationResult result{}; + result.matrix.prepare(5, 7); + result.matrix.set_value(1, 0, 0); + result.matrix.set_value(1, 0, 1); + result.matrix.set_value(1, 0, 2); + result.matrix.set_value(1, 1, 1); + result.matrix.set_value(1, 1, 2); + result.matrix.set_value(-1, 1, 3); + result.matrix.set_value(-1, 1, 6); + result.matrix.set_value(1, 2, 2); + result.matrix.set_value(-1, 2, 3); + result.matrix.set_value(-1, 2, 4); + result.matrix.set_value(-1, 2, 5); + result.matrix.set_value(-1, 2, 6); + result.matrix.set_value(1, 3, 5); + result.matrix.set_value(1, 3, 6); + + result.rhs = std::vector{{-1.0, -1.0}, {-1.0, -1.0}, {-2.0, -2.0}, {0.0, 0.0}, {0.0, 0.0}}; + + result.free_edge_indices = std::vector{3, 4, 6}; + result.pivot_edge_indices = std::vector{0, 1, 2, 5}; + + result.edges_history.resize(7); + result.edges_history[0].events = std::vector{Deleted}; + result.edges_history[0].rows = std::vector{0}; + result.edges_history[1].events = std::vector{Replaced, Deleted}; + result.edges_history[1].rows = std::vector{0, 1}; + result.edges_history[2].events = std::vector{Replaced, Replaced, Deleted}; + result.edges_history[2].rows = std::vector{0, 1, 2}; + result.edges_history[3].events = std::vector{Replaced, ContractedToPoint}; + result.edges_history[3].rows = std::vector{1, 2}; + result.edges_history[4].events = std::vector{ContractedToPoint}; + result.edges_history[4].rows = std::vector{2}; + result.edges_history[5].events = std::vector{Replaced, Deleted}; + result.edges_history[5].rows = std::vector{2, 3}; + result.edges_history[6].events = std::vector{Replaced, Replaced, ContractedToPoint}; + result.edges_history[6].rows = std::vector{1, 2, 3}; + + backward_substitution(result); + REQUIRE(result.matrix.data_map.size() == 11); + check_value(1, 0, 0, result.matrix); + check_value(1, 0, 3, result.matrix); + check_value(1, 0, 6, result.matrix); + check_value(1, 1, 1, result.matrix); + check_value(1, 1, 4, result.matrix); + check_value(-1, 1, 6, result.matrix); + check_value(1, 2, 2, result.matrix); + check_value(-1, 2, 3, result.matrix); + check_value(-1, 2, 4, result.matrix); + check_value(1, 3, 5, result.matrix); + check_value(1, 3, 6, result.matrix); + REQUIRE(result.rhs.size() == 5); + CHECK(result.rhs == + std::vector{{0.0, 0.0}, {1.0, 1.0}, {-2.0, -2.0}, {0.0, 0.0}, {0.0, 0.0}}); + } + } } } // namespace power_grid_model::link_solver From a1b591bfa91f40f3fd259897d253c95b27324d33 Mon Sep 17 00:00:00 2001 From: Santiago Figueroa Manrique Date: Wed, 18 Mar 2026 10:02:09 +0100 Subject: [PATCH 11/12] reduced echelon form function impl Signed-off-by: Santiago Figueroa Manrique --- .../include/power_grid_model/link_solver.hpp | 30 ++++++++++++------ tests/cpp_unit_tests/test_link_solver.cpp | 31 ++++++++++++++++--- 2 files changed, 47 insertions(+), 14 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp index 14f4b441e7..1c179d7506 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp @@ -131,18 +131,12 @@ inline void update_edge_info(uint64_t edge_idx, uint64_t matrix_row, std::vector // forward elimination is performed via a modified Gaussian elimination procedure // we name this procedure elimination-game // node_loads convention: caller passes negated external loads (as in RHS of power flow equations)) -[[nodiscard]] inline EliminationResult forward_elimination(std::vector edges, - std::vector node_loads) { - EliminationResult result{}; - uint64_t const edge_number{edges.size()}; - uint64_t const node_number{node_loads.size()}; - - result.edges_history.resize(edge_number); +inline void forward_elimination(EliminationResult& result, std::vector edges, + std::vector node_loads) { using enum EdgeEvent; using enum EdgeDirection; auto& [matrix, rhs, free_edge_indices, pivot_edge_indices, edges_history] = result; - matrix.prepare(edge_number, node_number); constexpr uint8_t starting_row{}; uint64_t matrix_row{starting_row}; @@ -186,7 +180,6 @@ inline void update_edge_info(uint64_t edge_idx, uint64_t matrix_row, std::vector ++matrix_row; } } - return result; } // backward substitution is performed in a sparse way @@ -220,4 +213,23 @@ inline void backward_substitution(EliminationResult& elimination_result) { } } } + +// reduced echelon form based in custom forward elimination and backward substitution procedures +[[nodiscard]] inline EliminationResult reduced_echelon_form(std::vector edges, + std::vector node_loads) { + EliminationResult result{}; + uint64_t const edge_number{edges.size()}; + uint64_t const node_number{node_loads.size()}; + + result.edges_history.resize(edge_number); + result.matrix.prepare(edge_number, node_number); + + // -1 because the loads represent the RHS + std::ranges::transform(node_loads, node_loads.begin(), [](auto load) { return -load; }); + + // both edges and node_loads are modified and consumed in the forward sweep + forward_elimination(result, std::move(edges), std::move(node_loads)); + backward_substitution(result); + return result; +} } // namespace power_grid_model::link_solver::detail diff --git a/tests/cpp_unit_tests/test_link_solver.cpp b/tests/cpp_unit_tests/test_link_solver.cpp index 0462eace76..a7f109b2bd 100644 --- a/tests/cpp_unit_tests/test_link_solver.cpp +++ b/tests/cpp_unit_tests/test_link_solver.cpp @@ -79,11 +79,16 @@ TEST_CASE("Test the link solver algorithm") { SUBCASE("Test forward elimination - elimination game") { using enum EdgeEvent; + EliminationResult result{}; SUBCASE("One edge, two nodes, two real loads") { auto const edges = std::vector{{0, 1}}; auto const node_loads = std::vector{{-1.0, 0.0}, {1.0, 0.0}}; - auto const result = forward_elimination(edges, node_loads); + uint64_t const edge_number{edges.size()}; + uint64_t const node_number{node_loads.size()}; + result.edges_history.resize(edge_number); + result.matrix.prepare(edge_number, node_number); + forward_elimination(result, edges, node_loads); REQUIRE(result.matrix.data_map.size() == 1); check_value(1, 0, 0, result.matrix); @@ -97,7 +102,11 @@ TEST_CASE("Test the link solver algorithm") { SUBCASE("Two edges, three nodes, two real loads") { auto const edges = std::vector{{1, 0}, {1, 2}}; auto const node_loads = std::vector{{-1.0, 0.0}, {1.0, 0.0}, {0.0, 0.0}}; - auto const result = forward_elimination(edges, node_loads); + uint64_t const edge_number{edges.size()}; + uint64_t const node_number{node_loads.size()}; + result.edges_history.resize(edge_number); + result.matrix.prepare(edge_number, node_number); + forward_elimination(result, edges, node_loads); REQUIRE(result.matrix.data_map.size() == 2); check_value(1, 0, 0, result.matrix); @@ -115,7 +124,11 @@ TEST_CASE("Test the link solver algorithm") { SUBCASE("Three edges, three nodes, two real loads") { auto const edges = std::vector{{0, 1}, {1, 2}, {2, 0}}; auto const node_loads = std::vector{{-1.0, 0.0}, {1.0, 0.0}, {0.0, 0.0}}; - auto const result = forward_elimination(edges, node_loads); + uint64_t const edge_number{edges.size()}; + uint64_t const node_number{node_loads.size()}; + result.edges_history.resize(edge_number); + result.matrix.prepare(edge_number, node_number); + forward_elimination(result, edges, node_loads); REQUIRE(result.matrix.data_map.size() == 4); check_value(1, 0, 0, result.matrix); @@ -138,7 +151,11 @@ TEST_CASE("Test the link solver algorithm") { SUBCASE("Two edges, two nodes, two real loads") { auto const edges = std::vector{{0, 1}, {0, 1}}; auto const node_loads = std::vector{{-1.0, 0.0}, {1.0, 0.0}}; - auto const result = forward_elimination(edges, node_loads); + uint64_t const edge_number{edges.size()}; + uint64_t const node_number{node_loads.size()}; + result.edges_history.resize(edge_number); + result.matrix.prepare(edge_number, node_number); + forward_elimination(result, edges, node_loads); REQUIRE(result.matrix.data_map.size() == 2); check_value(1, 0, 0, result.matrix); @@ -158,7 +175,11 @@ TEST_CASE("Test the link solver algorithm") { auto const edges = std::vector{{3, 0}, {1, 0}, {2, 0}, {3, 2}, {1, 2}, {1, 4}, {3, 4}}; auto const node_loads = std::vector{{-1.0, -1.0}, {-1.0, -1.0}, {2.0, 2.0}, {0.0, 0.0}, {0.0, 0.0}}; - auto const result = forward_elimination(edges, node_loads); + uint64_t const edge_number{edges.size()}; + uint64_t const node_number{node_loads.size()}; + result.edges_history.resize(edge_number); + result.matrix.prepare(edge_number, node_number); + forward_elimination(result, edges, node_loads); REQUIRE(result.matrix.data_map.size() == 14); check_value(1, 0, 0, result.matrix); From 66d5820300b2419ea6b0ea53830987025f99ecf1 Mon Sep 17 00:00:00 2001 From: Santiago Figueroa Manrique Date: Wed, 18 Mar 2026 15:52:09 +0100 Subject: [PATCH 12/12] more tests for backwards substitution Signed-off-by: Santiago Figueroa Manrique --- tests/cpp_unit_tests/test_link_solver.cpp | 90 +++++++++++++++++++++++ 1 file changed, 90 insertions(+) mode change 100644 => 100755 tests/cpp_unit_tests/test_link_solver.cpp diff --git a/tests/cpp_unit_tests/test_link_solver.cpp b/tests/cpp_unit_tests/test_link_solver.cpp old mode 100644 new mode 100755 index a7f109b2bd..aaabe13d49 --- a/tests/cpp_unit_tests/test_link_solver.cpp +++ b/tests/cpp_unit_tests/test_link_solver.cpp @@ -221,6 +221,96 @@ TEST_CASE("Test the link solver algorithm") { SUBCASE("Test backward substitution") { using enum EdgeEvent; + SUBCASE("One edge, two nodes, two real loads") { + EliminationResult result{}; + result.matrix.prepare(1, 2); + result.matrix.set_value(1, 0, 0); + result.rhs = std::vector{{1.0, 0.0}}; + result.free_edge_indices = {}; + result.pivot_edge_indices = {0}; + result.edges_history.resize(1); + result.edges_history[0].events = {Deleted}; + result.edges_history[0].rows = {0}; + + backward_substitution(result); + REQUIRE(result.matrix.data_map.size() == 1); + check_value(1, 0, 0, result.matrix); + REQUIRE(result.rhs.size() == 1); + CHECK(result.rhs == std::vector{{1.0, 0.0}}); + } + + SUBCASE("Two edges, three nodes, two real loads") { + EliminationResult result{}; + result.matrix.prepare(2, 3); + result.matrix.set_value(1, 0, 0); + result.matrix.set_value(1, 1, 1); + result.rhs = std::vector{{-1.0, 0.0}, {0.0, 0.0}}; + result.free_edge_indices = {}; + result.pivot_edge_indices = {0, 1}; + result.edges_history.resize(2); + result.edges_history[0].events = {Deleted}; + result.edges_history[0].rows = {0}; + result.edges_history[1].events = {Deleted}; + result.edges_history[1].rows = {1}; + + backward_substitution(result); + REQUIRE(result.matrix.data_map.size() == 2); + check_value(1, 0, 0, result.matrix); + check_value(1, 1, 1, result.matrix); + REQUIRE(result.rhs.size() == 2); + CHECK(result.rhs == std::vector{{-1.0, 0.0}, {0.0, 0.0}}); + } + + SUBCASE("Three edges, three nodes, two real loads") { + EliminationResult result{}; + result.matrix.prepare(4, 3); + result.matrix.set_value(1, 0, 0); + result.matrix.set_value(-1, 0, 1); + result.matrix.set_value(1, 1, 1); + result.matrix.set_value(-1, 1, 2); + result.rhs = std::vector{{1.0, 0.0}, {0.0, 0.0}}; + result.free_edge_indices = {2}; + result.pivot_edge_indices = {0, 1}; + result.edges_history.resize(3); + result.edges_history[0].events = {Deleted}; + result.edges_history[0].rows = {0}; + result.edges_history[1].events = {Replaced, Deleted}; + result.edges_history[1].rows = {0, 1}; + result.edges_history[2].events = {ContractedToPoint}; + result.edges_history[2].rows = {1}; + + backward_substitution(result); + REQUIRE(result.matrix.data_map.size() == 4); + check_value(1, 0, 0, result.matrix); + check_value(1, 1, 1, result.matrix); + check_value(-1, 1, 2, result.matrix); + check_value(-1, 0, 2, result.matrix); + REQUIRE(result.rhs.size() == 2); + CHECK(result.rhs == std::vector{{1.0, 0.0}, {0.0, 0.0}}); + } + + SUBCASE("Two edges, two nodes, two real loads") { + EliminationResult result{}; + result.matrix.prepare(2, 2); + result.matrix.set_value(1, 0, 0); + result.matrix.set_value(1, 0, 1); + result.rhs = std::vector{{1.0, 0.0}}; + result.free_edge_indices = {1}; + result.pivot_edge_indices = {0}; + result.edges_history.resize(2); + result.edges_history[0].events = {Deleted}; + result.edges_history[0].rows = {0}; + result.edges_history[1].events = {ContractedToPoint}; + result.edges_history[1].rows = {0}; + + backward_substitution(result); + REQUIRE(result.matrix.data_map.size() == 2); + check_value(1, 0, 0, result.matrix); + check_value(1, 0, 1, result.matrix); + REQUIRE(result.rhs.size() == 1); + CHECK(result.rhs == std::vector{{1.0, 0.0}}); + } + SUBCASE("Complex case with complex loads") { EliminationResult result{}; result.matrix.prepare(5, 7);