diff --git a/src/mir/method/MethodWeighted.cc b/src/mir/method/MethodWeighted.cc
index 51be149a9..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<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 = decltype(omp_orig)())
 
 namespace mir::method {
 
@@ -429,14 +433,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.outer()[r] == W.outer()[r + 1]) {
+            forceMissing.push_back(r);
         }
     }
 
diff --git a/src/mir/method/nonlinear/MissingIfHeaviestMissing.cc b/src/mir/method/nonlinear/MissingIfHeaviestMissing.cc
index f57eb1fe0..6bb54f8ef 100644
--- a/src/mir/method/nonlinear/MissingIfHeaviestMissing.cc
+++ b/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.outer();
+    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);
+        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 +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;