Repository: SparseLinearAlgebra/spla (language: python)

Issue:
Vector-matrix multiplication algo (coo format)
Implement vector-matrix multiplication algorithm to multiply vector-matrix blocks.

## Expected result

`w[j] = v[1] * m[1, j] + ....  + v[m] * m[m, j]`

## Tasks

- [x] Setup VxMCOO algo
- [x] Setup ParamsVxM structure
- [x] Implement algorithm
- [x] Support masking
- [x] Handle no-values case


Discussion / hints:

Fix patch:
diff --git a/include/spla-cpp/SplaExpression.hpp b/include/spla-cpp/SplaExpression.hpp
index 83ad1359e..43af7170e 100644
--- a/include/spla-cpp/SplaExpression.hpp
+++ b/include/spla-cpp/SplaExpression.hpp
@@ -253,14 +253,14 @@ namespace spla {
 
         /**
          * Make masked matrix-matrix multiplication expression node.
-         * Operation is evaluated as `w<mask> = a x(mult, add) b`.
+         * Operation is evaluated as `w<mask> = a mxm(mult, add) b`.
          * 
          * @note If mask is null, the whole result is saved
          * @note Mult must have signature `f: ta x tb -> tw`
          * @note Add must have signature `f: tw x tw -> tw`
          *
          * @param w Matrix to store result
-         * @param mask Mask to filter input matrices; may be null
+         * @param mask Mask to filter result; may be null
          * @param mult Binary function used to multiply a and b values
          * @param add Binary function used to accumulate multiplication results
          * @param a Input a matrix
@@ -277,6 +277,32 @@ namespace spla {
                                        const RefPtr<Matrix> &b,
                                        const RefPtr<Descriptor> &desc = nullptr);
 
+        /**
+         * Make masked vector-matrix multiplication expression node.
+         * Operation is evaluated as `w<mask> = a vxm(mult, add) b`.
+         *
+         * @note If mask is null, the whole result is saved
+         * @note Mult must have signature `f: ta x tb -> tw`
+         * @note Add must have signature `f: tw x tw -> tw`
+         *
+         * @param w Vector to store result
+         * @param mask Mask to filter result; may be null
+         * @param mult Binary function used to multiply a and b values
+         * @param add Binary function used to accumulate multiplication results
+         * @param a Input a vector
+         * @param b Input b matrix
+         * @param desc Operation descriptor
+         *
+         * @return Created expression node
+         */
+        RefPtr<ExpressionNode> MakeVxM(const RefPtr<Vector> &w,
+                                       const RefPtr<Vector> &mask,
+                                       const RefPtr<FunctionBinary> &mult,
+                                       const RefPtr<FunctionBinary> &add,
+                                       const RefPtr<Vector> &a,
+                                       const RefPtr<Matrix> &b,
+                                       const RefPtr<Descriptor> &desc = nullptr);
+
         /** @return Current expression state */
         State GetState() const;
 
diff --git a/sources/CMakeLists.txt b/sources/CMakeLists.txt
index 8fe08ec66..4832ed8a7 100644
--- a/sources/CMakeLists.txt
+++ b/sources/CMakeLists.txt
@@ -22,19 +22,24 @@ set(SPLA_ALGORITHM_SOURCES
         sources/algo/mxm/SplaMxMCOO.hpp
         sources/algo/vector/SplaVectorEWiseAddCOO.cpp
         sources/algo/vector/SplaVectorEWiseAddCOO.hpp
+        sources/algo/vxm/SplaVxMCOO.cpp
+        sources/algo/vxm/SplaVxMCOO.hpp
         sources/algo/SplaAlgorithm.hpp
         sources/algo/SplaAlgorithmManager.cpp
         sources/algo/SplaAlgorithmManager.hpp
         sources/algo/SplaAlgorithmParams.hpp)
 
 set(SPLA_COMPUTE_SOURCES
+        sources/compute/SplaApplyMask.hpp
         sources/compute/SplaGather.hpp
         sources/compute/SplaIndicesToRowOffsets.hpp
         sources/compute/SplaMaskByKey.hpp
         sources/compute/SplaMergeByKey.hpp
         sources/compute/SplaReduceByKey.hpp
         sources/compute/SplaReduceDuplicates.hpp
+        sources/compute/SplaSortByRow.hpp
         sources/compute/SplaSortByRowColumn.hpp
+        sources/compute/SplaTransformValues.hpp
         )
 
 set(SPLA_CORE_SOURCES
@@ -58,6 +63,8 @@ set(SPLA_EXPRESSION_SOURCES
         sources/expression/matrix/SplaMatrixEWiseAdd.hpp
         sources/expression/prod/SplaMxM.cpp
         sources/expression/prod/SplaMxM.hpp
+        sources/expression/prod/SplaVxM.cpp
+        sources/expression/prod/SplaVxM.hpp
         sources/expression/scalar/SplaScalarDataWrite.cpp
         sources/expression/scalar/SplaScalarDataWrite.hpp
         sources/expression/scalar/SplaScalarDataRead.cpp
diff --git a/sources/SplaExpression.cpp b/sources/SplaExpression.cpp
index b0244587e..5ca21375a 100644
--- a/sources/SplaExpression.cpp
+++ b/sources/SplaExpression.cpp
@@ -291,6 +291,37 @@ spla::Expression::MakeMxM(const spla::RefPtr<spla::Matrix> &w,
                     desc);
 }
 
+spla::RefPtr<spla::ExpressionNode>
+spla::Expression::MakeVxM(const spla::RefPtr<spla::Vector> &w,
+                          const spla::RefPtr<spla::Vector> &mask,
+                          const spla::RefPtr<spla::FunctionBinary> &mult,
+                          const spla::RefPtr<spla::FunctionBinary> &add,
+                          const spla::RefPtr<spla::Vector> &a,
+                          const spla::RefPtr<spla::Matrix> &b,
+                          const spla::RefPtr<spla::Descriptor> &desc) {
+    CHECK_RAISE_ERROR(w.IsNotNull(), NullPointer, "w can't be null");
+    CHECK_RAISE_ERROR(a.IsNotNull(), NullPointer, "a can't be null");
+    CHECK_RAISE_ERROR(b.IsNotNull(), NullPointer, "b can't be null");
+    CHECK_RAISE_ERROR(w->GetNrows() == b->GetNcols(), DimensionMismatch, "Incompatible size");
+    CHECK_RAISE_ERROR(a->GetNrows() == b->GetNrows(), DimensionMismatch, "Incompatible size");
+    CHECK_RAISE_ERROR(mask.IsNull() || w->GetNrows() == mask->GetNrows(), DimensionMismatch, "Incompatible size");
+    CHECK_RAISE_ERROR(mult.IsNull() || mult->CanApply(*a, *b, *w), InvalidType, "Cannot apply `mult` op to provided objects");
+    CHECK_RAISE_ERROR(add.IsNull() || add->CanApply(*w, *w, *w), InvalidType, "Cannot apply `add` op to provided objects");
+    CHECK_RAISE_ERROR(!w->GetType()->HasValues() || (mult.IsNotNull() && add.IsNotNull()), NullPointer, "If type has values, `mult` and `add` op must be provided");
+
+    std::vector<RefPtr<Object>> args = {
+            w.As<Object>(),
+            mask.As<Object>(),
+            mult.As<Object>(),
+            add.As<Object>(),
+            a.As<Object>(),
+            b.As<Object>()};
+
+    return MakeNode(ExpressionNode::Operation::VxM,
+                    std::move(args),
+                    desc);
+}
+
 void spla::Expression::SetState(State state) {
     mState.store(state);
 }
diff --git a/sources/algo/SplaAlgorithmManager.cpp b/sources/algo/SplaAlgorithmManager.cpp
index 4097df4dc..eecb3568d 100644
--- a/sources/algo/SplaAlgorithmManager.cpp
+++ b/sources/algo/SplaAlgorithmManager.cpp
@@ -29,6 +29,7 @@
 #include <algo/matrix/SplaMatrixEWiseAddCOO.hpp>
 #include <algo/mxm/SplaMxMCOO.hpp>
 #include <algo/vector/SplaVectorEWiseAddCOO.hpp>
+#include <algo/vxm/SplaVxMCOO.hpp>
 #include <cassert>
 #include <core/SplaError.hpp>
 
@@ -36,6 +37,7 @@ spla::AlgorithmManager::AlgorithmManager(Library &library) : mLibrary(library) {
     Register(new MatrixEWiseAddCOO());
     Register(new VectorEWiseAddCOO());
     Register(new MxMCOO());
+    Register(new VxMCOO());
 }
 
 void spla::AlgorithmManager::Register(const spla::RefPtr<spla::Algorithm> &algo) {
diff --git a/sources/algo/SplaAlgorithmParams.hpp b/sources/algo/SplaAlgorithmParams.hpp
index e1053fe8d..507429032 100644
--- a/sources/algo/SplaAlgorithmParams.hpp
+++ b/sources/algo/SplaAlgorithmParams.hpp
@@ -105,6 +105,23 @@ namespace spla {
         RefPtr<Type> tw;
     };
 
+    /** Blocked vector-matrix multiply params */
+    class ParamsVxM final : public AlgorithmParams {
+    public:
+        ~ParamsVxM() override = default;
+
+        bool hasMask = false;       // true if must apply mask
+        RefPtr<VectorBlock> w;      // tw
+        RefPtr<VectorBlock> mask;   // if has mask, must apply this
+        RefPtr<FunctionBinary> mult;// f: ta x tb -> tw
+        RefPtr<FunctionBinary> add; // f: tw x tw -> tw
+        RefPtr<VectorBlock> a;      // ta
+        RefPtr<MatrixBlock> b;      // tb
+        RefPtr<Type> ta;
+        RefPtr<Type> tb;
+        RefPtr<Type> tw;
+    };
+
     /**
      * @}
      */
diff --git a/sources/algo/vector/SplaVectorEWiseAddCOO.cpp b/sources/algo/vector/SplaVectorEWiseAddCOO.cpp
index c2dc1a8fe..eed10f389 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 &logger = library->GetLogger();
 
     auto device = library->GetDeviceManager().GetDevice(p->deviceId);
     compute::context ctx = library->GetContext();
@@ -218,6 +219,7 @@ void spla::VectorEWiseAddCOO::Process(spla::AlgorithmParams &params) {
                                        queue);
 
     p->w = VectorCOO::Make(blockA->GetNrows(), resultNvals, std::move(resultRows), std::move(resultVals)).As<VectorBlock>();
+    SPDLOG_LOGGER_TRACE(logger, "Merge vectors size={} nnz={}", blockA->GetNrows(), resultNvals);
 }
 
 spla::Algorithm::Type spla::VectorEWiseAddCOO::GetType() const {
diff --git a/sources/algo/vxm/SplaVxMCOO.cpp b/sources/algo/vxm/SplaVxMCOO.cpp
new file mode 100644
index 000000000..2a8a1cef2
--- /dev/null
+++ b/sources/algo/vxm/SplaVxMCOO.cpp
@@ -0,0 +1,184 @@
+/**********************************************************************************/
+/* This file is part of spla project                                              */
+/* https://github.com/JetBrains-Research/spla                                     */
+/**********************************************************************************/
+/* MIT License                                                                    */
+/*                                                                                */
+/* Copyright (c) 2021 JetBrains-Research                                          */
+/*                                                                                */
+/* Permission is hereby granted, free of charge, to any person obtaining a copy   */
+/* of this software and associated documentation files (the "Software"), to deal  */
+/* in the Software without restriction, including without limitation the rights   */
+/* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      */
+/* copies of the Software, and to permit persons to whom the Software is          */
+/* furnished to do so, subject to the following conditions:                       */
+/*                                                                                */
+/* The above copyright notice and this permission notice shall be included in all */
+/* copies or substantial portions of the Software.                                */
+/*                                                                                */
+/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     */
+/* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       */
+/* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    */
+/* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         */
+/* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  */
+/* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  */
+/* SOFTWARE.                                                                      */
+/**********************************************************************************/
+
+#include <algo/vxm/SplaVxMCOO.hpp>
+#include <boost/compute/algorithm.hpp>
+#include <boost/compute/algorithm/scatter_if.hpp>
+#include <compute/SplaApplyMask.hpp>
+#include <compute/SplaIndicesToRowOffsets.hpp>
+#include <compute/SplaReduceByKey.hpp>
+#include <compute/SplaReduceDuplicates.hpp>
+#include <compute/SplaSortByRow.hpp>
+#include <compute/SplaTransformValues.hpp>
+#include <core/SplaLibraryPrivate.hpp>
+#include <core/SplaQueueFinisher.hpp>
+#include <storage/block/SplaMatrixCOO.hpp>
+#include <storage/block/SplaVectorCOO.hpp>
+
+bool spla::VxMCOO::Select(const spla::AlgorithmParams &params) const {
+    auto p = dynamic_cast<const ParamsVxM *>(&params);
+
+    return p &&
+           p->w.Is<VectorCOO>() &&
+           p->mask.Is<VectorCOO>() &&
+           p->a.Is<VectorCOO>() &&
+           p->b.Is<MatrixCOO>();
+}
+
+void spla::VxMCOO::Process(spla::AlgorithmParams &params) {
+    using namespace boost;
+
+    auto p = dynamic_cast<ParamsVxM *>(&params);
+    auto w = p->w;
+    auto library = p->desc->GetLibrary().GetPrivatePtr();
+
+    auto device = library->GetDeviceManager().GetDevice(p->deviceId);
+    compute::context ctx = library->GetContext();
+    compute::command_queue queue(ctx, device);
+    QueueFinisher finisher(queue);
+
+    auto a = p->a.Cast<VectorCOO>();
+    auto b = p->b.Cast<MatrixCOO>();
+    auto mask = p->mask.Cast<VectorCOO>();
+
+    if (p->hasMask && mask.IsNull())
+        return;
+
+    if (a->GetNvals() == 0 || b->GetNvals() == 0)
+        return;
+
+    const auto &ta = p->ta;
+    const auto &tb = p->tb;
+    const auto &tw = p->tw;
+    auto hasValues = tw->HasValues();
+    auto M = a->GetNrows();
+    auto N = b->GetNcols();
+
+    // Compute rows offsets and rows lengths for matrix b
+    compute::vector<unsigned int> offsets(ctx);
+    compute::vector<unsigned int> lengths(ctx);
+    IndicesToRowOffsets(b->GetRows(), offsets, lengths, M, queue);
+
+    // 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);
+    compute::exclusive_scan(segmentLengths.begin(), segmentLengths.end(), outputPtr.begin(), 0u, queue);
+
+    // Number of products to count
+    std::size_t cooNnz = (outputPtr.end() - 1).read(queue);
+
+    // nothing to do, no a[i] * b[i,:] product
+    if (!cooNnz)
+        return;
+
+    // Determine location of a and b values (indices in coo arrays) to copy for products evaluation
+    compute::vector<unsigned int> aLocations(cooNnz, ctx);
+    compute::vector<unsigned int> bLocations(cooNnz, ctx);
+
+    compute::fill(aLocations.begin(), aLocations.end(), 0u, queue);
+    compute::scatter_if(compute::counting_iterator<unsigned int>(0),
+                        compute::counting_iterator<unsigned int>(a->GetNvals()),
+                        outputPtr.begin(),
+                        segmentLengths.begin(),
+                        aLocations.begin(),
+                        queue);
+    compute::inclusive_scan(aLocations.begin(), aLocations.end(), aLocations.begin(), compute::max<unsigned int>(), queue);
+
+    auto &aRows = a->GetRows();
+    BOOST_COMPUTE_CLOSURE(void, unfoldSegment, (unsigned int i), (outputPtr, offsets, aRows, aLocations, bLocations), {
+        uint locationOfRowIndex = aLocations[i];
+        uint rowIdx = aRows[locationOfRowIndex];
+        uint rowBaseOffset = offsets[rowIdx];
+        uint offsetOfRowSegment = outputPtr[locationOfRowIndex];
+        bLocations[i] = rowBaseOffset + (i - offsetOfRowSegment);
+    });
+    compute::for_each_n(compute::counting_iterator<unsigned int>(0), cooNnz, unfoldSegment, queue);
+
+    // Gather indices j for each product a[i] * b[i,j]
+    compute::vector<unsigned int> J(cooNnz, ctx);
+    compute::gather(bLocations.begin(), bLocations.end(), b->GetCols().begin(), J.begin(), queue);
+
+    // Store final result here
+    compute::vector<unsigned int> rows(ctx);
+    compute::vector<unsigned char> vals(ctx);
+
+    // Compute a[i] * b[i, j] for each value i and j
+    if (hasValues) {
+        compute::vector<unsigned char> V(cooNnz * tw->GetByteSize(), ctx);
+        TransformValues(aLocations, bLocations,
+                        a->GetVals(), b->GetVals(), V,
+                        ta->GetByteSize(),
+                        tb->GetByteSize(),
+                        tw->GetByteSize(),
+                        p->mult->GetSource(),
+                        queue);
+
+        // Sort a[i] * b[i, j] products, so all j products stored in sequence
+        SortByRow(J, V, tw->GetByteSize(), queue);
+
+        // 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) {
+            compute::vector<unsigned int> tmpRows(ctx);
+            compute::vector<unsigned char> tmpVals(ctx);
+            ApplyMask(mask->GetRows(), rows, vals, tmpRows, tmpVals, tw->GetByteSize(), queue);
+            std::swap(rows, tmpRows);
+            std::swap(vals, tmpVals);
+        }
+    } else {
+        // Sort result indices
+        compute::sort(J.begin(), J.end(), queue);
+
+        // Reduce duplicates (keep only first entry)
+        ReduceDuplicates(J, rows, queue);
+
+        // Apply mask to indices
+        if (p->hasMask) {
+            compute::vector<unsigned int> tmpRows(ctx);
+            MaskKeys(mask->GetRows(), rows, tmpRows, 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>();
+}
+
+spla::Algorithm::Type spla::VxMCOO::GetType() const {
+    return Type::VxM;
+}
+
+std::string spla::VxMCOO::GetName() const {
+    return "VxMCOO";
+}
\ No newline at end of file
diff --git a/sources/algo/vxm/SplaVxMCOO.hpp b/sources/algo/vxm/SplaVxMCOO.hpp
new file mode 100644
index 000000000..786a0df7a
--- /dev/null
+++ b/sources/algo/vxm/SplaVxMCOO.hpp
@@ -0,0 +1,46 @@
+/**********************************************************************************/
+/* This file is part of spla project                                              */
+/* https://github.com/JetBrains-Research/spla                                     */
+/**********************************************************************************/
+/* MIT License                                                                    */
+/*                                                                                */
+/* Copyright (c) 2021 JetBrains-Research                                          */
+/*                                                                                */
+/* Permission is hereby granted, free of charge, to any person obtaining a copy   */
+/* of this software and associated documentation files (the "Software"), to deal  */
+/* in the Software without restriction, including without limitation the rights   */
+/* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      */
+/* copies of the Software, and to permit persons to whom the Software is          */
+/* furnished to do so, subject to the following conditions:                       */
+/*                                                                                */
+/* The above copyright notice and this permission notice shall be included in all */
+/* copies or substantial portions of the Software.                                */
+/*                                                                                */
+/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     */
+/* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       */
+/* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    */
+/* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         */
+/* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  */
+/* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  */
+/* SOFTWARE.                                                                      */
+/**********************************************************************************/
+
+#ifndef SPLA_SPLAVXMCOO_HPP
+#define SPLA_SPLAVXMCOO_HPP
+
+#include <algo/SplaAlgorithm.hpp>
+
+namespace spla {
+
+    class VxMCOO final : public Algorithm {
+    public:
+        ~VxMCOO() override = default;
+        bool Select(const AlgorithmParams &params) const override;
+        void Process(AlgorithmParams &params) override;
+        Type GetType() const override;
+        std::string GetName() const override;
+    };
+
+}// namespace spla
+
+#endif//SPLA_SPLAVXMCOO_HPP
diff --git a/sources/compute/SplaApplyMask.hpp b/sources/compute/SplaApplyMask.hpp
new file mode 100644
index 000000000..ff32c55d2
--- /dev/null
+++ b/sources/compute/SplaApplyMask.hpp
@@ -0,0 +1,101 @@
+/**********************************************************************************/
+/* This file is part of spla project                                              */
+/* https://github.com/JetBrains-Research/spla                                     */
+/**********************************************************************************/
+/* MIT License                                                                    */
+/*                                                                                */
+/* Copyright (c) 2021 JetBrains-Research                                          */
+/*                                                                                */
+/* Permission is hereby granted, free of charge, to any person obtaining a copy   */
+/* of this software and associated documentation files (the "Software"), to deal  */
+/* in the Software without restriction, including without limitation the rights   */
+/* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      */
+/* copies of the Software, and to permit persons to whom the Software is          */
+/* furnished to do so, subject to the following conditions:                       */
+/*                                                                                */
+/* The above copyright notice and this permission notice shall be included in all */
+/* copies or substantial portions of the Software.                                */
+/*                                                                                */
+/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     */
+/* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       */
+/* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    */
+/* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         */
+/* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  */
+/* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  */
+/* SOFTWARE.                                                                      */
+/**********************************************************************************/
+
+#ifndef SPLA_SPLAAPPLYMASK_HPP
+#define SPLA_SPLAAPPLYMASK_HPP
+
+#include <boost/compute/command_queue.hpp>
+#include <boost/compute/container/vector.hpp>
+#include <compute/SplaGather.hpp>
+#include <compute/SplaMaskByKey.hpp>
+
+namespace spla {
+
+    /**
+     * @addtogroup Internal
+     * @{
+     */
+
+    /**
+     * @brief Apply mask to coo vector.
+     *
+     * @param mask Mask indices
+     * @param inputRows Row indices
+     * @param inputVals Values
+     * @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 queue Command queue to perform operation
+     */
+    inline void ApplyMask(const boost::compute::vector<unsigned int> &mask,
+                          const boost::compute::vector<unsigned int> &inputRows,
+                          const boost::compute::vector<unsigned char> &inputVals,
+                          boost::compute::vector<unsigned int> &outputRows,
+                          boost::compute::vector<unsigned char> &outputVals,
+                          std::size_t byteSize,
+                          boost::compute::command_queue &queue) {
+        using namespace boost;
+
+        if (mask.empty() || inputRows.empty())
+            return;
+
+        auto ctx = queue.get_context();
+        auto typeHasValues = byteSize != 0;
+        auto nnz = inputRows.size();
+        compute::vector<unsigned int> perm(ctx);
+
+        if (typeHasValues) {
+            // Fill permutation to link value to indices
+            perm.resize(nnz, queue);
+            compute::copy_n(compute::make_counting_iterator(0), nnz, perm.begin(), queue);
+        }
+
+        if (typeHasValues) {
+            // Mask with new permutation buffer
+            compute::vector<unsigned int> outputPerm(ctx);
+            MaskByKeys(mask,
+                       inputRows, perm,
+                       outputRows, outputPerm,
+                       queue);
+            std::swap(perm, outputPerm);
+        } else
+            MaskKeys(mask, inputRows, outputRows, queue);
+
+        if (typeHasValues) {
+            // Copy output values using indices buffer
+            outputVals.resize(outputRows.size() * byteSize, queue);
+            Gather(perm.begin(), perm.end(), inputVals.begin(), outputVals.begin(), byteSize, queue);
+        }
+    }
+
+    /**
+     * @}
+     */
+
+}// namespace spla
+
+#endif//SPLA_SPLAAPPLYMASK_HPP
diff --git a/sources/compute/SplaIndicesToRowOffsets.hpp b/sources/compute/SplaIndicesToRowOffsets.hpp
index 366f9c376..77d6bc6b2 100644
--- a/sources/compute/SplaIndicesToRowOffsets.hpp
+++ b/sources/compute/SplaIndicesToRowOffsets.hpp
@@ -53,11 +53,11 @@ namespace spla {
      * @param n Number of rows in matrix
      * @param queue Command queue to execute
      */
-    void IndicesToRowOffsets(const boost::compute::vector<unsigned int> &indices,
-                             boost::compute::vector<unsigned int> &offsets,
-                             boost::compute::vector<unsigned int> &lengths,
-                             std::size_t n,
-                             boost::compute::command_queue &queue) {
+    inline void IndicesToRowOffsets(const boost::compute::vector<unsigned int> &indices,
+                                    boost::compute::vector<unsigned int> &offsets,
+                                    boost::compute::vector<unsigned int> &lengths,
+                                    std::size_t n,
+                                    boost::compute::command_queue &queue) {
         using namespace boost;
 
         lengths.resize(n + 1, queue);
diff --git a/sources/compute/SplaSortByRow.hpp b/sources/compute/SplaSortByRow.hpp
new file mode 100644
index 000000000..953904395
--- /dev/null
+++ b/sources/compute/SplaSortByRow.hpp
@@ -0,0 +1,88 @@
+/**********************************************************************************/
+/* This file is part of spla project                                              */
+/* https://github.com/JetBrains-Research/spla                                     */
+/**********************************************************************************/
+/* MIT License                                                                    */
+/*                                                                                */
+/* Copyright (c) 2021 JetBrains-Research                                          */
+/*                                                                                */
+/* Permission is hereby granted, free of charge, to any person obtaining a copy   */
+/* of this software and associated documentation files (the "Software"), to deal  */
+/* in the Software without restriction, including without limitation the rights   */
+/* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      */
+/* copies of the Software, and to permit persons to whom the Software is          */
+/* furnished to do so, subject to the following conditions:                       */
+/*                                                                                */
+/* The above copyright notice and this permission notice shall be included in all */
+/* copies or substantial portions of the Software.                                */
+/*                                                                                */
+/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     */
+/* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       */
+/* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    */
+/* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         */
+/* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  */
+/* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  */
+/* SOFTWARE.                                                                      */
+/**********************************************************************************/
+
+#ifndef SPLA_SPLASORTBYROW_HPP
+#define SPLA_SPLASORTBYROW_HPP
+
+#include <boost/compute.hpp>
+#include <cassert>
+#include <compute/SplaGather.hpp>
+
+namespace spla {
+
+    /**
+     * @addtogroup Internal
+     * @{
+     */
+
+    /**
+     * @brief Sort vector data in coo format in row order.
+     * @note If elementsInSequence is 0 and vals is empty, sorts only vector indices.
+     *
+     * @param rows Vector with row indices to sort
+     * @param vals Vector with values byte data
+     * @param elementsInSequence Size in bytes of values in vals vector
+     * @param queue Queue to perform sort operation
+     */
+    inline void SortByRow(boost::compute::vector<unsigned int> &rows,
+                          boost::compute::vector<unsigned char> &vals,
+                          std::size_t elementsInSequence,
+                          boost::compute::command_queue &queue) {
+        using namespace boost;
+
+        compute::context ctx = queue.get_context();
+        std::size_t nvals = rows.size();
+        bool typeHasValues = elementsInSequence != 0;
+
+        assert(vals.size() == nvals * elementsInSequence);
+
+        // Use permutation buffer to synchronize position of indices of the same
+        // entries in all 2 buffers: rows and vals.
+        compute::vector<unsigned int> permutation(nvals, ctx);
+        compute::copy(compute::counting_iterator<unsigned int>(0),
+                      compute::counting_iterator<unsigned int>(nvals),
+                      permutation.begin(),
+                      queue);
+
+        // Sort in row order and then, using permutation, shuffle column indices
+        compute::sort_by_key(rows.begin(), rows.end(), permutation.begin(), queue);
+
+        // Copy values, using permutation
+        if (typeHasValues) {
+            compute::vector<unsigned char> valsTmp(nvals * elementsInSequence, ctx);
+            Gather(permutation.begin(), permutation.end(), vals.begin(), valsTmp.begin(), elementsInSequence, queue);
+            std::swap(vals, valsTmp);
+        }
+    }
+
+    /**
+     * @}
+     */
+
+}// namespace spla
+
+#endif//SPLA_SPLASORTBYROW_HPP
diff --git a/sources/compute/SplaTransformValues.hpp b/sources/compute/SplaTransformValues.hpp
new file mode 100644
index 000000000..07c157463
--- /dev/null
+++ b/sources/compute/SplaTransformValues.hpp
@@ -0,0 +1,153 @@
+/**********************************************************************************/
+/* This file is part of spla project                                              */
+/* https://github.com/JetBrains-Research/spla                                     */
+/**********************************************************************************/
+/* MIT License                                                                    */
+/*                                                                                */
+/* Copyright (c) 2021 JetBrains-Research                                          */
+/*                                                                                */
+/* Permission is hereby granted, free of charge, to any person obtaining a copy   */
+/* of this software and associated documentation files (the "Software"), to deal  */
+/* in the Software without restriction, including without limitation the rights   */
+/* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      */
+/* copies of the Software, and to permit persons to whom the Software is          */
+/* furnished to do so, subject to the following conditions:                       */
+/*                                                                                */
+/* The above copyright notice and this permission notice shall be included in all */
+/* copies or substantial portions of the Software.                                */
+/*                                                                                */
+/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     */
+/* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       */
+/* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    */
+/* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         */
+/* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  */
+/* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  */
+/* SOFTWARE.                                                                      */
+/**********************************************************************************/
+
+#ifndef SPLA_SPLATRANSFORMVALUES_HPP
+#define SPLA_SPLATRANSFORMVALUES_HPP
+
+#include <boost/compute/command_queue.hpp>
+#include <boost/compute/container/vector.hpp>
+
+namespace spla {
+
+    /**
+     * @addtogroup Internal
+     * @{
+     */
+
+    namespace detail {
+
+        class TransformValuesKernel : public boost::compute::detail::meta_kernel {
+        public:
+            TransformValuesKernel() : boost::compute::detail::meta_kernel("__spla_transform_values_kernel") {
+            }
+
+            template<typename InputMap,
+                     typename InputValues,
+                     typename OutputValues>
+            void SetRange(InputMap aMap,
+                          InputMap bMap,
+                          InputValues aValues,
+                          InputValues bValues,
+                          OutputValues outputValues,
+                          std::size_t aByteSize,
+                          std::size_t bByteSize,
+                          std::size_t resultByteSize,
+                          std::size_t count,
+                          const std::string &transformOp) {
+                using namespace boost;
+                mCount = count;
+
+                std::stringstream _spla_transform_op;
+                _spla_transform_op << "void _spla_transform_op(__global void* vp_a, __global void* vp_b, __global void* vp_c) {\n"
+                                   << "#define _ACCESS_A __global\n"
+                                   << "#define _ACCESS_B __global\n"
+                                   << "#define _ACCESS_C __global\n"
+                                   << "   " << transformOp << "\n"
+                                   << "#undef _ACCESS_A\n"
+                                   << "#undef _ACCESS_B\n"
+                                   << "#undef _ACCESS_C\n"
+                                   << "}";
+
+                add_function("_spla_transform_op", _spla_transform_op.str());
+
+                *this << "const uint i = get_global_id(0);\n"
+                      << "const uint a_idx = " << aMap[expr<compute::uint_>("i")] << ";\n"
+                      << "const uint b_idx = " << bMap[expr<compute::uint_>("i")] << ";\n"
+                      << "const uint a_offset = a_idx * " << aByteSize << ";\n"
+                      << "const uint b_offset = b_idx * " << bByteSize << ";\n"
+                      << "const uint result_offset = i * " << resultByteSize << ";\n"
+                      << "_spla_transform_op(&" << aValues[expr<compute::uint_>("a_offset")] << ", &" << bValues[expr<compute::uint_>("b_offset")] << ", &" << outputValues[expr<compute::uint_>("result_offset")] << ");\n";
+            }
+
+            boost::compute::event Exec(boost::compute::command_queue &queue) {
+                if (mCount == 0) {
+                    return boost::compute::event();
+                }
+
+                return exec_1d(queue, 0, mCount);
+            }
+
+        private:
+            std::size_t mCount = 0;
+        };
+
+    }// namespace detail
+
+    /**
+     * @brief Transforms a and b values using map and transfrom op
+     * Transforms sequence of a and b values using map to fetch actual
+     * a and b values. For each index i of computes result values as v[i] = transform(a[a_map[i]], b[b_map[i]]).
+     *
+     * @param aMap Map for a values
+     * @param bMap Map for b values
+     * @param aValues A values
+     * @param bValues B values
+     * @param values Where to store result of transformation; must be allocated by the user
+     * @param aByteSize Size of a values
+     * @param bByteSize Size of b values
+     * @param resultByteSize Size of result values
+     * @param transformOp Source code for transform function
+     * @param queue Command queue to perform operation
+     *
+     * @return Event to syn this operation
+     */
+    inline boost::compute::event TransformValues(const boost::compute::vector<unsigned int> &aMap,
+                                                 const boost::compute::vector<unsigned int> &bMap,
+                                                 const boost::compute::vector<unsigned char> &aValues,
+                                                 const boost::compute::vector<unsigned char> &bValues,
+                                                 boost::compute::vector<unsigned char> &values,
+                                                 std::size_t aByteSize,
+                                                 std::size_t bByteSize,
+                                                 std::size_t resultByteSize,
+                                                 const std::string &transformOp,
+                                                 boost::compute::command_queue &queue) {
+        using namespace boost;
+        auto count = aMap.size();
+
+        assert(aMap.size() == bMap.size());
+        assert(values.size() == count * resultByteSize);
+
+        detail::TransformValuesKernel kernel;
+        kernel.SetRange(aMap.begin(), bMap.begin(),
+                        aValues.begin(), bValues.begin(),
+                        values.begin(),
+                        aByteSize,
+                        bByteSize,
+                        resultByteSize,
+                        count,
+                        transformOp);
+
+        return kernel.Exec(queue);
+    }
+
+    /**
+     * @}
+     */
+
+}// namespace spla
+
+#endif//SPLA_SPLATRANSFORMVALUES_HPP
diff --git a/sources/expression/SplaExpressionManager.cpp b/sources/expression/SplaExpressionManager.cpp
index 4c34adf9c..655fb86a6 100644
--- a/sources/expression/SplaExpressionManager.cpp
+++ b/sources/expression/SplaExpressionManager.cpp
@@ -38,6 +38,7 @@
 #include <expression/matrix/SplaMatrixDataWrite.hpp>
 #include <expression/matrix/SplaMatrixEWiseAdd.hpp>
 #include <expression/prod/SplaMxM.hpp>
+#include <expression/prod/SplaVxM.hpp>
 #include <expression/scalar/SplaScalarDataRead.hpp>
 #include <expression/scalar/SplaScalarDataWrite.hpp>
 #include <expression/vector/SplaVectorDataRead.hpp>
@@ -56,6 +57,7 @@ spla::ExpressionManager::ExpressionManager(spla::Library &library) : mLibrary(li
     Register(new VectorDataRead());
     Register(new VectorEWiseAdd());
     Register(new MxM());
+    Register(new VxM());
 }
 
 void spla::ExpressionManager::Submit(const spla::RefPtr<spla::Expression> &expression) {
diff --git a/sources/expression/prod/SplaMxM.cpp b/sources/expression/prod/SplaMxM.cpp
index 82787e342..f250cbea7 100644
--- a/sources/expression/prod/SplaMxM.cpp
+++ b/sources/expression/prod/SplaMxM.cpp
@@ -232,7 +232,7 @@ void spla::MxM::Process(std::size_t nodeIdx, const spla::Expression &expression,
                     }
 
                     // Store final result
-                    MatrixStorage::Index index{static_cast<unsigned int>(i), static_cast<unsigned int>(j)};
+                    Index index{static_cast<unsigned int>(i), static_cast<unsigned int>(j)};
                     w->GetStorage()->SetBlock(index, block);
                 });
 
diff --git a/sources/expression/prod/SplaVxM.cpp b/sources/expression/prod/SplaVxM.cpp
new file mode 100644
index 000000000..ca0c26143
--- /dev/null
+++ b/sources/expression/prod/SplaVxM.cpp
@@ -0,0 +1,246 @@
+/**********************************************************************************/
+/* This file is part of spla project                                              */
+/* https://github.com/JetBrains-Research/spla                                     */
+/**********************************************************************************/
+/* MIT License                                                                    */
+/*                                                                                */
+/* Copyright (c) 2021 JetBrains-Research                                          */
+/*                                                                                */
+/* Permission is hereby granted, free of charge, to any person obtaining a copy   */
+/* of this software and associated documentation files (the "Software"), to deal  */
+/* in the Software without restriction, including without limitation the rights   */
+/* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      */
+/* copies of the Software, and to permit persons to whom the Software is          */
+/* furnished to do so, subject to the following conditions:                       */
+/*                                                                                */
+/* The above copyright notice and this permission notice shall be included in all */
+/* copies or substantial portions of the Software.                                */
+/*                                                                                */
+/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     */
+/* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       */
+/* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    */
+/* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         */
+/* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  */
+/* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  */
+/* SOFTWARE.                                                                      */
+/**********************************************************************************/
+
+#include <algo/SplaAlgorithmManager.hpp>
+#include <core/SplaLibraryPrivate.hpp>
+#include <expression/prod/SplaVxM.hpp>
+#include <storage/SplaMatrixStorage.hpp>
+#include <storage/SplaVectorStorage.hpp>
+
+namespace spla {
+    namespace {
+        /** Utility to fetch mask block from entry map */
+        inline RefPtr<VectorBlock> GetMaskBlock(VectorStorage::EntryMap &map, const VectorStorage::Index &idx) {
+            auto found = map.find(idx);
+            return found != map.end() ? found->second : RefPtr<VectorBlock>{};
+        }
+
+        /** Used to aggregate non empty results of a[i] x b[i,j] products for each (j) entry */
+        class ProductsResults {
+        public:
+            explicit ProductsResults(std::size_t n) : mBlocks(n) {}
+
+            void AddBlock(std::size_t j, const RefPtr<VectorBlock> &block) {
+                std::lock_guard<std::mutex> lockGuard(mMutex);
+                mBlocks[j].push_back(block);
+            }
+
+            void GetBlocks(std::size_t j, std::vector<RefPtr<VectorBlock>> &blocks) {
+                std::lock_guard<std::mutex> lockGuard(mMutex);
+                blocks = mBlocks[j];
+            }
+
+        private:
+            std::vector<std::vector<RefPtr<VectorBlock>>> mBlocks;
+            mutable std::mutex mMutex;
+        };
+    }// namespace
+}// namespace spla
+
+bool spla::VxM::Select(std::size_t nodeIdx, const spla::Expression &expression) {
+    return true;
+}
+
+void spla::VxM::Process(std::size_t nodeIdx, const spla::Expression &expression, spla::TaskBuilder &builder) {
+    auto &nodes = expression.GetNodes();
+    auto &node = nodes[nodeIdx];
+    auto library = node->GetLibrary().GetPrivatePtr();
+    auto logger = library->GetLogger();
+    auto &deviceMan = library->GetDeviceManager();
+
+    auto w = node->GetArg(0).Cast<Vector>();
+    auto mask = node->GetArg(1).Cast<Vector>();
+    auto mult = node->GetArg(2).Cast<FunctionBinary>();
+    auto add = node->GetArg(3).Cast<FunctionBinary>();
+    auto a = node->GetArg(4).Cast<Vector>();
+    auto b = node->GetArg(5).Cast<Matrix>();
+    auto desc = node->GetDescriptor();
+
+    assert(w.IsNotNull());
+    assert(a.IsNotNull());
+    assert(b.IsNotNull());
+    assert(desc.IsNotNull());
+
+    // Clear w, so by default empty result returned
+    w->GetStorage()->Clear();
+
+    auto ta = a->GetType();
+    auto tb = b->GetType();
+    auto tw = w->GetType();
+    auto hasMask = mask.IsNotNull();
+    auto &aStorage = a->GetStorage();
+    auto &bStorage = b->GetStorage();
+
+    // Query block storage dimensions for vectors and matrix
+    std::size_t nBlockM = aStorage->GetNblockRows(),
+                nBlockN = bStorage->GetNblockCols();
+
+    using IndexV = VectorStorage::Index;
+    using IndexM = MatrixStorage::Index;
+    struct ToProcess {
+        IndexV a;
+        IndexM b;
+    };
+
+    // Fetch blocks and store locally
+    VectorStorage::EntryMap aBlocks;
+    MatrixStorage::EntryMap bBlocks;
+    VectorStorage::EntryMap maskBlocks;
+    aStorage->GetBlocks(aBlocks);
+    bStorage->GetBlocks(bBlocks);
+
+    if (hasMask)
+        // If mask empty => does not apply mask at all
+        mask->GetStorage()->GetBlocks(maskBlocks);
+
+    // Determine number of block products and product pairs for each result block
+    std::size_t totalProducts = 0;
+    std::size_t totalBlocks = 0;
+    std::vector<std::vector<ToProcess>> blockProducts(nBlockN);
+    for (std::size_t j = 0; j < nBlockN; j++) {
+        auto &toProcess = blockProducts[j];
+        for (std::size_t i = 0; i < nBlockM; i++) {
+            IndexV aIndex{static_cast<unsigned int>(i)};
+            IndexM bIndex{static_cast<unsigned int>(i), static_cast<unsigned int>(j)};
+            auto aBlock = aBlocks.find(aIndex);
+            auto bBlock = bBlocks.find(bIndex);
+
+            // If has something to multiply in both a[i,k] and b[k,j] blocks
+            if (aBlock != aBlocks.end() && bBlock != bBlocks.end()) {
+                toProcess.push_back(ToProcess{aIndex, bIndex});
+                totalProducts += 1;
+            }
+        }
+        // Count number of potentially not empty final w[j] blocks
+        totalBlocks += toProcess.empty() ? 0 : 1;
+    }
+
+    // Edge case: if no products, return empty result
+    if (totalProducts == 0) {
+        return;
+    }
+
+    // Shared thread-safe storage to aggregate results of products
+    auto products = std::make_shared<ProductsResults>(nBlockN);
+
+    // Query required number of devices (strategy: device per product)
+    auto devicesForProducts = deviceMan.FetchDevices(totalProducts, node);
+
+    // Dispatch tasks to compute a.block[i] x b.block[i,j] products
+    std::size_t deviceToFetch = 0;
+    std::vector<std::vector<tf::Task>> blockProductsTasks(nBlockN);
+    for (std::size_t j = 0; j < nBlockN; j++) {
+        auto &tasks = blockProductsTasks[j];
+        for (auto &toProcess : blockProducts[j]) {
+            auto deviceId = devicesForProducts[deviceToFetch];
+            auto aIdx = toProcess.a;
+            auto bIdx = toProcess.b;
+            auto aBlock = aBlocks.find(aIdx)->second;
+            auto bBlock = bBlocks.find(bIdx)->second;
+            auto maskBlock = GetMaskBlock(maskBlocks, IndexV{bIdx.second});
+            auto task = builder.Emplace([=]() {
+                assert(aBlock->GetNrows() == bBlock->GetNrows());
+                ParamsVxM params;
+                params.desc = desc;
+                params.deviceId = deviceId;
+                params.hasMask = hasMask;
+                params.mask = maskBlock;
+                params.mult = mult;
+                params.add = add;
+                params.a = aBlock;
+                params.b = bBlock;
+                params.ta = ta;
+                params.tb = tb;
+                params.tw = tw;
+                library->GetAlgoManager()->Dispatch(Algorithm::Type::VxM, params);
+
+                if (params.w.IsNotNull()) {
+                    // If has not empty result, store it to sum later
+                    products->AddBlock(j, params.w);
+                    SPDLOG_LOGGER_TRACE(logger, "Blocks product ({})x({},{}) nnz={}",
+                                        aIdx, bIdx.first, bIdx.second, params.w->GetNvals());
+                }
+            });
+            deviceToFetch += 1;
+            tasks.push_back(std::move(task));
+        }
+    }
+
+    // Query required number of devices (strategy: device per not empty w[j])
+    auto devicesForFinalMerge = deviceMan.FetchDevices(totalBlocks, node);
+
+    // Finally, for each block w[i,j] we must aggregate intermediate
+    // blocks multiplications results as a series of element-wise additions of blocks
+    deviceToFetch = 0;
+    for (std::size_t j = 0; j < nBlockN; j++) {
+        auto &toProcess = blockProducts[j];
+        if (!toProcess.empty()) {
+            auto deviceId = devicesForFinalMerge[deviceToFetch];
+            auto task = builder.Emplace([=]() {
+                std::vector<RefPtr<VectorBlock>> blocks;
+                products->GetBlocks(j, blocks);
+
+                // Nothing to do, w[j] is empty
+                if (blocks.empty())
+                    return;
+
+                // Start to merge n blocks, number of merges n - 1
+                auto block = blocks[0];
+                for (std::size_t k = 1; k < blocks.size(); k++) {
+                    assert(block->GetNrows() == blocks[k]->GetNrows());
+                    ParamsVectorEWiseAdd params;
+                    params.desc = desc;
+                    params.deviceId = deviceId;
+                    params.op = add;
+                    params.a = block;
+                    params.b = blocks[k];
+                    params.type = tw;
+                    library->GetAlgoManager()->Dispatch(Algorithm::Type::VectorEWiseAdd, params);
+
+                    // Store result for next iteration
+                    block = params.w;
+                }
+
+                // Store final result
+                IndexV index{static_cast<unsigned int>(j)};
+                w->GetStorage()->SetBlock(index, block);
+            });
+
+            // Setup dependencies
+            // Start aggregation as soon as all sums are computed for w[j]
+            auto &deps = blockProductsTasks[j];
+            for (auto &parent : deps)
+                parent.precede(task);
+
+            deviceToFetch += 1;
+        }
+    }
+}
+
+spla::ExpressionNode::Operation spla::VxM::GetOperationType() const {
+    return ExpressionNode::Operation::VxM;
+}
diff --git a/sources/expression/prod/SplaVxM.hpp b/sources/expression/prod/SplaVxM.hpp
new file mode 100644
index 000000000..39dd0b0e7
--- /dev/null
+++ b/sources/expression/prod/SplaVxM.hpp
@@ -0,0 +1,45 @@
+/**********************************************************************************/
+/* This file is part of spla project                                              */
+/* https://github.com/JetBrains-Research/spla                                     */
+/**********************************************************************************/
+/* MIT License                                                                    */
+/*                                                                                */
+/* Copyright (c) 2021 JetBrains-Research                                          */
+/*                                                                                */
+/* Permission is hereby granted, free of charge, to any person obtaining a copy   */
+/* of this software and associated documentation files (the "Software"), to deal  */
+/* in the Software without restriction, including without limitation the rights   */
+/* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      */
+/* copies of the Software, and to permit persons to whom the Software is          */
+/* furnished to do so, subject to the following conditions:                       */
+/*                                                                                */
+/* The above copyright notice and this permission notice shall be included in all */
+/* copies or substantial portions of the Software.                                */
+/*                                                                                */
+/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     */
+/* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       */
+/* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    */
+/* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         */
+/* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  */
+/* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  */
+/* SOFTWARE.                                                                      */
+/**********************************************************************************/
+
+#ifndef SPLA_SPLAVXM_HPP
+#define SPLA_SPLAVXM_HPP
+
+#include <expression/SplaNodeProcessor.hpp>
+
+namespace spla {
+
+    class VxM final : public NodeProcessor {
+    public:
+        ~VxM() override = default;
+        bool Select(std::size_t nodeIdx, const Expression &expression) override;
+        void Process(std::size_t nodeIdx, const Expression &expression, TaskBuilder &builder) override;
+        ExpressionNode::Operation GetOperationType() const override;
+    };
+
+}// namespace spla
+
+#endif//SPLA_SPLAVXM_HPP
diff --git a/sources/storage/SplaMatrixStorage.cpp b/sources/storage/SplaMatrixStorage.cpp
index f99e92099..06c4dba4f 100644
--- a/sources/storage/SplaMatrixStorage.cpp
+++ b/sources/storage/SplaMatrixStorage.cpp
@@ -69,7 +69,6 @@ void spla::MatrixStorage::GetBlocks(spla::MatrixStorage::EntryMap &entryMap) con
     entryMap = mBlocks;
 }
 
-
 void spla::MatrixStorage::GetBlocks(spla::MatrixStorage::EntryRowList &entryList) const {
     std::lock_guard<std::mutex> lock(mMutex);
     entryList.clear();
diff --git a/sources/storage/SplaVectorStorage.cpp b/sources/storage/SplaVectorStorage.cpp
index b0617dfe4..b42cad8d5 100644
--- a/sources/storage/SplaVectorStorage.cpp
+++ b/sources/storage/SplaVectorStorage.cpp
@@ -29,6 +29,12 @@
 #include <core/SplaMath.hpp>
 #include <storage/SplaVectorStorage.hpp>
 
+void spla::VectorStorage::Clear() {
+    std::lock_guard<std::mutex> lock(mMutex);
+    mNvals = 0;
+    mBlocks.clear();
+}
+
 void spla::VectorStorage::SetBlock(const spla::VectorStorage::Index &index, const spla::RefPtr<spla::VectorBlock> &block) {
     assert(index < mNblockRows);
     assert(block.IsNotNull());
@@ -49,6 +55,11 @@ void spla::VectorStorage::GetBlocks(spla::VectorStorage::EntryList &entryList) c
     entryList.insert(entryList.begin(), mBlocks.begin(), mBlocks.end());
 }
 
+void spla::VectorStorage::GetBlocks(EntryMap &entryMap) const {
+    std::lock_guard<std::mutex> lock(mMutex);
+    entryMap = mBlocks;
+}
+
 void spla::VectorStorage::GetBlocksGrid(std::size_t &rows) const {
     rows = mNblockRows;
 }
diff --git a/sources/storage/SplaVectorStorage.hpp b/sources/storage/SplaVectorStorage.hpp
index d0f1fc630..ebc79a801 100644
--- a/sources/storage/SplaVectorStorage.hpp
+++ b/sources/storage/SplaVectorStorage.hpp
@@ -54,6 +54,9 @@ namespace spla {
 
         ~VectorStorage() override = default;
 
+        /** Remove all blocks from storage (empty matrix) */
+        void Clear();
+
         /** Set block at specified block index */
         void SetBlock(const Index &index, const RefPtr<VectorBlock> &block);
 
@@ -63,6 +66,9 @@ namespace spla {
         /** Get list of non-null presented blocks in storage */
         void GetBlocks(EntryList &entryList) const;
 
+        /** Get map of non-null presented blocks in storage */
+        void GetBlocks(EntryMap &entryMap) const;
+
         /** Get blocks grid (total number of blocks) */
         void GetBlocksGrid(std::size_t &rows) const;
 


Test patch:
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index dec9f7e7f..502a7f18d 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -18,3 +18,4 @@ spla_test_target(TestMxM)
 spla_test_target(TestReduceByKey)
 spla_test_target(TestReduceDuplicates)
 spla_test_target(TestVectorEWiseAdd)
+spla_test_target(TestVxM)
diff --git a/tests/TestVxM.cpp b/tests/TestVxM.cpp
new file mode 100644
index 000000000..34dfc0393
--- /dev/null
+++ b/tests/TestVxM.cpp
@@ -0,0 +1,209 @@
+/**********************************************************************************/
+/* This file is part of spla project                                              */
+/* https://github.com/JetBrains-Research/spla                                     */
+/**********************************************************************************/
+/* MIT License                                                                    */
+/*                                                                                */
+/* Copyright (c) 2021 JetBrains-Research                                          */
+/*                                                                                */
+/* Permission is hereby granted, free of charge, to any person obtaining a copy   */
+/* of this software and associated documentation files (the "Software"), to deal  */
+/* in the Software without restriction, including without limitation the rights   */
+/* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      */
+/* copies of the Software, and to permit persons to whom the Software is          */
+/* furnished to do so, subject to the following conditions:                       */
+/*                                                                                */
+/* The above copyright notice and this permission notice shall be included in all */
+/* copies or substantial portions of the Software.                                */
+/*                                                                                */
+/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     */
+/* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       */
+/* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    */
+/* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         */
+/* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  */
+/* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  */
+/* SOFTWARE.                                                                      */
+/**********************************************************************************/
+
+#include <Testing.hpp>
+
+template<typename Type, typename MultOp, typename AddOp>
+void testCommon(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();
+
+    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);
+
+    // 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 spExpr = spla::Expression::Make(library);
+    auto spWriteA = spExpr->MakeDataWrite(spA, a.GetData(library), spDesc);
+    auto spWriteB = spExpr->MakeDataWrite(spB, b.GetData(library), spDesc);
+    auto spVxM = spExpr->MakeVxM(spW, nullptr, spMult, spAdd, spA, spB);
+    spExpr->Dependency(spWriteA, spVxM);
+    spExpr->Dependency(spWriteB, spVxM);
+    spExpr->Submit();
+    spExpr->Wait();
+
+    ASSERT_EQ(spExpr->GetState(), spla::Expression::State::Evaluated);
+
+    utils::Vector<Type> c = utils::VxM(a, b, multOp, addOp);
+    ASSERT_TRUE(c.Equals(spW));
+}
+
+template<typename Type, typename MultOp, typename AddOp>
+void testMasked(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);
+
+    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);
+    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, a, b, multOp, addOp);
+    ASSERT_TRUE(c.Equals(spW));
+}
+
+void testNoValues(spla::Library &library, std::size_t M, std::size_t N, std::size_t nvals, std::size_t seed = 0) {
+    utils::Vector a = utils::Vector<unsigned char>::Generate(M, nvals, seed).SortReduceDuplicates();
+    utils::Matrix b = utils::Matrix<unsigned char>::Generate(M, N, nvals, seed + 1).SortReduceDuplicates();
+    utils::Vector mask = utils::Vector<unsigned char>::Generate(N, nvals, seed + 2).SortReduceDuplicates();
+
+    auto spT = spla::Types::Void(library);
+    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, spT, 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 spExpr = spla::Expression::Make(library);
+    auto spWriteA = spExpr->MakeDataWrite(spA, a.GetDataIndices(library), spDesc);
+    auto spWriteB = spExpr->MakeDataWrite(spB, b.GetDataIndices(library), spDesc);
+    auto spWriteMask = spExpr->MakeDataWrite(spMask, mask.GetDataIndices(library), spDesc);
+    auto spVxM = spExpr->MakeVxM(spW, spMask, nullptr, nullptr, spA, spB);
+    spExpr->Dependency(spWriteA, spVxM);
+    spExpr->Dependency(spWriteB, spVxM);
+    spExpr->Dependency(spWriteMask, spVxM);
+    spExpr->Submit();
+    spExpr->Wait();
+
+    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);
+    ASSERT_TRUE(c.EqualsStructure(spW));
+}
+
+void test(std::size_t M, std::size_t N, std::size_t base, std::size_t step, std::size_t iter, const std::vector<std::size_t> &blocksSizes) {
+    utils::testBlocks(blocksSizes, [=](spla::Library &library) {
+        using T = float;
+        auto spT = spla::Types::Float32(library);
+        auto spMult = spla::Functions::MultFloat32(library);
+        auto spAdd = spla::Functions::PlusFloat32(library);
+        auto mult = [](T a, T b) { return a * b; };
+        auto add = [](T a, T b) { return a + b; };
+
+        for (std::size_t i = 0; i < iter; i++) {
+            std::size_t nvals = base + i * step;
+            testCommon<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;
+            testMasked<T>(library, M, N, nvals, spT, spMult, spAdd, mult, add, i);
+        }
+    });
+
+    utils::testBlocks(blocksSizes, [=](spla::Library &library) {
+        using T = std::int32_t;
+        auto spT = spla::Types::Int32(library);
+        auto spMult = spla::Functions::MultInt32(library);
+        auto spAdd = spla::Functions::PlusInt32(library);
+        auto mult = [](T a, T b) { return a * b; };
+        auto add = [](T a, T b) { return a + b; };
+
+        for (std::size_t i = 0; i < iter; i++) {
+            std::size_t nvals = base + i * step;
+            testCommon<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;
+            testMasked<T>(library, M, N, nvals, spT, spMult, spAdd, mult, add, i);
+        }
+    });
+
+    utils::testBlocks(blocksSizes, [=](spla::Library &library) {
+        for (std::size_t i = 0; i < iter; i++) {
+            std::size_t nvals = base + i * step;
+            testNoValues(library, M, N, nvals, i);
+        }
+    });
+}
+
+TEST(VxM, Small) {
+    std::vector<std::size_t> blockSizes = {100, 1000};
+    std::size_t M = 120, N = 80;
+    test(M, N, M, M, 10, blockSizes);
+}
+
+TEST(MxM, Medium) {
+    std::vector<std::size_t> blockSizes = {1000, 10000};
+    std::size_t M = 1220, N = 880;
+    test(M, N, M, M, 10, blockSizes);
+}
+
+TEST(MxM, Large) {
+    std::vector<std::size_t> blockSizes = {10000, 100000};
+    std::size_t M = 12400, N = 8080;
+    test(M, N, M, M, 5, blockSizes);
+}
+
+SPLA_GTEST_MAIN
\ No newline at end of file
diff --git a/tests/Testing.hpp b/tests/Testing.hpp
index aa53645b5..2b21eb0c3 100644
--- a/tests/Testing.hpp
+++ b/tests/Testing.hpp
@@ -31,6 +31,7 @@
 #include <gtest/gtest.h>
 #include <utils/Compute.hpp>
 #include <utils/Matrix.hpp>
+#include <utils/Operations.hpp>
 #include <utils/Random.hpp>
 #include <utils/Setup.hpp>
 #include <utils/Typetraits.hpp>
diff --git a/tests/utils/Compute.hpp b/tests/utils/Compute.hpp
index ce0be6666..0938a72b2 100644
--- a/tests/utils/Compute.hpp
+++ b/tests/utils/Compute.hpp
@@ -34,7 +34,7 @@
 namespace utils {
 
     template<typename T>
-    std::vector<T> toOffsets(std::size_t n, const std::vector<T> &indices) {
+    std::vector<T> ToOffsets(std::size_t n, const std::vector<T> &indices) {
         std::vector<T> rowLengths(n + 1, 0);
         std::vector<T> offsets(n + 1);
         for (auto i : indices) rowLengths[i] += 1;
diff --git a/tests/utils/Matrix.hpp b/tests/utils/Matrix.hpp
index 23360996a..eca82e344 100644
--- a/tests/utils/Matrix.hpp
+++ b/tests/utils/Matrix.hpp
@@ -91,6 +91,18 @@ namespace utils {
             return mCols.data();
         }
 
+        [[nodiscard]] const std::vector<T> &GetValsVec() const {
+            return mVals;
+        }
+
+        [[nodiscard]] const std::vector<Index> &GetRowsVec() const {
+            return mRows;
+        }
+
+        [[nodiscard]] const std::vector<Index> &GetColsVec() const {
+            return mCols;
+        }
+
         [[nodiscard]] std::size_t GetNrows() const {
             return mNrows;
         }
@@ -221,7 +233,7 @@ namespace utils {
                 auto a = GetVals()[i];
                 auto b = spVals[i];
 
-                if ((useError && std::abs(a - b) > error) || (a != b))
+                if ((useError && std::abs(a - b) > error) || (!useError && a != b))
                     return false;
             }
 
@@ -414,8 +426,8 @@ namespace utils {
             std::vector<Index> cols;
             std::vector<W> vals;
 
-            auto aOffsets = toOffsets(GetNrows(), mRows);
-            auto bOffsets = toOffsets(b.GetNrows(), b.mRows);
+            auto aOffsets = ToOffsets(GetNrows(), mRows);
+            auto bOffsets = ToOffsets(b.GetNrows(), b.mRows);
 
             for (std::size_t i = 0; i < GetNrows(); i++) {
                 std::vector<Index> nnz;
@@ -454,6 +466,7 @@ 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);
         }
 
diff --git a/tests/utils/Operations.hpp b/tests/utils/Operations.hpp
new file mode 100644
index 000000000..99eea3261
--- /dev/null
+++ b/tests/utils/Operations.hpp
@@ -0,0 +1,88 @@
+/**********************************************************************************/
+/* This file is part of spla project                                              */
+/* https://github.com/JetBrains-Research/spla                                     */
+/**********************************************************************************/
+/* MIT License                                                                    */
+/*                                                                                */
+/* Copyright (c) 2021 JetBrains-Research                                          */
+/*                                                                                */
+/* Permission is hereby granted, free of charge, to any person obtaining a copy   */
+/* of this software and associated documentation files (the "Software"), to deal  */
+/* in the Software without restriction, including without limitation the rights   */
+/* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      */
+/* copies of the Software, and to permit persons to whom the Software is          */
+/* furnished to do so, subject to the following conditions:                       */
+/*                                                                                */
+/* The above copyright notice and this permission notice shall be included in all */
+/* copies or substantial portions of the Software.                                */
+/*                                                                                */
+/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     */
+/* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       */
+/* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    */
+/* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         */
+/* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  */
+/* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  */
+/* SOFTWARE.                                                                      */
+/**********************************************************************************/
+
+#ifndef SPLA_OPERATIONS_HPP
+#define SPLA_OPERATIONS_HPP
+
+#include <utils/Matrix.hpp>
+#include <utils/Vector.hpp>
+
+namespace utils {
+
+    template<typename T, typename Mult, typename Add>
+    inline Vector<T> VxM(const Vector<T> &a, const Matrix<T> &b, Mult mult, Add add) {
+        assert(a.GetNrows() == b.GetNrows());
+
+        using Index = typename Vector<T>::Index;
+
+        std::vector<Index> rows;
+        std::vector<T> vals;
+
+        std::size_t n = b.GetNcols();
+        std::vector<T> tmp(n);
+        std::vector<bool> mask(n, false);
+        std::vector<Index> nnz;
+
+        auto bOffsets = ToOffsets(b.GetNrows(), b.GetRowsVec());
+
+        for (std::size_t ak = 0; ak < a.GetNvals(); ak++) {
+            auto aRowId = a.GetRowsVec()[ak];
+            auto aVal = a.GetValsVec()[ak];
+
+            for (std::size_t bk = bOffsets[aRowId]; bk < bOffsets[aRowId + 1]; bk++) {
+                auto bColId = b.GetColsVec()[bk];
+                auto bVal = b.GetValsVec()[bk];
+
+                if (mask[bColId])
+                    tmp[bColId] = add(tmp[bColId], mult(aVal, bVal));
+                else {
+                    mask[bColId] = true;
+                    nnz.push_back(bColId);
+                    tmp[bColId] = mult(aVal, bVal);
+                }
+            }
+        }
+
+        if (!nnz.empty()) {
+            std::sort(nnz.begin(), nnz.end());
+            for (auto k : nnz) {
+                rows.push_back(k);
+                vals.push_back(tmp[k]);
+            }
+        }
+
+        return Vector<T>(n, std::move(rows), std::move(vals));
+    }
+
+    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);
+    }
+
+}// namespace utils
+
+#endif//SPLA_OPERATIONS_HPP
diff --git a/tests/utils/Typetraits.hpp b/tests/utils/Typetraits.hpp
index 5170b300c..77a1351a7 100644
--- a/tests/utils/Typetraits.hpp
+++ b/tests/utils/Typetraits.hpp
@@ -39,7 +39,7 @@ namespace utils {
     T GetError();
 
     template<>
-    bool UseError<float>() { return false; }
+    bool UseError<float>() { return true; }
 
     template<>
     float GetError<float>() { return 1e-5f; }
diff --git a/tests/utils/Vector.hpp b/tests/utils/Vector.hpp
index ac490ada5..febdf20eb 100644
--- a/tests/utils/Vector.hpp
+++ b/tests/utils/Vector.hpp
@@ -79,6 +79,14 @@ namespace utils {
             return mRows.data();
         }
 
+        [[nodiscard]] const std::vector<T> &GetValsVec() const {
+            return mVals;
+        }
+
+        [[nodiscard]] const std::vector<Index> &GetRowsVec() const {
+            return mRows;
+        }
+
         [[nodiscard]] std::size_t GetNrows() const {
             return mNrows;
         }
@@ -150,7 +158,9 @@ namespace utils {
             }
 
             if (v->GetNvals() != GetNvals()) {
-                std::cout << "Number of nnz mismatched" << std::endl;
+                std::cout << "Number of nnz mismatched "
+                          << "expected=" << GetNvals() << " "
+                          << "actual=" << v->GetNvals() << std::endl;
                 return false;
             }
 
@@ -192,7 +202,7 @@ namespace utils {
                 auto a = GetVals()[i];
                 auto b = spVals[i];
 
-                if ((useError && std::abs(a - b) > error) || (a != b))
+                if ((useError && std::abs(a - b) > error) || (!useError && a != b))
                     return false;
             }