Optimizations method weighted

src/mir/method/MethodWeighted.cc

@@ -45,7 +45,7 @@
 #include "mir/util/Mutex.h"
 #include "mir/util/Trace.h"
 #include "mir/util/Types.h"
-
+#include "mir/util/OmpReductions.h"
 
 namespace mir::method {
 
@@ -446,14 +446,10 @@ void MethodWeighted::execute(context::Context& ctx, const repres::Representation
     ASSERT(W.cols() == npts_inp);
 
     std::vector<size_t> 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;
+    #pragma omp parallel for reduction(vec_merge_sorted:forceMissing)
+    for (size_t r = 0; r < W.rows(); ++r) {
+        if (W.outerIndex()[r] == W.outerIndex()[r + 1]) {
+            forceMissing.push_back(r);
         }
     }
 

src/mir/method/nonlinear/MissingIfHeaviestMissing.cc

@@ -32,36 +32,34 @@ bool MissingIfHeaviestMissing::treatment(DenseMatrix& /*A*/, WeightMatrix& W, De
                                          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.outerIndex();
+    auto* inner = W.inner();
     auto* data = const_cast<WeightMatrix::Scalar*>(W.data());
     bool modif = false;
 
-    WeightMatrix::Size i = 0;
-    WeightMatrix::iterator it(W);
+    #pragma omp parallel for reduction(||:modif) 
     for (WeightMatrix::Size r = 0; r < W.rows(); ++r) {
-        const WeightMatrix::iterator end = W.end(r);
+        auto row_start = outer[r];
+        auto 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
+        auto i_missing         = row_start;
         size_t N_missing         = 0;
-        size_t N_entries         = 0;
+        auto 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 +68,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;

src/mir/util/OmpReductions.h

@@ -0,0 +1,23 @@
+/*
+ * (C) Copyright 1996- ECMWF.
+ *
+ * This software is licensed under the terms of the Apache Licence Version 2.0
+ * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
+ *
+ * In applying this licence, ECMWF does not waive the privileges and immunities
+ * granted to it by virtue of its status as an intergovernmental organisation nor
+ * does it submit to any jurisdiction.
+ */
+
+
+#pragma once
+
+#include <vector>
+#include <algorithm>
+#include <omp.h> 
+
+#pragma omp declare reduction(vec_merge_sorted : std::vector<size_t> : \
+    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 = std::vector<size_t>()) 
+