Repository: SparseLinearAlgebra/spla (language: python)

Issue:
Matrix-matrix masked multiplication (csr x csc)



Discussion / hints:

Fix patch:
diff --git a/examples/spla_tc.cpp b/examples/spla_tc.cpp
index a3118c84b..74d7a15e7 100644
--- a/examples/spla_tc.cpp
+++ b/examples/spla_tc.cpp
@@ -109,9 +109,13 @@ int main(int argc, const char *const *argv) {
     prepareData->SubmitWait();
     SPLA_ALGO_CHECK(prepareData);
 
-    // Release A to free used memory
+    // Release A to free unused in algo memory
     A.Reset();
 
+    // Set for U transposed decoration U^T (in case of undirected graph exactly L)
+    U->SetDecoration(spla::Decorated::Decoration::TransposedMatrix, L.As<spla::Object>());
+
+    // Descriptor to profile time if required
     spla::RefPtr<spla::Descriptor> desc = spla::Descriptor::Make(library);
     desc->SetParam(spla::Descriptor::Param::ProfileTime, debugTiming);
 
diff --git a/include/spla-algo/SplaAlgoTc.hpp b/include/spla-algo/SplaAlgoTc.hpp
index c60503f12..226a626b5 100644
--- a/include/spla-algo/SplaAlgoTc.hpp
+++ b/include/spla-algo/SplaAlgoTc.hpp
@@ -43,7 +43,7 @@ namespace spla {
      * @brief Triangle counting
      *
      * As an input algorithm expects undirected graph in form of adjacency matrix
-     * with values of type int32_t with 1 in the cells (i,j) when graph has eject (i,j).
+     * with values of type int32_t with 1 in the cells (i,j) when graph has edge (i,j).
      * Adjacency matrix must passed as lower and upper triangular parts @p spL and @p spU.
      *
      * As an output algorithm produces number of triangles @p ntrins in the graph
diff --git a/sources/CMakeLists.txt b/sources/CMakeLists.txt
index a8b394d44..55cf7c2e0 100644
--- a/sources/CMakeLists.txt
+++ b/sources/CMakeLists.txt
@@ -33,6 +33,8 @@ set(SPLA_ALGORITHM_SOURCES
         sources/algo/matrix/SplaMatrixTriaCOO.hpp
         sources/algo/mxm/SplaMxMCOO.cpp
         sources/algo/mxm/SplaMxMCOO.hpp
+        sources/algo/mxm/SplaMxMMaskedCSRCSC.cpp
+        sources/algo/mxm/SplaMxMMaskedCSRCSC.hpp
         sources/algo/vector/SplaVectorAssignCOO.cpp
         sources/algo/vector/SplaVectorAssignCOO.hpp
         sources/algo/vector/SplaVectorEWiseAddCOO.cpp
diff --git a/sources/algo/SplaAlgorithmManager.cpp b/sources/algo/SplaAlgorithmManager.cpp
index 7fb5eb19b..7370bede7 100644
--- a/sources/algo/SplaAlgorithmManager.cpp
+++ b/sources/algo/SplaAlgorithmManager.cpp
@@ -44,6 +44,7 @@
 #include <algo/vector/SplaVectorToDenseCOO.hpp>
 
 #include <algo/mxm/SplaMxMCOO.hpp>
+#include <algo/mxm/SplaMxMMaskedCSRCSC.hpp>
 #include <algo/vxm/SplaVxMCOO.hpp>
 #include <algo/vxm/SplaVxMCOOStructure.hpp>
 
@@ -59,6 +60,7 @@ spla::AlgorithmManager::AlgorithmManager(Library &library) : mLibrary(library) {
     Register(new VectorReduceCOO());
     Register(new VectorEWiseAddCOO());
     Register(new VectorEWiseAddDense());
+    Register(new MxMMaskedCSRCSC());
     Register(new MxMCOO());
     Register(new VxMCOOStructure());
     Register(new VxMCOO());
diff --git a/sources/algo/SplaAlgorithmParams.hpp b/sources/algo/SplaAlgorithmParams.hpp
index 45cbaf74b..e17299916 100644
--- a/sources/algo/SplaAlgorithmParams.hpp
+++ b/sources/algo/SplaAlgorithmParams.hpp
@@ -115,6 +115,7 @@ namespace spla {
         RefPtr<FunctionBinary> add; // f: tw x tw -> tw
         RefPtr<MatrixBlock> a;      // ta
         RefPtr<MatrixBlock> b;      // tb
+        RefPtr<MatrixBlock> bT;     // transposed, type tb
         RefPtr<Type> ta;
         RefPtr<Type> tb;
         RefPtr<Type> tw;
diff --git a/sources/algo/matrix/SplaMatrixTransposeCOO.cpp b/sources/algo/matrix/SplaMatrixTransposeCOO.cpp
index aa34544d9..469e69505 100644
--- a/sources/algo/matrix/SplaMatrixTransposeCOO.cpp
+++ b/sources/algo/matrix/SplaMatrixTransposeCOO.cpp
@@ -33,6 +33,7 @@
 #include <core/SplaQueueFinisher.hpp>
 #include <storage/SplaMatrixStorage.hpp>
 #include <storage/block/SplaMatrixCOO.hpp>
+#include <storage/block/SplaMatrixCSR.hpp>
 
 bool spla::MatrixTransposeCOO::Select(const spla::AlgorithmParams &params) const {
     auto p = dynamic_cast<const ParamsTranspose *>(&params);
@@ -105,7 +106,7 @@ void spla::MatrixTransposeCOO::Process(spla::AlgorithmParams &params) {
         auto nrowsT = a->GetNcols();
         auto ncolsT = a->GetNrows();
         auto nvalsT = rows.size();
-        p->w = MatrixCOO::Make(nrowsT, ncolsT, nvalsT, std::move(rows), std::move(cols), std::move(vals)).As<MatrixBlock>();
+        p->w = MatrixCSR::Make(nrowsT, ncolsT, nvalsT, std::move(rows), std::move(cols), std::move(vals), queue).As<MatrixBlock>();
     }
 }
 
diff --git a/sources/algo/mxm/SplaMxMCOO.cpp b/sources/algo/mxm/SplaMxMCOO.cpp
index 16f575e77..278486464 100644
--- a/sources/algo/mxm/SplaMxMCOO.cpp
+++ b/sources/algo/mxm/SplaMxMCOO.cpp
@@ -149,10 +149,13 @@ namespace spla::detail {
 bool spla::MxMCOO::Select(const spla::AlgorithmParams &params) const {
     auto p = dynamic_cast<const ParamsMxM *>(&params);
 
+    // Requires a and b be presented in csr format
     return p &&
            p->mask.Is<MatrixCSR>() &&
            p->a.Is<MatrixCSR>() &&
-           p->b.Is<MatrixCSR>();
+           p->b.Is<MatrixCSR>() &&
+           p->a.IsNotNull() &&
+           p->b.IsNotNull();
 }
 
 void spla::MxMCOO::Process(spla::AlgorithmParams &algoParams) {
diff --git a/sources/algo/mxm/SplaMxMMaskedCSRCSC.cpp b/sources/algo/mxm/SplaMxMMaskedCSRCSC.cpp
new file mode 100644
index 000000000..d1d05f9cd
--- /dev/null
+++ b/sources/algo/mxm/SplaMxMMaskedCSRCSC.cpp
@@ -0,0 +1,424 @@
+/**********************************************************************************/
+/* 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/mxm/SplaMxMMaskedCSRCSC.hpp>
+#include <compute/metautil/SplaMetaUtil.hpp>
+#include <core/SplaLibraryPrivate.hpp>
+#include <core/SplaQueueFinisher.hpp>
+#include <storage/block/SplaMatrixCOO.hpp>
+#include <storage/block/SplaMatrixCSR.hpp>
+#include <utils/SplaProfiling.hpp>
+
+#define HAS_VALUES_BEGIN if (hasValues) {
+#define HAS_VALUES_END }
+
+bool spla::MxMMaskedCSRCSC::Select(const spla::AlgorithmParams &params) const {
+    auto p = dynamic_cast<const ParamsMxM *>(&params);
+
+    // Requires a and bT be in csr format with mask in csr
+    // todo: support complement mask (its seems easy,
+    //  need to iterate over the row [0, N) and skip (i,j),
+    //  which belong to the mask)
+    return p &&
+           !p->desc->IsParamSet(Descriptor::Param::MaskComplement) &&
+           p->mask.Is<MatrixCSR>() &&
+           p->a.Is<MatrixCSR>() &&
+           p->bT.Is<MatrixCSR>() &&
+           p->mask.IsNotNull() &&
+           p->a.IsNotNull() &&
+           p->bT.IsNotNull();
+}
+
+void spla::MxMMaskedCSRCSC::Process(spla::AlgorithmParams &algoParams) {
+    using namespace boost;
+
+    PF_SCOPE(mxm, "-mxm-m-csr-csc-");
+
+    auto params = dynamic_cast<ParamsMxM *>(&algoParams);
+    auto library = params->desc->GetLibrary().GetPrivatePtr();
+    auto device = library->GetDeviceManager().GetDevice(params->deviceId);
+    auto hasValues = params->tw->HasValues();
+
+    compute::context ctx = library->GetContext();
+    compute::command_queue queue(ctx, device);
+    QueueFinisher finisher(queue, library->GetLogger());
+
+    // Algorithm strategy:
+    // 1. Allocate temporary storage of size nnz(mask), since mask defines upper bound on result size
+    // 2. For each mask value (i,j) evaluate product of A[i,:] x B^T[j,:]
+    // 3. Use a warp (32/64 threads) per row
+    // 4. Compact temporary storage into final storage (remove holes)
+
+    MatrixCSR &m = *params->mask.Cast<MatrixCSR>();
+    MatrixCSR &a = *params->a.Cast<MatrixCSR>();
+    MatrixCSR &bT = *params->bT.Cast<MatrixCSR>();
+
+    std::size_t M = a.GetNrows(),
+                N = bT.GetNrows();
+
+    std::size_t tpb = 32;
+    std::size_t workSize = M * tpb;
+
+    assert(m.GetNrows() == M);
+    assert(m.GetNcols() == N);
+
+    // Auxiliary binary search to find bCol in the bT row by aCol
+    std::stringstream binSearch;
+    binSearch << "uint binary_search(const uint x, __global const uint* buffer, const uint start, const uint range) {\n"
+              << "    int l = start, r = start + range - 1;\n"
+              << "    while (l <= r) {\n"
+              << "        int k = (l + r) / 2;\n"
+              << "        if (buffer[k] < x) l = k + 1;\n"
+              << "        else if (x < buffer[k]) r = k - 1;\n"
+              << "        else return k;\n"
+              << "    }\n"
+              << "    return 0xffffffff;\n"
+              << "}\n";
+
+    PF_SCOPE_MARK(mxm, "setup");
+
+    const std::size_t BUCKETS_COUNT = 5;
+    compute::vector<Index> buckets(BUCKETS_COUNT, ctx);
+    compute::vector<Index> bucketsOffsets(BUCKETS_COUNT, ctx);
+    compute::vector<Index> bucketsConfig(ctx);
+    std::vector<Index> bucketsHost(BUCKETS_COUNT);
+
+    auto &rowsLenM = m.GetRowLengths();
+    BOOST_COMPUTE_CLOSURE(void, countBuckets, (unsigned int rid), (rowsLenM, buckets), {
+        const uint length = rowsLenM[rid];
+        if (length <= 0) return;
+        const uint bucketId = (uint) (clamp(ceil(log2((float) length)) - 1.0f, 0.0f, 4.0f));
+        atomic_add(&buckets[bucketId], 1u);
+    });
+    compute::fill(buckets.begin(), buckets.end(), 0u, queue);
+    compute::for_each_n(compute::counting_iterator<unsigned int>(0), M, countBuckets, queue);
+    compute::exclusive_scan(buckets.begin(), buckets.end(), bucketsOffsets.begin(), queue);
+    compute::copy(buckets.begin(), buckets.end(), bucketsHost.begin(), queue);
+
+    std::size_t rowsToProcess = std::reduce(bucketsHost.begin(), bucketsHost.end(), 0u);
+    bucketsConfig.resize(rowsToProcess, queue);
+
+    BOOST_COMPUTE_CLOSURE(void, fillBuckets, (unsigned int rid), (rowsLenM, buckets, bucketsOffsets, bucketsConfig), {
+        const uint length = rowsLenM[rid];
+        if (length <= 0) return;
+        const uint bucketId = (uint) (clamp(ceil(log2((float) length)) - 1.0f, 0.0f, 4.0f));
+        const uint offset = atomic_add(&buckets[bucketId], 1u);
+        bucketsConfig[bucketsOffsets[bucketId] + offset] = rid;
+    });
+    compute::fill(buckets.begin(), buckets.end(), 0u, queue);
+    compute::for_each_n(compute::counting_iterator<unsigned int>(0), M, fillBuckets, queue);
+
+    PF_SCOPE_MARK(mxm, "buckets");
+
+    // Rows len to evaluate exact size of each row (can be less or equal to the same mask row)
+    compute::vector<unsigned int> rowsLenC(M + 1, ctx);
+    compute::fill(rowsLenC.begin(), rowsLenC.end(), 0u, queue);
+
+    // Allocate temporary buffers for result (use mask nnz)
+    compute::vector<unsigned int> rowsC(m.GetNvals(), ctx);
+    compute::vector<unsigned int> colsC(m.GetNvals(), ctx);
+    compute::vector<unsigned char> valsC(hasValues ? m.GetNvals() * params->tw->GetByteSize() : 0, ctx);
+    {
+        PF_SCOPE_MARK(mxm, "allocate tmp");
+
+        compute::detail::meta_kernel kernel("__spla_masked_csr_csc_prod");
+        kernel.add_function("binary_search", binSearch.str());
+
+        auto argRowsPtrM = kernel.add_arg<const cl_uint *>(compute::memory_object::global_memory, "rowsPtrM");
+        auto argColsM = kernel.add_arg<const cl_uint *>(compute::memory_object::global_memory, "colsM");
+        auto argRowsPtrA = kernel.add_arg<const cl_uint *>(compute::memory_object::global_memory, "rowsPtrA");
+        auto argColsA = kernel.add_arg<const cl_uint *>(compute::memory_object::global_memory, "colsA");
+        auto argRowsPtrB = kernel.add_arg<const cl_uint *>(compute::memory_object::global_memory, "rowsPtrB");
+        auto argColsB = kernel.add_arg<const cl_uint *>(compute::memory_object::global_memory, "colsB");
+        auto argRowsLenC = kernel.add_arg<cl_uint *>(compute::memory_object::global_memory, "rowsLenC");
+        auto argRowsC = kernel.add_arg<cl_uint *>(compute::memory_object::global_memory, "rowsC");
+        auto argColsC = kernel.add_arg<cl_uint *>(compute::memory_object::global_memory, "colsC");
+
+        auto argCount = kernel.add_arg<cl_uint>("count");
+        auto argRptb = kernel.add_arg<cl_uint>("rptb");
+        auto argBucketId = kernel.add_arg<cl_uint>("bucketId");
+        auto argBucketOffsets = kernel.add_arg<const cl_uint *>(compute::memory_object::global_memory, "bucketOffsets");
+        auto argBucketConfig = kernel.add_arg<const cl_uint *>(compute::memory_object::global_memory, "bucketConfig");
+
+        std::size_t argValsA;
+        std::size_t argValsB;
+        std::size_t argValsC;
+
+        HAS_VALUES_BEGIN
+        argValsA = kernel.add_arg<const cl_uchar *>(compute::memory_object::global_memory, "valsA");
+        argValsB = kernel.add_arg<const cl_uchar *>(compute::memory_object::global_memory, "valsB");
+        argValsC = kernel.add_arg<cl_uchar *>(compute::memory_object::global_memory, "valsC");
+
+        kernel.add_function("mult_op", detail::meta::MakeFunction("mult_op",
+                                                                  params->mult->GetSource(),
+                                                                  detail::meta::Visibility::Global,
+                                                                  detail::meta::Visibility::Global,
+                                                                  detail::meta::Visibility::Unspecified));
+
+        kernel.add_function("add_op", detail::meta::MakeFunction("add_op",
+                                                                 params->add->GetSource(),
+                                                                 detail::meta::Visibility::Unspecified,
+                                                                 detail::meta::Visibility::Unspecified,
+                                                                 detail::meta::Visibility::Unspecified));
+
+        kernel.add_function("reduce_op", detail::meta::MakeFunction("reduce_op",
+                                                                    params->add->GetSource(),
+                                                                    detail::meta::Visibility::Local,
+                                                                    detail::meta::Visibility::Local,
+                                                                    detail::meta::Visibility::Local));
+        HAS_VALUES_END
+
+        kernel << "#define TPB " << static_cast<cl_uint>(tpb) << "\n";
+        HAS_VALUES_BEGIN
+        kernel << "#define BYTE_SIZE_A " << static_cast<cl_uint>(params->ta->GetByteSize()) << "\n"
+               << "#define BYTE_SIZE_B " << static_cast<cl_uint>(params->tb->GetByteSize()) << "\n"
+               << "#define BYTE_SIZE_C " << static_cast<cl_uint>(params->tw->GetByteSize()) << "\n"
+               << "__local ulong buff_reduce[(TPB * BYTE_SIZE_C) /" << static_cast<cl_uint>(sizeof(cl_ulong)) << "];\n"
+               << "__local uchar* p_reduce = (__local uchar*)buff_reduce;\n"
+               << "__local uint p_reduce_flag[TPB];\n"
+               << "uchar p_accum[BYTE_SIZE_C];\n";
+        HAS_VALUES_END
+        kernel << "const uint gid = get_group_id(0) * rptb;\n"
+               << "const uint tid = get_local_id(0);\n"
+               << "__local uint hasValue;\n"
+               << "for (uint rtp = gid; rtp < count && rtp < gid + rptb; rtp++) {\n"
+               << "    const rid = bucketConfig[bucketOffsets[bucketId] + rtp];\n"
+               << "    const uint mRowStart = rowsPtrM[rid];\n"
+               << "    const uint mRowEnd = rowsPtrM[rid + 1];\n"
+               << "    uint writePos = mRowStart;\n"
+               << "    for (uint mIdx = mRowStart; mIdx < mRowEnd; mIdx++) {\n"
+               << "        if (tid == 0) hasValue = 0;\n"
+               << "        barrier(CLK_LOCAL_MEM_FENCE);\n"
+               << "        const uint cid = colsM[mIdx];\n"
+               << "        const uint aRowStart = rowsPtrA[rid];\n"
+               << "        const uint aRowEnd = rowsPtrA[rid + 1];\n"
+               << "        const uint bRowStart = rowsPtrB[cid];\n"
+               << "        const uint bRowEnd = rowsPtrB[cid + 1];\n"
+               << "        const uint bRowLen = bRowEnd - bRowStart;\n"
+               << "        if (bRowLen > 0) {\n"
+               << "            uint seenValue = 0;\n"
+               << "            for (uint aIdx = aRowStart + tid; aIdx < aRowEnd; aIdx += TPB) {\n"
+               << "                const uint aCol = colsA[aIdx];\n"
+               << "                const uint bIdx = binary_search(aCol, colsB, bRowStart, bRowLen);\n"
+               << "                if (bIdx != 0xffffffff) {\n ";
+        HAS_VALUES_BEGIN
+        kernel << "                    if (seenValue == 0) {\n"
+               << "                        mult_op(&valsA[aIdx * BYTE_SIZE_A], &valsB[bIdx * BYTE_SIZE_B], &p_accum[0]);\n"
+               << "                    } else {\n"
+               << "                        uchar tmp[BYTE_SIZE_C];\n"
+               << "                        mult_op(&valsA[aIdx * BYTE_SIZE_A], &valsB[bIdx * BYTE_SIZE_B], &tmp[0]);\n"
+               << "                        add_op(&p_accum[0], &tmp[0], &p_accum[0]);\n"
+               << "                    }\n";
+        HAS_VALUES_END
+        kernel << "                    seenValue = 1;\n"
+               << "                }\n"
+               << "            }\n"
+               << "            atomic_or(&hasValue, seenValue);\n"
+               << "            barrier(CLK_LOCAL_MEM_FENCE);\n"
+               << "            if (hasValue) {\n";
+        HAS_VALUES_BEGIN
+        kernel << "                p_reduce_flag[tid] = seenValue;\n"
+               << "                for (uint i = 0; i < BYTE_SIZE_C; i++)\n"
+               << "                    p_reduce[tid * BYTE_SIZE_C + i] = p_accum[i];\n"
+               << "                for (uint stride = TPB/2; stride > 0; stride = stride / 2) {\n"
+               << "                    barrier(CLK_LOCAL_MEM_FENCE);"
+               << "                    if (tid < stride) {\n"
+               << "                        if (p_reduce_flag[tid] && p_reduce_flag[tid + stride])\n"
+               << "                            reduce_op(&p_reduce[tid * BYTE_SIZE_C], &p_reduce[(tid + stride) * BYTE_SIZE_C], &p_reduce[tid * BYTE_SIZE_C]);\n"
+               << "                        else if (p_reduce_flag[tid + stride]) {\n"
+               << "                            p_reduce_flag[tid] = 1;\n"
+               << "                            for (uint i = 0; i < BYTE_SIZE_C; i++)\n"
+               << "                                p_reduce[tid * BYTE_SIZE_C + i] = p_reduce[(tid + stride) * BYTE_SIZE_C + i];\n"
+               << "                        }\n"
+               << "                    }\n"
+               << "                }\n";
+        HAS_VALUES_END
+        kernel << "                if (tid == 0) {\n"
+               << "                    rowsC[writePos] = rid;\n"
+               << "                    colsC[writePos] = cid;\n";
+        HAS_VALUES_BEGIN
+        kernel << "                    for (uint i = 0; i < BYTE_SIZE_C; i++)\n"
+               << "                        valsC[writePos * BYTE_SIZE_C + i] = p_reduce[0 * BYTE_SIZE_C + i];\n";
+        HAS_VALUES_END
+        kernel << "                    writePos += 1;\n"
+               << "                    atomic_inc(&rowsLenC[rid]);\n"
+               << "                }\n"
+               << "            }\n"
+               << "        }\n"
+               << "    }\n"
+               << "}\n";
+
+        auto compiledKernel = kernel.compile(ctx);
+        compiledKernel.set_arg(argRowsPtrM, m.GetRowsOffsets().get_buffer());
+        compiledKernel.set_arg(argColsM, m.GetCols().get_buffer());
+        compiledKernel.set_arg(argRowsPtrA, a.GetRowsOffsets().get_buffer());
+        compiledKernel.set_arg(argColsA, a.GetCols().get_buffer());
+        compiledKernel.set_arg(argRowsPtrB, bT.GetRowsOffsets().get_buffer());
+        compiledKernel.set_arg(argColsB, bT.GetCols().get_buffer());
+        compiledKernel.set_arg(argRowsLenC, rowsLenC.get_buffer());
+        compiledKernel.set_arg(argRowsC, rowsC.get_buffer());
+        compiledKernel.set_arg(argColsC, colsC.get_buffer());
+
+        HAS_VALUES_BEGIN
+        compiledKernel.set_arg(argValsA, a.GetVals().get_buffer());
+        compiledKernel.set_arg(argValsB, bT.GetVals().get_buffer());
+        compiledKernel.set_arg(argValsC, valsC.get_buffer());
+        HAS_VALUES_END
+
+        compiledKernel.set_arg(argBucketConfig, bucketsConfig.get_buffer());
+        compiledKernel.set_arg(argBucketOffsets, bucketsOffsets.get_buffer());
+
+        auto dispatchGroup = [&](compute::command_queue &q, std::size_t bucketId, std::size_t rptb) {
+            std::size_t rowsToProcess = bucketsHost[bucketId];
+            if (rowsToProcess > 0) {
+                std::size_t workSize = (rowsToProcess / rptb + (rowsToProcess % rptb ? 1 : 0)) * tpb;
+                PF_SCOPE_SHOW(mxm, tpb << " " << rptb << " " << rowsToProcess);
+                compiledKernel.set_arg(argCount, static_cast<cl_uint>(rowsToProcess));
+                compiledKernel.set_arg(argRptb, static_cast<cl_uint>(rptb));
+                compiledKernel.set_arg(argBucketId, static_cast<cl_uint>(bucketId));
+                q.enqueue_1d_range_kernel(compiledKernel, 0, workSize, tpb);
+            }
+        };
+
+        struct GroupInfo {
+            boost::compute::command_queue queue;
+            std::size_t rptb;
+        };
+
+        GroupInfo infos[] = {
+                {compute::command_queue(ctx, device), 16},
+                {compute::command_queue(ctx, device), 8},
+                {compute::command_queue(ctx, device), 4},
+                {compute::command_queue(ctx, device), 2},
+                {compute::command_queue(ctx, device), 1}};
+
+        auto before = queue.enqueue_marker();
+
+        for (std::size_t bucketId = 0; bucketId < BUCKETS_COUNT; bucketId++) {
+            auto &info = infos[bucketId];
+            info.queue.enqueue_barrier(before);
+            dispatchGroup(info.queue, bucketId, info.rptb);
+            queue.enqueue_barrier(info.queue.enqueue_marker());
+        }
+
+        PF_SCOPE_MARK(mxm, "evaluate");
+    }
+
+    // Evaluate offsets for each row of the result
+    compute::vector<unsigned int> rowsPtrC(M + 1, ctx);
+    compute::exclusive_scan(rowsLenC.begin(), rowsLenC.end(), rowsPtrC.begin(), 0u, queue);
+
+    PF_SCOPE_MARK(mxm, "offsets");
+
+    // Get nnz of the result
+    std::size_t nvalsC = (rowsPtrC.end() - 1).read(queue);
+
+    assert(nvalsC <= m.GetNvals());
+
+    // Early exit, nothing to do
+    if (nvalsC <= 0)
+        return;
+
+    // Need to compact result
+    if (nvalsC < m.GetNvals()) {
+        // Allocate final buffers for result (use nvalsC)
+        compute::vector<unsigned int> dstRowsC(nvalsC, ctx);
+        compute::vector<unsigned int> dstColsC(nvalsC, ctx);
+        compute::vector<unsigned char> dstValsC(hasValues ? nvalsC * params->tw->GetByteSize() : 0, ctx);
+
+        PF_SCOPE_MARK(mxm, "allocate final");
+
+        compute::detail::meta_kernel kernel("__spla_masked_csr_csc_compact");
+
+        auto argSrcRowsPtrC = kernel.add_arg<const cl_uint *>(compute::memory_object::global_memory, "srcRowsPtrC");
+        auto argSrcRowsC = kernel.add_arg<const cl_uint *>(compute::memory_object::global_memory, "srcRowsC");
+        auto argSrcColsC = kernel.add_arg<const cl_uint *>(compute::memory_object::global_memory, "srcColsC");
+        auto argDstRowsPtrC = kernel.add_arg<const cl_uint *>(compute::memory_object::global_memory, "dstRowsPtrC");
+        auto argDstRowsC = kernel.add_arg<cl_uint *>(compute::memory_object::global_memory, "dstRowsC");
+        auto argDstColsC = kernel.add_arg<cl_uint *>(compute::memory_object::global_memory, "dstColsC");
+
+        std::size_t argSrcValsC;
+        std::size_t argDstValsC;
+
+        HAS_VALUES_BEGIN
+        argSrcValsC = kernel.add_arg<const cl_uchar *>(compute::memory_object::global_memory, "srcValsC");
+        argDstValsC = kernel.add_arg<const cl_uchar *>(compute::memory_object::global_memory, "dstValsC");
+        kernel << "#define BYTE_SIZE " << static_cast<cl_uint>(params->tw->GetByteSize()) << "\n";
+        HAS_VALUES_END
+
+        kernel << "#define TPB " << static_cast<cl_uint>(tpb) << "\n"
+               << "const uint rid = get_group_id(0);\n"
+               << "const uint tid = get_local_id(0);\n"
+               << "const uint srcRowStart = srcRowsPtrC[rid];\n"
+               << "const uint dstRowStart = dstRowsPtrC[rid];\n"
+               << "const uint rowLength = dstRowsPtrC[rid + 1] - dstRowStart;\n"
+               << "for (uint idx = tid; idx < rowLength; idx += TPB) {\n"
+               << "    dstRowsC[dstRowStart + idx] = srcRowsC[srcRowStart + idx];\n"
+               << "    dstColsC[dstRowStart + idx] = srcColsC[srcRowStart + idx];\n";
+        HAS_VALUES_BEGIN
+        kernel << "    for (uint i = 0; i < BYTE_SIZE; i++)\n"
+               << "        dstValsC[(dstRowStart + idx) * BYTE_SIZE + i] = srcValsC[(srcRowStart + idx) * BYTE_SIZE + i];\n";
+        HAS_VALUES_END
+        kernel << "}\n";
+
+        auto compiledKernel = kernel.compile(ctx);
+        compiledKernel.set_arg(argSrcRowsPtrC, m.GetRowsOffsets().get_buffer());
+        compiledKernel.set_arg(argSrcRowsC, rowsC.get_buffer());
+        compiledKernel.set_arg(argSrcColsC, colsC.get_buffer());
+        compiledKernel.set_arg(argDstRowsPtrC, rowsPtrC.get_buffer());
+        compiledKernel.set_arg(argDstRowsC, dstRowsC.get_buffer());
+        compiledKernel.set_arg(argDstColsC, dstColsC.get_buffer());
+
+        HAS_VALUES_BEGIN
+        compiledKernel.set_arg(argSrcValsC, valsC.get_buffer());
+        compiledKernel.set_arg(argDstValsC, dstValsC.get_buffer());
+        HAS_VALUES_END
+
+        queue.enqueue_1d_range_kernel(compiledKernel, 0, workSize, tpb);
+        PF_SCOPE_MARK(mxm, "compact");
+
+        std::swap(rowsC, dstRowsC);
+        std::swap(colsC, dstColsC);
+        std::swap(valsC, dstValsC);
+    }
+
+    params->w = MatrixCSR::Make(M, N, nvalsC,
+                                std::move(rowsC), std::move(colsC), std::move(valsC),
+                                std::move(rowsPtrC), std::move(rowsLenC))
+                        .As<MatrixBlock>();
+}
+
+spla::Algorithm::Type spla::MxMMaskedCSRCSC::GetType() const {
+    return spla::Algorithm::Type::MxM;
+}
+
+std::string spla::MxMMaskedCSRCSC::GetName() const {
+    return "MxMMaskedCSRCSC";
+}
+
+#undef HAS_VALUES_BEGIN
+#undef HAS_VALUES_END
\ No newline at end of file
diff --git a/sources/algo/mxm/SplaMxMMaskedCSRCSC.hpp b/sources/algo/mxm/SplaMxMMaskedCSRCSC.hpp
new file mode 100644
index 000000000..764000112
--- /dev/null
+++ b/sources/algo/mxm/SplaMxMMaskedCSRCSC.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_SPLAMXMMASKEDCSRCSC_HPP
+#define SPLA_SPLAMXMMASKEDCSRCSC_HPP
+
+#include <algo/SplaAlgorithm.hpp>
+
+namespace spla {
+    class MxMMaskedCSRCSC final : public Algorithm {
+    public:
+        ~MxMMaskedCSRCSC() 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_SPLAMXMMASKEDCSRCSC_HPP
diff --git a/sources/algo/vxm/SplaVxMCOOStructure.cpp b/sources/algo/vxm/SplaVxMCOOStructure.cpp
index a029ae048..1d5233462 100644
--- a/sources/algo/vxm/SplaVxMCOOStructure.cpp
+++ b/sources/algo/vxm/SplaVxMCOOStructure.cpp
@@ -31,8 +31,6 @@
 #include <algo/vxm/SplaVxMCOOStructure.hpp>
 #include <compute/SplaApplyMask.hpp>
 #include <compute/SplaIndicesToRowOffsets.hpp>
-#include <compute/SplaReduceByKey.hpp>
-#include <compute/SplaTransformValues.hpp>
 #include <core/SplaError.hpp>
 #include <core/SplaLibraryPrivate.hpp>
 #include <core/SplaQueueFinisher.hpp>
@@ -91,7 +89,7 @@ void spla::VxMCOOStructure::Process(spla::AlgorithmParams &params) {
 
     compute::fill(buckets.begin(), buckets.end(), 0u, queue);
     BOOST_COMPUTE_CLOSURE(void, countBuckets, (unsigned int rowId), (lengths, buckets), {
-        const length = lengths[rowId];
+        const uint length = lengths[rowId];
         if (length <= 0) return;
         const uint bucketId = (uint) (clamp(ceil(log2((float) length)) - 3.0f, 0.0f, 6.0f));
         atomic_add(&buckets[bucketId], 1u);
@@ -101,11 +99,11 @@ void spla::VxMCOOStructure::Process(spla::AlgorithmParams &params) {
     compute::copy(buckets.begin(), buckets.end(), bucketsHost.begin(), queue);
 
     std::size_t rowsToProcess = std::reduce(bucketsHost.begin(), bucketsHost.end(), 0u);
-    bucketsConfig.resize(rowsToProcess);
+    bucketsConfig.resize(rowsToProcess, queue);
 
     compute::fill(buckets.begin(), buckets.end(), 0u, queue);
     BOOST_COMPUTE_CLOSURE(void, fillBuckets, (unsigned int rowId), (lengths, buckets, bucketsOffsets, bucketsConfig), {
-        const length = lengths[rowId];
+        const uint length = lengths[rowId];
         if (length <= 0) return;
         const uint bucketId = (uint) (clamp(ceil(log2((float) length)) - 3.0f, 0.0f, 6.0f));
         const uint offset = atomic_add(&buckets[bucketId], 1u);
@@ -160,7 +158,7 @@ void spla::VxMCOOStructure::Process(spla::AlgorithmParams &params) {
     };
 
     struct GroupInfo {
-        compute::command_queue queue;
+        boost::compute::command_queue queue;
         std::size_t vpt;
         std::size_t tpb;
     };
diff --git a/sources/expression/prod/SplaMxM.cpp b/sources/expression/prod/SplaMxM.cpp
index a0864e34c..a5e504a43 100644
--- a/sources/expression/prod/SplaMxM.cpp
+++ b/sources/expression/prod/SplaMxM.cpp
@@ -33,7 +33,7 @@
 namespace spla {
     namespace {
         /** Utility to fetch mask block from entry map */
-        inline RefPtr<MatrixBlock> GetMaskBlock(MatrixStorage::EntryMap &map, const MatrixStorage::Index &idx) {
+        inline RefPtr<MatrixBlock> GetBlock(MatrixStorage::EntryMap &map, const MatrixStorage::Index &idx) {
             auto found = map.find(idx);
             return found != map.end() ? found->second : RefPtr<MatrixBlock>{};
         }
@@ -88,6 +88,7 @@ void spla::MxM::Process(std::size_t nodeIdx, const spla::Expression &expression,
     // Clear w, so by default empty result returned
     w->GetStorage()->Clear();
 
+    auto bT = b->GetDecoration(Decorated::Decoration::TransposedMatrix).Cast<Matrix>();
     auto ta = a->GetType();
     auto tb = b->GetType();
     auto tw = w->GetType();
@@ -109,9 +110,13 @@ void spla::MxM::Process(std::size_t nodeIdx, const spla::Expression &expression,
     // Fetch blocks and store locally
     MatrixStorage::EntryMap aBlocks;
     MatrixStorage::EntryMap bBlocks;
+    MatrixStorage::EntryMap bTBlocks;
     MatrixStorage::EntryMap maskBlocks;
     aStorage->GetBlocks(aBlocks);
     bStorage->GetBlocks(bBlocks);
+    // If has b transposed, fetch its blocks
+    if (bT.IsNotNull())
+        bT->GetStorage()->GetBlocks(bTBlocks);
 
     if (hasMask)
         // If mask empty => does not apply mask at all
@@ -162,9 +167,11 @@ void spla::MxM::Process(std::size_t nodeIdx, const spla::Expression &expression,
                 auto deviceId = devicesForProducts[deviceToFetch];
                 auto aIdx = toProcess.a;
                 auto bIdx = toProcess.b;
+                auto bTIdx = SIndex{bIdx.second, bIdx.first};
                 auto aBlock = aBlocks.find(aIdx)->second;
                 auto bBlock = bBlocks.find(bIdx)->second;
-                auto maskBlock = GetMaskBlock(maskBlocks, SIndex{aIdx.first, bIdx.second});
+                auto bTBlock = GetBlock(bTBlocks, bTIdx);
+                auto maskBlock = GetBlock(maskBlocks, SIndex{aIdx.first, bIdx.second});
                 auto task = builder.Emplace("mxm", [=]() {
                     assert(aBlock->GetNcols() == bBlock->GetNrows());
                     ParamsMxM params;
@@ -176,6 +183,7 @@ void spla::MxM::Process(std::size_t nodeIdx, const spla::Expression &expression,
                     params.add = add;
                     params.a = aBlock;
                     params.b = bBlock;
+                    params.bT = bTBlock;
                     params.ta = ta;
                     params.tb = tb;
                     params.tw = tw;
diff --git a/sources/storage/block/SplaMatrixCSR.cpp b/sources/storage/block/SplaMatrixCSR.cpp
index 15b91ddd0..7dbe0fa5d 100644
--- a/sources/storage/block/SplaMatrixCSR.cpp
+++ b/sources/storage/block/SplaMatrixCSR.cpp
@@ -36,6 +36,10 @@ const spla::MatrixCSR::Indices &spla::MatrixCSR::GetRowLengths() const noexcept
     return mRowLengths;
 }
 
+spla::RefPtr<spla::MatrixCSR> spla::MatrixCSR::Make(std::size_t nrows, std::size_t ncols, std::size_t nvals, spla::MatrixCSR::Indices rows, spla::MatrixCSR::Indices cols, spla::MatrixCSR::Values vals, spla::MatrixCSR::Indices rowOffsets, spla::MatrixCSR::Indices rowLengths) {
+    return spla::RefPtr<spla::MatrixCSR>(new MatrixCSR(nrows, ncols, nvals, std::move(rows), std::move(cols), std::move(vals), std::move(rowOffsets), std::move(rowLengths)));
+}
+
 spla::RefPtr<spla::MatrixCSR> spla::MatrixCSR::Make(std::size_t nrows, std::size_t ncols, std::size_t nvals, spla::MatrixCSR::Indices rows, spla::MatrixCSR::Indices cols, spla::MatrixCSR::Values vals, boost::compute::command_queue &queue) {
     return spla::RefPtr<spla::MatrixCSR>(new MatrixCSR(nrows, ncols, nvals, std::move(rows), std::move(cols), std::move(vals), queue));
 }
@@ -54,3 +58,7 @@ spla::MatrixCSR::MatrixCSR(std::size_t nrows, std::size_t ncols, std::size_t nva
     // Compute and cache offsets and rows lengths (frequently used in vxm mxm operations)
     IndicesToRowOffsets(GetRows(), mRowOffsets, mRowLengths, GetNrows(), queue);
 }
+
+spla::MatrixCSR::MatrixCSR(std::size_t nrows, std::size_t ncols, std::size_t nvals, spla::MatrixCSR::Indices rows, spla::MatrixCSR::Indices cols, spla::MatrixCSR::Values vals, spla::MatrixCSR::Indices rowOffsets, spla::MatrixCSR::Indices rowLengths)
+    : MatrixCOO(nrows, ncols, nvals, std::move(rows), std::move(cols), std::move(vals)), mRowOffsets(std::move(rowOffsets)), mRowLengths(std::move(rowLengths)) {
+}
diff --git a/sources/storage/block/SplaMatrixCSR.hpp b/sources/storage/block/SplaMatrixCSR.hpp
index fd1ce6934..0026cbb9e 100644
--- a/sources/storage/block/SplaMatrixCSR.hpp
+++ b/sources/storage/block/SplaMatrixCSR.hpp
@@ -50,11 +50,13 @@ namespace spla {
 
         [[nodiscard]] const Indices &GetRowLengths() const noexcept;
 
+        static RefPtr<MatrixCSR> Make(std::size_t nrows, std::size_t ncols, std::size_t nvals, Indices rows, Indices cols, Values vals, Indices rowOffsets, Indices rowLengths);
         static RefPtr<MatrixCSR> Make(std::size_t nrows, std::size_t ncols, std::size_t nvals, Indices rows, Indices cols, Values vals, boost::compute::command_queue &queue);
         static RefPtr<MatrixCSR> Make(const RefPtr<MatrixCOO> &block, boost::compute::command_queue &queue);
 
     private:
         MatrixCSR(std::size_t nrows, std::size_t ncols, std::size_t nvals, Indices rows, Indices cols, Values vals, boost::compute::command_queue &queue);
+        MatrixCSR(std::size_t nrows, std::size_t ncols, std::size_t nvals, Indices rows, Indices cols, Values vals, Indices rowOffsets, Indices rowLengths);
 
         Indices mRowOffsets;
         Indices mRowLengths;


Test patch:
diff --git a/tests/TestAlgoTc.cpp b/tests/TestAlgoTc.cpp
index 927c1595f..6bcffb235 100644
--- a/tests/TestAlgoTc.cpp
+++ b/tests/TestAlgoTc.cpp
@@ -37,10 +37,15 @@ void testCase(spla::Library &library, std::size_t M, std::size_t nvals, std::siz
 
     auto sp_B = spla::Matrix::Make(M, M, sp_Int32, library);
     auto sp_B_tria = spla::Matrix::Make(M, M, sp_Int32, library);
+    auto sp_B_tria_t = spla::Matrix::Make(M, M, sp_Int32, library);
 
     auto sp_A = spla::Matrix::Make(M, M, sp_Int32, library);
     auto sp_L = spla::Matrix::Make(M, M, sp_Int32, library);
     auto sp_U = spla::Matrix::Make(M, M, sp_Int32, library);
+    auto sp_UT = spla::Matrix::Make(M, M, sp_Int32, library);
+
+    // UT has additional transposed decoration (matrix L)
+    sp_UT->SetDecoration(spla::Decorated::Decoration::TransposedMatrix, sp_L.As<spla::Object>());
 
     auto sp_setup = spla::Expression::Make(library);
     sp_setup->MakeDataWrite(sp_A, A.GetData(library));
@@ -50,9 +55,20 @@ void testCase(spla::Library &library, std::size_t M, std::size_t nvals, std::siz
     auto sp_setup_tria = spla::Expression::Make(library);
     sp_setup_tria->MakeTril(sp_L, sp_A);
     sp_setup_tria->MakeTriu(sp_U, sp_A);
+    sp_setup_tria->MakeTriu(sp_UT, sp_A);
     sp_setup_tria->SubmitWait();
     ASSERT_EQ(sp_setup_tria->GetState(), spla::Expression::State::Evaluated);
 
+    std::int32_t nTrinsCpu = 0;
+    auto host_A = A.ToHostMatrix();
+    auto host_B = spla::RefPtr<spla::HostMatrix>();
+
+    SPLA_TIME_BEGIN(tc_cpu);
+    spla::Tc(nTrinsCpu, host_B, host_A);
+    SPLA_TIME_END(tc_cpu, "cpu");
+
+    auto result = utils::Matrix<int>::FromHostMatrix(host_B);
+
     for (std::size_t k = 0; k < repeats; k++) {
         std::int32_t nTrinsSpla = 0;
         SPLA_TIME_BEGIN(tc_spla);
@@ -64,16 +80,10 @@ void testCase(spla::Library &library, std::size_t M, std::size_t nvals, std::siz
         spla::Tc(nTrinsSplaTria, sp_B_tria, sp_L, sp_U, nullptr);
         SPLA_TIME_END(tc_spla_tria, "spla-tria");
 
-        auto host_A = A.ToHostMatrix();
-        auto host_B = spla::RefPtr<spla::HostMatrix>();
-
-        std::int32_t nTrinsCpu = 0;
-
-        SPLA_TIME_BEGIN(tc_cpu);
-        spla::Tc(nTrinsCpu, host_B, host_A);
-        SPLA_TIME_END(tc_cpu, "cpu");
-
-        auto result = utils::Matrix<int>::FromHostMatrix(host_B);
+        std::int32_t nTrinsSplaTriaT = 0;
+        SPLA_TIME_BEGIN(tc_spla_tria_t);
+        spla::Tc(nTrinsSplaTriaT, sp_B_tria_t, sp_L, sp_UT, nullptr);
+        SPLA_TIME_END(tc_spla_tria_t, "spla-tria-t");
 
         EXPECT_TRUE(result.Equals(sp_B));
         EXPECT_EQ(nTrinsCpu, nTrinsSpla);
@@ -81,6 +91,7 @@ void testCase(spla::Library &library, std::size_t M, std::size_t nvals, std::siz
         // Note: cannot check B, since it has another structure
         // EXPECT_TRUE(result.Equals(sp_B_tria));
         EXPECT_EQ(nTrinsCpu, nTrinsSplaTria);
+        EXPECT_EQ(nTrinsCpu, nTrinsSplaTriaT);
     }
 }
 
diff --git a/tests/TestMxM.cpp b/tests/TestMxM.cpp
index 8adfede9f..9778d23b7 100644
--- a/tests/TestMxM.cpp
+++ b/tests/TestMxM.cpp
@@ -151,61 +151,170 @@ void testNoValues(spla::Library &library, std::size_t M, std::size_t K, std::siz
     ASSERT_TRUE(c.EqualsStructure(spW));
 }
 
+template<typename Type, typename MultOp, typename AddOp, typename Random>
+void testMaskedCSRCSC(spla::Library &library,
+                      std::size_t M, std::size_t K, 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 indexGenerationSeed,
+                      Random valueGenerator,
+                      bool maskComplement) {
+    utils::Matrix a = utils::Matrix<Type>::Generate(M, K, nvals, indexGenerationSeed).SortReduceDuplicates();
+    utils::Matrix b = utils::Matrix<Type>::Generate(K, N, nvals, indexGenerationSeed + 2).SortReduceDuplicates();
+    utils::Matrix mask = utils::Matrix<unsigned char>::Generate(M, N, nvals, indexGenerationSeed + 3).SortReduceDuplicates();
+
+    a.Fill(std::ref(valueGenerator));
+    b.Fill(std::ref(valueGenerator));
+
+    auto spA = spla::Matrix::Make(M, K, spT, library);
+    auto spB = spla::Matrix::Make(K, N, spT, library);
+    auto spBT = spla::Matrix::Make(N, K, spT, library);
+    auto spW = spla::Matrix::Make(M, N, spT, library);
+    auto spMask = spla::Matrix::Make(M, N, spla::Types::Void(library), library);
+
+    // Specify transposed B version
+    spB->SetDecoration(spla::Decorated::Decoration::TransposedMatrix, spBT.As<spla::Object>());
+
+    // 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 if required
+    auto spOpDesc = spla::Descriptor::Make(library);
+    if (maskComplement) {
+        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 spTransposeB = spExpr->MakeTranspose(spBT, nullptr, nullptr, spB);
+    auto spWriteMask = spExpr->MakeDataWrite(spMask, mask.GetDataIndices(library), spDesc);
+    auto spMxM = spExpr->MakeMxM(spW, spMask, spMult, spAdd, spA, spB, spOpDesc);
+    spExpr->Dependency(spWriteA, spMxM);
+    spExpr->Dependency(spWriteB, spTransposeB);
+    spExpr->Dependency(spTransposeB, spMxM);
+    spExpr->Dependency(spWriteMask, spMxM);
+    spExpr->Submit();
+    spExpr->Wait();
+    ASSERT_EQ(spExpr->GetState(), spla::Expression::State::Evaluated);
+
+    utils::Matrix<Type> c = a.template MxM<Type>(mask, maskComplement, b, multOp, addOp);
+    ASSERT_TRUE(c.Equals(spW));
+}
+
+void testNoValuesMaskedCSRCSC(spla::Library &library, std::size_t M, std::size_t K, std::size_t N, std::size_t nvals, std::size_t seed = 0) {
+    utils::Matrix a = utils::Matrix<unsigned char>::Generate(M, K, nvals, seed).SortReduceDuplicates();
+    utils::Matrix b = utils::Matrix<unsigned char>::Generate(K, N, nvals, seed + 1).SortReduceDuplicates();
+    utils::Matrix mask = utils::Matrix<unsigned char>::Generate(M, N, nvals, seed + 2).SortReduceDuplicates();
+
+    auto spT = spla::Types::Void(library);
+    auto spA = spla::Matrix::Make(M, K, spT, library);
+    auto spB = spla::Matrix::Make(K, N, spT, library);
+    auto spBT = spla::Matrix::Make(N, K, spT, library);
+    auto spW = spla::Matrix::Make(M, N, spT, library);
+    auto spMask = spla::Matrix::Make(M, N, spT, library);
+
+    // Specify transposed B version
+    spB->SetDecoration(spla::Decorated::Decoration::TransposedMatrix, spBT.As<spla::Object>());
+
+    // 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 spTransposeB = spExpr->MakeTranspose(spBT, nullptr, nullptr, spB);
+    auto spWriteMask = spExpr->MakeDataWrite(spMask, mask.GetDataIndices(library), spDesc);
+    auto spMxM = spExpr->MakeMxM(spW, spMask, nullptr, nullptr, spA, spB);
+    spExpr->Dependency(spWriteA, spMxM);
+    spExpr->Dependency(spWriteB, spTransposeB);
+    spExpr->Dependency(spTransposeB, spMxM);
+    spExpr->Dependency(spWriteMask, spMxM);
+    spExpr->Submit();
+    spExpr->Wait();
+    ASSERT_EQ(spExpr->GetState(), spla::Expression::State::Evaluated);
+
+    auto dummy = [](unsigned char, unsigned char) { return 0; };
+
+    utils::Matrix<unsigned char> c = a.template MxM<unsigned char>(mask, false, b, dummy, dummy);
+    ASSERT_TRUE(c.EqualsStructure(spW));
+}
+
 void test(std::size_t M, std::size_t K, std::size_t N, std::size_t base, std::size_t step, std::size_t iter, const std::vector<std::size_t> &blocksSizes, utils::UniformIntGenerator<std::int32_t> intGen = utils::UniformIntGenerator<std::int32_t>()) {
     utils::testBlocks(blocksSizes, [=](spla::Library &library) {
+        using Type = float;
         auto spT = spla::Types::Float32(library);
         auto spMult = spla::Functions::MultFloat32(library);
         auto spAdd = spla::Functions::PlusFloat32(library);
-        auto mult = [](float a, float b) { return a * b; };
-        auto add = [](float a, float b) { return a + b; };
+        auto mult = [](Type a, Type b) { return a * b; };
+        auto add = [](Type a, Type b) { return a + b; };
 
         for (std::size_t i = 0; i < iter; i++) {
             std::size_t nvals = base + i * step;
-            testCommon<float>(library, M, K, N, nvals, spT, spMult, spAdd, mult, add, i, utils::UniformRealGenerator<float>(i));
+            testCommon<Type>(library, M, K, N, nvals, spT, spMult, spAdd, mult, add, i, utils::UniformRealGenerator<float>(i));
         }
 
         for (std::size_t i = 0; i < iter; i++) {
             std::size_t nvals = base + i * step;
-            testMasked<float>(library, M, K, N, nvals, spT, spMult, spAdd, mult, add, i, utils::UniformRealGenerator<float>(i), false);
+            testMasked<Type>(library, M, K, N, nvals, spT, spMult, spAdd, mult, add, i, utils::UniformRealGenerator<float>(i), false);
         }
 
         for (std::size_t i = 0; i < iter; i++) {
             std::size_t nvals = base + i * step;
-            testMasked<float>(library, M, K, N, nvals, spT, spMult, spAdd, mult, add, i, utils::UniformRealGenerator<float>(i), true);
+            testMasked<Type>(library, M, K, N, nvals, spT, spMult, spAdd, mult, add, i, utils::UniformRealGenerator<float>(i), true);
         }
 
         for (std::size_t i = 0; i < iter; i++) {
             std::size_t nvals = base + i * step;
-            testNoValues(library, M, K, N, nvals, i * i);
+            testMaskedCSRCSC<Type>(library, M, K, N, nvals, spT, spMult, spAdd, mult, add, i, utils::UniformRealGenerator<float>(i), false);
         }
     });
 
     utils::testBlocks(blocksSizes, [=](spla::Library &library) {
+        using Type = std::int32_t;
         auto spT = spla::Types::Int32(library);
         auto spMult = spla::Functions::MultInt32(library);
         auto spAdd = spla::Functions::PlusInt32(library);
-        auto mult = [](int a, int b) { return a * b; };
-        auto add = [](int a, int b) { return a + b; };
+        auto mult = [](Type a, Type b) { return a * b; };
+        auto add = [](Type a, Type b) { return a + b; };
 
         for (std::size_t i = 0; i < iter; i++) {
             std::size_t nvals = base + i * step;
-            testCommon<std::int32_t>(library, M, K, N, nvals, spT, spMult, spAdd, mult, add, i, intGen);
+            testCommon<Type>(library, M, K, N, nvals, spT, spMult, spAdd, mult, add, i, intGen);
         }
 
         for (std::size_t i = 0; i < iter; i++) {
             std::size_t nvals = base + i * step;
-            testMasked<std::int32_t>(library, M, K, N, nvals, spT, spMult, spAdd, mult, add, i, intGen, false);
+            testMasked<Type>(library, M, K, N, nvals, spT, spMult, spAdd, mult, add, i, intGen, false);
         }
 
         for (std::size_t i = 0; i < iter; i++) {
             std::size_t nvals = base + i * step;
-            testMasked<std::int32_t>(library, M, K, N, nvals, spT, spMult, spAdd, mult, add, i, intGen, true);
+            testMasked<Type>(library, M, K, N, nvals, spT, spMult, spAdd, mult, add, i, intGen, true);
         }
 
+        for (std::size_t i = 0; i < iter; i++) {
+            std::size_t nvals = base + i * step;
+            testMaskedCSRCSC<Type>(library, M, K, N, nvals, spT, spMult, spAdd, mult, add, i, intGen, false);
+        }
+    });
+
+    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, K, N, nvals, i * i);
         }
+
+        for (std::size_t i = 0; i < iter; i++) {
+            std::size_t nvals = base + i * step;
+            testNoValuesMaskedCSRCSC(library, M, K, N, nvals, i * i);
+        }
     });
 }
 
diff --git a/tests/Testing.hpp b/tests/Testing.hpp
index 94e93158a..bfa72e7fd 100644
--- a/tests/Testing.hpp
+++ b/tests/Testing.hpp
@@ -45,7 +45,7 @@
     auto timer = std::chrono::steady_clock::now();
 
 #define SPLA_TIME_END(timer, message)                                           \
-    std::cout << std::setw(10)                                                  \
+    std::cout << std::setw(12)                                                  \
               << message << " "                                                 \
               << std::setw(10)                                                  \
               << std::setprecision(9)                                           \
diff --git a/tests/utils/Matrix.hpp b/tests/utils/Matrix.hpp
index 0cf28d51f..9e5d6689f 100644
--- a/tests/utils/Matrix.hpp
+++ b/tests/utils/Matrix.hpp
@@ -236,6 +236,7 @@ namespace utils {
                 auto b = spVals[i];
 
                 if (!utils::EqWithError(a, b)) {
+                    std::cout << "Values not equal at " << i << " a=" << a << " b=" << b << std::endl;
                     return false;
                 }
             }