Skip to content

Commit a9940b8

Browse files
figueroa1395figueroa1395
committed
code quality
Signed-off-by: Santiago Figueroa Manrique <figueroa1395@gmail.com>
1 parent d809098 commit a9940b8

2 files changed

Lines changed: 55 additions & 53 deletions

File tree

‎power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp‎

Lines changed: 15 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -319,46 +319,47 @@ inline std::vector<std::vector<DoubleComplex>> set_projection_system(Idx free_in
319319
return projection_system;
320320
};
321321

322-
inline void gauss_elimination(std::vector<std::vector<DoubleComplex>>& system) {
322+
inline void naive_gauss_elimination(std::vector<std::vector<DoubleComplex>>& system) {
323323
auto const system_size = narrow_cast<Idx>(std::ssize(system));
324324

325325
// we skip pivoting since the matrix system is mostly diagonally dominant
326326
// and given the real power systems topology of super nodes, this should not introduce
327327
// any numerical instabilities in practice
328328

329329
// forward elimination
330-
for (Idx column : IdxRange{system_size}) {
331-
for (Idx row : IdxRange{column + Idx{1}, system_size}) {
332-
system[row][column] = -system[row][column] / system[column][column];
333-
334-
for (Idx column_part : IdxRange{column + Idx{1}, system_size + Idx{1}}) {
335-
system[row][column_part] += system[row][column] * system[column][column_part];
330+
for (Idx const column : IdxRange{system_size}) {
331+
auto const pivot = system[column][column];
332+
for (Idx const row : IdxRange{column + Idx{1}, system_size}) {
333+
auto multiplier = -system[row][column] / pivot;
334+
system[row][column] = multiplier;
335+
for (Idx const column_part : IdxRange{column + Idx{1}, system_size + Idx{1}}) {
336+
system[row][column_part] += multiplier * system[column][column_part];
336337
}
337338
}
338339
}
339340

340341
// backward substitution
341342
system[system_size - 1][system_size] /= system[system_size - 1][system_size - 1];
342-
for (Idx row = system_size - 1; row-- > 0;) {
343+
for (Idx const row : IdxRange{system_size - 1} | std::views::reverse) {
343344
auto element_sum = DoubleComplex{};
344-
for (Idx column : IdxRange{row + Idx{1}, system_size}) {
345+
for (Idx const column : IdxRange{row + Idx{1}, system_size}) {
345346
element_sum -= system[row][column] * system[column][system_size];
346347
}
347348
system[row][system_size] = (system[row][system_size] + element_sum) / system[row][row];
348349
}
349350
};
350351

351-
inline std::vector<DoubleComplex> compute_internal_loads(SolutionSet& solution_set,
352-
std::vector<std::vector<DoubleComplex>>& system) {
352+
inline std::vector<DoubleComplex> compute_internal_loads(SolutionSet const& solution_set,
353+
std::vector<std::vector<DoubleComplex>> const& system) {
353354
auto const number_of_rows = narrow_cast<Idx>(solution_set.extended_rhs.size());
354355
auto const number_of_columns = narrow_cast<Idx>(system.size());
355356
std::vector<DoubleComplex> internal_loads(number_of_rows);
356357

357-
for (auto row : IdxRange{number_of_rows}) {
358+
for (auto const row : IdxRange{number_of_rows}) {
358359
internal_loads[row] = solution_set.extended_rhs[row];
359360
auto sum_value = DoubleComplex{};
360-
for (auto column : IdxRange{number_of_columns}) {
361-
solution_set.dfs_matrix.get_value(row, column).transform([&sum_value, &system, &column](IntS value) {
361+
for (auto const column : IdxRange{number_of_columns}) {
362+
solution_set.dfs_matrix.get_value(row, column).transform([&sum_value, &system, column](IntS value) {
362363
sum_value += static_cast<DoubleComplex>(value) * system[column].back();
363364
return value;
364365
});

‎tests/cpp_unit_tests/test_link_solver.cpp‎

Lines changed: 40 additions & 39 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44

55
#include <power_grid_model/calculation_parameters.hpp>
66
#include <power_grid_model/common/common.hpp>
7+
#include <power_grid_model/common/counting_iterator.hpp>
78
#include <power_grid_model/common/typing.hpp>
89
#include <power_grid_model/link_solver.hpp>
910

@@ -523,45 +524,40 @@ TEST_CASE("Test the link solver algorithm") {
523524
}
524525

525526
SUBCASE("Testing the gauss elimination routine") {
526-
auto compare_systems = [](std::vector<std::vector<DoubleComplex>>& system,
527-
std::vector<std::vector<DoubleComplex>>& test_system) {
528-
auto const system_size = narrow_cast<Idx>(system.size());
529-
double element_sum{};
530-
531-
for (Idx column = 0; column < system_size + 1; column++) {
532-
for (Idx row = 0; row < system_size; row++) {
533-
element_sum += abs(system[row][column] - test_system[row][column]);
527+
auto compare_systems = [](std::vector<std::vector<DoubleComplex>> const& solution,
528+
std::vector<std::vector<DoubleComplex>> const& reference, Idx col_number,
529+
Idx row_number, double tolerance) {
530+
for (Idx const col : IdxRange{col_number}) {
531+
for (Idx const row : IdxRange{row_number}) {
532+
CHECK(solution[row][col].real() == doctest::Approx(reference[row][col].real()).epsilon(tolerance));
533+
CHECK(solution[row][col].imag() == doctest::Approx(reference[row][col].imag()).epsilon(tolerance));
534534
}
535535
}
536-
537-
return element_sum;
538536
};
539537

540-
auto compare_vectors = [](std::vector<std::vector<DoubleComplex>>& system,
541-
std::vector<DoubleComplex>& test_solution) {
542-
auto const system_size = narrow_cast<Idx>(system.size());
543-
double element_sum{};
544-
545-
for (Idx row = 0; row < system_size; row++) {
546-
element_sum += abs(system[row][system_size] - test_solution[row]);
538+
auto compare_vectors = [](std::vector<std::vector<DoubleComplex>>& solution,
539+
std::vector<DoubleComplex>& reference, double tolerance) {
540+
auto const solution_size = narrow_cast<Idx>(solution.size());
541+
// we only test against the last column of the solution as this lambda is intended for large edge test cases
542+
// this makes our lifes easier
543+
for (Idx const row : IdxRange{solution_size}) {
544+
CHECK(solution[row][solution_size].real() == doctest::Approx(reference[row].real()).epsilon(tolerance));
545+
CHECK(solution[row][solution_size].imag() == doctest::Approx(reference[row].imag()).epsilon(tolerance));
547546
}
548-
549-
return element_sum;
550547
};
551548

552549
SUBCASE("Linear system of the complex case") {
553550
std::vector<std::vector<DoubleComplex>> system = {{{3, 0}, {1, 0}, {1, 0}, {2, 2}},
554551
{{1, 0}, {3, 0}, {-1, 0}, {3, 3}},
555552
{{1, 0}, {-1, 0}, {4, 0}, {-1, -1}}};
556-
557-
gauss_elimination(system);
558-
553+
naive_gauss_elimination(system);
559554
std::vector<std::vector<DoubleComplex>> test_system = {
560555
{{3, 0}, {1, 0}, {1, 0}, {0.458333, 0.458333}},
561556
{{-0.333333, 0}, {2.66667, 0}, {-1.33333, 0}, {0.791667, 0.791667}},
562557
{{-0.333333, 0}, {0.5, -0}, {3, 0}, {-0.166667, -0.166667}}};
563558

564-
CHECK(compare_systems(system, test_system) < 1.e-5);
559+
// the tolerance is set to 1e-5 because that's the test system tolerance
560+
compare_systems(system, test_system, Idx{4}, Idx{3}, 1e-5);
565561
}
566562

567563
SUBCASE("A system that consisng of a 15 X 15 matrix with externally randomly generated elements") {
@@ -811,33 +807,37 @@ TEST_CASE("Test the link solver algorithm") {
811807
{-0.03845678, 0}, {-0.00652489, 0}, {-0.08356931, 0}, {0.05730963, 0}, {0.01390954, 0},
812808
{0.026622, 0}, {0.04469859, 0}, {-0.00946348, 0}, {0.08945877, 0}, {0.00377452, 0}};
813809

814-
gauss_elimination(system);
810+
naive_gauss_elimination(system);
815811

816-
CHECK(compare_vectors(system, test_solution) < 1.e-7);
812+
// error tolerance is increased as this is a stress test
813+
// the test system solution now includes for significant digits for this reason
814+
// this is important because we are using naive_gauss_elimination, which skips pivoting
815+
// so this test makes sure that assumption holds
816+
compare_vectors(system, test_solution, 1e-7);
817817
}
818818
}
819819

820820
SUBCASE("Testing the compute_internal_loads_routine") {
821-
822-
auto generate_input_solution_set = [](std::span<const IntS> data, std::span<const Idx> row,
823-
std::span<const Idx> col, Idx row_number, Idx col_number) {
821+
auto generate_input_solution_set = [](std::span<IntS const> data, std::span<Idx const> rows,
822+
std::span<Idx const> cols, Idx row_number, Idx col_number) {
824823
SolutionSet solution_set{};
825824
solution_set.dfs_matrix.prepare(row_number, col_number);
826825
for (auto idx = size_t{0}; idx < data.size(); ++idx) {
827-
solution_set.dfs_matrix.set_value(data[idx], row[idx], col[idx]);
826+
solution_set.dfs_matrix.set_value(data[idx], rows[idx], cols[idx]);
828827
}
829828
return solution_set;
830829
};
831830

832-
auto compare_vectors = [](std::vector<DoubleComplex>& load, std::vector<DoubleComplex>& test_load) {
833-
auto const load_size = narrow_cast<Idx>(load.size());
834-
double element_sum{};
831+
auto compare_vectors = [](std::vector<DoubleComplex>& loads, std::vector<DoubleComplex>& test_loads,
832+
double tolerance) {
833+
auto const loads_size = narrow_cast<Idx>(loads.size());
834+
auto const test_loads_size = narrow_cast<Idx>(test_loads.size());
835+
REQUIRE(loads_size == test_loads_size);
835836

836-
for (Idx idx = 0; idx < load_size; idx++) {
837-
element_sum += abs(load[idx] - test_load[idx]);
837+
for (Idx const idx : IdxRange(loads_size)) {
838+
CHECK(loads[idx].real() == doctest::Approx(test_loads[idx].real()).epsilon(tolerance));
839+
CHECK(loads[idx].imag() == doctest::Approx(test_loads[idx].imag()).epsilon(tolerance));
838840
}
839-
840-
return element_sum;
841841
};
842842

843843
SUBCASE("Complex case with complex loads") {
@@ -859,11 +859,11 @@ TEST_CASE("Test the link solver algorithm") {
859859
std::vector<DoubleComplex> test_loads = {
860860
{-0.291667, -0.291667}, {0.0416667, 0.0416667}, {-0.75, -0.75}, {0.458333, 0.458333},
861861
{0.791667, 0.791667}, {0.166667, 0.166667}, {-0.166667, -0.166667}};
862-
CHECK(compare_vectors(internal_loads, test_loads) < 1.e-5);
862+
863+
compare_vectors(internal_loads, test_loads, 1e-5);
863864
}
864865

865866
SUBCASE("Four edges, four nodes, two real loads") {
866-
867867
std::vector<IntS> data = {1, 1, -1};
868868
std::vector<Idx> row = {1, 2, 3};
869869
std::vector<Idx> col = {0, 0, 0};
@@ -876,7 +876,8 @@ TEST_CASE("Test the link solver algorithm") {
876876
std::vector<DoubleComplex> internal_loads = compute_internal_loads(solution_set, test_system);
877877

878878
std::vector<DoubleComplex> test_loads = {{1, 0}, {-0.333333, -0}, {-0.333333, -0}, {-0.666667, 0}};
879-
CHECK(compare_vectors(internal_loads, test_loads) < 1.e-5);
879+
880+
compare_vectors(internal_loads, test_loads, 1.e-5);
880881
}
881882
}
882883
}

0 commit comments

Comments
 (0)