From 2c7e430258d1e300609f6ed72b53f09885330875 Mon Sep 17 00:00:00 2001 From: Santiago Figueroa Manrique Date: Tue, 16 Sep 2025 11:28:14 +0200 Subject: [PATCH 1/6] working implementation Signed-off-by: Santiago Figueroa Manrique --- .../iterative_linear_se_solver.hpp | 4 +- .../math_solver/newton_raphson_se_solver.hpp | 5 ++- .../math_solver/observability.hpp | 18 +++++---- .../math_solver/sparse_lu_solver.hpp | 40 ++++++++++--------- .../cpp_unit_tests/test_sparse_lu_solver.cpp | 39 +++++++++++++++--- 5 files changed, 72 insertions(+), 34 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp index 121f6006fe..e4c1dbdcc4 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp @@ -100,12 +100,14 @@ template class IterativeLinearSESolver { MeasuredValues const measured_values{y_bus.shared_topology(), input}; auto const observability_result = observability_check(measured_values, y_bus.math_topology(), y_bus.y_bus_structure()); + std::pair> const possibly_ill_conditioned_pivots{ + observability_result.use_perturbation(), observability_result.row_is_possibly_ill_conditioned}; // prepare matrix sub_timer = Timer{log, LogEvent::prepare_matrix_including_prefactorization}; prepare_matrix(y_bus, measured_values); // prefactorize - sparse_solver_.prefactorize(data_gain_, perm_, observability_result.use_perturbation()); + sparse_solver_.prefactorize(data_gain_, perm_, possibly_ill_conditioned_pivots); // initialize voltage with initial angle sub_timer = Timer{log, LogEvent::initialize_voltages}; // TODO(mgovers): make scoped subtimers diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp index 54b9ba9c48..e2993fd228 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp @@ -166,8 +166,11 @@ template class NewtonRaphsonSESolver { // preprocess measured value sub_timer = Timer{log, LogEvent::preprocess_measured_value}; MeasuredValues const measured_values{y_bus.shared_topology(), input}; + auto const observability_result = observability_check(measured_values, y_bus.math_topology(), y_bus.y_bus_structure()); + std::pair> const possibly_ill_conditioned_pivots{ + observability_result.use_perturbation(), observability_result.row_is_possibly_ill_conditioned}; // initialize voltage with initial angle sub_timer = Timer{log, LogEvent::initialize_voltages}; @@ -184,7 +187,7 @@ template class NewtonRaphsonSESolver { // solve with prefactorization sub_timer = Timer{log, LogEvent::solve_sparse_linear_equation}; sparse_solver_.prefactorize_and_solve(data_gain_, perm_, delta_x_rhs_, delta_x_rhs_, - observability_result.use_perturbation()); + possibly_ill_conditioned_pivots); sub_timer = Timer{log, LogEvent::iterate_unknown}; max_dev = iterate_unknown(output.u, measured_values); }; diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp index 0890f39627..d7068401a2 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp @@ -15,7 +15,7 @@ namespace detail { struct ObservabilitySensorsResult { std::vector flow_sensors; std::vector voltage_phasor_sensors; - bool is_possibly_ill_conditioned{false}; + std::vector row_is_possibly_ill_conditioned; }; // count flow and voltage phasor sensors for the observability check @@ -26,7 +26,7 @@ struct ObservabilitySensorsResult { // return a ObservabilitySensorsResult struct with // a vector of flow sensor count // a vector of voltage phasor sensor count -// a boolean indicating if the system is possibly ill-conditioned +// a vector of boolean indicating if each row is possibly ill-conditioned template ObservabilitySensorsResult count_observability_sensors(MeasuredValues const& measured_values, MathModelTopology const& topo, @@ -35,7 +35,7 @@ ObservabilitySensorsResult count_observability_sensors(MeasuredValues const ObservabilitySensorsResult result{.flow_sensors = std::vector(y_bus_structure.row_indptr.back(), 0), .voltage_phasor_sensors = std::vector(n_bus, 0), - .is_possibly_ill_conditioned = false}; + .row_is_possibly_ill_conditioned = std::vector(n_bus, 0)}; auto has_flow_sensor = [&measured_values](Idx branch) { return measured_values.has_branch_from_power(branch) || measured_values.has_branch_to_power(branch) || @@ -81,7 +81,7 @@ ObservabilitySensorsResult count_observability_sensors(MeasuredValues const } // the system could be ill-conditioned if there is no flow sensor for one bus, except the last bus if (!has_at_least_one_sensor && row != n_bus - 1) { - result.is_possibly_ill_conditioned = true; + result.row_is_possibly_ill_conditioned[row] = 1; } } return result; @@ -195,8 +195,11 @@ inline bool sufficient_observability_condition(YBusStructure const& y_bus_struct struct ObservabilityResult { bool is_observable{false}; - bool is_possibly_ill_conditioned{false}; - constexpr bool use_perturbation() const { return is_possibly_ill_conditioned && is_observable; } + std::vector row_is_possibly_ill_conditioned{}; + bool use_perturbation() const { + return std::ranges::any_of(row_is_possibly_ill_conditioned, [](int8_t element) { return element == 1; }) && + is_observable; + } }; template @@ -227,7 +230,8 @@ inline ObservabilityResult observability_check(MeasuredValues const& measur } // ToDo(JGuo): meshed network will require a different treatment return ObservabilityResult{.is_observable = is_necessary_condition_met && is_sufficient_condition_met, - .is_possibly_ill_conditioned = observability_sensors.is_possibly_ill_conditioned}; + .row_is_possibly_ill_conditioned = + observability_sensors.row_is_possibly_ill_conditioned}; } } // namespace power_grid_model::math_solver diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp index 15b9d83033..6373261851 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp @@ -24,9 +24,9 @@ constexpr double epsilon_sqrt = 1.49011745e-8; // sqrt(epsilon), sqrt does // it will get modified if perturbation happens // also has_pivot_perturbation will be updated if perturbation happens template -inline void perturb_pivot_if_needed(double perturb_threshold, Scalar& value, double& abs_value, +inline void perturb_pivot_if_needed(int8_t pertub_pivot, double perturb_threshold, Scalar& value, double& abs_value, bool& has_pivot_perturbation) { - if (abs_value < perturb_threshold) { + if (pertub_pivot == 1 && abs_value < perturb_threshold) { Scalar const scale = (abs_value == 0.0) ? Scalar{1.0} : (value / abs_value); value = scale * perturb_threshold; has_pivot_perturbation = true; @@ -60,8 +60,9 @@ template class DenseLUFactor { // put pivot perturbation in has_pivot_perturbation in place template static void factorize_block_in_place(Eigen::MatrixBase&& matrix, BlockPerm& block_perm, - double perturb_threshold, bool use_pivot_perturbation, - bool& has_pivot_perturbation) + double perturb_threshold, + std::pair> possibly_ill_conditioned_pivots, + Idx block_node_idx, bool& has_pivot_perturbation) requires(std::same_as && rk2_tensor && (Derived::RowsAtCompileTime == size) && (Derived::ColsAtCompileTime == size)) { @@ -84,7 +85,7 @@ template class DenseLUFactor { assert(col_biggest_eigen + pivot < size); // check absolute singular matrix - if (biggest_score == 0.0 && !use_pivot_perturbation) { + if (biggest_score == 0.0 && !possibly_ill_conditioned_pivots.first) { // pivot perturbation not possible, cannot proceed // set identity permutation and break the loop for (int8_t remaining_rows_cols = pivot; remaining_rows_cols != size; ++remaining_rows_cols) { @@ -96,8 +97,10 @@ template class DenseLUFactor { // perturb pivot if needed double abs_pivot = sqrt(biggest_score); - perturb_pivot_if_needed(perturb_threshold, matrix(row_biggest, col_biggest), abs_pivot, - has_pivot_perturbation); + int8_t const possibly_ill_conditioned_pivot = + possibly_ill_conditioned_pivots.first ? possibly_ill_conditioned_pivots.second[block_node_idx] : 0; + perturb_pivot_if_needed(possibly_ill_conditioned_pivot, perturb_threshold, matrix(row_biggest, col_biggest), + abs_pivot, has_pivot_perturbation); max_pivot = std::max(max_pivot, abs_pivot); // swap rows and columns @@ -214,8 +217,8 @@ template class SparseLUSolver { prefactorize_and_solve(std::vector& data, // matrix data, factorize in-place BlockPermArray& block_perm_array, // pre-allocated permutation array, will be overwritten std::vector const& rhs, std::vector& x, - bool use_pivot_perturbation = false) { - prefactorize(data, block_perm_array, use_pivot_perturbation); + std::pair> possibly_ill_conditioned_pivots = {false, {}}) { + prefactorize(data, block_perm_array, possibly_ill_conditioned_pivots); // call solve with const method solve_with_prefactorized_matrix((std::vector const&)data, block_perm_array, rhs, x); } @@ -239,9 +242,10 @@ template class SparseLUSolver { // fill-ins should be pre-allocated with zero // block permutation array should be pre-allocated void prefactorize(std::vector& data, BlockPermArray& block_perm_array, - bool use_pivot_perturbation = false) { + std::pair> possibly_ill_conditioned_pivots = {false, {}}) { reset_matrix_cache(); - if (use_pivot_perturbation) { + if (possibly_ill_conditioned_pivots.first) { + assert(static_cast(possibly_ill_conditioned_pivots.second.size()) == size_); initialize_pivot_perturbation(data); } double const perturb_threshold = epsilon_perturbation * matrix_norm_; @@ -268,15 +272,17 @@ template class SparseLUSolver { // use machine precision by default // record block permutation LUFactor::factorize_block_in_place(lu_matrix[pivot_idx].matrix(), block_perm_array[pivot_row_col], - perturb_threshold, use_pivot_perturbation, + perturb_threshold, possibly_ill_conditioned_pivots, + pivot_row_col, // Pass the node index has_pivot_perturbation_); return block_perm_array[pivot_row_col]; } else { - if (use_pivot_perturbation) { + if (possibly_ill_conditioned_pivots.first) { // use machine precision by default // record pivot perturbation double abs_pivot = cabs(lu_matrix[pivot_idx]); - perturb_pivot_if_needed(perturb_threshold, lu_matrix[pivot_idx], abs_pivot, + perturb_pivot_if_needed(possibly_ill_conditioned_pivots.second[pivot_row_col], + perturb_threshold, lu_matrix[pivot_idx], abs_pivot, has_pivot_perturbation_); } if (!is_normal(lu_matrix[pivot_idx])) { @@ -543,15 +549,11 @@ template class SparseLUSolver { // 3. take the maximum of all the sums matrix_norm_ = 0.0; auto const& row_indptr = *row_indptr_; - auto const& col_indices = *col_indices_; + // auto const& col_indices = *col_indices_; for (Idx row = 0; row != size_; ++row) { // calculate the sum of the norms of the blocks in the row double row_norm = 0.0; for (Idx idx = row_indptr[row]; idx != row_indptr[row + 1]; ++idx) { - // skip diagonal - if (col_indices[idx] == row) { - continue; - } if constexpr (is_block) { row_norm += cabs(data[idx]).rowwise().sum().maxCoeff(); } else { diff --git a/tests/cpp_unit_tests/test_sparse_lu_solver.cpp b/tests/cpp_unit_tests/test_sparse_lu_solver.cpp index d1915de718..568c6f14a1 100644 --- a/tests/cpp_unit_tests/test_sparse_lu_solver.cpp +++ b/tests/cpp_unit_tests/test_sparse_lu_solver.cpp @@ -158,11 +158,32 @@ TEST_CASE("LU solver with ill-conditioned system") { Idx perm{}; SUBCASE("Error without perturbation") { - CHECK_THROWS_AS(solver.prefactorize(data, perm, false), SparseMatrixError); + CHECK_THROWS_AS(solver.prefactorize(data, perm, {false, {}}), SparseMatrixError); } - SUBCASE("Success with perturbation") { - CHECK_NOTHROW(solver.prefactorize(data, perm, true)); + // first row is ill conditioned, others are fine: the cases where the first row is not perturbed should fail + SUBCASE("Error with wrong perturbation") { + // generate all possible combinations with 4 possibly ill conditioned rows (16 total) + for (int8_t combination = 0; combination < 16; ++combination) { + std::string binary_possibly_ill_conditioned_rows = std::bitset<4>(combination).to_string(); + // we test combinations where the first row is not perturbed + if (binary_possibly_ill_conditioned_rows[0] == '1') { + continue; + } + std::vector possibly_ill_conditioned_rows(4); + std::ranges::transform(binary_possibly_ill_conditioned_rows, possibly_ill_conditioned_rows.begin(), + [](char c) -> int8_t { return c == '1' ? 1 : 0; }); + CAPTURE(combination); + CAPTURE(possibly_ill_conditioned_rows); + // first row is not perturbed: should fail + CHECK_THROWS_AS(solver.prefactorize(data, perm, {true, possibly_ill_conditioned_rows}), + SparseMatrixError); + } + } + + SUBCASE("Success with correct perturbation") { + // Only the first row is ill conditioned: possibly trigger perturbation on the correct row + CHECK_NOTHROW(solver.prefactorize(data, perm, {true, {1, 0, 0, 0}})); solver.solve_with_prefactorized_matrix(data, perm, rhs, x); check_result(x, x_ref); } @@ -186,11 +207,17 @@ TEST_CASE("LU solver with ill-conditioned system") { SparseLUSolver solver{row_indptr, col_indices, diag_lu}; SUBCASE("Error without perturbation") { - CHECK_THROWS_AS(solver.prefactorize(data, block_perm, false), SparseMatrixError); + CHECK_THROWS_AS(solver.prefactorize(data, block_perm, {false, {}}), SparseMatrixError); + } + + SUBCASE("Error with wrong perturbation") { + // first block is ill conditioned, second is fine: trigger perturbation on the wrong block + CHECK_THROWS_AS(solver.prefactorize(data, block_perm, {true, {0, 1}}), SparseMatrixError); } - SUBCASE("Success with perturbation") { - CHECK_NOTHROW(solver.prefactorize(data, block_perm, true)); + SUBCASE("Success with correct perturbation") { + // first block is ill conditioned, second is fine: trigger perturbation on the correct block + CHECK_NOTHROW(solver.prefactorize(data, block_perm, {true, {1, 0}})); solver.solve_with_prefactorized_matrix(data, block_perm, rhs, x); check_result(x, x_ref); } From 6a24417c7153ecfc51346ef3ce32006f49f46f64 Mon Sep 17 00:00:00 2001 From: Santiago Figueroa Manrique Date: Tue, 16 Sep 2025 12:50:27 +0200 Subject: [PATCH 2/6] minor Signed-off-by: Santiago Figueroa Manrique --- tests/cpp_unit_tests/test_sparse_lu_solver.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/cpp_unit_tests/test_sparse_lu_solver.cpp b/tests/cpp_unit_tests/test_sparse_lu_solver.cpp index 568c6f14a1..978a74abd6 100644 --- a/tests/cpp_unit_tests/test_sparse_lu_solver.cpp +++ b/tests/cpp_unit_tests/test_sparse_lu_solver.cpp @@ -7,6 +7,8 @@ #include +#include + namespace power_grid_model::math_solver { namespace { using lu_trait_double = math_solver::sparse_lu_entry_trait; From 9c49eaaa027f3a2a273598b30ad081730b6ffc91 Mon Sep 17 00:00:00 2001 From: Santiago Figueroa Manrique Date: Tue, 16 Sep 2025 13:35:37 +0200 Subject: [PATCH 3/6] clang tidy Signed-off-by: Santiago Figueroa Manrique --- .../power_grid_model/math_solver/sparse_lu_solver.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp index 6373261851..76ceed780f 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp @@ -97,8 +97,9 @@ template class DenseLUFactor { // perturb pivot if needed double abs_pivot = sqrt(biggest_score); - int8_t const possibly_ill_conditioned_pivot = - possibly_ill_conditioned_pivots.first ? possibly_ill_conditioned_pivots.second[block_node_idx] : 0; + int8_t const possibly_ill_conditioned_pivot = possibly_ill_conditioned_pivots.first + ? possibly_ill_conditioned_pivots.second[block_node_idx] + : int8_t{0}; perturb_pivot_if_needed(possibly_ill_conditioned_pivot, perturb_threshold, matrix(row_biggest, col_biggest), abs_pivot, has_pivot_perturbation); max_pivot = std::max(max_pivot, abs_pivot); From c1fa814ce36dd8868e8a85c5d9688433991a1b63 Mon Sep 17 00:00:00 2001 From: Santiago Figueroa Manrique Date: Wed, 17 Sep 2025 16:42:40 +0200 Subject: [PATCH 4/6] dirty implementation of max_pivot based in original matrix Signed-off-by: Santiago Figueroa Manrique --- .../power_grid_model/math_solver/sparse_lu_solver.hpp | 10 +++++++--- .../ill-conditioned-by-line-meshed/params.json | 4 ---- .../ill-conditioned-by-link-meshed/params.json | 4 ---- 3 files changed, 7 insertions(+), 11 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp index 76ceed780f..87686cd5be 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp @@ -62,7 +62,7 @@ template class DenseLUFactor { static void factorize_block_in_place(Eigen::MatrixBase&& matrix, BlockPerm& block_perm, double perturb_threshold, std::pair> possibly_ill_conditioned_pivots, - Idx block_node_idx, bool& has_pivot_perturbation) + Idx block_node_idx, bool& has_pivot_perturbation, double max_pivot_original) requires(std::same_as && rk2_tensor && (Derived::RowsAtCompileTime == size) && (Derived::ColsAtCompileTime == size)) { @@ -136,7 +136,7 @@ template class DenseLUFactor { // throw SparseMatrixError if the matrix is ill-conditioned // only check condition number if pivot perturbation is not used - double const pivot_threshold = has_pivot_perturbation ? 0.0 : epsilon * max_pivot; + double const pivot_threshold = has_pivot_perturbation ? 0.0 : epsilon * max_pivot_original; for (int8_t pivot = 0; pivot != size; ++pivot) { if (cabs(matrix(pivot, pivot)) < pivot_threshold || !is_normal(matrix(pivot, pivot))) { throw SparseMatrixError{}; // can not specify error code @@ -275,7 +275,7 @@ template class SparseLUSolver { LUFactor::factorize_block_in_place(lu_matrix[pivot_idx].matrix(), block_perm_array[pivot_row_col], perturb_threshold, possibly_ill_conditioned_pivots, pivot_row_col, // Pass the node index - has_pivot_perturbation_); + has_pivot_perturbation_, max_pivot_original_); return block_perm_array[pivot_row_col]; } else { if (possibly_ill_conditioned_pivots.first) { @@ -420,6 +420,7 @@ template class SparseLUSolver { // cache value for pivot perturbation for the factorize step bool has_pivot_perturbation_{false}; double matrix_norm_{}; + double max_pivot_original_{}; std::optional> original_matrix_; // cache value for iterative refinement for the solve step std::optional> dx_; @@ -549,6 +550,7 @@ template class SparseLUSolver { // 2. sum all norms of the blocks per row, except the diagonal block // 3. take the maximum of all the sums matrix_norm_ = 0.0; + max_pivot_original_ = 0.0; auto const& row_indptr = *row_indptr_; // auto const& col_indices = *col_indices_; for (Idx row = 0; row != size_; ++row) { @@ -557,8 +559,10 @@ template class SparseLUSolver { for (Idx idx = row_indptr[row]; idx != row_indptr[row + 1]; ++idx) { if constexpr (is_block) { row_norm += cabs(data[idx]).rowwise().sum().maxCoeff(); + max_pivot_original_ = std::max(max_pivot_original_, cabs(data[idx]).maxCoeff()); } else { row_norm += cabs(data[idx]); + max_pivot_original_ = std::max(max_pivot_original_, cabs(data[idx])); } } matrix_norm_ = std::max(matrix_norm_, row_norm); diff --git a/tests/data/state_estimation/ill-conditioned-by-line-meshed/params.json b/tests/data/state_estimation/ill-conditioned-by-line-meshed/params.json index 3358a01ad8..9e558925ef 100644 --- a/tests/data/state_estimation/ill-conditioned-by-line-meshed/params.json +++ b/tests/data/state_estimation/ill-conditioned-by-line-meshed/params.json @@ -3,9 +3,5 @@ "rtol": 1e-8, "atol": { "default": 1e-8 - }, - "raises": { - "raises": "SparseMatrixError", - "reason": "Observability check for meshed grids is not yet implemented. Also, pivot perturbation wouldn't be triggered as Matrix norm is smaller than the max pivot. See https://github.com/PowerGridModel/power-grid-model/pull/1118#issuecomment-3269339573" } } \ No newline at end of file diff --git a/tests/data/state_estimation/ill-conditioned-by-link-meshed/params.json b/tests/data/state_estimation/ill-conditioned-by-link-meshed/params.json index 6dd5b9a8dc..88cb11fee3 100644 --- a/tests/data/state_estimation/ill-conditioned-by-link-meshed/params.json +++ b/tests/data/state_estimation/ill-conditioned-by-link-meshed/params.json @@ -6,9 +6,5 @@ "rtol": 1e-8, "atol": { "default": 1e-8 - }, - "raises": { - "raises": "SparseMatrixError", - "reason": "Observability check for meshed grids is not yet implemented. Also, pivot perturbation wouldn't be triggered as Matrix norm is smaller than the max pivot. See https://github.com/PowerGridModel/power-grid-model/pull/1118#issuecomment-3269339573" } } \ No newline at end of file From a9f0158e14dd6f6abbf1ae439865a4c098bb4e50 Mon Sep 17 00:00:00 2001 From: Santiago Figueroa Manrique Date: Wed, 17 Sep 2025 16:56:05 +0200 Subject: [PATCH 5/6] fix quick implementation Signed-off-by: Santiago Figueroa Manrique --- .../math_solver/sparse_lu_solver.hpp | 21 ++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp index 87686cd5be..e6430fa565 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp @@ -245,6 +245,7 @@ template class SparseLUSolver { void prefactorize(std::vector& data, BlockPermArray& block_perm_array, std::pair> possibly_ill_conditioned_pivots = {false, {}}) { reset_matrix_cache(); + calculate_max_pivot(data); if (possibly_ill_conditioned_pivots.first) { assert(static_cast(possibly_ill_conditioned_pivots.second.size()) == size_); initialize_pivot_perturbation(data); @@ -547,28 +548,38 @@ template class SparseLUSolver { // calculate the block-wise non-diagonal infinite norm of the matrix // that is: // 1. calculate the infinite norm of each individual block - // 2. sum all norms of the blocks per row, except the diagonal block + // 2. sum all norms of the blocks per row, // 3. take the maximum of all the sums matrix_norm_ = 0.0; - max_pivot_original_ = 0.0; auto const& row_indptr = *row_indptr_; - // auto const& col_indices = *col_indices_; for (Idx row = 0; row != size_; ++row) { // calculate the sum of the norms of the blocks in the row double row_norm = 0.0; for (Idx idx = row_indptr[row]; idx != row_indptr[row + 1]; ++idx) { if constexpr (is_block) { row_norm += cabs(data[idx]).rowwise().sum().maxCoeff(); - max_pivot_original_ = std::max(max_pivot_original_, cabs(data[idx]).maxCoeff()); } else { row_norm += cabs(data[idx]); - max_pivot_original_ = std::max(max_pivot_original_, cabs(data[idx])); } } matrix_norm_ = std::max(matrix_norm_, row_norm); } } + void calculate_max_pivot(std::vector const& data) { + max_pivot_original_ = 0.0; + auto const& row_indptr = *row_indptr_; + for (Idx row = 0; row != size_; ++row) { + for (Idx idx = row_indptr[row]; idx != row_indptr[row + 1]; ++idx) { + if constexpr (is_block) { + max_pivot_original_ = std::max(max_pivot_original_, cabs(data[idx]).maxCoeff()); + } else { + max_pivot_original_ = std::max(max_pivot_original_, cabs(data[idx])); + } + } + } + } + void reset_matrix_cache() { has_pivot_perturbation_ = false; matrix_norm_ = 0.0; From 5d379a908e71eb30d08ecb0f0419d96a93abf23e Mon Sep 17 00:00:00 2001 From: Santiago Figueroa Manrique Date: Fri, 19 Sep 2025 15:42:42 +0200 Subject: [PATCH 6/6] abs and rel threshold Signed-off-by: Santiago Figueroa Manrique --- .../math_solver/sparse_lu_solver.hpp | 27 +++++-------------- 1 file changed, 7 insertions(+), 20 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp index e6430fa565..b39a809723 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp @@ -62,7 +62,7 @@ template class DenseLUFactor { static void factorize_block_in_place(Eigen::MatrixBase&& matrix, BlockPerm& block_perm, double perturb_threshold, std::pair> possibly_ill_conditioned_pivots, - Idx block_node_idx, bool& has_pivot_perturbation, double max_pivot_original) + Idx block_node_idx, bool& has_pivot_perturbation, double matrix_norm) requires(std::same_as && rk2_tensor && (Derived::RowsAtCompileTime == size) && (Derived::ColsAtCompileTime == size)) { @@ -136,9 +136,12 @@ template class DenseLUFactor { // throw SparseMatrixError if the matrix is ill-conditioned // only check condition number if pivot perturbation is not used - double const pivot_threshold = has_pivot_perturbation ? 0.0 : epsilon * max_pivot_original; + double const absolute_threshold = epsilon; + double const relative_threshold = has_pivot_perturbation ? 0.0 : epsilon * matrix_norm; for (int8_t pivot = 0; pivot != size; ++pivot) { - if (cabs(matrix(pivot, pivot)) < pivot_threshold || !is_normal(matrix(pivot, pivot))) { + auto const pivot_abs_v = cabs(matrix(pivot, pivot)); + if (pivot_abs_v < relative_threshold || pivot_abs_v < absolute_threshold || + !is_normal(matrix(pivot, pivot))) { throw SparseMatrixError{}; // can not specify error code } } @@ -245,7 +248,6 @@ template class SparseLUSolver { void prefactorize(std::vector& data, BlockPermArray& block_perm_array, std::pair> possibly_ill_conditioned_pivots = {false, {}}) { reset_matrix_cache(); - calculate_max_pivot(data); if (possibly_ill_conditioned_pivots.first) { assert(static_cast(possibly_ill_conditioned_pivots.second.size()) == size_); initialize_pivot_perturbation(data); @@ -276,7 +278,7 @@ template class SparseLUSolver { LUFactor::factorize_block_in_place(lu_matrix[pivot_idx].matrix(), block_perm_array[pivot_row_col], perturb_threshold, possibly_ill_conditioned_pivots, pivot_row_col, // Pass the node index - has_pivot_perturbation_, max_pivot_original_); + has_pivot_perturbation_, matrix_norm_); return block_perm_array[pivot_row_col]; } else { if (possibly_ill_conditioned_pivots.first) { @@ -421,7 +423,6 @@ template class SparseLUSolver { // cache value for pivot perturbation for the factorize step bool has_pivot_perturbation_{false}; double matrix_norm_{}; - double max_pivot_original_{}; std::optional> original_matrix_; // cache value for iterative refinement for the solve step std::optional> dx_; @@ -566,20 +567,6 @@ template class SparseLUSolver { } } - void calculate_max_pivot(std::vector const& data) { - max_pivot_original_ = 0.0; - auto const& row_indptr = *row_indptr_; - for (Idx row = 0; row != size_; ++row) { - for (Idx idx = row_indptr[row]; idx != row_indptr[row + 1]; ++idx) { - if constexpr (is_block) { - max_pivot_original_ = std::max(max_pivot_original_, cabs(data[idx]).maxCoeff()); - } else { - max_pivot_original_ = std::max(max_pivot_original_, cabs(data[idx])); - } - } - } - } - void reset_matrix_cache() { has_pivot_perturbation_ = false; matrix_norm_ = 0.0;