From f23b471cdf27595c6ca3aa119c1c3a97dc2db081 Mon Sep 17 00:00:00 2001 From: Carlos Date: Tue, 12 Nov 2024 08:52:26 +0100 Subject: [PATCH 1/6] Avoid iterator creation in forceMissing creation --- src/mir/method/MethodWeighted.cc | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/src/mir/method/MethodWeighted.cc b/src/mir/method/MethodWeighted.cc index 51be149a9..69ef9456b 100644 --- a/src/mir/method/MethodWeighted.cc +++ b/src/mir/method/MethodWeighted.cc @@ -429,14 +429,9 @@ void MethodWeighted::execute(context::Context& ctx, const repres::Representation ASSERT(W.cols() == npts_inp); std::vector forceMissing; // reserving size unnecessary (not the general case) - { - auto begin = W.begin(0); - auto end(begin); - for (size_t r = 0; r < W.rows(); r++) { - if (begin == (end = W.end(r))) { - forceMissing.push_back(r); - } - begin = end; + for (size_t r = 0; r < W.rows(); ++r) { + if (W.outer()[r] == W.outer()[r + 1]) { + forceMissing.push_back(r); } } From 7503bdd5fd785b7c13c8262325fa35da7bed9868 Mon Sep 17 00:00:00 2001 From: Carlos Date: Wed, 20 Nov 2024 13:08:34 +0100 Subject: [PATCH 2/6] Added OMP to forceMissing creation --- src/mir/method/MethodWeighted.cc | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/mir/method/MethodWeighted.cc b/src/mir/method/MethodWeighted.cc index 69ef9456b..7f1e122d6 100644 --- a/src/mir/method/MethodWeighted.cc +++ b/src/mir/method/MethodWeighted.cc @@ -44,6 +44,10 @@ #include "mir/util/Trace.h" #include "mir/util/Types.h" +#pragma omp declare reduction(vec_merge_sorted : std::vector : \ + omp_out.insert(omp_out.end(), omp_in.begin(), omp_in.end()); \ + std::inplace_merge(omp_out.begin(), omp_out.begin() + omp_out.size() - omp_in.size(), omp_out.end())) \ + initializer(omp_priv = decltype(omp_orig)()) namespace mir::method { @@ -429,6 +433,7 @@ void MethodWeighted::execute(context::Context& ctx, const repres::Representation ASSERT(W.cols() == npts_inp); std::vector forceMissing; // reserving size unnecessary (not the general case) + #pragma omp parallel for reduction(vec_merge_sorted:forceMissing) for (size_t r = 0; r < W.rows(); ++r) { if (W.outer()[r] == W.outer()[r + 1]) { forceMissing.push_back(r); From 0267e161613f1b1581d0b29c8b174f361061a66b Mon Sep 17 00:00:00 2001 From: Carlos Date: Fri, 15 Nov 2024 14:40:20 +0100 Subject: [PATCH 3/6] Modified MissingifHeaviestMissing to improve access to sparse matrix --- .../nonlinear/MissingIfHeaviestMissing.cc | 48 +++++++++---------- 1 file changed, 22 insertions(+), 26 deletions(-) diff --git a/src/mir/method/nonlinear/MissingIfHeaviestMissing.cc b/src/mir/method/nonlinear/MissingIfHeaviestMissing.cc index f57eb1fe0..5bb0aa285 100644 --- a/src/mir/method/nonlinear/MissingIfHeaviestMissing.cc +++ b/src/mir/method/nonlinear/MissingIfHeaviestMissing.cc @@ -28,40 +28,39 @@ namespace mir::method::nonlinear { MissingIfHeaviestMissing::MissingIfHeaviestMissing(const param::MIRParametrisation& param) : NonLinear(param) {} -bool MissingIfHeaviestMissing::treatment(DenseMatrix& /*A*/, WeightMatrix& W, DenseMatrix& /*B*/, - const MIRValuesVector& values, const double& missingValue) const { +bool MissingIfHeaviestMissing::treatment(MethodWeighted::Matrix& /*A*/, MethodWeighted::WeightMatrix& W, + MethodWeighted::Matrix& /*B*/, const MIRValuesVector& values, + const double& missingValue) const { // correct matrix weigths for the missing values - ASSERT(W.cols() == values.size()); + ASSERT(W.cols() == values.size()); + auto* outer = W.outer(); + auto* inner = W.inner(); auto* data = const_cast(W.data()); bool modif = false; - WeightMatrix::Size i = 0; - WeightMatrix::iterator it(W); + for (WeightMatrix::Size r = 0; r < W.rows(); ++r) { - const WeightMatrix::iterator end = W.end(r); + WeightMatrix::Size row_start = outer[r]; + WeightMatrix::Size row_end = outer[r + 1]; // Marks the end of this row - // count missing values, accumulate weights (disregarding missing values) and find maximum weight in row - size_t i_missing = i; + // Initialize variables for tracking missing values and weights in the row + size_t i_missing = row_start; size_t N_missing = 0; - size_t N_entries = 0; + size_t N_entries = row_end - row_start; double sum = 0.; double heaviest = -1.; bool heaviest_is_missing = false; - WeightMatrix::iterator kt(it); - WeightMatrix::Size k = i; - for (; it != end; ++it, ++i, ++N_entries) { - - const bool miss = values[it.col()] == missingValue; + for (WeightMatrix::Size i = row_start; i < row_end; ++i) { + const bool miss = values[inner[i]] == missingValue; if (miss) { ++N_missing; i_missing = i; - } - else { - sum += *it; + } else { + sum += data[i]; } if (heaviest < data[i]) { @@ -70,21 +69,18 @@ bool MissingIfHeaviestMissing::treatment(DenseMatrix& /*A*/, WeightMatrix& W, De } } - // weights redistribution: zero-weight all missing values, linear re-weighting for the others; - // if all values are missing, or the closest value is missing, force missing value if (N_missing > 0) { if (N_missing == N_entries || heaviest_is_missing || eckit::types::is_approximately_equal(sum, 0.)) { - for (WeightMatrix::Size j = k; j < k + N_entries; ++j) { - data[j] = j == i_missing ? 1. : 0.; + for (WeightMatrix::Size i = row_start; i < row_end; ++i) { + data[i] = (i == i_missing) ? 1. : 0.; } - } - else { + } else { const double factor = 1. / sum; - for (WeightMatrix::Size j = k; j < k + N_entries; ++j, ++kt) { - const bool miss = values[kt.col()] == missingValue; - data[j] = miss ? 0. : (factor * data[j]); + for (WeightMatrix::Size i = row_start; i < row_end; ++i) { + const bool miss = values[inner[i]] == missingValue; + data[i] = miss ? 0. : (factor * data[i]); } } modif = true; From d7931c89332e68dc15f2032d711ceb20cd3a2dfe Mon Sep 17 00:00:00 2001 From: Carlos Date: Wed, 20 Nov 2024 13:04:07 +0100 Subject: [PATCH 4/6] Added OpenMP to compute rows of MissingIfHeaviestMissing --- src/mir/method/nonlinear/MissingIfHeaviestMissing.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mir/method/nonlinear/MissingIfHeaviestMissing.cc b/src/mir/method/nonlinear/MissingIfHeaviestMissing.cc index 5bb0aa285..d6dca7221 100644 --- a/src/mir/method/nonlinear/MissingIfHeaviestMissing.cc +++ b/src/mir/method/nonlinear/MissingIfHeaviestMissing.cc @@ -40,7 +40,7 @@ bool MissingIfHeaviestMissing::treatment(MethodWeighted::Matrix& /*A*/, MethodWe auto* data = const_cast(W.data()); bool modif = false; - + #pragma omp parallel for reduction(||:modif) shared(outer, inner, data, values, missingValue, W) for (WeightMatrix::Size r = 0; r < W.rows(); ++r) { WeightMatrix::Size row_start = outer[r]; WeightMatrix::Size row_end = outer[r + 1]; // Marks the end of this row From 61866b6e8e596af4370a973ae9413e7e0cb6e0a9 Mon Sep 17 00:00:00 2001 From: Carlos Date: Wed, 20 Nov 2024 13:31:12 +0100 Subject: [PATCH 5/6] Rebased file to follow develop --- src/mir/method/nonlinear/MissingIfHeaviestMissing.cc | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/mir/method/nonlinear/MissingIfHeaviestMissing.cc b/src/mir/method/nonlinear/MissingIfHeaviestMissing.cc index d6dca7221..ffd44d4c3 100644 --- a/src/mir/method/nonlinear/MissingIfHeaviestMissing.cc +++ b/src/mir/method/nonlinear/MissingIfHeaviestMissing.cc @@ -28,9 +28,8 @@ namespace mir::method::nonlinear { MissingIfHeaviestMissing::MissingIfHeaviestMissing(const param::MIRParametrisation& param) : NonLinear(param) {} -bool MissingIfHeaviestMissing::treatment(MethodWeighted::Matrix& /*A*/, MethodWeighted::WeightMatrix& W, - MethodWeighted::Matrix& /*B*/, const MIRValuesVector& values, - const double& missingValue) const { +bool MissingIfHeaviestMissing::treatment(DenseMatrix& /*A*/, WeightMatrix& W, DenseMatrix& /*B*/, + const MIRValuesVector& values, const double& missingValue) const { // correct matrix weigths for the missing values From d08fa5f1fbf3168d0770273c9f7e45f0ef1e9d60 Mon Sep 17 00:00:00 2001 From: Carlos Date: Thu, 21 Nov 2024 08:39:12 +0100 Subject: [PATCH 6/6] Simplified and fixed OMP pragma --- src/mir/method/nonlinear/MissingIfHeaviestMissing.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mir/method/nonlinear/MissingIfHeaviestMissing.cc b/src/mir/method/nonlinear/MissingIfHeaviestMissing.cc index ffd44d4c3..6bb54f8ef 100644 --- a/src/mir/method/nonlinear/MissingIfHeaviestMissing.cc +++ b/src/mir/method/nonlinear/MissingIfHeaviestMissing.cc @@ -39,7 +39,7 @@ bool MissingIfHeaviestMissing::treatment(DenseMatrix& /*A*/, WeightMatrix& W, De auto* data = const_cast(W.data()); bool modif = false; - #pragma omp parallel for reduction(||:modif) shared(outer, inner, data, values, missingValue, W) + #pragma omp parallel for reduction(||:modif) for (WeightMatrix::Size r = 0; r < W.rows(); ++r) { WeightMatrix::Size row_start = outer[r]; WeightMatrix::Size row_end = outer[r + 1]; // Marks the end of this row