Repository: SparseLinearAlgebra/spla (language: python)

Issue:
Triangular matrix support (tc)
Support expression nodes and backend algos for making lower/upper triangular matrix.


Discussion / hints:

Fix patch:
diff --git a/examples/spla_tc.cpp b/examples/spla_tc.cpp
index 00db7f13f..a3118c84b 100644
--- a/examples/spla_tc.cpp
+++ b/examples/spla_tc.cpp
@@ -37,7 +37,6 @@ int main(int argc, const char *const *argv) {
     options.add_option("", cxxopts::Option("niters", "number of iterations to run", cxxopts::value<int>()->default_value("4")));
     options.add_option("", cxxopts::Option("bsize", "size of block to store matrix/vector", cxxopts::value<int>()->default_value("10000000")));
     options.add_option("", cxxopts::Option("devices-count", "amount of devices for execution", cxxopts::value<int>()->default_value("1")));
-    options.add_option("", cxxopts::Option("undirected", "force graph to be undirected", cxxopts::value<bool>()->default_value("false")));
     options.add_option("", cxxopts::Option("verbose", "verbose std output", cxxopts::value<bool>()->default_value("true")));
     options.add_option("", cxxopts::Option("debug-timing", "timing for each iteration of algorithm", cxxopts::value<bool>()->default_value("false")));
     auto args = options.parse(argc, argv);
@@ -51,20 +50,20 @@ int main(int argc, const char *const *argv) {
     int niters;
     int bsize;
     int devicesCount;
-    bool undirected;
     bool removeLoops = true;
     bool ignoreValues = true;
     bool verbose;
+    bool debugTiming;
 
     try {
         mtxpath = args["mtxpath"].as<std::string>();
         niters = args["niters"].as<int>();
         bsize = args["bsize"].as<int>();
         devicesCount = args["devices-count"].as<int>();
-        undirected = args["undirected"].as<bool>();
         verbose = args["verbose"].as<bool>();
+        debugTiming = args["debug-timing"].as<bool>();
     } catch (const std::exception &e) {
-        std::cerr << "Invalid input arguments: " << e.what() << std::endl;
+        std::cerr << "Invalid input arguments: " << e.what();
         return 1;
     }
 
@@ -76,7 +75,7 @@ int main(int argc, const char *const *argv) {
     spla::MatrixLoader<std::int32_t> loader;
 
     try {
-        loader.Load(mtxpath, undirected, removeLoops, ignoreValues, verbose);
+        loader.Load(mtxpath, true, removeLoops, ignoreValues, verbose);
     } catch (const std::exception &e) {
         std::cerr << "Failed load matrix: " << e.what();
         return 1;
@@ -87,36 +86,56 @@ int main(int argc, const char *const *argv) {
     // Fill matrix uniformly with 1
     loader.Fill(1);
 
+    // Configure library
     spla::Library::Config config;
+    config.SetDeviceType(spla::Library::Config::DeviceType::GPU);
     config.SetBlockSize(bsize);
     config.SetWorkersCount(devicesCount * devicesCount);
     config.LimitAmount(devicesCount);
     spla::Library library(config);
 
     // A and B tc args
-    spla::RefPtr<spla::Matrix> A = spla::Matrix::Make(loader.GetNrows(), loader.GetNcols(), spla::Types::Int32(library), library);
-    spla::RefPtr<spla::Matrix> B = spla::Matrix::Make(loader.GetNrows(), loader.GetNcols(), spla::Types::Int32(library), library);
+    auto M = loader.GetNrows(), N = loader.GetNcols();
+    auto T = spla::Types::Int32(library);
+    spla::RefPtr<spla::Matrix> A = spla::Matrix::Make(M, N, T, library);
+    spla::RefPtr<spla::Matrix> L = spla::Matrix::Make(M, N, T, library);
+    spla::RefPtr<spla::Matrix> U = spla::Matrix::Make(M, N, T, library);
 
-    // Prepare data and fill A
+    // Prepare data, fill A, and get L and U matrices
     spla::RefPtr<spla::Expression> prepareData = spla::Expression::Make(library);
-    prepareData->MakeDataWrite(A, spla::DataMatrix::Make(loader.GetRowIndices().data(), loader.GetColIndices().data(), loader.GetValues().data(), loader.GetNvals(), library));
+    auto writeA = prepareData->MakeDataWrite(A, spla::DataMatrix::Make(loader.GetRowIndices().data(), loader.GetColIndices().data(), loader.GetValues().data(), loader.GetNvals(), library));
+    prepareData->Dependency(writeA, prepareData->MakeTril(L, A));
+    prepareData->Dependency(writeA, prepareData->MakeTriu(U, A));
     prepareData->SubmitWait();
     SPLA_ALGO_CHECK(prepareData);
 
+    // Release A to free used memory
+    A.Reset();
+
+    spla::RefPtr<spla::Descriptor> desc = spla::Descriptor::Make(library);
+    desc->SetParam(spla::Descriptor::Param::ProfileTime, debugTiming);
+
     std::int32_t nTrins = 0;
 
     // Warm up phase
     spla::CpuTimer tWarmUp;
-    tWarmUp.Start();
-    spla::Tc(nTrins, B, A, !undirected);
-    tWarmUp.Stop();
+    {
+        spla::RefPtr<spla::Matrix> B = spla::Matrix::Make(M, N, T, library);
+        tWarmUp.Start();
+        spla::Tc(nTrins, B, L, U, desc);
+        tWarmUp.Stop();
+    }
 
     // Main phase, measure iterations
     std::vector<spla::CpuTimer> tIters(niters);
     for (int i = 0; i < niters; i++) {
-        tIters[i].Start();
-        spla::Tc(nTrins, B, A, !undirected);
-        tIters[i].Stop();
+        spla::CpuTimer &tIter = tIters[i];
+        {
+            spla::RefPtr<spla::Matrix> B = spla::Matrix::Make(M, N, T, library);
+            tIter.Start();
+            spla::Tc(nTrins, B, L, U, desc);
+            tIter.Stop();
+        }
     }
 
     spla::OutputMeasurements(tWarmUp, tIters);
diff --git a/include/spla-algo/SplaAlgoTc.hpp b/include/spla-algo/SplaAlgoTc.hpp
index b6a930c51..c60503f12 100644
--- a/include/spla-algo/SplaAlgoTc.hpp
+++ b/include/spla-algo/SplaAlgoTc.hpp
@@ -28,6 +28,7 @@
 #define SPLA_SPLAALGOTC_HPP
 
 #include <spla-algo/SplaAlgoCommon.hpp>
+#include <spla-cpp/SplaDescriptor.hpp>
 #include <spla-cpp/SplaMatrix.hpp>
 #include <spla-cpp/SplaVector.hpp>
 
@@ -41,66 +42,64 @@ namespace spla {
     /**
      * @brief Triangle counting
      *
-     * @note If graph is
-     *  - Directed, then a triangle is a set of edges
-     *      - a -> b
-     *      - a -> c
-     *      - b -> c
-     *   and in the resulting matrix @p spB in the position [a, b]
-     *   will be the number of such triangles
-     *   And in the @p ntrins will be total number of triangles.
+     * 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).
+     * Adjacency matrix must passed as lower and upper triangular parts @p spL and @p spU.
      *
-     *  - Undirected, than a triangle is a set of edges
+     * As an output algorithm produces number of triangles @p ntrins in the graph
+     * as well as buffer matrix @p spB with triangles per vertex.
+     *
+     * @param[out] ntris Number of triangles in the graph
+     * @param[out] spB  Matrix b with number of triangles per vertex
+     * @param[in] spL Lower triangular part of adjacency matrix
+     * @param[in] spU Upper triangular part of adjacency matrix
+     * @param[in] descriptor Algorithm descriptor
+     */
+    SPLA_API void Tc(std::int32_t &ntrins, RefPtr<Matrix> &spB, const RefPtr<Matrix> &spL, const RefPtr<Matrix> &spU, const RefPtr<Descriptor> &descriptor);
+
+    /**
+     * @brief Triangle counting
+     *
+     * As an input expects undirected graph.
+     * Triangle is a set of edges:
      *      - a <-> b
      *      - b <-> c
      *      - c <-> a
-     *   and in the resulting matrix in the position [a, b]
-     *   will be the number of such triangles. But
-     *   total number of triangles @p ntrins will not be the sum of
-     *   elements of @p spB. Instead, it will be the sum divided by 6,
-     *   because each triangle was counted 6 times.
-     *
+     * and in the resulting matrix in the position [a, b]
+     * will be the number of such triangles. But
+     * total number of triangles @p ntrins will not be the sum of
+     * elements of @p spB. Instead, it will be the sum divided by 6,
+     * because each triangle was counted 6 times.
      *
-     * @param [out] ntris - number of triangles in the graph
-     * @param [out] spB - matrix b with number of triangles per vertex
-     * @param [in]  spA - input adjacency matrix of directed graph with values of type int32_t,
-     *                   with 1 stored where has edge between i and j
-     * @param [in]  dir - is input graph directed or not
+     * @param[out] ntris Number of triangles in the graph
+     * @param[out] spB  Matrix b with number of triangles per vertex
+     * @param[in] spA Input adjacency matrix of directed graph with values of type int32_t,
+     *                with 1 stored where has edge between i and j
      */
-    SPLA_API void Tc(std::int32_t &ntrins, RefPtr<Matrix> &spB, const RefPtr<Matrix> &spA, bool dir);
+    SPLA_API void Tc(std::int32_t &ntrins, RefPtr<Matrix> &spB, const RefPtr<Matrix> &spA);
 
     /**
      * @brief Triangle counting
      *
-     * @note If graph is
-     *  - Directed, then a triangle is a set of edges
-     *      - a -> b
-     *      - a -> c
-     *      - b -> c
-     *   and in the resulting matrix @p spB in the position [a, b]
-     *   will be the number of such triangles
-     *   And in the @p ntrins will be total number of triangles.
-     *
-     *  - Undirected, than a triangle is a set of edges
+     * As an input expects undirected graph.
+     * Triangle is a set of edges:
      *      - a <-> b
      *      - b <-> c
      *      - c <-> a
-     *   and in the resulting matrix in the position [a, b]
-     *   will be the number of such triangles. But
-     *   total number of triangles @p ntrins will not be the sum of
-     *   elements of @p spB. Instead, it will be the sum divided by 6,
-     *   because each triangle was counted 6 times.
-     *
+     * and in the resulting matrix in the position [a, b]
+     * will be the number of such triangles. But
+     * total number of triangles @p ntrins will not be the sum of
+     * elements of @p spB. Instead, it will be the sum divided by 6,
+     * because each triangle was counted 6 times.
      *
      * @note Naive cpu reference algo implementation for correctness only
      *
-     * @param [out] ntris - number of triangles in the graph
-     * @param [out] B - matrix b with number of triangles per vertex
-     * @param [in]  A - input adjacency matrix of directed graph with values of type int32_t,
-     *                 with 1 stored where has edge between i and j
-     * @param [in]  dir - is input graph directed or not
+     * @param[out] ntris Number of triangles in the graph
+     * @param[out] B Matrix b with number of triangles per vertex
+     * @param[in] A Input adjacency matrix of directed graph with values of type int32_t,
+     *              with 1 stored where has edge between i and j
      */
-    SPLA_API void Tc(std::int32_t &ntrins, RefPtr<HostMatrix> &B, const RefPtr<HostMatrix> &A, bool dir);
+    SPLA_API void Tc(std::int32_t &ntrins, RefPtr<HostMatrix> &B, const RefPtr<HostMatrix> &A);
 
     /**
      * @}
diff --git a/include/spla-cpp/SplaExpression.hpp b/include/spla-cpp/SplaExpression.hpp
index 9e0e49510..a453802e4 100644
--- a/include/spla-cpp/SplaExpression.hpp
+++ b/include/spla-cpp/SplaExpression.hpp
@@ -424,6 +424,32 @@ namespace spla {
                                              const RefPtr<Matrix> &a,
                                              const RefPtr<Descriptor> &desc = nullptr);
 
+        /**
+         * @brief Make lower-triangular matrix expression node
+         *
+         * @param w Matrix to store result
+         * @param a Matrix source to take lower triangular part
+         * @param desc Operation descriptor
+         *
+         * @return Created expression node
+         */
+        RefPtr<ExpressionNode> MakeTril(const RefPtr<Matrix> &w,
+                                        const RefPtr<Matrix> &a,
+                                        const RefPtr<Descriptor> &desc = nullptr);
+
+        /**
+         * @brief Make upper-triangular matrix expression node
+         *
+         * @param w Matrix to store result
+         * @param a Matrix source to take upper triangular part
+         * @param desc Operation descriptor
+         *
+         * @return Created expression node
+         */
+        RefPtr<ExpressionNode> MakeTriu(const RefPtr<Matrix> &w,
+                                        const RefPtr<Matrix> &a,
+                                        const RefPtr<Descriptor> &desc = nullptr);
+
         /** @return Current expression state */
         State GetState() const;
 
diff --git a/include/spla-cpp/SplaExpressionNode.hpp b/include/spla-cpp/SplaExpressionNode.hpp
index 361a170d6..75489cc83 100644
--- a/include/spla-cpp/SplaExpressionNode.hpp
+++ b/include/spla-cpp/SplaExpressionNode.hpp
@@ -86,7 +86,11 @@ namespace spla {
             /** Matrix-vector multiplication */
             MxV,
             /** Transpose matrix */
-            Transpose
+            Transpose,
+            /** Make lower triangular matrix */
+            Tril,
+            /** Make upper triangular matrix */
+            Triu
         };
 
         /** @return Node argument at specified index */
@@ -170,6 +174,10 @@ namespace spla {
                     return "MxV";
                 case ExpressionNode::Operation::Transpose:
                     return "Transpose";
+                case ExpressionNode::Operation::Tril:
+                    return "Tril";
+                case ExpressionNode::Operation::Triu:
+                    return "Triu";
                 default:
                     return "Unknown";
             }
diff --git a/sources/CMakeLists.txt b/sources/CMakeLists.txt
index dd5eaa661..a8b394d44 100644
--- a/sources/CMakeLists.txt
+++ b/sources/CMakeLists.txt
@@ -25,10 +25,12 @@ set(SPLA_ALGO_SOURCES
 set(SPLA_ALGORITHM_SOURCES
         sources/algo/matrix/SplaMatrixEWiseAddCOO.cpp
         sources/algo/matrix/SplaMatrixEWiseAddCOO.hpp
-        sources/algo/matrix/SplaMatrixTransposeCOO.cpp
-        sources/algo/matrix/SplaMatrixTransposeCOO.hpp
         sources/algo/matrix/SplaMatrixReduceScalarCOO.cpp
         sources/algo/matrix/SplaMatrixReduceScalarCOO.hpp
+        sources/algo/matrix/SplaMatrixTransposeCOO.cpp
+        sources/algo/matrix/SplaMatrixTransposeCOO.hpp
+        sources/algo/matrix/SplaMatrixTriaCOO.cpp
+        sources/algo/matrix/SplaMatrixTriaCOO.hpp
         sources/algo/mxm/SplaMxMCOO.cpp
         sources/algo/mxm/SplaMxMCOO.hpp
         sources/algo/vector/SplaVectorAssignCOO.cpp
@@ -93,10 +95,14 @@ set(SPLA_EXPRESSION_SOURCES
         sources/expression/matrix/SplaMatrixDataWrite.hpp
         sources/expression/matrix/SplaMatrixEWiseAdd.cpp
         sources/expression/matrix/SplaMatrixEWiseAdd.hpp
-        sources/expression/matrix/SplaMatrixTranspose.cpp
-        sources/expression/matrix/SplaMatrixTranspose.hpp
         sources/expression/matrix/SplaMatrixReduceScalar.hpp
         sources/expression/matrix/SplaMatrixReduceScalar.cpp
+        sources/expression/matrix/SplaMatrixTranspose.cpp
+        sources/expression/matrix/SplaMatrixTranspose.hpp
+        sources/expression/matrix/SplaMatrixTril.cpp
+        sources/expression/matrix/SplaMatrixTril.hpp
+        sources/expression/matrix/SplaMatrixTriu.cpp
+        sources/expression/matrix/SplaMatrixTriu.hpp
         sources/expression/prod/SplaMxM.cpp
         sources/expression/prod/SplaMxM.hpp
         sources/expression/prod/SplaVxM.cpp
diff --git a/sources/SplaAlgoTc.cpp b/sources/SplaAlgoTc.cpp
index 29f1cc06b..52f1030fc 100644
--- a/sources/SplaAlgoTc.cpp
+++ b/sources/SplaAlgoTc.cpp
@@ -35,7 +35,54 @@
 #include <unordered_set>
 #include <vector>
 
-void spla::Tc(std::int32_t &ntrins, RefPtr<Matrix> &spB, const RefPtr<Matrix> &spA, bool dir) {
+void spla::Tc(std::int32_t &ntrins, RefPtr<Matrix> &spB, const RefPtr<Matrix> &spL, const RefPtr<Matrix> &spU, const RefPtr<Descriptor> &descriptor) {
+    CHECK_RAISE_ERROR(spB.IsNotNull(), NullPointer, "Passed null argument");
+    CHECK_RAISE_ERROR(spL.IsNotNull(), NullPointer, "Passed null argument");
+    CHECK_RAISE_ERROR(spU.IsNotNull(), NullPointer, "Passed null argument");
+    CHECK_RAISE_ERROR(spL->GetNrows() == spL->GetNcols(), DimensionMismatch, "Matrix must be nxn");
+    CHECK_RAISE_ERROR(spU->GetNrows() == spU->GetNcols(), DimensionMismatch, "Matrix must be nxn");
+    CHECK_RAISE_ERROR(spL->GetDim() == spU->GetDim(), DimensionMismatch, "Must be of the same size");
+
+    bool timing = false;
+
+    if (descriptor.IsNotNull())
+        timing = descriptor->IsParamSet(Descriptor::Param::ProfileTime);
+
+    ntrins = 0;
+
+    auto &library = spL->GetLibrary();
+
+    CHECK_RAISE_ERROR(spL->GetType() == Types::Int32(library), InvalidArgument, "Must be int32 type");
+    CHECK_RAISE_ERROR(spU->GetType() == Types::Int32(library), InvalidArgument, "Must be int32 type");
+
+    auto multOp = Functions::MultInt32(library);
+    auto addOp = Functions::PlusInt32(library);
+
+    auto spNTriangles = Scalar::Make(Types::Int32(library), library);
+    auto spNTrianglesData = DataScalar::Make(&ntrins, library);
+
+    CpuTimer tightTimer;// Tight timer to measure iterations
+    tightTimer.Start();
+
+    auto spTcExpr = Expression::Make(library);
+    auto spNodeMxM = spTcExpr->MakeMxM(spB, spL, multOp, addOp, spL, spU, nullptr);
+    auto spNodeAccum = spTcExpr->MakeReduceScalar(spNTriangles, nullptr, nullptr, addOp, spB, nullptr);
+    auto spNodeReadNTriangles = spTcExpr->MakeDataRead(spNTriangles, spNTrianglesData);
+    spTcExpr->Dependency(spNodeMxM, spNodeAccum);
+    spTcExpr->Dependency(spNodeAccum, spNodeReadNTriangles);
+    spTcExpr->SubmitWait();
+    SPLA_ALGO_CHECK(spTcExpr);
+
+    tightTimer.Stop();
+
+    if (timing) {
+        std::cout << "Result: " << ntrins << " ntrins\n";
+        std::cout << " run tight: " << tightTimer.GetElapsedMs() << "\n";
+    }
+}
+
+void spla::Tc(std::int32_t &ntrins, RefPtr<Matrix> &spB, const RefPtr<Matrix> &spA) {
+    CHECK_RAISE_ERROR(spB.IsNotNull(), NullPointer, "Passed null argument");
     CHECK_RAISE_ERROR(spA.IsNotNull(), NullPointer, "Passed null argument");
     CHECK_RAISE_ERROR(spA->GetNrows() == spA->GetNcols(), DimensionMismatch, "Matrix must be nxn");
 
@@ -43,42 +90,33 @@ void spla::Tc(std::int32_t &ntrins, RefPtr<Matrix> &spB, const RefPtr<Matrix> &s
 
     auto &library = spA->GetLibrary();
     auto type = spA->GetType();
-    auto n = spA->GetNrows();
 
     CHECK_RAISE_ERROR(type->IsBuiltIn(), InvalidArgument, "Must be built-in type");
     CHECK_RAISE_ERROR(type == Types::Int32(library), InvalidArgument, "Must be int32 type");
 
-    auto spTcExpr = Expression::Make(library);
-
     auto multOp = Functions::MultInt32(library);
     auto addOp = Functions::PlusInt32(library);
 
-    auto spATransposed = dir ? Matrix::Make(n, n, type, library) : spA;
+    const auto &spATransposed = spA;
     auto spNTriangles = Scalar::Make(type, library);
     auto spNTrianglesData = spla::DataScalar::Make(library);
-
     spNTrianglesData->SetValue(&ntrins);
 
+    auto spTcExpr = Expression::Make(library);
     auto spNodeMxM = spTcExpr->MakeMxM(spB, spA, multOp, addOp, spA, spATransposed, nullptr);
     auto spNodeAccum = spTcExpr->MakeReduceScalar(spNTriangles, nullptr, nullptr, addOp, spB, nullptr);
     auto spNodeReadNTriangles = spTcExpr->MakeDataRead(spNTriangles, spNTrianglesData);
 
-    if (dir) {
-        auto spNodeTransposeA = spTcExpr->MakeTranspose(spATransposed, nullptr, nullptr, spA);
-        spTcExpr->Dependency(spNodeTransposeA, spNodeMxM);
-    }
     spTcExpr->Dependency(spNodeMxM, spNodeAccum);
     spTcExpr->Dependency(spNodeAccum, spNodeReadNTriangles);
 
     spTcExpr->SubmitWait();
     SPLA_ALGO_CHECK(spTcExpr);
 
-    if (!dir) {
-        ntrins /= 6;
-    }
+    ntrins /= 6;
 }
 
-void spla::Tc(std::int32_t &ntrins, RefPtr<HostMatrix> &B, const RefPtr<HostMatrix> &A, bool dir) {
+void spla::Tc(std::int32_t &ntrins, RefPtr<HostMatrix> &B, const RefPtr<HostMatrix> &A) {
     CHECK_RAISE_ERROR(A.IsNotNull(), NullPointer, "Passed null argument");
     CHECK_RAISE_ERROR(A->GetNrows() == A->GetNcols(), DimensionMismatch, "Matrix must be nxn");
     CHECK_RAISE_ERROR(A->GetElementSize() == sizeof(std::int32_t), InvalidType, "Matrix A must have int32 values of size 4");
@@ -136,13 +174,7 @@ void spla::Tc(std::int32_t &ntrins, RefPtr<HostMatrix> &B, const RefPtr<HostMatr
     std::vector<unsigned char> data(nTrins.size() * sizeof(std::int32_t));
     std::memcpy(data.data(), nTrins.data(), nTrins.size() * sizeof(std::int32_t));
 
-    if (!trinsFrom.empty()) {
-        B = RefPtr<HostMatrix>(new HostMatrix(n, n, std::move(trinsFrom), std::move(trinsTo), std::move(data)));
-    } else {
-        B = nullptr;
-    }
+    B = RefPtr<HostMatrix>(new HostMatrix(n, n, std::move(trinsFrom), std::move(trinsTo), std::move(data)));
 
-    if (!dir) {
-        ntrins /= 6;
-    }
+    ntrins /= 6;
 }
diff --git a/sources/SplaExpression.cpp b/sources/SplaExpression.cpp
index 05ae627b0..c76fd498e 100644
--- a/sources/SplaExpression.cpp
+++ b/sources/SplaExpression.cpp
@@ -307,6 +307,80 @@ spla::Expression::MakeEWiseAdd(const spla::RefPtr<spla::Vector> &w,
                     desc);
 }
 
+spla::RefPtr<spla::ExpressionNode>
+spla::Expression::MakeReduce(const spla::RefPtr<spla::Scalar> &s,
+                             const spla::RefPtr<spla::FunctionBinary> &op,
+                             const spla::RefPtr<spla::Vector> &v,
+                             const spla::RefPtr<spla::Descriptor> &desc) {
+    CHECK_RAISE_ERROR(s.IsNotNull(), NullPointer, "s op can't be null");
+    CHECK_RAISE_ERROR(op.IsNotNull(), NullPointer, "op can't be null");
+    CHECK_RAISE_ERROR(v.IsNotNull(), NullPointer, "v can't be null");
+    CHECK_RAISE_ERROR(op->CanApply(*v, *v, *s), InvalidType, "Can't apply reduce op to provided objects");
+
+    std::vector<RefPtr<Object>> args = {
+            s.As<Object>(),
+            op.As<Object>(),
+            v.As<Object>()};
+
+    return MakeNode(ExpressionNode::Operation::VectorReduce,
+                    std::move(args),
+                    desc);
+}
+
+
+spla::RefPtr<spla::ExpressionNode>
+spla::Expression::MakeReduceScalar(const spla::RefPtr<spla::Scalar> &w,
+                                   const spla::RefPtr<spla::Matrix> &mask,
+                                   const spla::RefPtr<spla::FunctionBinary> &accum,
+                                   const spla::RefPtr<spla::FunctionBinary> &op,
+                                   const spla::RefPtr<spla::Matrix> &a,
+                                   const spla::RefPtr<spla::Descriptor> &desc) {
+    CHECK_RAISE_ERROR(w.IsNotNull(), NullPointer, "scalar w can't be null");
+    CHECK_RAISE_ERROR(op.IsNotNull(), NullPointer, "reduce op can't be null");
+    CHECK_RAISE_ERROR(a.IsNotNull(), NullPointer, "matrix a can't be null");
+    CHECK_RAISE_ERROR(op->CanApply(*a, *a, *w), InvalidType, "Can't apply op to provided objects");
+    if (accum.IsNotNull()) {
+        CHECK_RAISE_ERROR(accum->CanApply(*w, *w, *w), InvalidType, "Can't apply accum to provided objects");
+    }
+    if (mask.IsNotNull()) {
+        CHECK_RAISE_ERROR(mask->GetDim() == a->GetDim(), DimensionMismatch, "Mask has incompatible size");
+    }
+
+    std::vector<RefPtr<Object>> args = {
+            w.As<Object>(),
+            mask.As<Object>(),
+            accum.As<Object>(),
+            op.As<Object>(),
+            a.As<Object>()};
+
+    return MakeNode(ExpressionNode::Operation::MatrixReduceScalar,
+                    std::move(args),
+                    desc);
+}
+
+spla::RefPtr<spla::ExpressionNode>
+spla::Expression::MakeEWiseAdd(const spla::RefPtr<spla::Scalar> &w,
+                               const spla::RefPtr<spla::FunctionBinary> &op,
+                               const spla::RefPtr<spla::Scalar> &a,
+                               const spla::RefPtr<spla::Scalar> &b,
+                               const spla::RefPtr<spla::Descriptor> &desc) {
+    CHECK_RAISE_ERROR(w.IsNotNull(), NullPointer, "Scalar w can't be null");
+    CHECK_RAISE_ERROR(a.IsNotNull(), NullPointer, "Scalar a can't be null");
+    CHECK_RAISE_ERROR(b.IsNotNull(), NullPointer, "Scalar b can't be null");
+    CHECK_RAISE_ERROR(op.IsNotNull(), NullPointer, "Function op can't be null");
+    CHECK_RAISE_ERROR(op->CanApply(*a, *b, *w), InvalidType, "Can't apply op ot provided objects");
+
+    std::vector<RefPtr<Object>> args = {
+            w.As<Object>(),
+            op.As<Object>(),
+            a.As<Object>(),
+            b.As<Object>()};
+
+    return MakeNode(ExpressionNode::Operation::ScalarEWiseAdd,
+                    std::move(args),
+                    desc);
+}
+
 spla::RefPtr<spla::ExpressionNode>
 spla::Expression::MakeMxM(const spla::RefPtr<spla::Matrix> &w,
                           const spla::RefPtr<spla::Matrix> &mask,
@@ -397,85 +471,50 @@ spla::Expression::MakeTranspose(const spla::RefPtr<spla::Matrix> &w,
                     desc);
 }
 
-void spla::Expression::SetState(State state) {
-    mState.store(state);
-}
-
-void spla::Expression::SetFuture(std::unique_ptr<class ExpressionFuture> &&future) {
-    mFuture = std::move(future);
-}
-
-void spla::Expression::SetTasks(std::unique_ptr<class ExpressionTasks> &&tasks) {
-    mTasks = std::move(tasks);
-}
-
-spla::RefPtr<spla::ExpressionNode> spla::Expression::MakeReduce(const spla::RefPtr<spla::Scalar> &s,
-                                                                const spla::RefPtr<spla::FunctionBinary> &op,
-                                                                const spla::RefPtr<spla::Vector> &v,
-                                                                const spla::RefPtr<spla::Descriptor> &desc) {
-    CHECK_RAISE_ERROR(s.IsNotNull(), NullPointer, "s op can't be null");
-    CHECK_RAISE_ERROR(op.IsNotNull(), NullPointer, "op can't be null");
-    CHECK_RAISE_ERROR(v.IsNotNull(), NullPointer, "v can't be null");
-    CHECK_RAISE_ERROR(op->CanApply(*v, *v, *s), InvalidType, "Can't apply reduce op to provided objects");
+spla::RefPtr<spla::ExpressionNode>
+spla::Expression::MakeTril(const spla::RefPtr<spla::Matrix> &w,
+                           const spla::RefPtr<spla::Matrix> &a,
+                           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(w->IsCompatible(*a), InvalidType, "w and a must have the same type");
+    CHECK_RAISE_ERROR(w->GetDim() == a->GetDim(), DimensionMismatch, "Incompatible size");
 
     std::vector<RefPtr<Object>> args = {
-            s.As<Object>(),
-            op.As<Object>(),
-            v.As<Object>()};
+            w.As<Object>(),
+            a.As<Object>()};
 
-    return MakeNode(ExpressionNode::Operation::VectorReduce,
+    return MakeNode(ExpressionNode::Operation::Tril,
                     std::move(args),
                     desc);
 }
 
-
-spla::RefPtr<spla::ExpressionNode> spla::Expression::MakeReduceScalar(const spla::RefPtr<spla::Scalar> &w,
-                                                                      const spla::RefPtr<spla::Matrix> &mask,
-                                                                      const spla::RefPtr<spla::FunctionBinary> &accum,
-                                                                      const spla::RefPtr<spla::FunctionBinary> &op,
-                                                                      const spla::RefPtr<spla::Matrix> &a,
-                                                                      const spla::RefPtr<spla::Descriptor> &desc) {
-    CHECK_RAISE_ERROR(w.IsNotNull(), NullPointer, "scalar w can't be null");
-    CHECK_RAISE_ERROR(op.IsNotNull(), NullPointer, "reduce op can't be null");
-    CHECK_RAISE_ERROR(a.IsNotNull(), NullPointer, "matrix a can't be null");
-    CHECK_RAISE_ERROR(op->CanApply(*a, *a, *w), InvalidType, "Can't apply op to provided objects");
-    if (accum.IsNotNull()) {
-        CHECK_RAISE_ERROR(accum->CanApply(*w, *w, *w), InvalidType, "Can't apply accum to provided objects");
-    }
-    if (mask.IsNotNull()) {
-        CHECK_RAISE_ERROR(mask->GetDim() == a->GetDim(), DimensionMismatch, "Mask has incompatible size");
-    }
+spla::RefPtr<spla::ExpressionNode>
+spla::Expression::MakeTriu(const spla::RefPtr<spla::Matrix> &w,
+                           const spla::RefPtr<spla::Matrix> &a,
+                           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(w->IsCompatible(*a), InvalidType, "w and a must have the same type");
+    CHECK_RAISE_ERROR(w->GetDim() == a->GetDim(), DimensionMismatch, "Incompatible size");
 
     std::vector<RefPtr<Object>> args = {
             w.As<Object>(),
-            mask.As<Object>(),
-            accum.As<Object>(),
-            op.As<Object>(),
             a.As<Object>()};
 
-    return MakeNode(ExpressionNode::Operation::MatrixReduceScalar,
+    return MakeNode(ExpressionNode::Operation::Triu,
                     std::move(args),
                     desc);
 }
 
-spla::RefPtr<spla::ExpressionNode> spla::Expression::MakeEWiseAdd(const spla::RefPtr<spla::Scalar> &w,
-                                                                  const spla::RefPtr<spla::FunctionBinary> &op,
-                                                                  const spla::RefPtr<spla::Scalar> &a,
-                                                                  const spla::RefPtr<spla::Scalar> &b,
-                                                                  const spla::RefPtr<spla::Descriptor> &desc) {
-    CHECK_RAISE_ERROR(w.IsNotNull(), NullPointer, "Scalar w can't be null");
-    CHECK_RAISE_ERROR(a.IsNotNull(), NullPointer, "Scalar a can't be null");
-    CHECK_RAISE_ERROR(b.IsNotNull(), NullPointer, "Scalar b can't be null");
-    CHECK_RAISE_ERROR(op.IsNotNull(), NullPointer, "Function op can't be null");
-    CHECK_RAISE_ERROR(op->CanApply(*a, *b, *w), InvalidType, "Can't apply op ot provided objects");
-
-    std::vector<RefPtr<Object>> args = {
-            w.As<Object>(),
-            op.As<Object>(),
-            a.As<Object>(),
-            b.As<Object>()};
+void spla::Expression::SetState(State state) {
+    mState.store(state);
+}
 
-    return MakeNode(ExpressionNode::Operation::ScalarEWiseAdd,
-                    std::move(args),
-                    desc);
+void spla::Expression::SetFuture(std::unique_ptr<class ExpressionFuture> &&future) {
+    mFuture = std::move(future);
 }
+
+void spla::Expression::SetTasks(std::unique_ptr<class ExpressionTasks> &&tasks) {
+    mTasks = std::move(tasks);
+}
\ No newline at end of file
diff --git a/sources/algo/SplaAlgorithm.hpp b/sources/algo/SplaAlgorithm.hpp
index fd3e5760c..990bdfb57 100644
--- a/sources/algo/SplaAlgorithm.hpp
+++ b/sources/algo/SplaAlgorithm.hpp
@@ -83,7 +83,9 @@ namespace spla {
             /** Vector-matrix multiplication */
             VxM,
             /** Matrix block transpose */
-            Transpose
+            Transpose,
+            /** Matrix block lower triangular matrix */
+            Tria
         };
 
         /**
@@ -143,6 +145,8 @@ namespace spla {
                     return "VxM";
                 case Algorithm::Type::Transpose:
                     return "Transpose";
+                case Algorithm::Type::Tria:
+                    return "Tria";
                 default:
                     return "Unknown";
             }
diff --git a/sources/algo/SplaAlgorithmManager.cpp b/sources/algo/SplaAlgorithmManager.cpp
index d55934f46..7fb5eb19b 100644
--- a/sources/algo/SplaAlgorithmManager.cpp
+++ b/sources/algo/SplaAlgorithmManager.cpp
@@ -33,6 +33,7 @@
 #include <algo/matrix/SplaMatrixEWiseAddCOO.hpp>
 #include <algo/matrix/SplaMatrixReduceScalarCOO.hpp>
 #include <algo/matrix/SplaMatrixTransposeCOO.hpp>
+#include <algo/matrix/SplaMatrixTriaCOO.hpp>
 
 #include <algo/vector/SplaVectorAssignCOO.hpp>
 #include <algo/vector/SplaVectorEWiseAddCOO.hpp>
@@ -49,6 +50,7 @@
 spla::AlgorithmManager::AlgorithmManager(Library &library) : mLibrary(library) {
     Register(new MatrixEWiseAddCOO());
     Register(new MatrixTransposeCOO());
+    Register(new MatrixTriaCOO());
     Register(new MatrixReduceScalarCOO());
     Register(new VectorToDenseCOO());
     Register(new VectorAssignCOO());
diff --git a/sources/algo/SplaAlgorithmParams.hpp b/sources/algo/SplaAlgorithmParams.hpp
index 7099dc163..45cbaf74b 100644
--- a/sources/algo/SplaAlgorithmParams.hpp
+++ b/sources/algo/SplaAlgorithmParams.hpp
@@ -184,6 +184,24 @@ namespace spla {
         RefPtr<Type> type;       // t
     };
 
+    /** Lower triangular matrix block */
+    class ParamsTria final : public AlgorithmParams {
+    public:
+        ~ParamsTria() override = default;
+
+        enum class Mode {
+            Lower,
+            Upper
+        };
+
+        Mode mode = Mode::Lower;// Select desired mode
+        Size firstI;            // Offset row indices of block
+        Size firstJ;            // Offset col indices of block
+        RefPtr<MatrixBlock> w;  // of type t
+        RefPtr<MatrixBlock> a;  // of type t
+        RefPtr<Type> type;      // t
+    };
+
     /** Vector reduce params */
     class ParamsVectorReduce final : public AlgorithmParams {
     public:
diff --git a/sources/algo/matrix/SplaMatrixTriaCOO.cpp b/sources/algo/matrix/SplaMatrixTriaCOO.cpp
new file mode 100644
index 000000000..9a5026b0a
--- /dev/null
+++ b/sources/algo/matrix/SplaMatrixTriaCOO.cpp
@@ -0,0 +1,160 @@
+/**********************************************************************************/
+/* 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/matrix/SplaMatrixTriaCOO.hpp>
+#include <compute/SplaGather.hpp>
+#include <core/SplaLibraryPrivate.hpp>
+#include <core/SplaQueueFinisher.hpp>
+#include <storage/SplaMatrixStorage.hpp>
+#include <storage/block/SplaMatrixCOO.hpp>
+#include <storage/block/SplaMatrixCSR.hpp>
+#include <utils/SplaProfiling.hpp>
+
+bool spla::MatrixTriaCOO::Select(const spla::AlgorithmParams &params) const {
+    auto p = dynamic_cast<const ParamsTria *>(&params);
+
+    return p &&
+           p->a.Is<MatrixCOO>();
+}
+
+void spla::MatrixTriaCOO::Process(spla::AlgorithmParams &params) {
+    using namespace boost;
+
+    PF_SCOPE(tria, "-tria-");
+
+    auto p = dynamic_cast<ParamsTria *>(&params);
+    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, library->GetLogger());
+
+    auto &type = p->type;
+    auto byteSize = type->GetByteSize();
+    auto typeHasValues = type->HasValues();
+
+    auto a = p->a.Cast<MatrixCOO>();
+
+    // Do not process empty block
+    if (a.IsNull())
+        return;
+
+    assert(a->GetNvals() > 0);
+
+    const auto sourceNvals = a->GetNvals();
+    const auto &sourceRows = a->GetRows();
+    const auto &sourceCols = a->GetCols();
+
+    PF_SCOPE_MARK(tria, "setup");
+
+    // Define indices which below main diagonal
+    compute::vector<unsigned int> indicesCount(1, ctx);
+    compute::vector<unsigned int> indices(sourceNvals, ctx);
+    compute::fill_n(indicesCount.begin(), 1, 0u, queue);
+    {
+        compute::detail::meta_kernel kernel("__spla_select_tria");
+        auto argOffsetI = kernel.add_arg<cl_uint>("offsetI");
+        auto argOffsetJ = kernel.add_arg<cl_uint>("offsetJ");
+        auto argCount = kernel.add_arg<cl_uint>("count");
+        auto argRows = kernel.add_arg<const cl_uint *>(compute::memory_object::global_memory, "rows");
+        auto argCols = kernel.add_arg<const cl_uint *>(compute::memory_object::global_memory, "cols");
+        auto argIndices = kernel.add_arg<cl_uint *>(compute::memory_object::global_memory, "indices");
+        auto argIndicesCount = kernel.add_arg<cl_uint *>(compute::memory_object::global_memory, "indicesCount");
+
+        kernel << "const uint i = get_global_id(0);\n"
+               << "if (i < count) {\n"
+               << "    const uint rowId = offsetI + rows[i];\n"
+               << "    const uint colId = offsetJ + cols[i];\n"
+               << "    if (" << (p->mode == ParamsTria::Mode::Lower ? "rowId > colId" : "rowId < colId") << ") {\n"
+               << "        indices[atomic_add(&indicesCount[0], 1u)] = i;\n"
+               << "    }\n"
+               << "}\n";
+
+        auto compiledKernel = kernel.compile(ctx);
+        compiledKernel.set_arg(argOffsetI, static_cast<cl_uint>(p->firstI));
+        compiledKernel.set_arg(argOffsetJ, static_cast<cl_uint>(p->firstJ));
+        compiledKernel.set_arg(argCount, static_cast<cl_uint>(sourceNvals));
+        compiledKernel.set_arg(argRows, sourceRows.get_buffer());
+        compiledKernel.set_arg(argCols, sourceCols.get_buffer());
+        compiledKernel.set_arg(argIndices, indices.get_buffer());
+        compiledKernel.set_arg(argIndicesCount, indicesCount.get_buffer());
+
+        std::size_t tpb = 64;
+        std::size_t workSize = compute::detail::calculate_work_size(sourceNvals, 1, tpb);
+
+        queue.enqueue_1d_range_kernel(compiledKernel, 0, workSize, tpb);
+    }
+
+    PF_SCOPE_MARK(tria, "select tria");
+
+    // Get result number of values
+    auto resultNvals = (indicesCount.begin()).read(queue);
+
+    // No value, nothing to do
+    if (!resultNvals)
+        return;
+
+    PF_SCOPE_MARK(tria, "read nnz");
+
+    // Resize indices and sort before copy
+    indices.resize(resultNvals, queue);
+    compute::sort(indices.begin(), indices.end(), queue);
+
+    PF_SCOPE_MARK(tria, "sort indices");
+
+    compute::vector<unsigned int> resultRows(resultNvals, ctx);
+    compute::vector<unsigned int> resultCols(resultNvals, ctx);
+    compute::vector<unsigned char> resultVals(ctx);
+
+    // Copy indices
+    compute::gather(indices.begin(), indices.end(), sourceRows.begin(), resultRows.begin(), queue);
+    compute::gather(indices.begin(), indices.end(), sourceCols.begin(), resultCols.begin(), queue);
+
+    PF_SCOPE_MARK(tria, "gather indices");
+
+    // Copy values if present
+    if (typeHasValues) {
+        resultVals.resize(resultNvals * byteSize, queue);
+        Gather(indices.begin(), indices.end(), a->GetVals().begin(), resultVals.begin(), byteSize, queue);
+    }
+
+    PF_SCOPE_MARK(tria, "gather values");
+
+    p->w = MatrixCSR::Make(a->GetNrows(), a->GetNcols(), resultNvals,
+                           std::move(resultRows), std::move(resultCols), std::move(resultVals),
+                           queue)
+                   .As<MatrixBlock>();
+}
+
+spla::Algorithm::Type spla::MatrixTriaCOO::GetType() const {
+    return spla::Algorithm::Type::Tria;
+}
+
+std::string spla::MatrixTriaCOO::GetName() const {
+    return "MatrixTriaCOO";
+}
diff --git a/sources/algo/matrix/SplaMatrixTriaCOO.hpp b/sources/algo/matrix/SplaMatrixTriaCOO.hpp
new file mode 100644
index 000000000..a24cca0f7
--- /dev/null
+++ b/sources/algo/matrix/SplaMatrixTriaCOO.hpp
@@ -0,0 +1,44 @@
+/**********************************************************************************/
+/* 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_SPLAMATRIXTRIACOO_HPP
+#define SPLA_SPLAMATRIXTRIACOO_HPP
+
+#include <algo/SplaAlgorithm.hpp>
+
+namespace spla {
+    class MatrixTriaCOO final : public Algorithm {
+    public:
+        ~MatrixTriaCOO() 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_SPLAMATRIXTRIACOO_HPP
diff --git a/sources/algo/mxm/SplaMxMCOO.cpp b/sources/algo/mxm/SplaMxMCOO.cpp
index 555b5e073..16f575e77 100644
--- a/sources/algo/mxm/SplaMxMCOO.cpp
+++ b/sources/algo/mxm/SplaMxMCOO.cpp
@@ -26,7 +26,6 @@
 /**********************************************************************************/
 
 #include <algo/mxm/SplaMxMCOO.hpp>
-#include <boost/compute/algorithm/scatter_if.hpp>
 #include <compute/SplaApplyMask.hpp>
 #include <compute/SplaIndicesToRowOffsets.hpp>
 #include <compute/SplaReduceByKey.hpp>
@@ -37,6 +36,7 @@
 #include <core/SplaLibraryPrivate.hpp>
 #include <core/SplaQueueFinisher.hpp>
 #include <storage/block/SplaMatrixCOO.hpp>
+#include <storage/block/SplaMatrixCSR.hpp>
 #include <utils/SplaProfiling.hpp>
 
 using IndeciesVector = boost::compute::vector<unsigned int>;
@@ -150,10 +150,9 @@ bool spla::MxMCOO::Select(const spla::AlgorithmParams &params) const {
     auto p = dynamic_cast<const ParamsMxM *>(&params);
 
     return p &&
-           p->w.Is<MatrixCOO>() &&
-           p->mask.Is<MatrixCOO>() &&
-           p->a.Is<MatrixCOO>() &&
-           p->b.Is<MatrixCOO>();
+           p->mask.Is<MatrixCSR>() &&
+           p->a.Is<MatrixCSR>() &&
+           p->b.Is<MatrixCSR>();
 }
 
 void spla::MxMCOO::Process(spla::AlgorithmParams &algoParams) {
@@ -168,8 +167,8 @@ void spla::MxMCOO::Process(spla::AlgorithmParams &algoParams) {
     compute::command_queue queue(ctx, device);
     QueueFinisher finisher(queue, library->GetLogger());
 
-    MatrixCOO &a = *params->a.Cast<MatrixCOO>();
-    MatrixCOO &b = *params->b.Cast<MatrixCOO>();
+    MatrixCSR &a = *params->a.Cast<MatrixCSR>();
+    MatrixCSR &b = *params->b.Cast<MatrixCSR>();
 
     const std::size_t wValueByteSize = params->tw->GetByteSize();
     const bool maskIsComplement = params->desc->IsParamSet(Descriptor::Param::MaskComplement);
@@ -185,8 +184,6 @@ void spla::MxMCOO::Process(spla::AlgorithmParams &algoParams) {
         return;
     }
 
-    auto &logger = library->GetLogger();
-
     const auto &typeW = params->tw;
     const bool hasValues = typeW->HasValues();
     const std::size_t valueByteSize = typeW->GetByteSize();
@@ -195,9 +192,8 @@ void spla::MxMCOO::Process(spla::AlgorithmParams &algoParams) {
     PF_SCOPE_MARK(mxm, "setup");
 
     // compute row offsets and row lengths for B
-    compute::vector<unsigned int> bRowOffsets(ctx);
-    compute::vector<unsigned int> bRowLengths(ctx);
-    IndicesToRowOffsets(b.GetRows(), bRowOffsets, bRowLengths, b.GetNrows(), queue);
+    const compute::vector<unsigned int> &bRowOffsets = b.GetRowsOffsets();
+    const compute::vector<unsigned int> &bRowLengths = b.GetRowLengths();
 
     PF_SCOPE_MARK(mxm, "to row offsets");
 
@@ -222,25 +218,27 @@ void spla::MxMCOO::Process(spla::AlgorithmParams &algoParams) {
     std::size_t cooNumNonZeros = (outputPtr.end() - 1).read(queue);
     std::size_t workspaceCapacity = cooNumNonZeros;
     {
-        const auto maxGlobalMem = device.global_memory_size();
-        const auto maxAllocSize = device.get_info<cl_ulong>(CL_DEVICE_MAX_MEM_ALLOC_SIZE);
-
         // See issue #97 for more info https://github.com/JetBrains-Research/spla/issues/97
         // - CL_DEVICE_MAX_MEM_ALLOC_SIZE (max single allocation, nearly max single buffer size, CL_DEVICE_MAX_MEM_ALLOC_SIZE <= CL_DEVICE_GLOBAL_MEM_SIZE)
         // - CL_DEVICE_GLOBAL_MEM_SIZE (total device memory, might be virtualized)
-        const std::size_t factor = std::max<std::size_t>(maxGlobalMem / maxAllocSize, 3);
-        const std::size_t free = maxGlobalMem;
-        const std::size_t maxWorkspaceCapacity = free / (6 * sizeof(unsigned int) + valueByteSize);
-        const std::size_t maxWorkspaceCapacityToSelect = maxWorkspaceCapacity / factor;
 
-        // use at most one third of the remaining capacity
-        workspaceCapacity = std::min(maxWorkspaceCapacityToSelect, workspaceCapacity);
+        // query info from OpenCL
+        const auto maxGlobalMem = device.global_memory_size();
+        const auto maxAllocSize = device.get_info<cl_ulong>(CL_DEVICE_MAX_MEM_ALLOC_SIZE);
+
+        // use at most one fourth of the remaining capacity
+        const std::size_t factor = 4;
+        const std::size_t maxWorkspaceByGlobalMem = (maxGlobalMem / factor) / (6 * sizeof(unsigned int) + valueByteSize);
+        const std::size_t maxWorkspaceByAllocSize = maxAllocSize / std::max(valueByteSize, sizeof(unsigned int));
+        const std::size_t maxWorkspace = std::min(maxWorkspaceByGlobalMem, maxWorkspaceByAllocSize);
+
+        workspaceCapacity = std::min(workspaceCapacity, maxWorkspace);
 
         // Log for info only
-        SPDLOG_LOGGER_TRACE(logger, "Global mem={} KiB alloc={} KiB ({}%) required={} selected={} available={}",
+        SPDLOG_LOGGER_TRACE(library->GetLogger(), "Global mem={} KiB alloc={} KiB ({}%) required={} selected={} available={}",
                             maxGlobalMem / 1024, maxAllocSize / 1024,
                             static_cast<double>(maxAllocSize) / static_cast<double>(maxGlobalMem) * 100.0f,
-                            cooNumNonZeros, workspaceCapacity, maxWorkspaceCapacityToSelect);
+                            cooNumNonZeros, workspaceCapacity, maxWorkspace);
     }
 
     compute::vector<unsigned int> aGatherLocations(ctx), bGatherLocations(ctx);
@@ -267,7 +265,7 @@ void spla::MxMCOO::Process(spla::AlgorithmParams &algoParams) {
                                         aGatherLocations, bGatherLocations,
                                         I, J, V,
                                         params->mult, params->add,
-                                        queue, logger);
+                                        queue, library->GetLogger());
 
         PF_SCOPE_MARK(mxm, "single step");
     } else {
@@ -292,8 +290,7 @@ void spla::MxMCOO::Process(spla::AlgorithmParams &algoParams) {
         std::deque<MatrixSlice> slices;
 
         // compute row offsets for A
-        compute::vector<unsigned int> aRowOffsets(ctx);
-        IndicesToRowOffsets(a.GetRows(), aRowOffsets, a.GetNrows(), queue);
+        const compute::vector<unsigned int> &aRowOffsets = a.GetRowsOffsets();
 
         // compute workspace requirements for each row
         compute::vector<unsigned int> cumulativeRowWorkspace(a.GetNrows(), ctx);
@@ -329,11 +326,16 @@ void spla::MxMCOO::Process(spla::AlgorithmParams &algoParams) {
                                              aGatherLocations, bGatherLocations,
                                              I, J, V,
                                              params->mult, params->add,
-                                             queue, logger);
+                                             queue, library->GetLogger());
             slices.emplace_back(std::move(wSliceRows),
                                 std::move(wSliceCols),
                                 std::move(wSliceVals));
             beginRow = endRow;
+
+            PF_SCOPE_MARK(mxm, "slice");
+
+            SPDLOG_LOGGER_TRACE(library->GetLogger(), "Slice {}/{} nnz={}",
+                                workspaceSize, cooNumNonZeros, wTmpNnz);
         }
 
         PF_SCOPE_MARK(mxm, "slices eval");
diff --git a/sources/algo/vxm/SplaVxMCOO.cpp b/sources/algo/vxm/SplaVxMCOO.cpp
index bfac50fb0..50d6a083e 100644
--- a/sources/algo/vxm/SplaVxMCOO.cpp
+++ b/sources/algo/vxm/SplaVxMCOO.cpp
@@ -45,10 +45,9 @@ 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>();
+           p->b.Is<MatrixCSR>();
 }
 
 void spla::VxMCOO::Process(spla::AlgorithmParams &params) {
@@ -64,7 +63,7 @@ void spla::VxMCOO::Process(spla::AlgorithmParams &params) {
     QueueFinisher finisher(queue, library->GetLogger());
 
     auto a = p->a.Cast<VectorCOO>();
-    auto b = p->b.Cast<MatrixCOO>();
+    auto b = p->b.Cast<MatrixCSR>();
     auto mask = p->mask.Cast<VectorCOO>();
     auto complementMask = desc->IsParamSet(Descriptor::Param::MaskComplement);
 
@@ -78,33 +77,17 @@ void spla::VxMCOO::Process(spla::AlgorithmParams &params) {
     const auto &tb = p->tb;
     const auto &tw = p->tw;
     auto hasValues = tw->HasValues();
-    auto M = a->GetNrows();
     auto N = b->GetNcols();
 
-    // todo: beautify and avoid pinters usage?
-    const compute::vector<Index> *offsets;
-    const compute::vector<Index> *lengths;
-
-    compute::vector<Index> offsetsBuffer(ctx);
-    compute::vector<Index> lengthsBuffer(ctx);
-
-    if (b.Is<MatrixCSR>()) {
-        // Get rows offsets and rows lengths for matrix b
-        auto bCSR = b.Cast<MatrixCSR>();
-        offsets = &bCSR->GetRowsOffsets();
-        lengths = &bCSR->GetRowLengths();
-    } else {
-        // Compute new rows offsets and rows lengths for matrix b
-        offsets = &offsetsBuffer;
-        lengths = &lengthsBuffer;
-        IndicesToRowOffsets(b->GetRows(), offsetsBuffer, lengthsBuffer, M, queue);
-    }
+    // Get row lengths and offsets for B
+    const compute::vector<Index> &offsets = b->GetRowsOffsets();
+    const compute::vector<Index> &lengths = b->GetRowLengths();
 
     PF_SCOPE(vxm, "-vxm-");
 
     // 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);
+    compute::gather(a->GetRows().begin(), a->GetRows().end(), lengths.begin(), segmentLengths.begin(), queue);
 
     PF_SCOPE_MARK(vxm, "segment lengths");
 
@@ -138,11 +121,10 @@ void spla::VxMCOO::Process(spla::AlgorithmParams &params) {
     PF_SCOPE_MARK(vxm, "a-loc scan");
 
     auto &aRows = a->GetRows();
-    auto &offsetsRef = *offsets;
-    BOOST_COMPUTE_CLOSURE(unsigned int, unfoldSegment, (unsigned int i), (outputPtr, offsetsRef, aRows, aLocations), {
+    BOOST_COMPUTE_CLOSURE(unsigned int, unfoldSegment, (unsigned int i), (outputPtr, offsets, aRows, aLocations), {
         uint locationOfRowIndex = aLocations[i];
         uint rowIdx = aRows[locationOfRowIndex];
-        uint rowBaseOffset = offsetsRef[rowIdx];
+        uint rowBaseOffset = offsets[rowIdx];
         uint offsetOfRowSegment = outputPtr[locationOfRowIndex];
         return rowBaseOffset + (i - offsetOfRowSegment);
     });
diff --git a/sources/algo/vxm/SplaVxMCOOStructure.cpp b/sources/algo/vxm/SplaVxMCOOStructure.cpp
index 77a9cb670..a029ae048 100644
--- a/sources/algo/vxm/SplaVxMCOOStructure.cpp
+++ b/sources/algo/vxm/SplaVxMCOOStructure.cpp
@@ -48,9 +48,8 @@ bool spla::VxMCOOStructure::Select(const spla::AlgorithmParams &params) const {
 
     return p &&
            !p->tw->HasValues() &&
-           p->w.Is<VectorCOO>() &&
            p->a.Is<VectorCOO>() &&
-           p->b.Is<MatrixCOO>() &&
+           p->b.Is<MatrixCSR>() &&
            (p->mask.Is<VectorCOO>() || p->mask.Is<VectorDense>());
 }
 
@@ -67,7 +66,7 @@ void spla::VxMCOOStructure::Process(spla::AlgorithmParams &params) {
     QueueFinisher finisher(queue, library->GetLogger());
 
     auto a = p->a.Cast<VectorCOO>();
-    auto b = p->b.Cast<MatrixCOO>();
+    auto b = p->b.Cast<MatrixCSR>();
     auto complementMask = desc->IsParamSet(Descriptor::Param::MaskComplement);
 
     if (p->hasMask && !complementMask && p->mask.IsNull())
@@ -76,29 +75,11 @@ void spla::VxMCOOStructure::Process(spla::AlgorithmParams &params) {
     if (a->GetNvals() == 0 || b->GetNvals() == 0)
         return;
 
-    auto M = a->GetNrows();
     auto N = b->GetNcols();
 
-    assert(!p->tw->HasValues());
-
     // Get row lengths and offsets for B
-    const compute::vector<Index> *offsetsB;
-    const compute::vector<Index> *lengthsB;
-
-    compute::vector<Index> offsetsBuffer(ctx);
-    compute::vector<Index> lengthsBuffer(ctx);
-
-    if (b.Is<MatrixCSR>()) {
-        // Get rows offsets and rows lengths for matrix b
-        auto bCSR = b.Cast<MatrixCSR>();
-        offsetsB = &bCSR->GetRowsOffsets();
-        lengthsB = &bCSR->GetRowLengths();
-    } else {
-        // Compute new rows offsets and rows lengths for matrix b
-        offsetsB = &offsetsBuffer;
-        lengthsB = &lengthsBuffer;
-        IndicesToRowOffsets(b->GetRows(), offsetsBuffer, lengthsBuffer, M, queue);
-    }
+    const compute::vector<Index> &offsets = b->GetRowsOffsets();
+    const compute::vector<Index> &lengths = b->GetRowLengths();
 
     const std::size_t BUCKETS_COUNT = 7;
     compute::vector<Index> buckets(BUCKETS_COUNT, ctx);
@@ -108,10 +89,9 @@ void spla::VxMCOOStructure::Process(spla::AlgorithmParams &params) {
 
     PF_SCOPE(vxm, "-vxm-");
 
-    auto &bLengths = *lengthsB;
     compute::fill(buckets.begin(), buckets.end(), 0u, queue);
-    BOOST_COMPUTE_CLOSURE(void, countBuckets, (unsigned int rowId), (bLengths, buckets), {
-        const length = bLengths[rowId];
+    BOOST_COMPUTE_CLOSURE(void, countBuckets, (unsigned int rowId), (lengths, buckets), {
+        const 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);
@@ -124,8 +104,8 @@ void spla::VxMCOOStructure::Process(spla::AlgorithmParams &params) {
     bucketsConfig.resize(rowsToProcess);
 
     compute::fill(buckets.begin(), buckets.end(), 0u, queue);
-    BOOST_COMPUTE_CLOSURE(void, fillBuckets, (unsigned int rowId), (bLengths, buckets, bucketsOffsets, bucketsConfig), {
-        const length = bLengths[rowId];
+    BOOST_COMPUTE_CLOSURE(void, fillBuckets, (unsigned int rowId), (lengths, buckets, bucketsOffsets, bucketsConfig), {
+        const 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);
@@ -159,7 +139,7 @@ void spla::VxMCOOStructure::Process(spla::AlgorithmParams &params) {
     compute::fill_n(resultStructure.begin(), N, 0u, queue);
 
     auto kernelCompiled = kernel.compile(ctx);
-    kernelCompiled.set_arg(argRowOffsets, offsetsB->get_buffer());
+    kernelCompiled.set_arg(argRowOffsets, offsets.get_buffer());
     kernelCompiled.set_arg(argColIndices, b->GetCols().get_buffer());
     kernelCompiled.set_arg(argResult, resultStructure.get_buffer());
     kernelCompiled.set_arg(argBucketOffsets, bucketsOffsets.get_buffer());
diff --git a/sources/expression/SplaExpressionManager.cpp b/sources/expression/SplaExpressionManager.cpp
index 3300c03f9..c8a5ee57b 100644
--- a/sources/expression/SplaExpressionManager.cpp
+++ b/sources/expression/SplaExpressionManager.cpp
@@ -25,21 +25,29 @@
 /* SOFTWARE.                                                                      */
 /**********************************************************************************/
 
+#include <algorithm>
+#include <vector>
+
 #include <core/SplaError.hpp>
 #include <core/SplaLibraryPrivate.hpp>
-
 #include <expression/SplaExpressionFuture.hpp>
 #include <expression/SplaExpressionTasks.hpp>
+
 #include <expression/matrix/SplaMatrixDataRead.hpp>
 #include <expression/matrix/SplaMatrixDataWrite.hpp>
 #include <expression/matrix/SplaMatrixEWiseAdd.hpp>
 #include <expression/matrix/SplaMatrixReduceScalar.hpp>
 #include <expression/matrix/SplaMatrixTranspose.hpp>
+#include <expression/matrix/SplaMatrixTril.hpp>
+#include <expression/matrix/SplaMatrixTriu.hpp>
+
 #include <expression/prod/SplaMxM.hpp>
 #include <expression/prod/SplaVxM.hpp>
+
 #include <expression/scalar/SplaScalarDataRead.hpp>
 #include <expression/scalar/SplaScalarDataWrite.hpp>
 #include <expression/scalar/SplaScalarEWiseAdd.hpp>
+
 #include <expression/vector/SplaVectorAssign.hpp>
 #include <expression/vector/SplaVectorDataRead.hpp>
 #include <expression/vector/SplaVectorDataWrite.hpp>
@@ -47,9 +55,6 @@
 #include <expression/vector/SplaVectorReduce.hpp>
 #include <expression/vector/SplaVectorToDense.hpp>
 
-#include <algorithm>
-#include <vector>
-
 spla::ExpressionManager::ExpressionManager(spla::Library &library) : mLibrary(library) {
     // Here we can register built-in processors, one by on
     // NOTE: order of registration matters. First one has priority among others for the same op.
@@ -57,6 +62,8 @@ spla::ExpressionManager::ExpressionManager(spla::Library &library) : mLibrary(li
     Register(new MatrixDataWrite());
     Register(new MatrixEWiseAdd());
     Register(new MatrixTranspose());
+    Register(new MatrixTril());
+    Register(new MatrixTriu());
     Register(new MatrixReduceScalar());
     Register(new ScalarDataRead());
     Register(new ScalarDataWrite());
diff --git a/sources/expression/matrix/SplaMatrixDataRead.cpp b/sources/expression/matrix/SplaMatrixDataRead.cpp
index 092727b85..896af6d30 100644
--- a/sources/expression/matrix/SplaMatrixDataRead.cpp
+++ b/sources/expression/matrix/SplaMatrixDataRead.cpp
@@ -72,7 +72,7 @@ void spla::MatrixDataRead::Process(std::size_t nodeIdx, const spla::Expression &
     storage->GetBlocks(shared->entries);
 
     for (const auto &entry : shared->entries) {
-        CHECK_RAISE_ERROR(entry.second->GetFormat() == MatrixBlock::Format::COO, NotImplemented,
+        CHECK_RAISE_ERROR(entry.second.Cast<MatrixCOO>(), NotImplemented,
                           "Supported only COO matrix block format");
     }
 
diff --git a/sources/expression/matrix/SplaMatrixTril.cpp b/sources/expression/matrix/SplaMatrixTril.cpp
new file mode 100644
index 000000000..39dedfcbc
--- /dev/null
+++ b/sources/expression/matrix/SplaMatrixTril.cpp
@@ -0,0 +1,84 @@
+/**********************************************************************************/
+/* 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/matrix/SplaMatrixTril.hpp>
+#include <storage/SplaMatrixStorage.hpp>
+#include <storage/SplaVectorStorage.hpp>
+
+bool spla::MatrixTril::Select(std::size_t, const spla::Expression &) {
+    return true;
+}
+
+void spla::MatrixTril::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 w = node->GetArg(0).Cast<Matrix>();
+    auto a = node->GetArg(1).Cast<Matrix>();
+    auto desc = node->GetDescriptor();
+
+    assert(w.IsNotNull());
+    assert(a.IsNotNull());
+    assert(desc.IsNotNull());
+
+    w->GetStorage()->Clear();
+    auto blockSize = a->GetStorage()->GetBlockSize();
+
+    std::size_t requiredDeviceCount = w->GetStorage()->GetNblockRows() * w->GetStorage()->GetNblockCols();
+    auto deviceIds = library->GetDeviceManager().FetchDevices(requiredDeviceCount, node);
+
+    for (std::size_t i = 0; i < w->GetStorage()->GetNblockRows(); i++) {
+        for (std::size_t j = 0; j < w->GetStorage()->GetNblockCols(); j++) {
+            auto deviceId = deviceIds[i * w->GetStorage()->GetNblockCols() + j];
+            builder.Emplace("tril", [=]() {
+                MatrixStorage::Index index{static_cast<unsigned int>(i), static_cast<unsigned int>(j)};
+
+                ParamsTria params;
+                params.desc = desc;
+                params.deviceId = deviceId;
+                params.firstI = blockSize * i;
+                params.firstJ = blockSize * j;
+                params.a = a->GetStorage()->GetBlock(index);
+                params.type = w->GetType();
+                params.mode = ParamsTria::Mode::Lower;
+                library->GetAlgoManager()->Dispatch(Algorithm::Type::Tria, params);
+
+                if (params.w.IsNotNull()) {
+                    w->GetStorage()->SetBlock(index, params.w);
+                    SPDLOG_LOGGER_TRACE(library->GetLogger(), "Tria block=({},{})", index.first, index.second);
+                }
+            });
+        }
+    }
+}
+
+spla::ExpressionNode::Operation spla::MatrixTril::GetOperationType() const {
+    return spla::ExpressionNode::Operation::Tril;
+}
diff --git a/sources/expression/matrix/SplaMatrixTril.hpp b/sources/expression/matrix/SplaMatrixTril.hpp
new file mode 100644
index 000000000..eb31970f3
--- /dev/null
+++ b/sources/expression/matrix/SplaMatrixTril.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_SPLAMATRIXTRIL_HPP
+#define SPLA_SPLAMATRIXTRIL_HPP
+
+#include <expression/SplaNodeProcessor.hpp>
+
+namespace spla {
+
+    class MatrixTril final : public NodeProcessor {
+    public:
+        ~MatrixTril() 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_SPLAMATRIXTRIL_HPP
diff --git a/sources/expression/matrix/SplaMatrixTriu.cpp b/sources/expression/matrix/SplaMatrixTriu.cpp
new file mode 100644
index 000000000..19505446c
--- /dev/null
+++ b/sources/expression/matrix/SplaMatrixTriu.cpp
@@ -0,0 +1,84 @@
+/**********************************************************************************/
+/* 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/matrix/SplaMatrixTriu.hpp>
+#include <storage/SplaMatrixStorage.hpp>
+#include <storage/SplaVectorStorage.hpp>
+
+bool spla::MatrixTriu::Select(std::size_t, const spla::Expression &) {
+    return true;
+}
+
+void spla::MatrixTriu::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 w = node->GetArg(0).Cast<Matrix>();
+    auto a = node->GetArg(1).Cast<Matrix>();
+    auto desc = node->GetDescriptor();
+
+    assert(w.IsNotNull());
+    assert(a.IsNotNull());
+    assert(desc.IsNotNull());
+
+    w->GetStorage()->Clear();
+    auto blockSize = a->GetStorage()->GetBlockSize();
+
+    std::size_t requiredDeviceCount = w->GetStorage()->GetNblockRows() * w->GetStorage()->GetNblockCols();
+    auto deviceIds = library->GetDeviceManager().FetchDevices(requiredDeviceCount, node);
+
+    for (std::size_t i = 0; i < w->GetStorage()->GetNblockRows(); i++) {
+        for (std::size_t j = 0; j < w->GetStorage()->GetNblockCols(); j++) {
+            auto deviceId = deviceIds[i * w->GetStorage()->GetNblockCols() + j];
+            builder.Emplace("triu", [=]() {
+                MatrixStorage::Index index{static_cast<unsigned int>(i), static_cast<unsigned int>(j)};
+
+                ParamsTria params;
+                params.desc = desc;
+                params.deviceId = deviceId;
+                params.firstI = blockSize * i;
+                params.firstJ = blockSize * j;
+                params.a = a->GetStorage()->GetBlock(index);
+                params.type = w->GetType();
+                params.mode = ParamsTria::Mode::Upper;
+                library->GetAlgoManager()->Dispatch(Algorithm::Type::Tria, params);
+
+                if (params.w.IsNotNull()) {
+                    w->GetStorage()->SetBlock(index, params.w);
+                    SPDLOG_LOGGER_TRACE(library->GetLogger(), "Tria block=({},{})", index.first, index.second);
+                }
+            });
+        }
+    }
+}
+
+spla::ExpressionNode::Operation spla::MatrixTriu::GetOperationType() const {
+    return spla::ExpressionNode::Operation::Triu;
+}
diff --git a/sources/expression/matrix/SplaMatrixTriu.hpp b/sources/expression/matrix/SplaMatrixTriu.hpp
new file mode 100644
index 000000000..d75780c07
--- /dev/null
+++ b/sources/expression/matrix/SplaMatrixTriu.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_SPLAMATRIXTRIU_HPP
+#define SPLA_SPLAMATRIXTRIU_HPP
+
+#include <expression/SplaNodeProcessor.hpp>
+
+namespace spla {
+
+    class MatrixTriu final : public NodeProcessor {
+    public:
+        ~MatrixTriu() 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_SPLAMATRIXTRIU_HPP
diff --git a/sources/storage/SplaMatrixStorage.cpp b/sources/storage/SplaMatrixStorage.cpp
index a2cccacd3..4c462b3c4 100644
--- a/sources/storage/SplaMatrixStorage.cpp
+++ b/sources/storage/SplaMatrixStorage.cpp
@@ -117,6 +117,10 @@ std::size_t spla::MatrixStorage::GetNblockCols() const noexcept {
     return mNblockCols;
 }
 
+std::size_t spla::MatrixStorage::GetBlockSize() const noexcept {
+    return mBlockSize;
+}
+
 void spla::MatrixStorage::Dump(std::ostream &stream) const {
     std::lock_guard<std::mutex> lock(mMutex);
 
diff --git a/sources/storage/SplaMatrixStorage.hpp b/sources/storage/SplaMatrixStorage.hpp
index f1b26ef4c..848c9a045 100644
--- a/sources/storage/SplaMatrixStorage.hpp
+++ b/sources/storage/SplaMatrixStorage.hpp
@@ -125,6 +125,9 @@ namespace spla {
         /** @return Number of cols of blocks */
         [[nodiscard]] std::size_t GetNblockCols() const noexcept;
 
+        /** @return Size of single block */
+        [[nodiscard]] std::size_t GetBlockSize() const noexcept;
+
         /** Dump matrix content to provided stream */
         void Dump(std::ostream &stream) const;
 


Test patch:
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 523f7466b..52f855679 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -26,6 +26,8 @@ spla_test_target(TestReduceDuplicates)
 spla_test_target(TestScalarEWiseAdd)
 spla_test_target(TestStableSortByColumn)
 spla_test_target(TestTranspose)
+spla_test_target(TestTril)
+spla_test_target(TestTriu)
 spla_test_target(TestVectorAssign)
 spla_test_target(TestVectorEWiseAdd)
 spla_test_target(TestVectorReduce)
diff --git a/tests/TestAlgoTc.cpp b/tests/TestAlgoTc.cpp
index e2959a7cd..927c1595f 100644
--- a/tests/TestAlgoTc.cpp
+++ b/tests/TestAlgoTc.cpp
@@ -28,43 +28,59 @@
 #include <Testing.hpp>
 #include <spla-algo/SplaAlgoTc.hpp>
 
-void testCase(spla::Library &library, std::size_t M, std::size_t nvals, bool dir, std::size_t seed = 0, std::size_t repeats = 4) {
-    utils::Matrix A = utils::Matrix<int>::Generate(M, M, nvals, seed).SortReduceDuplicates();
+void testCase(spla::Library &library, std::size_t M, std::size_t nvals, std::size_t seed = 0, std::size_t repeats = 4) {
+    utils::Matrix A = utils::Matrix<int>::Generate(M, M, nvals, seed).SortReduceDuplicates().RemoveLoops();
     A.Fill([]() { return 1; });
-
-    if (!dir) {
-        A = A.EWiseAdd(A.Transpose(), [](int a, int b) { return a | b; });
-    }
+    A = A.EWiseAdd(A.Transpose(), [](int a, int b) { return a | b; });
 
     auto sp_Int32 = spla::Types::Int32(library);
 
     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_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_setup = spla::Expression::Make(library);
     sp_setup->MakeDataWrite(sp_A, A.GetData(library));
     sp_setup->SubmitWait();
     ASSERT_EQ(sp_setup->GetState(), spla::Expression::State::Evaluated);
 
+    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->SubmitWait();
+    ASSERT_EQ(sp_setup_tria->GetState(), spla::Expression::State::Evaluated);
+
     for (std::size_t k = 0; k < repeats; k++) {
         std::int32_t nTrinsSpla = 0;
-
         SPLA_TIME_BEGIN(tc_spla);
-        spla::Tc(nTrinsSpla, sp_B, sp_A, dir);
+        spla::Tc(nTrinsSpla, sp_B, sp_A);
         SPLA_TIME_END(tc_spla, "spla");
 
+        std::int32_t nTrinsSplaTria = 0;
+        SPLA_TIME_BEGIN(tc_spla_tria);
+        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, dir);
+        spla::Tc(nTrinsCpu, host_B, host_A);
         SPLA_TIME_END(tc_cpu, "cpu");
 
         auto result = utils::Matrix<int>::FromHostMatrix(host_B);
+
         EXPECT_TRUE(result.Equals(sp_B));
-        EXPECT_EQ(nTrinsSpla, nTrinsCpu);
+        EXPECT_EQ(nTrinsCpu, nTrinsSpla);
+
+        // Note: cannot check B, since it has another structure
+        // EXPECT_TRUE(result.Equals(sp_B_tria));
+        EXPECT_EQ(nTrinsCpu, nTrinsSplaTria);
     }
 }
 
@@ -72,8 +88,7 @@ void test(std::size_t M, std::size_t base, std::size_t step, std::size_t iter, c
     utils::testBlocks(blocksSizes, "", 1, [=](spla::Library &library) {
         for (std::size_t i = 0; i < iter; i++) {
             std::size_t nvals = base + i * step;
-            testCase(library, M, nvals, true, i);
-            testCase(library, M, nvals, false, i);
+            testCase(library, M, nvals, i, repeats);
         }
     });
 }
@@ -81,11 +96,13 @@ void test(std::size_t M, std::size_t base, std::size_t step, std::size_t iter, c
 void testCpu(std::size_t n,
              std::vector<spla::Index> from,
              std::vector<spla::Index> to,
-             bool dir,
              std::int32_t nTrins,
              const std::vector<std::int32_t> &trinsPerEdge) {
     std::int32_t actualNTrins = 0;
 
+    assert(to.size() == from.size());
+    assert(from.size() == trinsPerEdge.size());
+
     const std::size_t nVals = from.size();
     assert(to.size() == nVals);
 
@@ -95,13 +112,9 @@ void testCpu(std::size_t n,
     std::memcpy(values.data(), valuesData.data(), nVals * sizeof(std::int32_t));
 
     auto trinsActual = spla::RefPtr<spla::HostMatrix>();
-    auto A = spla::RefPtr<spla::HostMatrix>(new spla::HostMatrix(n,
-                                                                 n,
-                                                                 from,
-                                                                 to,
-                                                                 valuesData));
+    auto A = spla::RefPtr<spla::HostMatrix>(new spla::HostMatrix(n, n, from, to, valuesData));
 
-    spla::Tc(actualNTrins, trinsActual, A, dir);
+    spla::Tc(actualNTrins, trinsActual, A);
 
     if (trinsActual.IsNull()) {
         EXPECT_EQ(nTrins, 0);
@@ -126,23 +139,12 @@ void testCpu(std::size_t n,
     }
 }
 
-TEST(TC, CpuShallowDirectedDir) {
-    testCpu(3, {1, 1, 2}, {0, 2, 0}, true, 1, {0, 1, 0});
-    testCpu(3, {0, 1, 2}, {1, 2, 0}, true, 0, {0, 0, 0});
-    testCpu(5,
-            {2, 2, 2, 3, 3, 3, 4},
-            {0, 1, 3, 0, 1, 4, 1}, true, 3,
-            {0, 0, 2, 0, 0, 1, 0});
-}
-
-// todo: fix CpuShallowDirectedNonDir fails on third testCpu-call case
-TEST(TC, CpuShallowDirectedNonDir) {
-    testCpu(3, {0, 0, 1, 1, 2, 2}, {1, 2, 0, 2, 0, 1}, false, 1, {1, 1, 1});
-    testCpu(3, {0, 1, 1, 2}, {1, 0, 2, 1}, false, 0, {0, 0, 0});
+TEST(TC, CpuShallow) {
+    testCpu(3, {0, 0, 1, 1, 2, 2}, {1, 2, 0, 2, 0, 1}, 1, {1, 1, 1, 1, 1, 1});
+    testCpu(3, {0, 1, 1, 2}, {1, 0, 2, 1}, 0, {0, 0, 0, 0});
     testCpu(6,
             {0, 0, 1, 1, 1, 2, 2, 3, 3, 3, 4, 4, 5, 5},
             {1, 2, 0, 2, 3, 0, 1, 1, 4, 5, 3, 5, 3, 4},
-            false,
             2,
             {1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1});
 }
@@ -150,13 +152,13 @@ TEST(TC, CpuShallowDirectedNonDir) {
 TEST(TC, Small) {
     std::vector<std::size_t> blockSizes = {100, 1000};
     std::size_t M = 120;
-    test(M, M, M, 10, blockSizes);
+    test(M, M, M, 10, blockSizes, 4);
 }
 
 TEST(TC, Medium) {
     std::vector<std::size_t> blockSizes = {1000, 10000};
     std::size_t M = 1220;
-    test(M, M, M, 10, blockSizes, 3);
+    test(M, M, M, 10, blockSizes, 4);
 }
 
 TEST(TC, Large) {
@@ -168,7 +170,7 @@ TEST(TC, Large) {
 TEST(TC, MegaLarge) {
     std::vector<std::size_t> blockSizes = {302400};
     std::size_t M = 102400;
-    test(M, M, M, 5, blockSizes, 3);
+    test(M, M, M, 3, blockSizes, 3);
 }
 
 SPLA_GTEST_MAIN
\ No newline at end of file
diff --git a/tests/TestTril.cpp b/tests/TestTril.cpp
new file mode 100644
index 000000000..f83215192
--- /dev/null
+++ b/tests/TestTril.cpp
@@ -0,0 +1,138 @@
+/**********************************************************************************/
+/* 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>
+void testCommon(spla::Library &library, std::size_t M, std::size_t N, std::size_t nvals,
+                const spla::RefPtr<spla::Type> &spT, std::size_t seed = 0) {
+    utils::Matrix a = utils::Matrix<Type>::Generate(M, N, nvals, seed).SortReduceDuplicates();
+
+    a.Fill(utils::UniformGenerator<Type>());
+
+    auto spA = spla::Matrix::Make(M, 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 spTria = spExpr->MakeTril(spA, spA);
+    spExpr->Dependency(spWriteA, spTria);
+    spExpr->SubmitWait();
+
+    ASSERT_EQ(spExpr->GetState(), spla::Expression::State::Evaluated);
+
+    utils::Matrix result = a.Tril();
+    EXPECT_TRUE(result.Equals(spA));
+}
+
+void testNoValues(spla::Library &library, std::size_t M, std::size_t N, std::size_t nvals, std::size_t seed = 0) {
+    using Type = char;
+
+    utils::Matrix a = utils::Matrix<Type>::Generate(M, N, nvals, seed).SortReduceDuplicates();
+
+    auto spT = spla::Types::Void(library);
+    auto spA = spla::Matrix::Make(M, 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 spTria = spExpr->MakeTril(spA, spA);
+    spExpr->Dependency(spWriteA, spTria);
+    spExpr->SubmitWait();
+
+    ASSERT_EQ(spExpr->GetState(), spla::Expression::State::Evaluated);
+
+    utils::Matrix result = a.Tril();
+    EXPECT_TRUE(result.EqualsStructure(spA));
+}
+
+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);
+
+        for (std::size_t i = 0; i < iter; i++) {
+            std::size_t nvals = base + i * step;
+            testCommon<T>(library, M, N, nvals, spT, i);
+        }
+    });
+
+    utils::testBlocks(blocksSizes, [=](spla::Library &library) {
+        using T = std::int32_t;
+        auto spT = spla::Types::Int32(library);
+
+        for (std::size_t i = 0; i < iter; i++) {
+            std::size_t nvals = base + i * step;
+            testCommon<T>(library, M, N, nvals, spT, 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(Tril, Small) {
+    std::vector<std::size_t> blocksSizes{100, 1000};
+    std::size_t M = 100;
+    std::size_t N = 120;
+    test(M, N, M, M, 10, blocksSizes);
+}
+
+TEST(Tril, Medium) {
+    std::vector<std::size_t> blocksSizes{1000, 10000};
+    std::size_t M = 900;
+    std::size_t N = 1200;
+    test(M, N, M, M, 10, blocksSizes);
+}
+
+TEST(Tril, Large) {
+    std::vector<std::size_t> blocksSizes{10000, 100000};
+    std::size_t M = 8800;
+    std::size_t N = 12700;
+    test(M, N, M, M, 5, blocksSizes);
+}
+
+TEST(Tril, MegaLarge) {
+    std::vector<std::size_t> blocksSizes{1000000};
+    std::size_t M = 900100;
+    std::size_t N = 970500;
+    test(M, N, M, M, 2, blocksSizes);
+}
+
+SPLA_GTEST_MAIN
\ No newline at end of file
diff --git a/tests/TestTriu.cpp b/tests/TestTriu.cpp
new file mode 100644
index 000000000..a41f17c3e
--- /dev/null
+++ b/tests/TestTriu.cpp
@@ -0,0 +1,138 @@
+/**********************************************************************************/
+/* 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>
+void testCommon(spla::Library &library, std::size_t M, std::size_t N, std::size_t nvals,
+                const spla::RefPtr<spla::Type> &spT, std::size_t seed = 0) {
+    utils::Matrix a = utils::Matrix<Type>::Generate(M, N, nvals, seed).SortReduceDuplicates();
+
+    a.Fill(utils::UniformGenerator<Type>());
+
+    auto spA = spla::Matrix::Make(M, 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 spTria = spExpr->MakeTriu(spA, spA);
+    spExpr->Dependency(spWriteA, spTria);
+    spExpr->SubmitWait();
+
+    ASSERT_EQ(spExpr->GetState(), spla::Expression::State::Evaluated);
+
+    utils::Matrix result = a.Triu();
+    EXPECT_TRUE(result.Equals(spA));
+}
+
+void testNoValues(spla::Library &library, std::size_t M, std::size_t N, std::size_t nvals, std::size_t seed = 0) {
+    using Type = char;
+
+    utils::Matrix a = utils::Matrix<Type>::Generate(M, N, nvals, seed).SortReduceDuplicates();
+
+    auto spT = spla::Types::Void(library);
+    auto spA = spla::Matrix::Make(M, 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 spTria = spExpr->MakeTriu(spA, spA);
+    spExpr->Dependency(spWriteA, spTria);
+    spExpr->SubmitWait();
+
+    ASSERT_EQ(spExpr->GetState(), spla::Expression::State::Evaluated);
+
+    utils::Matrix result = a.Triu();
+    EXPECT_TRUE(result.EqualsStructure(spA));
+}
+
+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);
+
+        for (std::size_t i = 0; i < iter; i++) {
+            std::size_t nvals = base + i * step;
+            testCommon<T>(library, M, N, nvals, spT, i);
+        }
+    });
+
+    utils::testBlocks(blocksSizes, [=](spla::Library &library) {
+        using T = std::int32_t;
+        auto spT = spla::Types::Int32(library);
+
+        for (std::size_t i = 0; i < iter; i++) {
+            std::size_t nvals = base + i * step;
+            testCommon<T>(library, M, N, nvals, spT, 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(Triu, Small) {
+    std::vector<std::size_t> blocksSizes{100, 1000};
+    std::size_t M = 100;
+    std::size_t N = 120;
+    test(M, N, M, M, 10, blocksSizes);
+}
+
+TEST(Triu, Medium) {
+    std::vector<std::size_t> blocksSizes{1000, 10000};
+    std::size_t M = 900;
+    std::size_t N = 1200;
+    test(M, N, M, M, 10, blocksSizes);
+}
+
+TEST(Triu, Large) {
+    std::vector<std::size_t> blocksSizes{10000, 100000};
+    std::size_t M = 8800;
+    std::size_t N = 12700;
+    test(M, N, M, M, 5, blocksSizes);
+}
+
+TEST(Triu, MegaLarge) {
+    std::vector<std::size_t> blocksSizes{1000000};
+    std::size_t M = 900100;
+    std::size_t N = 970500;
+    test(M, N, M, M, 2, blocksSizes);
+}
+
+SPLA_GTEST_MAIN
\ No newline at end of file
diff --git a/tests/utils/Matrix.hpp b/tests/utils/Matrix.hpp
index 63edebe4a..0cf28d51f 100644
--- a/tests/utils/Matrix.hpp
+++ b/tests/utils/Matrix.hpp
@@ -287,6 +287,22 @@ namespace utils {
             return Matrix<T>(GetNrows(), GetNcols(), std::move(rows), std::move(cols), std::move(vals));
         }
 
+        [[nodiscard]] Matrix RemoveLoops() const {
+            std::vector<Index> rows;
+            std::vector<Index> cols;
+            std::vector<T> vals;
+
+            for (std::size_t k = 0; k < GetNvals(); k++) {
+                if (GetRows()[k] != GetCols()[k]) {
+                    rows.push_back(GetRows()[k]);
+                    cols.push_back(GetCols()[k]);
+                    vals.push_back(GetVals()[k]);
+                }
+            }
+
+            return Matrix<T>(GetNrows(), GetNcols(), std::move(rows), std::move(cols), std::move(vals));
+        }
+
         static Matrix Generate(size_t nrows, size_t ncols, size_t nvals, size_t seed = 0, const T &value = T()) {
             std::vector<Index> rows;
             std::vector<Index> cols;
@@ -533,6 +549,35 @@ namespace utils {
             return Matrix<T>(N, M, std::move(rowsT), std::move(colsT), std::move(valsT));
         }
 
+        template<typename Selector>
+        Matrix<T> Select(Selector selector) const {
+            std::size_t M = GetNrows();
+            std::size_t N = GetNcols();
+            std::size_t nvals = GetNvals();
+
+            std::vector<Index> rowsTria;
+            std::vector<Index> colsTria;
+            std::vector<T> valsTria;
+
+            for (std::size_t k = 0; k < nvals; k++) {
+                if (selector(GetRows()[k], GetCols()[k])) {
+                    rowsTria.push_back(GetRows()[k]);
+                    colsTria.push_back(GetCols()[k]);
+                    valsTria.push_back(GetVals()[k]);
+                }
+            }
+
+            return Matrix<T>(M, N, std::move(rowsTria), std::move(colsTria), std::move(valsTria));
+        }
+
+        Matrix<T> Tril() const {
+            return Select([](Index i, Index j) { return i > j; });
+        }
+
+        Matrix<T> Triu() const {
+            return Select([](Index i, Index j) { return i < j; });
+        }
+
         template<typename ReduceT, typename R = std::invoke_result_t<ReduceT, T, T>>
         [[nodiscard]] R Reduce(ReduceT reduce) {
             if (mVals.empty()) {