resolve comments

Signed-off-by: Santiago Figueroa Manrique <figueroa1395@gmail.com>

Author: figueroa1395

power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp

@@ -320,24 +320,28 @@ inline std::vector<std::vector<DoubleComplex>> set_projection_system(Idx free_in
 };
 
 inline void gauss_elimination(std::vector<std::vector<DoubleComplex>>& system) {
+    auto const system_size = narrow_cast<Idx>(std::ssize(system));
 
-    auto const system_size = std::ssize(system);
-
-    for (Idx column = 0; column < system_size; column++) {
-        for (Idx row = column + 1; row < system_size; row++) {
+    // we skip pivoting since the matrix system is mostly diagonally dominant
+    // and given the real power systems topology of super nodes, this should not introduce
+    // any numerical instabilities in practice
 
+    // forward elimination
+    for (Idx column : IdxRange{system_size}) {
+        for (Idx row : IdxRange{column + Idx{1}, system_size}) {
             system[row][column] = -system[row][column] / system[column][column];
 
-            for (Idx column_part = column + 1; column_part < system_size + 1; column_part++) {
+            for (Idx column_part : IdxRange{column + Idx{1}, system_size + Idx{1}}) {
                 system[row][column_part] += system[row][column] * system[column][column_part];
             }
         }
     }
 
+    // backward substitution
     system[system_size - 1][system_size] /= system[system_size - 1][system_size - 1];
     for (Idx row = system_size - 1; row-- > 0;) {
         auto element_sum = DoubleComplex{};
-        for (Idx column = row + 1; column < system_size; column++) {
+        for (Idx column : IdxRange{row + Idx{1}, system_size}) {
             element_sum -= system[row][column] * system[column][system_size];
         }
         system[row][system_size] = (system[row][system_size] + element_sum) / system[row][row];
@@ -346,23 +350,18 @@ inline void gauss_elimination(std::vector<std::vector<DoubleComplex>>& system) {
 
 inline std::vector<DoubleComplex> compute_internal_loads(SolutionSet& solution_set,
                                                          std::vector<std::vector<DoubleComplex>>& system) {
-
-    std::vector<DoubleComplex> internal_loads{};
-
     auto const number_of_rows = narrow_cast<Idx>(solution_set.extended_rhs.size());
     auto const number_of_columns = narrow_cast<Idx>(system.size());
-
-    internal_loads.resize(number_of_rows);
+    std::vector<DoubleComplex> internal_loads(number_of_rows);
 
     for (auto row : IdxRange{number_of_rows}) {
-
         internal_loads[row] = solution_set.extended_rhs[row];
         auto sum_value = DoubleComplex{};
         for (auto column : IdxRange{number_of_columns}) {
-            auto const value = solution_set.dfs_matrix.get_value(row, column);
-            if (value.has_value()) {
-                sum_value += static_cast<DoubleComplex>(value.value()) * system[column].back();
-            }
+            solution_set.dfs_matrix.get_value(row, column).transform([&sum_value, &system, &column](IntS value) {
+                sum_value += static_cast<DoubleComplex>(value) * system[column].back();
+                return value;
+            });
         }
         internal_loads[row] -= sum_value;
     }