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..1c179d7506
--- /dev/null
+++ b/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp
@@ -0,0 +1,235 @@
+// SPDX-FileCopyrightText: Contributors to the Power Grid Model project <powergridmodel@lfenergy.org>
+//
+// SPDX-License-Identifier: MPL-2.0
+
+#pragma once
+
+#include "calculation_parameters.hpp"
+#include "common/common.hpp"
+#include "common/counting_iterator.hpp"
+
+#include <array>
+#include <cstdint>
+#include <ranges>
+#include <span>
+#include <unordered_map>
+#include <unordered_set>
+#include <utility>
+#include <vector>
+
+namespace power_grid_model::link_solver::detail {
+
+using AdjacencyMap = std::unordered_map<Idx, std::unordered_set<uint64_t>>;
+using ContractedEdgesSet = std::unordered_set<uint64_t>;
+
+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 {
+    uint64_t row_number{};
+    std::unordered_map<uint64_t, IntS> 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;
+    }
+    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;
+            return true;
+        }
+        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 = static_cast<IntS>(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 {
+    std::vector<uint64_t> rows{};
+    std::vector<EdgeEvent> events{};
+};
+
+struct EliminationResult {
+    COOSparseMatrix matrix{};
+    std::vector<DoubleComplex> rhs{};           // RHS value at each pivot row
+    std::vector<uint64_t> free_edge_indices{};  // index of degrees of freedom (self loop edges)
+    std::vector<uint64_t> pivot_edge_indices{}; // index of pivot edges
+    std::vector<EdgeHistory> edges_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]] inline auto build_adjacency_map(std::span<BranchIdx> edges) -> AdjacencyMap {
+    AdjacencyMap adjacency_map{};
+    for (auto const& [raw_index, edge] : enumerate(std::as_const(edges))) {
+        auto const index = static_cast<uint64_t>(raw_index);
+        auto const [from_node, to_node] = edge;
+        adjacency_map[from_node].insert(index);
+        adjacency_map[to_node].insert(index);
+    }
+    return adjacency_map;
+}
+
+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 replace_and_write(uint64_t edge_idx, Idx from_node_idx, Idx to_node_idx, uint64_t matrix_row,
+                              std::vector<BranchIdx>& edges, COOSparseMatrix& matrix) {
+    using enum EdgeDirection;
+
+    if (from_node(edges[edge_idx]) == to_node_idx) {
+        from_node(edges[edge_idx]) = from_node_idx; // re-attachment step
+        matrix.set_value(std::to_underlying(Away), matrix_row, edge_idx);
+    } else {
+        to_node(edges[edge_idx]) = from_node_idx; // re-attachment step
+        matrix.set_value(std::to_underlying(Towards), matrix_row, edge_idx);
+    }
+}
+
+inline void update_edge_info(uint64_t edge_idx, uint64_t matrix_row, std::vector<BranchIdx>& edges,
+                             std::vector<EdgeHistory>& edges_history, ContractedEdgesSet& edges_contracted_to_point) {
+    using enum EdgeEvent;
+
+    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);
+    }
+}
+
+// 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))
+inline void forward_elimination(EliminationResult& result, std::vector<BranchIdx> edges,
+                                std::vector<DoubleComplex> node_loads) {
+    using enum EdgeEvent;
+    using enum EdgeDirection;
+
+    auto& [matrix, rhs, free_edge_indices, pivot_edge_indices, edges_history] = result;
+
+    constexpr uint8_t starting_row{};
+    uint64_t matrix_row{starting_row};
+    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<uint64_t>(raw_index);
+        if (edges_contracted_to_point.contains(index)) {
+            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);
+            Idx const to_node_idx = to_node(edge);
+
+            // update adjacency list for deleted edge
+            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<uint64_t>(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_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);
+            }
+            ++matrix_row;
+        }
+    }
+}
+
+// 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(
+                static_cast<IntS>(-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(static_cast<IntS>(-multiplier_value * pivot_value), row_idx,
+                                                           backward_col_idx);
+                }
+            }
+
+            elimination_result.rhs[row_idx] +=
+                -static_cast<DoubleComplex>(multiplier_value) * elimination_result.rhs[pivot_row_idx];
+        }
+    }
+}
+
+// reduced echelon form based in custom forward elimination and backward substitution procedures
+[[nodiscard]] inline EliminationResult reduced_echelon_form(std::vector<BranchIdx> edges,
+                                                            std::vector<DoubleComplex> 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/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(
diff --git a/tests/cpp_unit_tests/test_link_solver.cpp b/tests/cpp_unit_tests/test_link_solver.cpp
new file mode 100755
index 0000000000..aaabe13d49
--- /dev/null
+++ b/tests/cpp_unit_tests/test_link_solver.cpp
@@ -0,0 +1,373 @@
+// SPDX-FileCopyrightText: Contributors to the Power Grid Model project <powergridmodel@lfenergy.org>
+//
+// SPDX-License-Identifier: MPL-2.0
+
+#include <power_grid_model/calculation_parameters.hpp>
+#include <power_grid_model/common/common.hpp>
+#include <power_grid_model/link_solver.hpp>
+
+#include <doctest/doctest.h>
+
+#include <cstdint>
+#include <unordered_map>
+#include <unordered_set>
+#include <vector>
+
+namespace power_grid_model::link_solver {
+namespace {
+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") {
+    using namespace detail;
+
+    SUBCASE("Test build adjacency list") {
+
+        SUBCASE("One edge, two nodes") {
+            auto edges = std::vector<BranchIdx>{{0, 1}};
+            AdjacencyMap const adjacency_map = build_adjacency_map(edges);
+
+            REQUIRE(adjacency_map.size() == 2);
+            CHECK(adjacency_map.at(0) == std::unordered_set<uint64_t>{0});
+            CHECK(adjacency_map.at(1) == std::unordered_set<uint64_t>{0});
+        }
+
+        SUBCASE("Two edges, three nodes") {
+            auto edges = std::vector<BranchIdx>{{1, 0}, {1, 2}};
+            AdjacencyMap const adjacency_map = build_adjacency_map(edges);
+
+            REQUIRE(adjacency_map.size() == 3);
+            CHECK(adjacency_map.at(0) == std::unordered_set<uint64_t>{0});
+            CHECK(adjacency_map.at(1) == std::unordered_set<uint64_t>{0, 1});
+            CHECK(adjacency_map.at(2) == std::unordered_set<uint64_t>{1});
+        }
+
+        SUBCASE("Three edges, three nodes") {
+            auto edges = std::vector<BranchIdx>{{0, 1}, {1, 2}, {2, 0}};
+            AdjacencyMap const adjacency_map = build_adjacency_map(edges);
+
+            REQUIRE(adjacency_map.size() == 3);
+            CHECK(adjacency_map.at(0) == std::unordered_set<uint64_t>{0, 2});
+            CHECK(adjacency_map.at(1) == std::unordered_set<uint64_t>{0, 1});
+            CHECK(adjacency_map.at(2) == std::unordered_set<uint64_t>{1, 2});
+        }
+
+        SUBCASE("Two edges, two nodes") {
+            auto edges = std::vector<BranchIdx>{{0, 1}, {0, 1}};
+            AdjacencyMap const adjacency_map = build_adjacency_map(edges);
+
+            REQUIRE(adjacency_map.size() == 2);
+            CHECK(adjacency_map.at(0) == std::unordered_set<uint64_t>{0, 1});
+            CHECK(adjacency_map.at(1) == std::unordered_set<uint64_t>{0, 1});
+        }
+
+        SUBCASE("Seven edges, five nodes") {
+            auto edges = std::vector<BranchIdx>{{3, 0}, {1, 0}, {2, 0}, {3, 2}, {1, 2}, {1, 4}, {3, 4}};
+            AdjacencyMap const adjacency_map = build_adjacency_map(edges);
+
+            REQUIRE(adjacency_map.size() == 5);
+            CHECK(adjacency_map.at(0) == std::unordered_set<uint64_t>{0, 1, 2});
+            CHECK(adjacency_map.at(1) == std::unordered_set<uint64_t>{1, 4, 5});
+            CHECK(adjacency_map.at(2) == std::unordered_set<uint64_t>{2, 3, 4});
+            CHECK(adjacency_map.at(3) == std::unordered_set<uint64_t>{0, 3, 6});
+            CHECK(adjacency_map.at(4) == std::unordered_set<uint64_t>{5, 6});
+        }
+    }
+
+    SUBCASE("Test forward elimination - elimination game") {
+        using enum EdgeEvent;
+        EliminationResult result{};
+
+        SUBCASE("One edge, two nodes, two real loads") {
+            auto const edges = std::vector<BranchIdx>{{0, 1}};
+            auto const node_loads = std::vector<DoubleComplex>{{-1.0, 0.0}, {1.0, 0.0}};
+            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);
+            CHECK(result.rhs == std::vector<DoubleComplex>{{1.0, 0.0}});
+            CHECK(result.free_edge_indices.empty());
+            REQUIRE(result.edges_history.size() == 1);
+            CHECK(result.edges_history[0].events == std::vector<EdgeEvent>{Deleted});
+            CHECK(result.edges_history[0].rows == std::vector<uint64_t>{0});
+        }
+
+        SUBCASE("Two edges, three nodes, two real loads") {
+            auto const edges = std::vector<BranchIdx>{{1, 0}, {1, 2}};
+            auto const node_loads = std::vector<DoubleComplex>{{-1.0, 0.0}, {1.0, 0.0}, {0.0, 0.0}};
+            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);
+            check_value(1, 1, 1, result.matrix);
+            REQUIRE(result.rhs.size() == 2);
+            CHECK(result.rhs == std::vector<DoubleComplex>{{-1.0, 0.0}, {0.0, 0.0}});
+            CHECK(result.free_edge_indices.empty());
+            REQUIRE(result.edges_history.size() == 2);
+            CHECK(result.edges_history[0].events == std::vector<EdgeEvent>{Deleted});
+            CHECK(result.edges_history[0].rows == std::vector<uint64_t>{0});
+            CHECK(result.edges_history[1].events == std::vector<EdgeEvent>{Deleted});
+            CHECK(result.edges_history[1].rows == std::vector<uint64_t>{1});
+        }
+
+        SUBCASE("Three edges, three nodes, two real loads") {
+            auto const edges = std::vector<BranchIdx>{{0, 1}, {1, 2}, {2, 0}};
+            auto const node_loads = std::vector<DoubleComplex>{{-1.0, 0.0}, {1.0, 0.0}, {0.0, 0.0}};
+            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);
+            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<DoubleComplex>{{1.0, 0.0}, {0.0, 0.0}});
+            REQUIRE(result.free_edge_indices.size() == 1);
+            CHECK(result.free_edge_indices == std::vector<uint64_t>{2});
+            REQUIRE(result.edges_history.size() == 3);
+            CHECK(result.edges_history[0].events == std::vector<EdgeEvent>{Deleted});
+            CHECK(result.edges_history[0].rows == std::vector<uint64_t>{0});
+            CHECK(result.edges_history[1].events == std::vector<EdgeEvent>{Replaced, Deleted});
+            CHECK(result.edges_history[1].rows == std::vector<uint64_t>{0, 1});
+            CHECK(result.edges_history[2].events == std::vector<EdgeEvent>{ContractedToPoint});
+            CHECK(result.edges_history[2].rows == std::vector<uint64_t>{1});
+        }
+
+        SUBCASE("Two edges, two nodes, two real loads") {
+            auto const edges = std::vector<BranchIdx>{{0, 1}, {0, 1}};
+            auto const node_loads = std::vector<DoubleComplex>{{-1.0, 0.0}, {1.0, 0.0}};
+            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);
+            check_value(1, 0, 1, result.matrix);
+            REQUIRE(result.rhs.size() == 1);
+            CHECK(result.rhs == std::vector<DoubleComplex>{{1.0, 0.0}});
+            REQUIRE(result.free_edge_indices.size() == 1);
+            CHECK(result.free_edge_indices == std::vector<uint64_t>{1});
+            REQUIRE(result.edges_history.size() == 2);
+            CHECK(result.edges_history[0].events == std::vector<EdgeEvent>{Deleted});
+            CHECK(result.edges_history[0].rows == std::vector<uint64_t>{0});
+            CHECK(result.edges_history[1].events == std::vector<EdgeEvent>{ContractedToPoint});
+            CHECK(result.edges_history[1].rows == std::vector<uint64_t>{0});
+        }
+
+        SUBCASE("Complex case with complex loads") {
+            auto const edges = std::vector<BranchIdx>{{3, 0}, {1, 0}, {2, 0}, {3, 2}, {1, 2}, {1, 4}, {3, 4}};
+            auto const node_loads =
+                std::vector<DoubleComplex>{{-1.0, -1.0}, {-1.0, -1.0}, {2.0, 2.0}, {0.0, 0.0}, {0.0, 0.0}};
+            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);
+            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<DoubleComplex>{{-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<uint64_t>{3, 4, 6});
+            REQUIRE(result.edges_history.size() == 7);
+            CHECK(result.edges_history[0].events == std::vector<EdgeEvent>{Deleted});
+            CHECK(result.edges_history[0].rows == std::vector<uint64_t>{0});
+            CHECK(result.edges_history[1].events == std::vector<EdgeEvent>{Replaced, Deleted});
+            CHECK(result.edges_history[1].rows == std::vector<uint64_t>{0, 1});
+            CHECK(result.edges_history[2].events == std::vector<EdgeEvent>{Replaced, Replaced, Deleted});
+            CHECK(result.edges_history[2].rows == std::vector<uint64_t>{0, 1, 2});
+            CHECK(result.edges_history[3].events == std::vector<EdgeEvent>{Replaced, ContractedToPoint});
+            CHECK(result.edges_history[3].rows == std::vector<uint64_t>{1, 2});
+            CHECK(result.edges_history[4].events == std::vector<EdgeEvent>{ContractedToPoint});
+            CHECK(result.edges_history[4].rows == std::vector<uint64_t>{2});
+            CHECK(result.edges_history[5].events == std::vector<EdgeEvent>{Replaced, Deleted});
+            CHECK(result.edges_history[5].rows == std::vector<uint64_t>{2, 3});
+            CHECK(result.edges_history[6].events == std::vector<EdgeEvent>{Replaced, Replaced, ContractedToPoint});
+            CHECK(result.edges_history[6].rows == std::vector<uint64_t>{1, 2, 3});
+        }
+    }
+
+    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<DoubleComplex>{{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<DoubleComplex>{{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<DoubleComplex>{{-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<DoubleComplex>{{-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<DoubleComplex>{{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<DoubleComplex>{{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<DoubleComplex>{{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<DoubleComplex>{{1.0, 0.0}});
+        }
+
+        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<DoubleComplex>{{-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<uint64_t>{3, 4, 6};
+            result.pivot_edge_indices = std::vector<uint64_t>{0, 1, 2, 5};
+
+            result.edges_history.resize(7);
+            result.edges_history[0].events = std::vector<EdgeEvent>{Deleted};
+            result.edges_history[0].rows = std::vector<uint64_t>{0};
+            result.edges_history[1].events = std::vector<EdgeEvent>{Replaced, Deleted};
+            result.edges_history[1].rows = std::vector<uint64_t>{0, 1};
+            result.edges_history[2].events = std::vector<EdgeEvent>{Replaced, Replaced, Deleted};
+            result.edges_history[2].rows = std::vector<uint64_t>{0, 1, 2};
+            result.edges_history[3].events = std::vector<EdgeEvent>{Replaced, ContractedToPoint};
+            result.edges_history[3].rows = std::vector<uint64_t>{1, 2};
+            result.edges_history[4].events = std::vector<EdgeEvent>{ContractedToPoint};
+            result.edges_history[4].rows = std::vector<uint64_t>{2};
+            result.edges_history[5].events = std::vector<EdgeEvent>{Replaced, Deleted};
+            result.edges_history[5].rows = std::vector<uint64_t>{2, 3};
+            result.edges_history[6].events = std::vector<EdgeEvent>{Replaced, Replaced, ContractedToPoint};
+            result.edges_history[6].rows = std::vector<uint64_t>{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<DoubleComplex>{{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