Repository: SparseLinearAlgebra/spla (language: python)

Issue:
Inverse mask apply
Support `!mask` apply to vector and matrix inoperations.

## Exprected behavior

For primitive v of size N:
- `v[!mask] = m <=> v[i] = m[i] for i in range(0, N) i not in mask`
 
## Tasks

- [ ] Add descriptor `MaskComplement` value
- [ ] Add complement mask algo
- [ ] Wrap mask apply for `bool complement` variant
- [ ] Use new mask apply in algorithms 


Discussion / hints:

Fix patch:
diff --git a/include/spla-cpp/SplaDescriptor.hpp b/include/spla-cpp/SplaDescriptor.hpp
index 309b3f2d8..8ac4a5866 100644
--- a/include/spla-cpp/SplaDescriptor.hpp
+++ b/include/spla-cpp/SplaDescriptor.hpp
@@ -57,6 +57,8 @@ namespace spla {
             ValuesSorted,
             /** Provided matrix/vector values has no duplicated */
             NoDuplicates,
+            /** Apply !mask (complementary mask) to result or input arguments */
+            MaskComplement,
             /** Profiles time of the operation and outputs result to the log */
             ProfileTime,
             /** Transpose operation arg 1 matrix before operation */
diff --git a/sources/algo/matrix/SplaMatrixEWiseAddCOO.cpp b/sources/algo/matrix/SplaMatrixEWiseAddCOO.cpp
index 0ec6cbd52..b6c5dc49c 100644
--- a/sources/algo/matrix/SplaMatrixEWiseAddCOO.cpp
+++ b/sources/algo/matrix/SplaMatrixEWiseAddCOO.cpp
@@ -51,6 +51,7 @@ void spla::MatrixEWiseAddCOO::Process(spla::AlgorithmParams &params) {
     auto p = dynamic_cast<ParamsMatrixEWiseAdd *>(&params);
     auto w = p->w;
     auto library = p->desc->GetLibrary().GetPrivatePtr();
+    auto &desc = p->desc;
 
     auto device = library->GetDeviceManager().GetDevice(p->deviceId);
     compute::context ctx = library->GetContext();
@@ -90,22 +91,27 @@ void spla::MatrixEWiseAddCOO::Process(spla::AlgorithmParams &params) {
     compute::vector<unsigned int> tmpColsB(ctx);
 
     auto maskBlock = p->mask.Cast<MatrixCOO>();
+    auto complementMask = desc->IsParamSet(Descriptor::Param::MaskComplement);
+
     auto applyMask = [&](RefPtr<MatrixCOO> &block,
                          compute::vector<unsigned int> &tmpRows,
                          compute::vector<unsigned int> &tmpCols,
                          compute::vector<unsigned int> &perm,
                          const compute::vector<unsigned int> *&outRows,
                          const compute::vector<unsigned int> *&outCols) {
+        // Nothing to do
         if (block.IsNull())
             return;
 
-        if (!p->hasMask) {
+        // No mask or must apply inverse mask (but !null = all)
+        if (!p->hasMask || (complementMask && maskBlock.IsNull())) {
             outRows = &block->GetRows();
             outCols = &block->GetCols();
             return;
         }
 
-        if (maskBlock.IsNull())
+        // No inverse and no block - null result
+        if (!complementMask && maskBlock.IsNull())
             return;
 
         if (typeHasValues) {
@@ -113,12 +119,14 @@ void spla::MatrixEWiseAddCOO::Process(spla::AlgorithmParams &params) {
             MaskByPairKeys(maskBlock->GetRows(), maskBlock->GetCols(),
                            block->GetRows(), block->GetCols(), perm,
                            tmpRows, tmpCols, tmpPerm,
+                           complementMask,
                            queue);
             std::swap(perm, tmpPerm);
         } else
             MaskPairKeys(maskBlock->GetRows(), maskBlock->GetCols(),
                          block->GetRows(), block->GetCols(),
                          tmpRows, tmpCols,
+                         complementMask,
                          queue);
 
         outRows = &tmpRows;
diff --git a/sources/algo/vector/SplaVectorEWiseAddCOO.cpp b/sources/algo/vector/SplaVectorEWiseAddCOO.cpp
index eed10f389..2b764b6ed 100644
--- a/sources/algo/vector/SplaVectorEWiseAddCOO.cpp
+++ b/sources/algo/vector/SplaVectorEWiseAddCOO.cpp
@@ -51,6 +51,7 @@ void spla::VectorEWiseAddCOO::Process(spla::AlgorithmParams &params) {
     auto p = dynamic_cast<ParamsVectorEWiseAdd *>(&params);
     auto w = p->w;
     auto library = p->desc->GetLibrary().GetPrivatePtr();
+    auto &desc = p->desc;
     auto &logger = library->GetLogger();
 
     auto device = library->GetDeviceManager().GetDevice(p->deviceId);
@@ -87,16 +88,21 @@ void spla::VectorEWiseAddCOO::Process(spla::AlgorithmParams &params) {
     compute::vector<unsigned int> tmpRowsB(ctx);
 
     auto maskBlock = p->mask.Cast<VectorCOO>();
+    auto complementMask = desc->IsParamSet(Descriptor::Param::MaskComplement);
+
     auto applyMask = [&](RefPtr<VectorCOO> &block, compute::vector<unsigned int> &tmpRows, compute::vector<unsigned int> &perm, const compute::vector<unsigned int> *&out) {
+        // Nothing to do
         if (block.IsNull())
             return;
 
-        if (!p->hasMask) {
+        // No mask or must apply inverse mask (but !null = all)
+        if (!p->hasMask || (complementMask && maskBlock.IsNull())) {
             out = &block->GetRows();
             return;
         }
 
-        if (maskBlock.IsNull())
+        // No inverse and no block - null result
+        if (!complementMask && maskBlock.IsNull())
             return;
 
         if (typeHasValues) {
@@ -104,12 +110,14 @@ void spla::VectorEWiseAddCOO::Process(spla::AlgorithmParams &params) {
             MaskByKeys(maskBlock->GetRows(),
                        block->GetRows(), perm,
                        tmpRows, tmpPerm,
+                       complementMask,
                        queue);
             std::swap(perm, tmpPerm);
         } else
             MaskKeys(maskBlock->GetRows(),
                      block->GetRows(),
                      tmpRows,
+                     complementMask,
                      queue);
 
         out = &tmpRows;
diff --git a/sources/algo/vxm/SplaVxMCOO.cpp b/sources/algo/vxm/SplaVxMCOO.cpp
index 2a8a1cef2..59e93b56c 100644
--- a/sources/algo/vxm/SplaVxMCOO.cpp
+++ b/sources/algo/vxm/SplaVxMCOO.cpp
@@ -55,6 +55,7 @@ void spla::VxMCOO::Process(spla::AlgorithmParams &params) {
     auto p = dynamic_cast<ParamsVxM *>(&params);
     auto w = p->w;
     auto library = p->desc->GetLibrary().GetPrivatePtr();
+    auto &desc = p->desc;
 
     auto device = library->GetDeviceManager().GetDevice(p->deviceId);
     compute::context ctx = library->GetContext();
@@ -64,8 +65,9 @@ void spla::VxMCOO::Process(spla::AlgorithmParams &params) {
     auto a = p->a.Cast<VectorCOO>();
     auto b = p->b.Cast<MatrixCOO>();
     auto mask = p->mask.Cast<VectorCOO>();
+    auto complementMask = desc->IsParamSet(Descriptor::Param::MaskComplement);
 
-    if (p->hasMask && mask.IsNull())
+    if (p->hasMask && (complementMask && mask.IsNull()))
         return;
 
     if (a->GetNvals() == 0 || b->GetNvals() == 0)
@@ -86,7 +88,6 @@ void spla::VxMCOO::Process(spla::AlgorithmParams &params) {
     // Compute number of products for each a[i] x b[i,:]
     compute::vector<unsigned int> segmentLengths(a->GetNvals() + 1, ctx);
     compute::gather(a->GetRows().begin(), a->GetRows().end(), lengths.begin(), segmentLengths.begin(), queue);
-    (segmentLengths.end() - 1).write(0, queue);
 
     // Compute offsets between each a[i] x b[i,:] products
     compute::vector<unsigned int> outputPtr(a->GetNvals() + 1, ctx);
@@ -147,11 +148,11 @@ void spla::VxMCOO::Process(spla::AlgorithmParams &params) {
         // Reduce all produces a[i] * b[i, j] for j using provided add op
         ReduceByKey(J, V, rows, vals, tw->GetByteSize(), p->add->GetSource(), queue);
 
-        // Apply mask if required (at this point we know, that mask is not null)
-        if (p->hasMask) {
+        // Apply mask if required
+        if (p->hasMask && mask.IsNotNull()) {
             compute::vector<unsigned int> tmpRows(ctx);
             compute::vector<unsigned char> tmpVals(ctx);
-            ApplyMask(mask->GetRows(), rows, vals, tmpRows, tmpVals, tw->GetByteSize(), queue);
+            ApplyMask(mask->GetRows(), rows, vals, tmpRows, tmpVals, tw->GetByteSize(), complementMask, queue);
             std::swap(rows, tmpRows);
             std::swap(vals, tmpVals);
         }
@@ -163,16 +164,19 @@ void spla::VxMCOO::Process(spla::AlgorithmParams &params) {
         ReduceDuplicates(J, rows, queue);
 
         // Apply mask to indices
-        if (p->hasMask) {
+        if (p->hasMask && mask.IsNotNull()) {
             compute::vector<unsigned int> tmpRows(ctx);
-            MaskKeys(mask->GetRows(), rows, tmpRows, queue);
+            MaskKeys(mask->GetRows(), rows, tmpRows, complementMask, queue);
             std::swap(rows, tmpRows);
         }
     }
 
     // Store result
-    auto nvals = rows.size();
-    p->w = VectorCOO::Make(N, nvals, std::move(rows), std::move(vals)).As<VectorBlock>();
+    // NOTE: can be empty due mask application as the last step
+    if (!rows.empty()) {
+        auto nvals = rows.size();
+        p->w = VectorCOO::Make(N, nvals, std::move(rows), std::move(vals)).As<VectorBlock>();
+    }
 }
 
 spla::Algorithm::Type spla::VxMCOO::GetType() const {
diff --git a/sources/compute/SplaApplyMask.hpp b/sources/compute/SplaApplyMask.hpp
index ff32c55d2..b2114501f 100644
--- a/sources/compute/SplaApplyMask.hpp
+++ b/sources/compute/SplaApplyMask.hpp
@@ -49,6 +49,7 @@ namespace spla {
      * @param outputRows Result row indices; resized automatically
      * @param outputVals Result values; resized automatically
      * @param byteSize Size of values; if 0, apply only indices mask
+     * @param complement Pass true to apply inverse (complementary mask)
      * @param queue Command queue to perform operation
      */
     inline void ApplyMask(const boost::compute::vector<unsigned int> &mask,
@@ -57,6 +58,7 @@ namespace spla {
                           boost::compute::vector<unsigned int> &outputRows,
                           boost::compute::vector<unsigned char> &outputVals,
                           std::size_t byteSize,
+                          bool complement,
                           boost::compute::command_queue &queue) {
         using namespace boost;
 
@@ -80,10 +82,11 @@ namespace spla {
             MaskByKeys(mask,
                        inputRows, perm,
                        outputRows, outputPerm,
+                       complement,
                        queue);
             std::swap(perm, outputPerm);
         } else
-            MaskKeys(mask, inputRows, outputRows, queue);
+            MaskKeys(mask, inputRows, outputRows, complement, queue);
 
         if (typeHasValues) {
             // Copy output values using indices buffer
diff --git a/sources/compute/SplaMaskByKey.hpp b/sources/compute/SplaMaskByKey.hpp
index a9d61d8a3..f4ce8f2b1 100644
--- a/sources/compute/SplaMaskByKey.hpp
+++ b/sources/compute/SplaMaskByKey.hpp
@@ -168,7 +168,8 @@ namespace spla {
                            InputIterator4 tile_first2,
                            OutputIterator counts,
                            Compare compare,
-                           Equals equals) {
+                           Equals equals,
+                           bool complement) {
                 using uint_ = boost::compute::uint_;
                 m_count = boost::compute::detail::iterator_range_size(tile_first1, tile_last1) - 1;
 
@@ -179,17 +180,36 @@ namespace spla {
                       << "uint end2 = " << tile_first2[expr<uint_>("i+1")] << ";\n"
                       << "uint count = 0;\n"
                       << "while(start1<end1 && start2<end2)\n"
-                      << "{\n"
-                      << "   if(" << equals(maskFirst[expr<uint_>("start1")], keyFirsts[expr<uint_>("start2")]) << ")\n"
-                      << "   {\n"
-                      << "       count++;\n"
-                      << "       start1++; start2++;\n"
-                      << "   }\n"
-                      << "   else if(" << compare(maskFirst[expr<uint_>("start1")], keyFirsts[expr<uint_>("start2")]) << ")\n"
-                      << "       start1++;\n"
-                      << "   else start2++;\n"
-                      << "}\n"
-                      << counts[expr<uint_>("i")] << " = count;\n";
+                      << "{\n";
+
+                if (complement) {
+                    *this << "   if(" << equals(maskFirst[expr<uint_>("start1")], keyFirsts[expr<uint_>("start2")]) << ")\n"
+                          << "   {\n"
+                          << "       start1++; start2++;\n"
+                          << "   }\n"
+                          << "   else if(" << compare(maskFirst[expr<uint_>("start1")], keyFirsts[expr<uint_>("start2")]) << ")\n"
+                          << "       start1++;\n"
+                          << "   else\n"
+                          << "   {\n"
+                          << "       count++; start2++;\n"
+                          << "   }\n";
+                } else {
+                    *this << "   if(" << equals(maskFirst[expr<uint_>("start1")], keyFirsts[expr<uint_>("start2")]) << ")\n"
+                          << "   {\n"
+                          << "       count++;\n"
+                          << "       start1++; start2++;\n"
+                          << "   }\n"
+                          << "   else if(" << compare(maskFirst[expr<uint_>("start1")], keyFirsts[expr<uint_>("start2")]) << ")\n"
+                          << "       start1++;\n"
+                          << "   else start2++;\n";
+                }
+
+                *this << "}\n";
+
+                if (complement)
+                    *this << counts[expr<uint_>("i")] << " = count + (end2 - start2);\n";
+                else
+                    *this << counts[expr<uint_>("i")] << " = count;\n";
             }
 
             boost::compute::event exec(boost::compute::command_queue &queue) {
@@ -225,7 +245,8 @@ namespace spla {
                            InputIterator5 counts,
                            Compare compare,
                            Equals equals,
-                           AssignResult assignResult) {
+                           AssignResult assignResult,
+                           bool complement) {
                 using uint_ = boost::compute::uint_;
                 m_count = boost::compute::detail::iterator_range_size(tile_first1, tile_last1) - 1;
 
@@ -236,17 +257,41 @@ namespace spla {
                       << "uint end2 = " << tile_first2[expr<uint_>("i+1")] << ";\n"
                       << "uint count = " << counts[expr<uint_>("i")] << ";\n"
                       << "while(start1<end1 && start2<end2)\n"
-                      << "{\n"
-                      << "   if(" << equals(maskFirst[expr<uint_>("start1")], keyFirsts[expr<uint_>("start2")]) << ")\n"
-                      << "   {\n";
-                assignResult(*this, expr<uint_>("count"), expr<uint_>("start2"));
-                *this << "       count++;\n"
-                      << "       start1++; start2++;\n"
-                      << "   }\n"
-                      << "   else if(" << compare(maskFirst[expr<uint_>("start1")], keyFirsts[expr<uint_>("start2")]) << ")\n"
-                      << "       start1++;\n"
-                      << "   else start2++;\n"
-                      << "}\n";
+                      << "{\n";
+
+                if (complement) {
+                    *this << "   if(" << equals(maskFirst[expr<uint_>("start1")], keyFirsts[expr<uint_>("start2")]) << ")\n"
+                          << "   {\n"
+                          << "       start1++; start2++;\n"
+                          << "   }\n"
+                          << "   else if(" << compare(maskFirst[expr<uint_>("start1")], keyFirsts[expr<uint_>("start2")]) << ")\n"
+                          << "       start1++;\n"
+                          << "   else\n"
+                          << "   {\n";
+                    assignResult(*this, expr<uint_>("count"), expr<uint_>("start2"));
+                    *this << "       count++; start2++;\n"
+                          << "   }\n";
+                } else {
+                    *this << "   if(" << equals(maskFirst[expr<uint_>("start1")], keyFirsts[expr<uint_>("start2")]) << ")\n"
+                          << "   {\n";
+                    assignResult(*this, expr<uint_>("count"), expr<uint_>("start2"));
+                    *this << "       count++;\n"
+                          << "       start1++; start2++;\n"
+                          << "   }\n"
+                          << "   else if(" << compare(maskFirst[expr<uint_>("start1")], keyFirsts[expr<uint_>("start2")]) << ")\n"
+                          << "       start1++;\n"
+                          << "   else start2++;\n";
+                }
+
+                *this << "}\n";
+
+                if (complement) {
+                    *this << "while (start2 < end2)\n"
+                          << "{\n";
+                    assignResult(*this, expr<uint_>("count"), expr<uint_>("start2"));
+                    *this << "    count++; start2++;"
+                          << "}\n";
+                }
             }
 
             boost::compute::event exec(boost::compute::command_queue &queue) {
@@ -276,6 +321,7 @@ namespace spla {
                               AssignResult assignResult,
                               std::size_t maskCount,
                               std::size_t keyCount,
+                              bool complement,
                               boost::compute::command_queue &queue) {
             using namespace boost;
             typedef typename std::iterator_traits<InputMask>::value_type key_type;
@@ -308,7 +354,7 @@ namespace spla {
                                               tileA.begin(), tileA.end(),
                                               tileB.begin(),
                                               counts.begin(),
-                                              compare, equals);
+                                              compare, equals, complement);
             intersectionCountKernel.exec(queue);
 
             // Compute actual counts offsets
@@ -323,7 +369,7 @@ namespace spla {
             intersectionKernel.set_range(maskFirst, keyFirst, tileA.begin(), tileA.end(),
                                          tileB.begin(),
                                          counts.begin(),
-                                         compare, equals, assignResult);
+                                         compare, equals, assignResult, complement);
             intersectionKernel.exec(queue);
 
             return resultCount;
@@ -338,10 +384,12 @@ namespace spla {
      * keys range and stores it in range starting at resultKeys.
      *
      * @note Automatically resizes result containers to result count size.
+     * @note Use complement flag to apply direct or inverse (complement) mask
      *
      * @param mask Mask elements
      * @param keys Keys elements
      * @param resultKeys Result keys elements
+     * @param complement Pass true to apply !mask (complementary mask)
      * @param queue Command queue to perform operations on
      *
      * @return Count of values in intersected region
@@ -349,6 +397,7 @@ namespace spla {
     inline std::size_t MaskKeys(const boost::compute::vector<unsigned int> &mask,
                                 const boost::compute::vector<unsigned int> &keys,
                                 boost::compute::vector<unsigned int> &resultKeys,
+                                bool complement,
                                 boost::compute::command_queue &queue) {
         using namespace boost;
         using MetaKernel = boost::compute::detail::meta_kernel;
@@ -378,6 +427,7 @@ namespace spla {
                                  assignResult,
                                  mask.size(),
                                  keys.size(),
+                                 complement,
                                  queue);
     }
 
@@ -389,12 +439,14 @@ namespace spla {
      *
      * @note Manages associated values with keys.
      * @note Automatically resizes result containers to result count size.
+     * @note Use complement flag to apply direct or inverse (complement) mask
      *
      * @param mask Mask elements
      * @param keys Keys elements
      * @param values Associated with keys values
      * @param resultKeys Result keys elements
      * @param resultValues Result values associated with result keys.
+     * @param complement Pass true to apply !mask (complementary mask)
      * @param queue Command queue to perform operations on
      *
      * @return Count of values in intersected region
@@ -404,6 +456,7 @@ namespace spla {
                                   const boost::compute::vector<unsigned int> &values,
                                   boost::compute::vector<unsigned int> &resultKeys,
                                   boost::compute::vector<unsigned int> &resultValues,
+                                  bool complement,
                                   boost::compute::command_queue &queue) {
         using namespace boost;
         using MetaKernel = boost::compute::detail::meta_kernel;
@@ -437,6 +490,7 @@ namespace spla {
                                  assignResult,
                                  mask.size(),
                                  keys.size(),
+                                 complement,
                                  queue);
     }
 
@@ -448,6 +502,7 @@ namespace spla {
      *
      * @note Interprets keys as pairs, where first and second elements stored in separate arrays.
      * @note Automatically resizes result containers to result count size.
+     * @note Use complement flag to apply direct or inverse (complement) mask
      *
      * @param mask1 Mask first elements
      * @param mask2 Mask second elements
@@ -455,6 +510,7 @@ namespace spla {
      * @param keys2 Keys second elements
      * @param resultKeys1 Result keys first elements
      * @param resultKeys2 Result keys second elements
+     * @param complement Pass true to apply !mask (complementary mask)
      * @param queue Command queue to perform operations on
      *
      * @return Count of values in intersected region
@@ -465,6 +521,7 @@ namespace spla {
                                     const boost::compute::vector<unsigned int> &keys2,
                                     boost::compute::vector<unsigned int> &resultKeys1,
                                     boost::compute::vector<unsigned int> &resultKeys2,
+                                    bool complement,
                                     boost::compute::command_queue &queue) {
         using namespace boost;
         using MetaKernel = boost::compute::detail::meta_kernel;
@@ -500,6 +557,7 @@ namespace spla {
                                  assignResult,
                                  mask1.size(),
                                  keys1.size(),
+                                 complement,
                                  queue);
     }
 
@@ -512,6 +570,7 @@ namespace spla {
      * @note Manages associated values with keys.
      * @note Interprets keys as pairs, where first and second elements stored in separate arrays.
      * @note Automatically resizes result containers to result count size.
+     * @note Use complement flag to apply direct or inverse (complement) mask
      *
      * @param mask1 Mask first elements
      * @param mask2 Mask second elements
@@ -521,6 +580,7 @@ namespace spla {
      * @param resultKeys1 Result keys first elements
      * @param resultKeys2 Result keys second elements
      * @param resultValues Result values associated with result keys.
+     * @param complement Pass true to apply !mask (complementary mask)
      * @param queue Command queue to perform operations on
      *
      * @return Count of values in intersected region
@@ -533,6 +593,7 @@ namespace spla {
                                       boost::compute::vector<unsigned int> &resultKeys1,
                                       boost::compute::vector<unsigned int> &resultKeys2,
                                       boost::compute::vector<unsigned int> &resultValues,
+                                      bool complement,
                                       boost::compute::command_queue &queue) {
         using namespace boost;
         using MetaKernel = boost::compute::detail::meta_kernel;
@@ -571,6 +632,7 @@ namespace spla {
                                  assignResult,
                                  mask1.size(),
                                  keys1.size(),
+                                 complement,
                                  queue);
     }
 


Test patch:
diff --git a/tests/TestMaskByKey.cpp b/tests/TestMaskByKey.cpp
index 587725e68..4c5e867cf 100644
--- a/tests/TestMaskByKey.cpp
+++ b/tests/TestMaskByKey.cpp
@@ -61,6 +61,53 @@ TEST(Compute, MaskByKeys) {
                                   deviceValues,
                                   resultKeys,
                                   resultValues,
+                                  false,
+                                  queue);
+
+    ASSERT_EQ(count, expectedKeys.size());
+    ASSERT_EQ(count, expectedValues.size());
+
+    for (std::size_t i = 0; i < count; i++)
+        EXPECT_EQ(expectedKeys[i], (resultKeys.begin() + i).read(queue));
+
+    for (std::size_t i = 0; i < count; i++)
+        EXPECT_EQ(expectedValues[i], (resultValues.begin() + i).read(queue));
+}
+
+TEST(Compute, MaskByKeyComp) {
+    namespace compute = boost::compute;
+
+    // get the default compute device
+    compute::device gpu = compute::system::default_device();
+
+    // create a compute context and command queue
+    compute::context ctx(gpu);
+    compute::command_queue queue(ctx, gpu);
+
+    std::vector<unsigned int> mask = {0, 1, 3, 8, 10, 27};
+    std::vector<unsigned int> keys = {0, 2, 3, 4, 7, 9, 10, 15, 19, 27};
+    std::vector<unsigned int> values = {66, 24, 13, 0, 7, 119, 1, 0, 1, 200};
+
+    std::vector<unsigned int> expectedKeys = {2, 4, 7, 9, 15, 19};
+    std::vector<unsigned int> expectedValues = {24, 0, 7, 119, 0, 1};
+
+    compute::vector<unsigned int> deviceMask(mask.size(), ctx);
+    compute::vector<unsigned int> deviceKeys(keys.size(), ctx);
+    compute::vector<unsigned int> deviceValues(values.size(), ctx);
+
+    compute::vector<unsigned int> resultKeys(ctx);
+    compute::vector<unsigned int> resultValues(ctx);
+
+    compute::copy(mask.begin(), mask.end(), deviceMask.begin(), queue);
+    compute::copy(keys.begin(), keys.end(), deviceKeys.begin(), queue);
+    compute::copy(values.begin(), values.end(), deviceValues.begin(), queue);
+
+    auto count = spla::MaskByKeys(deviceMask,
+                                  deviceKeys,
+                                  deviceValues,
+                                  resultKeys,
+                                  resultValues,
+                                  true,
                                   queue);
 
     ASSERT_EQ(count, expectedKeys.size());
@@ -100,6 +147,43 @@ TEST(Compute, MaskKeys) {
     auto count = spla::MaskKeys(deviceMask,
                                 deviceKeys,
                                 resultKeys,
+                                false,
+                                queue);
+
+    ASSERT_EQ(count, expectedKeys.size());
+
+    for (std::size_t i = 0; i < count; i++)
+        EXPECT_EQ(expectedKeys[i], (resultKeys.begin() + i).read(queue));
+}
+
+TEST(Compute, MaskKeysComp) {
+    namespace compute = boost::compute;
+
+    // get the default compute device
+    compute::device gpu = compute::system::default_device();
+
+    // create a compute context and command queue
+    compute::context ctx(gpu);
+    compute::command_queue queue(ctx, gpu);
+
+    std::vector<unsigned int> mask = {0, 1, 3, 8, 10, 27};
+    std::vector<unsigned int> keys = {0, 2, 3, 4, 7, 9, 10, 15, 19, 27};
+
+    std::vector<unsigned int> expectedKeys = {2, 4, 7, 9, 15, 19};
+
+    compute::vector<unsigned int> deviceMask(mask.size(), ctx);
+    compute::vector<unsigned int> deviceKeys(keys.size(), ctx);
+
+    compute::vector<unsigned int> resultKeys(ctx);
+    compute::vector<unsigned int> resultValues(ctx);
+
+    compute::copy(mask.begin(), mask.end(), deviceMask.begin(), queue);
+    compute::copy(keys.begin(), keys.end(), deviceKeys.begin(), queue);
+
+    auto count = spla::MaskKeys(deviceMask,
+                                deviceKeys,
+                                resultKeys,
+                                true,
                                 queue);
 
     ASSERT_EQ(count, expectedKeys.size());
@@ -147,6 +231,63 @@ TEST(Compute, MaskByPairKeys) {
     auto count = spla::MaskByPairKeys(deviceMask1, deviceMask2,
                                       deviceKeys1, deviceKeys2, deviceValues,
                                       resultKeys1, resultKeys2, resultValues,
+                                      false,
+                                      queue);
+
+    ASSERT_EQ(count, expectedKeys1.size());
+    ASSERT_EQ(count, expectedKeys2.size());
+    ASSERT_EQ(count, expectedValues.size());
+
+    for (std::size_t i = 0; i < count; i++)
+        EXPECT_EQ(expectedKeys1[i], (resultKeys1.begin() + i).read(queue));
+
+    for (std::size_t i = 0; i < count; i++)
+        EXPECT_EQ(expectedKeys2[i], (resultKeys2.begin() + i).read(queue));
+
+    for (std::size_t i = 0; i < count; i++)
+        EXPECT_EQ(expectedValues[i], (resultValues.begin() + i).read(queue));
+}
+
+TEST(Compute, MaskByPairKeysComp) {
+    namespace compute = boost::compute;
+
+    // get the default compute device
+    compute::device gpu = compute::system::default_device();
+
+    // create a compute context and command queue
+    compute::context ctx(gpu);
+    compute::command_queue queue(ctx, gpu);
+
+    std::vector<unsigned int> mask1 = {0, 1, 1, 8, 10, 27};
+    std::vector<unsigned int> mask2 = {10, 0, 2, 0, 0, 100};
+    std::vector<unsigned int> keys1 = {0, 0, 1, 7, 10, 15, 19, 27};
+    std::vector<unsigned int> keys2 = {10, 20, 2, 0, 0, 10, 40, 100};
+    std::vector<unsigned int> values = {66, 24, 13, 7, 1, 0, 1, 200};
+
+    std::vector<unsigned int> expectedKeys1 = {0, 7, 15, 19};
+    std::vector<unsigned int> expectedKeys2 = {20, 0, 10, 40};
+    std::vector<unsigned int> expectedValues = {24, 7, 0, 1};
+
+    compute::vector<unsigned int> deviceMask1(mask1.size(), ctx);
+    compute::vector<unsigned int> deviceMask2(mask2.size(), ctx);
+    compute::vector<unsigned int> deviceKeys1(keys1.size(), ctx);
+    compute::vector<unsigned int> deviceKeys2(keys2.size(), ctx);
+    compute::vector<unsigned int> deviceValues(values.size(), ctx);
+
+    compute::vector<unsigned int> resultKeys1(ctx);
+    compute::vector<unsigned int> resultKeys2(ctx);
+    compute::vector<unsigned int> resultValues(ctx);
+
+    compute::copy(mask1.begin(), mask1.end(), deviceMask1.begin(), queue);
+    compute::copy(mask2.begin(), mask2.end(), deviceMask2.begin(), queue);
+    compute::copy(keys1.begin(), keys1.end(), deviceKeys1.begin(), queue);
+    compute::copy(keys2.begin(), keys2.end(), deviceKeys2.begin(), queue);
+    compute::copy(values.begin(), values.end(), deviceValues.begin(), queue);
+
+    auto count = spla::MaskByPairKeys(deviceMask1, deviceMask2,
+                                      deviceKeys1, deviceKeys2, deviceValues,
+                                      resultKeys1, resultKeys2, resultValues,
+                                      true,
                                       queue);
 
     ASSERT_EQ(count, expectedKeys1.size());
@@ -197,6 +338,54 @@ TEST(Compute, MaskPairKeys) {
     auto count = spla::MaskPairKeys(deviceMask1, deviceMask2,
                                     deviceKeys1, deviceKeys2,
                                     resultKeys1, resultKeys2,
+                                    false,
+                                    queue);
+
+    ASSERT_EQ(count, expectedKeys1.size());
+    ASSERT_EQ(count, expectedKeys2.size());
+
+    for (std::size_t i = 0; i < count; i++)
+        EXPECT_EQ(expectedKeys1[i], (resultKeys1.begin() + i).read(queue));
+
+    for (std::size_t i = 0; i < count; i++)
+        EXPECT_EQ(expectedKeys2[i], (resultKeys2.begin() + i).read(queue));
+}
+
+TEST(Compute, MaskPairKeysComp) {
+    namespace compute = boost::compute;
+
+    // get the default compute device
+    compute::device gpu = compute::system::default_device();
+
+    // create a compute context and command queue
+    compute::context ctx(gpu);
+    compute::command_queue queue(ctx, gpu);
+
+    std::vector<unsigned int> mask1 = {0, 1, 1, 8, 10, 27};
+    std::vector<unsigned int> mask2 = {10, 0, 2, 0, 0, 100};
+    std::vector<unsigned int> keys1 = {0, 0, 1, 7, 10, 15, 19, 27};
+    std::vector<unsigned int> keys2 = {10, 20, 2, 0, 0, 10, 40, 100};
+
+    std::vector<unsigned int> expectedKeys1 = {0, 7, 15, 19};
+    std::vector<unsigned int> expectedKeys2 = {20, 0, 10, 40};
+
+    compute::vector<unsigned int> deviceMask1(mask1.size(), ctx);
+    compute::vector<unsigned int> deviceMask2(mask2.size(), ctx);
+    compute::vector<unsigned int> deviceKeys1(keys1.size(), ctx);
+    compute::vector<unsigned int> deviceKeys2(keys2.size(), ctx);
+
+    compute::vector<unsigned int> resultKeys1(ctx);
+    compute::vector<unsigned int> resultKeys2(ctx);
+
+    compute::copy(mask1.begin(), mask1.end(), deviceMask1.begin(), queue);
+    compute::copy(mask2.begin(), mask2.end(), deviceMask2.begin(), queue);
+    compute::copy(keys1.begin(), keys1.end(), deviceKeys1.begin(), queue);
+    compute::copy(keys2.begin(), keys2.end(), deviceKeys2.begin(), queue);
+
+    auto count = spla::MaskPairKeys(deviceMask1, deviceMask2,
+                                    deviceKeys1, deviceKeys2,
+                                    resultKeys1, resultKeys2,
+                                    true,
                                     queue);
 
     ASSERT_EQ(count, expectedKeys1.size());
diff --git a/tests/TestMatrixEWiseAdd.cpp b/tests/TestMatrixEWiseAdd.cpp
index 3d109dfe8..502bbb5a6 100644
--- a/tests/TestMatrixEWiseAdd.cpp
+++ b/tests/TestMatrixEWiseAdd.cpp
@@ -98,7 +98,48 @@ void testMasked(spla::Library &library, std::size_t M, std::size_t N, std::size_
     spExpr->Wait();
     ASSERT_EQ(spExpr->GetState(), spla::Expression::State::Evaluated);
 
-    utils::Matrix<Type> c = a.EWiseAdd(mask, b, op);
+    utils::Matrix<Type> c = a.EWiseAdd(mask, false, b, op);
+    ASSERT_TRUE(c.Equals(spW));
+}
+
+template<typename Type, typename BinaryOp>
+void testMaskedComplement(spla::Library &library, std::size_t M, std::size_t N, std::size_t nvals,
+                          const spla::RefPtr<spla::Type> &spT,
+                          const spla::RefPtr<spla::FunctionBinary> &spOp,
+                          BinaryOp op, std::size_t seed = 0) {
+    utils::Matrix a = utils::Matrix<Type>::Generate(M, N, nvals, seed).SortReduceDuplicates();
+    utils::Matrix b = utils::Matrix<Type>::Generate(M, N, nvals, seed + 1).SortReduceDuplicates();
+    utils::Matrix mask = utils::Matrix<unsigned short>::Generate(M, N, nvals, seed + 2).SortReduceDuplicates();
+
+    a.Fill(utils::UniformGenerator<Type>());
+    b.Fill(utils::UniformGenerator<Type>());
+
+    auto spA = spla::Matrix::Make(M, N, spT, library);
+    auto spB = spla::Matrix::Make(M, N, spT, library);
+    auto spW = spla::Matrix::Make(M, N, spT, library);
+    auto spMask = spla::Matrix::Make(M, N, spla::Types::Void(library), library);
+
+    // Specify, that values already in row order + no duplicates
+    auto spDesc = spla::Descriptor::Make(library);
+    spDesc->SetParam(spla::Descriptor::Param::ValuesSorted);
+    spDesc->SetParam(spla::Descriptor::Param::NoDuplicates);
+
+    auto spOpDesc = spla::Descriptor::Make(library);
+    spOpDesc->SetParam(spla::Descriptor::Param::MaskComplement);
+
+    auto spExpr = spla::Expression::Make(library);
+    auto spWriteA = spExpr->MakeDataWrite(spA, a.GetData(library), spDesc);
+    auto spWriteB = spExpr->MakeDataWrite(spB, b.GetData(library), spDesc);
+    auto spWriteMask = spExpr->MakeDataWrite(spMask, mask.GetDataIndices(library), spDesc);
+    auto spEAddAB = spExpr->MakeEWiseAdd(spW, spMask, spOp, spA, spB, spOpDesc);
+    spExpr->Dependency(spWriteA, spEAddAB);
+    spExpr->Dependency(spWriteB, spEAddAB);
+    spExpr->Dependency(spWriteMask, spEAddAB);
+    spExpr->Submit();
+    spExpr->Wait();
+    ASSERT_EQ(spExpr->GetState(), spla::Expression::State::Evaluated);
+
+    utils::Matrix<Type> c = a.EWiseAdd(mask, true, b, op);
     ASSERT_TRUE(c.Equals(spW));
 }
 
@@ -137,7 +178,7 @@ void testOneIsEmpty(spla::Library &library, std::size_t M, std::size_t N, std::s
     spExpr->Wait();
     ASSERT_EQ(spExpr->GetState(), spla::Expression::State::Evaluated);
 
-    utils::Matrix<Type> c = a.EWiseAdd(mask, b, op);
+    utils::Matrix<Type> c = a.EWiseAdd(mask, false, b, op);
     ASSERT_TRUE(c.Equals(spW1));
     ASSERT_TRUE(c.Equals(spW2));
 }
@@ -170,7 +211,7 @@ void testNoValues(spla::Library &library, std::size_t M, std::size_t N, std::siz
     spExpr->Wait();
     ASSERT_EQ(spExpr->GetState(), spla::Expression::State::Evaluated);
 
-    utils::Matrix<unsigned char> c = a.EWiseAdd(mask, b, [](unsigned char a, unsigned char b) { return 0; });
+    utils::Matrix<unsigned char> c = a.EWiseAdd(mask, false, b, [](unsigned char a, unsigned char b) { return 0; });
     ASSERT_TRUE(c.EqualsStructure(spW));
 }
 
@@ -190,6 +231,11 @@ void test(std::size_t M, std::size_t N, std::size_t base, std::size_t step, std:
             testMasked<float>(library, M, N, nvals, spT, spOp, op, i);
         }
 
+        for (std::size_t i = 0; i < iter; i++) {
+            std::size_t nvals = base + i * step;
+            testMaskedComplement<float>(library, M, N, nvals, spT, spOp, op, i);
+        }
+
         for (std::size_t i = 0; i < iter; i++) {
             std::size_t nvals = base + i * step;
             testOneIsEmpty<float>(library, M, N, nvals, spT, spOp, op, i);
@@ -211,6 +257,11 @@ void test(std::size_t M, std::size_t N, std::size_t base, std::size_t step, std:
             testMasked<std::int32_t>(library, M, N, nvals, spT, spOp, op, i);
         }
 
+        for (std::size_t i = 0; i < iter; i++) {
+            std::size_t nvals = base + i * step;
+            testMaskedComplement<std::int32_t>(library, M, N, nvals, spT, spOp, op, i);
+        }
+
         for (std::size_t i = 0; i < iter; i++) {
             std::size_t nvals = base + i * step;
             testOneIsEmpty<std::int32_t>(library, M, N, nvals, spT, spOp, op, i);
diff --git a/tests/TestMxM.cpp b/tests/TestMxM.cpp
index 750f9eca3..5ae934ac1 100644
--- a/tests/TestMxM.cpp
+++ b/tests/TestMxM.cpp
@@ -103,7 +103,7 @@ void testMasked(spla::Library &library,
 
     ASSERT_EQ(spExpr->GetState(), spla::Expression::State::Evaluated);
 
-    utils::Matrix<Type> c = a.template MxM<Type>(mask, b, multOp, addOp);
+    utils::Matrix<Type> c = a.template MxM<Type>(mask, false, b, multOp, addOp);
     // ASSERT_TRUE(c.Equals(spW));
 }
 
diff --git a/tests/TestVectorEWiseAdd.cpp b/tests/TestVectorEWiseAdd.cpp
index e586feb35..41ca6bd1f 100644
--- a/tests/TestVectorEWiseAdd.cpp
+++ b/tests/TestVectorEWiseAdd.cpp
@@ -95,7 +95,48 @@ void testMasked(spla::Library &library, std::size_t M, std::size_t nvals,
     spExpr->Wait();
     ASSERT_EQ(spExpr->GetState(), spla::Expression::State::Evaluated);
 
-    utils::Vector<Type> c = a.EWiseAdd(mask, b, op);
+    utils::Vector<Type> c = a.EWiseAdd(mask, false, b, op);
+    ASSERT_TRUE(c.Equals(spW));
+}
+
+template<typename Type, typename BinaryOp>
+void testMaskedComplement(spla::Library &library, std::size_t M, std::size_t nvals,
+                          const spla::RefPtr<spla::Type> &spT,
+                          const spla::RefPtr<spla::FunctionBinary> &spOp,
+                          BinaryOp op, std::size_t seed = 0) {
+    utils::Vector a = utils::Vector<Type>::Generate(M, nvals, seed).SortReduceDuplicates();
+    utils::Vector b = utils::Vector<Type>::Generate(M, nvals, seed + 1).SortReduceDuplicates();
+    utils::Vector mask = utils::Vector<unsigned short>::Generate(M, nvals, seed + 2).SortReduceDuplicates();
+
+    a.Fill(utils::UniformGenerator<Type>());
+    b.Fill(utils::UniformGenerator<Type>());
+
+    auto spA = spla::Vector::Make(M, spT, library);
+    auto spB = spla::Vector::Make(M, spT, library);
+    auto spW = spla::Vector::Make(M, spT, library);
+    auto spMask = spla::Vector::Make(M, spla::Types::Void(library), library);
+
+    // Specify, that values already in row order + no duplicates
+    auto spDesc = spla::Descriptor::Make(library);
+    spDesc->SetParam(spla::Descriptor::Param::ValuesSorted);
+    spDesc->SetParam(spla::Descriptor::Param::NoDuplicates);
+
+    auto spOpDesc = spla::Descriptor::Make(library);
+    spOpDesc->SetParam(spla::Descriptor::Param::MaskComplement);
+
+    auto spExpr = spla::Expression::Make(library);
+    auto spWriteA = spExpr->MakeDataWrite(spA, a.GetData(library), spDesc);
+    auto spWriteB = spExpr->MakeDataWrite(spB, b.GetData(library), spDesc);
+    auto spWriteMask = spExpr->MakeDataWrite(spMask, mask.GetDataIndices(library), spDesc);
+    auto spEAddAB = spExpr->MakeEWiseAdd(spW, spMask, spOp, spA, spB, spOpDesc);
+    spExpr->Dependency(spWriteA, spEAddAB);
+    spExpr->Dependency(spWriteB, spEAddAB);
+    spExpr->Dependency(spWriteMask, spEAddAB);
+    spExpr->Submit();
+    spExpr->Wait();
+    ASSERT_EQ(spExpr->GetState(), spla::Expression::State::Evaluated);
+
+    utils::Vector<Type> c = a.EWiseAdd(mask, true, b, op);
     ASSERT_TRUE(c.Equals(spW));
 }
 
@@ -134,7 +175,7 @@ void testOneIsEmpty(spla::Library &library, std::size_t M, std::size_t nvals,
     spExpr->Wait();
     ASSERT_EQ(spExpr->GetState(), spla::Expression::State::Evaluated);
 
-    utils::Vector<Type> c = a.EWiseAdd(mask, b, op);
+    utils::Vector<Type> c = a.EWiseAdd(mask, false, b, op);
     ASSERT_TRUE(c.Equals(spW1));
     ASSERT_TRUE(c.Equals(spW2));
 }
@@ -167,7 +208,7 @@ void testNoValues(spla::Library &library, std::size_t M, std::size_t nvals, std:
     spExpr->Wait();
     ASSERT_EQ(spExpr->GetState(), spla::Expression::State::Evaluated);
 
-    utils::Vector<unsigned char> c = a.EWiseAdd(mask, b, [](unsigned char a, unsigned char b) { return 0; });
+    utils::Vector<unsigned char> c = a.EWiseAdd(mask, false, b, [](unsigned char a, unsigned char b) { return 0; });
     ASSERT_TRUE(c.EqualsStructure(spW));
 }
 
@@ -187,6 +228,11 @@ void test(std::size_t M, std::size_t base, std::size_t step, std::size_t iter, c
             testMasked<float>(library, M, nvals, spT, spOp, op, i);
         }
 
+        for (std::size_t i = 0; i < iter; i++) {
+            std::size_t nvals = base + i * step;
+            testMaskedComplement<float>(library, M, nvals, spT, spOp, op, i);
+        }
+
         for (std::size_t i = 0; i < iter; i++) {
             std::size_t nvals = base + i * step;
             testOneIsEmpty<float>(library, M, nvals, spT, spOp, op, i);
@@ -208,6 +254,11 @@ void test(std::size_t M, std::size_t base, std::size_t step, std::size_t iter, c
             testMasked<std::int32_t>(library, M, nvals, spT, spOp, op, i);
         }
 
+        for (std::size_t i = 0; i < iter; i++) {
+            std::size_t nvals = base + i * step;
+            testMaskedComplement<std::int32_t>(library, M, nvals, spT, spOp, op, i);
+        }
+
         for (std::size_t i = 0; i < iter; i++) {
             std::size_t nvals = base + i * step;
             testOneIsEmpty<std::int32_t>(library, M, nvals, spT, spOp, op, i);
diff --git a/tests/TestVxM.cpp b/tests/TestVxM.cpp
index 34dfc0393..c13447777 100644
--- a/tests/TestVxM.cpp
+++ b/tests/TestVxM.cpp
@@ -103,7 +103,53 @@ void testMasked(spla::Library &library,
 
     ASSERT_EQ(spExpr->GetState(), spla::Expression::State::Evaluated);
 
-    utils::Vector<Type> c = utils::VxM(mask, a, b, multOp, addOp);
+    utils::Vector<Type> c = utils::VxM(mask, false, a, b, multOp, addOp);
+    ASSERT_TRUE(c.Equals(spW));
+}
+
+template<typename Type, typename MultOp, typename AddOp>
+void testMaskedComplement(spla::Library &library,
+                          std::size_t M, std::size_t N, std::size_t nvals,
+                          const spla::RefPtr<spla::Type> spT,
+                          const spla::RefPtr<spla::FunctionBinary> spMult,
+                          const spla::RefPtr<spla::FunctionBinary> spAdd,
+                          MultOp multOp, AddOp addOp,
+                          std::size_t seed = 0) {
+    utils::Vector a = utils::Vector<Type>::Generate(M, nvals, seed).SortReduceDuplicates();
+    utils::Matrix b = utils::Matrix<Type>::Generate(M, N, nvals, seed + 1).SortReduceDuplicates();
+    utils::Vector mask = utils::Vector<unsigned char>::Generate(N, nvals, seed + 2).SortReduceDuplicates();
+
+    a.Fill(utils::UniformGenerator<Type>());
+    b.Fill(utils::UniformGenerator<Type>());
+
+    auto spA = spla::Vector::Make(M, spT, library);
+    auto spB = spla::Matrix::Make(M, N, spT, library);
+    auto spW = spla::Vector::Make(N, spT, library);
+    auto spMask = spla::Vector::Make(N, spla::Types::Void(library), library);
+
+    // Specify, that values already in row order + no duplicates
+    auto spDesc = spla::Descriptor::Make(library);
+    spDesc->SetParam(spla::Descriptor::Param::ValuesSorted);
+    spDesc->SetParam(spla::Descriptor::Param::NoDuplicates);
+
+    // Use complementary mask
+    auto spOpDesc = spla::Descriptor::Make(library);
+    spOpDesc->SetParam(spla::Descriptor::Param::MaskComplement);
+
+    auto spExpr = spla::Expression::Make(library);
+    auto spWriteA = spExpr->MakeDataWrite(spA, a.GetData(library), spDesc);
+    auto spWriteB = spExpr->MakeDataWrite(spB, b.GetData(library), spDesc);
+    auto spWriteMask = spExpr->MakeDataWrite(spMask, mask.GetDataIndices(library), spDesc);
+    auto spVxM = spExpr->MakeVxM(spW, spMask, spMult, spAdd, spA, spB, spOpDesc);
+    spExpr->Dependency(spWriteA, spVxM);
+    spExpr->Dependency(spWriteB, spVxM);
+    spExpr->Dependency(spWriteMask, spVxM);
+    spExpr->Submit();
+    spExpr->Wait();
+
+    ASSERT_EQ(spExpr->GetState(), spla::Expression::State::Evaluated);
+
+    utils::Vector<Type> c = utils::VxM(mask, true, a, b, multOp, addOp);
     ASSERT_TRUE(c.Equals(spW));
 }
 
@@ -137,7 +183,7 @@ void testNoValues(spla::Library &library, std::size_t M, std::size_t N, std::siz
     ASSERT_EQ(spExpr->GetState(), spla::Expression::State::Evaluated);
 
     auto dummy = [](unsigned char a, unsigned char b) { return 0; };
-    utils::Vector<unsigned char> c = utils::VxM(mask, a, b, dummy, dummy);
+    utils::Vector<unsigned char> c = utils::VxM(mask, false, a, b, dummy, dummy);
     ASSERT_TRUE(c.EqualsStructure(spW));
 }
 
@@ -159,6 +205,11 @@ void test(std::size_t M, std::size_t N, std::size_t base, std::size_t step, std:
             std::size_t nvals = base + i * step;
             testMasked<T>(library, M, N, nvals, spT, spMult, spAdd, mult, add, i);
         }
+
+        for (std::size_t i = 0; i < iter; i++) {
+            std::size_t nvals = base + i * step;
+            testMaskedComplement<T>(library, M, N, nvals, spT, spMult, spAdd, mult, add, i);
+        }
     });
 
     utils::testBlocks(blocksSizes, [=](spla::Library &library) {
@@ -178,6 +229,11 @@ void test(std::size_t M, std::size_t N, std::size_t base, std::size_t step, std:
             std::size_t nvals = base + i * step;
             testMasked<T>(library, M, N, nvals, spT, spMult, spAdd, mult, add, i);
         }
+
+        for (std::size_t i = 0; i < iter; i++) {
+            std::size_t nvals = base + i * step;
+            testMaskedComplement<T>(library, M, N, nvals, spT, spMult, spAdd, mult, add, i);
+        }
     });
 
     utils::testBlocks(blocksSizes, [=](spla::Library &library) {
diff --git a/tests/utils/Matrix.hpp b/tests/utils/Matrix.hpp
index eca82e344..4fe7430ea 100644
--- a/tests/utils/Matrix.hpp
+++ b/tests/utils/Matrix.hpp
@@ -326,7 +326,7 @@ namespace utils {
         }
 
         template<typename M>
-        [[nodiscard]] Matrix Mask(const Matrix<M> &mask) const {
+        [[nodiscard]] Matrix Mask(const Matrix<M> &mask, bool complement) const {
             assert(GetNrows() == mask.GetNrows());
             assert(GetNcols() == mask.GetNcols());
 
@@ -339,20 +339,46 @@ namespace utils {
             std::size_t endA = GetNvals();
             std::size_t endB = mask.GetNvals();
 
-            while (a != endA && b != endB) {
-                auto thisCrd = std::pair{GetRows()[a], GetCols()[a]};
-                auto maskCrd = std::pair{mask.GetRows()[b], mask.GetCols()[b]};
+            if (complement) {
+                while (a != endA && b != endB) {
+                    auto thisCrd = std::pair{GetRows()[a], GetCols()[a]};
+                    auto maskCrd = std::pair{mask.GetRows()[b], mask.GetCols()[b]};
+
+                    if (thisCrd == maskCrd) {
+                        a += 1;
+                        b += 1;
+                    } else if (thisCrd < maskCrd) {
+                        rows.push_back(GetRows()[a]);
+                        cols.push_back(GetCols()[a]);
+                        vals.push_back(GetVals()[a]);
+                        a += 1;
+                    } else {
+                        b += 1;
+                    }
+                }
 
-                if (thisCrd == maskCrd) {
+                while (a != endA) {
                     rows.push_back(GetRows()[a]);
                     cols.push_back(GetCols()[a]);
                     vals.push_back(GetVals()[a]);
                     a += 1;
-                    b += 1;
-                } else if (thisCrd < maskCrd) {
-                    a += 1;
-                } else {
-                    b += 1;
+                }
+            } else {
+                while (a != endA && b != endB) {
+                    auto thisCrd = std::pair{GetRows()[a], GetCols()[a]};
+                    auto maskCrd = std::pair{mask.GetRows()[b], mask.GetCols()[b]};
+
+                    if (thisCrd == maskCrd) {
+                        rows.push_back(GetRows()[a]);
+                        cols.push_back(GetCols()[a]);
+                        vals.push_back(GetVals()[a]);
+                        a += 1;
+                        b += 1;
+                    } else if (thisCrd < maskCrd) {
+                        a += 1;
+                    } else {
+                        b += 1;
+                    }
                 }
             }
 
@@ -414,8 +440,8 @@ namespace utils {
         }
 
         template<typename M, typename BinaryOp>
-        [[nodiscard]] Matrix EWiseAdd(const Matrix<M> &mask, const Matrix<T> &other, BinaryOp op) const {
-            return Mask(mask).EWiseAdd(other.Mask(mask), op);
+        [[nodiscard]] Matrix EWiseAdd(const Matrix<M> &mask, bool complement, const Matrix<T> &other, BinaryOp op) const {
+            return Mask(mask, complement).EWiseAdd(other.Mask(mask, complement), op);
         }
 
         template<typename W, typename B, typename MultOp, typename AddOp>
@@ -465,9 +491,8 @@ namespace utils {
         }
 
         template<typename W, typename M, typename B, typename MultOp, typename AddOp>
-        [[nodiscard]] Matrix<W> MxM(const Matrix<M> &mask, const Matrix<B> &b, MultOp multOp, AddOp addOp) {
-
-            return MxM<W>(b, multOp, addOp).Mask(mask);
+        [[nodiscard]] Matrix<W> MxM(const Matrix<M> &mask, bool complement, const Matrix<B> &b, MultOp multOp, AddOp addOp) {
+            return MxM<W>(b, multOp, addOp).Mask(mask, complement);
         }
 
         static Matrix Empty(std::size_t nrows, std::size_t ncols) {
diff --git a/tests/utils/Operations.hpp b/tests/utils/Operations.hpp
index 99eea3261..7da0e6497 100644
--- a/tests/utils/Operations.hpp
+++ b/tests/utils/Operations.hpp
@@ -79,8 +79,8 @@ namespace utils {
     }
 
     template<typename T, typename M, typename Mult, typename Add>
-    inline Vector<T> VxM(const Vector<M> &mask, const Vector<T> &a, const Matrix<T> &b, Mult mult, Add add) {
-        return VxM(a, b, mult, add).Mask(mask);
+    inline Vector<T> VxM(const Vector<M> &mask, bool complement, const Vector<T> &a, const Matrix<T> &b, Mult mult, Add add) {
+        return VxM(a, b, mult, add).Mask(mask, complement);
     }
 
 }// namespace utils
diff --git a/tests/utils/Vector.hpp b/tests/utils/Vector.hpp
index febdf20eb..992c10965 100644
--- a/tests/utils/Vector.hpp
+++ b/tests/utils/Vector.hpp
@@ -245,7 +245,7 @@ namespace utils {
         }
 
         template<typename M>
-        [[nodiscard]] Vector Mask(const Vector<M> &mask) const {
+        [[nodiscard]] Vector Mask(const Vector<M> &mask, bool complement) const {
             assert(GetNrows() == mask.GetNrows());
 
             std::vector<Index> rows;
@@ -256,16 +256,37 @@ namespace utils {
             std::size_t endA = GetNvals();
             std::size_t endB = mask.GetNvals();
 
-            while (a != endA && b != endB) {
-                if (GetRows()[a] == mask.GetRows()[b]) {
+            if (complement) {
+                while (a != endA && b != endB) {
+                    if (GetRows()[a] == mask.GetRows()[b]) {
+                        a += 1;
+                        b += 1;
+                    } else if (GetRows()[a] < mask.GetRows()[b]) {
+                        rows.push_back(GetRows()[a]);
+                        vals.push_back(GetVals()[a]);
+                        a += 1;
+                    } else {
+                        b += 1;
+                    }
+                }
+
+                while (a != endA) {
                     rows.push_back(GetRows()[a]);
                     vals.push_back(GetVals()[a]);
                     a += 1;
-                    b += 1;
-                } else if (GetRows()[a] < mask.GetRows()[b]) {
-                    a += 1;
-                } else {
-                    b += 1;
+                }
+            } else {
+                while (a != endA && b != endB) {
+                    if (GetRows()[a] == mask.GetRows()[b]) {
+                        rows.push_back(GetRows()[a]);
+                        vals.push_back(GetVals()[a]);
+                        a += 1;
+                        b += 1;
+                    } else if (GetRows()[a] < mask.GetRows()[b]) {
+                        a += 1;
+                    } else {
+                        b += 1;
+                    }
                 }
             }
 
@@ -317,8 +338,8 @@ namespace utils {
         }
 
         template<typename M, typename BinaryOp>
-        [[nodiscard]] Vector EWiseAdd(const Vector<M> &mask, const Vector<T> &other, BinaryOp op) const {
-            return Mask(mask).EWiseAdd(other.Mask(mask), op);
+        [[nodiscard]] Vector EWiseAdd(const Vector<M> &mask, bool complement, const Vector<T> &other, BinaryOp op) const {
+            return Mask(mask, complement).EWiseAdd(other.Mask(mask, complement), op);
         }
 
         static Vector Generate(std::size_t nrows, std::size_t nvals, std::size_t seed = 0, const T &value = T()) {