From 11c8a4b0079fe81b1cc003b55a56014f9d4891c3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Tom=C3=A1=C5=A1=20Oberhuber?= <oberhuber.tomas@gmail.com>
Date: Thu, 22 Apr 2021 19:59:02 +0200
Subject: [PATCH] Added sanbox sparse matrix.

---
 .../Algorithms/Segments/SegmentsPrinting.h    |    3 +
 .../Matrices/Sandbox/SparseSandboxMatrix.h    | 1176 ++++++++++++++++
 .../Matrices/Sandbox/SparseSandboxMatrix.hpp  | 1197 +++++++++++++++++
 .../Sandbox/SparseSandboxMatrixRowView.h      |  282 ++++
 .../Sandbox/SparseSandboxMatrixRowView.hpp    |  240 ++++
 .../Sandbox/SparseSandboxMatrixView.h         |  871 ++++++++++++
 .../Sandbox/SparseSandboxMatrixView.hpp       | 1026 ++++++++++++++
 src/TNL/Matrices/SparseMatrix.hpp             |    1 +
 src/UnitTests/Matrices/CMakeLists.txt         |    2 +
 .../SparseMatrixTest_SandboxMatrix.cpp        |   11 +
 .../SparseMatrixTest_SandboxMatrix.cu         |    1 +
 .../Matrices/SparseMatrixTest_SandboxMatrix.h |   45 +
 ...eMatrixVectorProductTest_SandboxMatrix.cpp |   11 +
 ...seMatrixVectorProductTest_SandboxMatrix.cu |    1 +
 ...rseMatrixVectorProductTest_SandboxMatrix.h |   45 +
 15 files changed, 4912 insertions(+)
 create mode 100644 src/TNL/Matrices/Sandbox/SparseSandboxMatrix.h
 create mode 100644 src/TNL/Matrices/Sandbox/SparseSandboxMatrix.hpp
 create mode 100644 src/TNL/Matrices/Sandbox/SparseSandboxMatrixRowView.h
 create mode 100644 src/TNL/Matrices/Sandbox/SparseSandboxMatrixRowView.hpp
 create mode 100644 src/TNL/Matrices/Sandbox/SparseSandboxMatrixView.h
 create mode 100644 src/TNL/Matrices/Sandbox/SparseSandboxMatrixView.hpp
 create mode 100644 src/UnitTests/Matrices/SparseMatrixTest_SandboxMatrix.cpp
 create mode 120000 src/UnitTests/Matrices/SparseMatrixTest_SandboxMatrix.cu
 create mode 100644 src/UnitTests/Matrices/SparseMatrixTest_SandboxMatrix.h
 create mode 100644 src/UnitTests/Matrices/SparseMatrixVectorProductTest_SandboxMatrix.cpp
 create mode 120000 src/UnitTests/Matrices/SparseMatrixVectorProductTest_SandboxMatrix.cu
 create mode 100644 src/UnitTests/Matrices/SparseMatrixVectorProductTest_SandboxMatrix.h

diff --git a/src/TNL/Algorithms/Segments/SegmentsPrinting.h b/src/TNL/Algorithms/Segments/SegmentsPrinting.h
index f8fd7412ed..491f059fa6 100644
--- a/src/TNL/Algorithms/Segments/SegmentsPrinting.h
+++ b/src/TNL/Algorithms/Segments/SegmentsPrinting.h
@@ -49,6 +49,8 @@ std::ostream& printSegments( const Segments& segments, std::ostream& str )
    return str;
 }
 
+/// This is to prevent from appearing in Doxygen documentation.
+/// \cond HIDDEN_CLASS
 template< typename Segments,
           typename Fetch >
 struct SegmentsPrinter
@@ -130,6 +132,7 @@ std::ostream& printSegments( const Segments& segments, Fetch&& fetch, std::ostre
    }
    return str;
 }
+/// \endcond
 
       } // namespace Segments
    } // namespace Algorithms
diff --git a/src/TNL/Matrices/Sandbox/SparseSandboxMatrix.h b/src/TNL/Matrices/Sandbox/SparseSandboxMatrix.h
new file mode 100644
index 0000000000..6a2a6a5652
--- /dev/null
+++ b/src/TNL/Matrices/Sandbox/SparseSandboxMatrix.h
@@ -0,0 +1,1176 @@
+/***************************************************************************
+                          SparseSandboxMatrix.h -  description
+                             -------------------
+    begin                : Apr 19, 2021
+    copyright            : (C) 2021 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#pragma once
+
+#include <map>
+#include <TNL/Matrices/Matrix.h>
+#include <TNL/Matrices/MatrixType.h>
+#include <TNL/Allocators/Default.h>
+#include <TNL/Algorithms/Segments/CSR.h>
+#include <TNL/Matrices/Sandbox/SparseSandboxMatrixRowView.h>
+#include <TNL/Matrices/Sandbox/SparseSandboxMatrixView.h>
+#include <TNL/Matrices/DenseMatrix.h>
+
+namespace TNL {
+   namespace Matrices {
+      namespace Sandbox {
+
+/**
+ * \brief Template of a sparse matrix that can be used for testing of new sparse-matrix formats.
+ *
+ * \tparam Real is a type of matrix elements. If \e Real equals \e bool the matrix is treated
+ *    as binary and so the matrix elements values are not stored in the memory since we need
+ *    to remember only coordinates of non-zero elements( which equal one).
+ * \tparam Device is a device where the matrix is allocated.
+ * \tparam Index is a type for indexing of the matrix elements.
+ * \tparam MatrixType specifies a symmetry of matrix. See \ref MatrixType. Symmetric
+ *    matrices store only lower part of the matrix and its diagonal. The upper part is reconstructed on the fly.
+ *    GeneralMatrix with no symmetry is used by default.
+ * \tparam RealAllocator is allocator for the matrix elements values.
+ * \tparam IndexAllocator is allocator for the matrix elements column indexes.
+ *
+ * This class can be used for rapid testing and development of new formats for sparse matrices. One may profit from
+ * several TNL tools compatible with interface of this templated class like:
+ *
+ * 1. Large set of existing unit tests.
+ * 3. Matrix reading from MTX files - to use \ref TNL::Matrices::MatrixReader, the following methods must be functional
+ *    a. \ref TNL::Matrices::SandboxSparseMatrix::setRowCapacities
+ *    b. \ref TNL::Matrices::SandboxSparseMatrix::setElement
+ *    c. \ref TNL::Matrices::SandboxSparseMatrix::operator= between different devices
+ * 4. Matrix benchmarks - the following methods must be functional
+ *    a. \ref TNL::Matrices::SandboxSparseMatrix::vectorProduct - for SpMV benchmark
+ * 5. Linear solvers
+ * 6. Simple comparison of performance with other matrix formats
+ *
+ * In the core of this class there is:
+ *
+ * 1. Vector 'values` (\ref TNL::Matrices::Matrix::values) which is inheritted from \ref TNL::Matrices::Matrix. This vector is used for storing
+ *    of matrix elements values.
+ * 2. Vector `columnIndexes` (\ref TNL::Matrices::SendboxMatrix::columnIndexes). This vector is used for storing of matrix elements column indexes.
+ *
+ * This class contains fully functional implementation of CSR format and so the user have to replace just what he needs to. Once you have
+ * successfully implemented the sparse matrix format in this form, you may consider to extract it into a form of segments to make it accessible
+ * even for other algorithms then SpMV.
+ *
+ * Parts of the code, that need to be modified are marked by SANDBOX_TODO tag. The whole implementation consits of the following classes:
+ *
+ * 1. \ref TNL::Matrices::Sandbox::SparseSandboxMatrix - this class, it serves for matrix setup and performing of the main operations.
+ * 2. \ref TNL::Matrices::Sandbox::SparseSandboxMatrixView - view class which is necessary mainly for passing the matrix to GPU kernels. Most methods of `SparseSandboxMatrix` are common
+ *    with `SparseSandboxMatrixView` and in this case they are implemented in the view class (and there is just redirection from this class). For this reason, `SparseSandboxMatrix` contains instance of the view class
+ *    (\ref TNL::Matrices::Sandbox::SparseSandboxMatrix::view) which needs to be regularly updated each time when metadata are changed. This is usually done by the means of
+ *    method \ref TNL::Matrices::Sandbox::SparseSandboxMatrix::getView.
+ * 3. \ref TNL::Matrices::Sandbox::SparseSandboxMatrixRowView - is a class for accessing particular matrix rows. It will, likely, require some changes as well.
+ *
+ * We suggest the following way of implementation of the new sparse matrix format:
+ *
+ * 1. Add metadata required by your format next to \ref TNL::Matrices::Sandbox::SparseSandboxMatrix::rowPointers but do not replace the row pointers. It will allow you
+ *    to implement your new format next to the original CSR and to check/compare with the valid CSR implementation any time you get into troubles. The adventage is that all
+ *    unit tests are working properly and you may just focus on modifying one method after another. The unit tests are called from
+ *    `src/UnitTests/Matrices/SparseMatrixTests_SandboxMatrix.h` and `src/UnitTests/Matrices/SparseMatrixVectorProductTests_SandboxMatrix.h`
+ * 2. Modify first the method \ref TNL::Matrices::Sandbox::SparseSandboxMatrix::setRowCapacities which is responsible for the setup of the format metadata.
+ * 3. Continue with modification of constructors, view class, \ref TNL::Matrices::Sandbox::SparseSandoxMatrix::getView and \ref TNL::Matrices::Sandbox::SparseSandoxMatrix::getConstView.
+ * 4. Next you need to modify \ref TNL::Matrices::Sandbox::SparseSandboxMatrix::setElement and \ref TNL::Matrices::Sandbox::SparseSandboxMatrix::getElement methods and assignment operator
+ *    at least for copying the matrix across different devices (i.e. from CPU to GPU). It will allow you to use \ref TNL::Matrices::MatrixReader. We recommend to have the same data layout
+ *    on both CPU and GPU so that the transfer of the matrix from CPU to GPU is trivial.
+ * 5. Finally proceed to \ref TNL::Matrices::Sandbox::SparseSandboxMatrix::vectorProduct to implement SpMV operation. We recommend to implement first the CPU version which is easier to
+ *     debug. Next proceed to GPU version.
+ * 6. When SpMV works it is time to delete the original CSR implementation, i.e. everything around `rowPointers`.
+ * 7. Optimize your implementation to the best performance and test with `tnl-benchmark-spmv` - you need to include your new matrix to `src/Benchmarks/SpMV/spmv.h` and modify this file
+ *    accordingly.
+ * 8. If you want, you may now generalize SpMV to \ref TNL::Matrices::Sandbox::SparseSandboxMatrix::reduceRows method.
+ * 9. If you have `reduceRows` implemented, you may use the original implementation of SpMV based just on the `reduceRows` method.
+ * 10. You may implement \ref TNL::Matrices::Sandbox::SparseSandboxMatrix::forRows and \ref TNL::Matrices::Sandbox::SparseSandboxMatrix::forElements.
+ * 11. Now you have complete implementation of new sparse matrix format. You may turn it into new type of segments (\ref TNL::Algorithms::Segments).
+ *
+ * During the implementation some unit tests may crash. If you do not need them at the moment, you may comment them in files
+ * `src/UnitTests/Matrices/SparseMatrixTests.h` and `src/UnitTests/Matrices/SparseMatrixVectorProductTests.h`
+ */
+template< typename Real =  double,
+          typename Device = Devices::Host,
+          typename Index = int,
+          typename MatrixType = GeneralMatrix,
+          typename RealAllocator = typename Allocators::Default< Device >::template Allocator< Real >,
+          typename IndexAllocator = typename Allocators::Default< Device >::template Allocator< Index > >
+class SparseSandboxMatrix : public Matrix< Real, Device, Index, RealAllocator >
+{
+   static_assert(
+         ! MatrixType::isSymmetric() ||
+         ! std::is_same< Device, Devices::Cuda >::value ||
+         ( std::is_same< Real, float >::value || std::is_same< Real, double >::value || std::is_same< Real, int >::value || std::is_same< Real, long long int >::value ),
+         "Given Real type is not supported by atomic operations on GPU which are necessary for symmetric operations." );
+
+   public:
+
+      // Supporting types - they are not important for the user
+      using BaseType = Matrix< Real, Device, Index, RealAllocator >;
+      using ValuesVectorType = typename Matrix< Real, Device, Index, RealAllocator >::ValuesType;
+      using ValuesViewType = typename ValuesVectorType::ViewType;
+      using ConstValuesViewType = typename ValuesViewType::ConstViewType;
+      using ColumnsIndexesVectorType = Containers::Vector< typename TNL::copy_const< Index >::template from< Real >::type, Device, Index, IndexAllocator >;
+      using ColumnsIndexesViewType = typename ColumnsIndexesVectorType::ViewType;
+      using ConstColumnsIndexesViewType = typename ColumnsIndexesViewType::ConstViewType;
+      using RowsCapacitiesType = Containers::Vector< std::remove_const_t< Index >, Device, Index, IndexAllocator >;
+      using RowsCapacitiesView = Containers::VectorView< std::remove_const_t< Index >, Device, Index >;
+      using ConstRowsCapacitiesView = typename RowsCapacitiesView::ConstViewType;
+
+      /**
+       * \brief Test of symmetric matrix type.
+       *
+       * \return \e true if the matrix is stored as symmetric and \e false otherwise.
+       */
+      static constexpr bool isSymmetric() { return MatrixType::isSymmetric(); };
+
+      /**
+       * \brief Test of binary matrix type.
+       *
+       * \return \e true if the matrix is stored as binary and \e false otherwise.
+       */
+      static constexpr bool isBinary() { return std::is_same< Real, bool >::value; };
+
+      /**
+       * \brief The type of matrix elements.
+       */
+      using RealType = std::remove_const_t< Real >;
+
+      /**
+       * \brief The device where the matrix is allocated.
+       */
+      using DeviceType = Device;
+
+      /**
+       * \brief The type used for matrix elements indexing.
+       */
+      using IndexType = Index;
+
+      /**
+       * \brief The allocator for matrix elements values.
+       */
+      using RealAllocatorType = RealAllocator;
+
+      /**
+       * \brief The allocator for matrix elements column indexes.
+       */
+      using IndexAllocatorType = IndexAllocator;
+
+      /**
+       * \brief Type of related matrix view.
+       *
+       * See \ref SparseSandboxMatrixView.
+       */
+      using ViewType = SparseSandboxMatrixView< Real, Device, Index, MatrixType >;
+
+      /**
+       * \brief Matrix view type for constant instances.
+       *
+       * See \ref SparseSandboxMatrixView.
+       */
+      using ConstViewType = SparseSandboxMatrixView< std::add_const_t< Real >, Device, Index, MatrixType >;
+
+      /**
+       * \brief Type for accessing matrix rows.
+       */
+      using RowView = SparseSandboxMatrixRowView< ValuesViewType, ColumnsIndexesViewType, isBinary() >;
+
+      /**
+       * \brief Type for accessing constant matrix rows.
+       */
+      using ConstRowView = SparseSandboxMatrixRowView< ConstValuesViewType, ConstColumnsIndexesViewType, isBinary() >;;
+
+      /**
+       * \brief Helper type for getting self type or its modifications.
+       */
+      template< typename _Real = Real,
+                typename _Device = Device,
+                typename _Index = Index,
+                typename _MatrixType = MatrixType,
+                typename _RealAllocator = typename Allocators::Default< _Device >::template Allocator< _Real >,
+                typename _IndexAllocator = typename Allocators::Default< _Device >::template Allocator< _Index > >
+      using Self = SparseSandboxMatrix< _Real, _Device, _Index, _MatrixType, _RealAllocator, _IndexAllocator >;
+
+      /**
+       * \brief Type of container for CSR row pointers.
+       *
+       * SANDBOX_TODO: You may replace it with containers for metadata of your format.
+       */
+      using RowPointers = TNL::Containers::Vector< IndexType, DeviceType, IndexType >;
+
+      /**
+       * \brief Constructor only with values and column indexes allocators.
+       *
+       * \param realAllocator is used for allocation of matrix elements values.
+       * \param indexAllocator is used for allocation of matrix elements column indexes.
+       */
+      SparseSandboxMatrix( const RealAllocatorType& realAllocator = RealAllocatorType(),
+                           const IndexAllocatorType& indexAllocator = IndexAllocatorType() );
+
+      /**
+       * \brief Copy constructor.
+       *
+       * \param matrix is the source matrix
+       */
+      SparseSandboxMatrix( const SparseSandboxMatrix& matrix1 ) = default;
+
+      /**
+       * \brief Move constructor.
+       *
+       * \param matrix is the source matrix
+       */
+      SparseSandboxMatrix( SparseSandboxMatrix&& matrix ) = default;
+
+      /**
+       * \brief Constructor with matrix dimensions.
+       *
+       * \param rows is number of matrix rows.
+       * \param columns is number of matrix columns.
+       * \param realAllocator is used for allocation of matrix elements values.
+       * \param indexAllocator is used for allocation of matrix elements column indexes.
+       */
+      template< typename Index_t, std::enable_if_t< std::is_integral< Index_t >::value, int > = 0 >
+      SparseSandboxMatrix( const Index_t rows,
+                           const Index_t columns,
+                           const RealAllocatorType& realAllocator = RealAllocatorType(),
+                           const IndexAllocatorType& indexAllocator = IndexAllocatorType() );
+
+      /**
+       * \brief Constructor with matrix rows capacities and number of columns.
+       *
+       * The number of matrix rows is given by the size of \e rowCapacities list.
+       *
+       * \tparam ListIndex is the initializer list values type.
+       * \param rowCapacities is a list telling how many matrix elements must be
+       *    allocated in each row.
+       * \param columns is the number of matrix columns.
+       * \param realAllocator is used for allocation of matrix elements values.
+       * \param indexAllocator is used for allocation of matrix elements column indexes.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_Constructor_init_list_1.cpp
+       * \par Output
+       * \include SparseMatrixExample_Constructor_init_list_1.out
+       */
+      template< typename ListIndex >
+      explicit SparseSandboxMatrix( const std::initializer_list< ListIndex >& rowCapacities,
+                                    const IndexType columns,
+                                    const RealAllocatorType& realAllocator = RealAllocatorType(),
+                                    const IndexAllocatorType& indexAllocator = IndexAllocatorType() );
+
+      /**
+       * \brief Constructor with matrix rows capacities given as a vector and number of columns.
+       *
+       * The number of matrix rows is given by the size of \e rowCapacities vector.
+       *
+       * \tparam RowCapacitiesVector is the row capacities vector type. Usually it is some of
+       *    \ref TNL::Containers::Array, \ref TNL::Containers::ArrayView, \ref TNL::Containers::Vector or
+       *    \ref TNL::Containers::VectorView.
+       * \param rowCapacities is a vector telling how many matrix elements must be
+       *    allocated in each row.
+       * \param columns is the number of matrix columns.
+       * \param realAllocator is used for allocation of matrix elements values.
+       * \param indexAllocator is used for allocation of matrix elements column indexes.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_Constructor_rowCapacities_vector.cpp
+       * \par Output
+       * \include SparseMatrixExample_Constructor_rowCapacities_vector.out
+       */
+      template< typename RowCapacitiesVector, std::enable_if_t< TNL::IsArrayType< RowCapacitiesVector >::value, int > = 0 >
+      explicit SparseSandboxMatrix( const RowCapacitiesVector& rowCapacities,
+                                    const IndexType columns,
+                                    const RealAllocatorType& realAllocator = RealAllocatorType(),
+                                    const IndexAllocatorType& indexAllocator = IndexAllocatorType() );
+
+      /**
+       * \brief Constructor with matrix dimensions and data in initializer list.
+       *
+       * The matrix elements values are given as a list \e data of triples:
+       * { { row1, column1, value1 },
+       *   { row2, column2, value2 },
+       * ... }.
+       *
+       * \param rows is number of matrix rows.
+       * \param columns is number of matrix columns.
+       * \param data is a list of matrix elements values.
+       * \param realAllocator is used for allocation of matrix elements values.
+       * \param indexAllocator is used for allocation of matrix elements column indexes.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_Constructor_init_list_2.cpp
+       * \par Output
+       * \include SparseMatrixExample_Constructor_init_list_2.out
+       */
+      explicit SparseSandboxMatrix( const IndexType rows,
+                                    const IndexType columns,
+                                    const std::initializer_list< std::tuple< IndexType, IndexType, RealType > >& data,
+                                    const RealAllocatorType& realAllocator = RealAllocatorType(),
+                                    const IndexAllocatorType& indexAllocator = IndexAllocatorType() );
+
+      /**
+       * \brief Constructor with matrix dimensions and data in std::map.
+       *
+       * The matrix elements values are given as a map \e data where keys are
+       * std::pair of matrix coordinates ( {row, column} ) and value is the
+       * matrix element value.
+       *
+       * \tparam MapIndex is a type for indexing rows and columns.
+       * \tparam MapValue is a type for matrix elements values in the map.
+       *
+       * \param rows is number of matrix rows.
+       * \param columns is number of matrix columns.
+       * \param map is std::map containing matrix elements.
+       * \param realAllocator is used for allocation of matrix elements values.
+       * \param indexAllocator is used for allocation of matrix elements column indexes.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_Constructor_std_map.cpp
+       * \par Output
+       * \include SparseMatrixExample_Constructor_std_map.out
+       */
+      template< typename MapIndex,
+                typename MapValue >
+      explicit SparseSandboxMatrix( const IndexType rows,
+                                    const IndexType columns,
+                                    const std::map< std::pair< MapIndex, MapIndex >, MapValue >& map,
+                                    const RealAllocatorType& realAllocator = RealAllocatorType(),
+                                    const IndexAllocatorType& indexAllocator = IndexAllocatorType() );
+
+      /**
+       * \brief Returns a modifiable view of the sparse matrix.
+       *
+       * See \ref SparseSandboxMatrixView.
+       *
+       * \return sparse matrix view.
+       */
+      ViewType getView() const; // TODO: remove const
+
+      /**
+       * \brief Returns a non-modifiable view of the sparse matrix.
+       *
+       * See \ref SparseSandboxMatrixView.
+       *
+       * \return sparse matrix view.
+       */
+      ConstViewType getConstView() const;
+
+      /**
+       * \brief Returns string with serialization type.
+       *
+       * The string has a form `Matrices::SparseSandboxMatrix< RealType,  [any_device], IndexType, General/Symmetric, Format, [any_allocator] >`.
+       *
+       * \return \ref String with the serialization type.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_getSerializationType.cpp
+       * \par Output
+       * \include SparseMatrixExample_getSerializationType.out
+       */
+      static String getSerializationType();
+
+      /**
+       * \brief Returns string with serialization type.
+       *
+       * See \ref SparseSandboxMatrix::getSerializationType.
+       *
+       * \return \e String with the serialization type.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_getSerializationType.cpp
+       * \par Output
+       * \include SparseMatrixExample_getSerializationType.out
+       */
+      virtual String getSerializationTypeVirtual() const override;
+
+      /**
+       * \brief Set number of rows and columns of this matrix.
+       *
+       * \param rows is the number of matrix rows.
+       * \param columns is the number of matrix columns.
+       */
+      virtual void setDimensions( const IndexType rows,
+                                  const IndexType columns ) override;
+
+      /**
+       * \brief Set the number of matrix rows and columns by the given matrix.
+       *
+       * \tparam Matrix is matrix type. This can be any matrix having methods
+       *  \ref getRows and \ref getColumns.
+       *
+       * \param matrix in the input matrix dimensions of which are to be adopted.
+       */
+      template< typename Matrix >
+      void setLike( const Matrix& matrix );
+
+      /**
+       * \brief Allocates memory for non-zero matrix elements.
+       *
+       * The size of the input vector must be equal to the number of matrix rows.
+       * The number of allocated matrix elements for each matrix row depends on
+       * the sparse matrix format. Some formats may allocate more elements than
+       * required.
+       *
+       * \tparam RowsCapacitiesVector is a type of vector/array used for row
+       *    capacities setting.
+       *
+       * \param rowCapacities is a vector telling the number of required non-zero
+       *    matrix elements in each row.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_setRowCapacities.cpp
+       * \par Output
+       * \include SparseMatrixExample_setRowCapacities.out
+       */
+      template< typename RowsCapacitiesVector >
+      void setRowCapacities( const RowsCapacitiesVector& rowCapacities );
+
+      /**
+       * \brief Compute capacities of all rows.
+       *
+       * The row capacities are not stored explicitly and must be computed.
+       *
+       * \param rowCapacities is a vector where the row capacities will be stored.
+       */
+      template< typename Vector >
+      void getRowCapacities( Vector& rowCapacities ) const;
+
+      /**
+       * \brief This method sets the sparse matrix elements from initializer list.
+       *
+       * The number of matrix rows and columns must be set already.
+       * The matrix elements values are given as a list \e data of triples:
+       * { { row1, column1, value1 },
+       *   { row2, column2, value2 },
+       * ... }.
+       *
+       * \param data is a initializer list of initializer lists representing
+       * list of matrix rows.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_setElements.cpp
+       * \par Output
+       * \include SparseMatrixExample_setElements.out
+       */
+      void setElements( const std::initializer_list< std::tuple< IndexType, IndexType, RealType > >& data );
+
+      /**
+       * \brief This method sets the sparse matrix elements from std::map.
+       *
+       * The matrix elements values are given as a map \e data where keys are
+       * std::pair of matrix coordinates ( {row, column} ) and value is the
+       * matrix element value.
+       *
+       * \tparam MapIndex is a type for indexing rows and columns.
+       * \tparam MapValue is a type for matrix elements values in the map.
+       *
+       * \param map is std::map containing matrix elements.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_setElements_map.cpp
+       * \par Output
+       * \include SparseMatrixExample_setElements_map.out
+       */
+      template< typename MapIndex,
+                typename MapValue >
+      void setElements( const std::map< std::pair< MapIndex, MapIndex > , MapValue >& map );
+
+      /**
+       * \brief Computes number of non-zeros in each row.
+       *
+       * \param rowLengths is a vector into which the number of non-zeros in each row
+       * will be stored.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_getCompressedRowLengths.cpp
+       * \par Output
+       * \include SparseMatrixExample_getCompressedRowLengths.out
+       */
+      template< typename Vector >
+      void getCompressedRowLengths( Vector& rowLengths ) const;
+
+
+      /**
+       * \brief Returns capacity of given matrix row.
+       *
+       * \param row index of matrix row.
+       * \return number of matrix elements allocated for the row.
+       */
+      __cuda_callable__
+      IndexType getRowCapacity( const IndexType row ) const;
+
+      /**
+       * \brief Returns number of non-zero matrix elements.
+       *
+       * This method really counts the non-zero matrix elements and so
+       * it returns zero for matrix having all allocated elements set to zero.
+       *
+       * \return number of non-zero matrix elements.
+       */
+      IndexType getNonzeroElementsCount() const;
+
+      /**
+       * \brief Resets the matrix to zero dimensions.
+       */
+      void reset();
+
+      /**
+       * \brief Constant getter of simple structure for accessing given matrix row.
+       *
+       * \param rowIdx is matrix row index.
+       *
+       * \return RowView for accessing given matrix row.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_getConstRow.cpp
+       * \par Output
+       * \include SparseMatrixExample_getConstRow.out
+       *
+       * See \ref SparseMatrixRowView.
+       */
+      __cuda_callable__
+      const ConstRowView getRow( const IndexType& rowIdx ) const;
+
+      /**
+       * \brief Non-constant getter of simple structure for accessing given matrix row.
+       *
+       * \param rowIdx is matrix row index.
+       *
+       * \return RowView for accessing given matrix row.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_getRow.cpp
+       * \par Output
+       * \include SparseMatrixExample_getRow.out
+       *
+       * See \ref SparseMatrixRowView.
+       */
+      __cuda_callable__
+      RowView getRow( const IndexType& rowIdx );
+
+      /**
+       * \brief Sets element at given \e row and \e column to given \e value.
+       *
+       * This method can be called from the host system (CPU) no matter
+       * where the matrix is allocated. If the matrix is allocated on GPU this method
+       * can be called even from device kernels. If the matrix is allocated in GPU device
+       * this method is called from CPU, it transfers values of each matrix element separately and so the
+       * performance is very low. For higher performance see. \ref SparseMatrix::getRow
+       * or \ref SparseMatrix::forElements and \ref SparseMatrix::forAllElements.
+       * The call may fail if the matrix row capacity is exhausted.
+       *
+       * \param row is row index of the element.
+       * \param column is columns index of the element.
+       * \param value is the value the element will be set to.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_setElement.cpp
+       * \par Output
+       * \include SparseMatrixExample_setElement.out
+       */
+      __cuda_callable__
+      void setElement( const IndexType row,
+                       const IndexType column,
+                       const RealType& value );
+
+      /**
+       * \brief Add element at given \e row and \e column to given \e value.
+       *
+       * This method can be called from the host system (CPU) no matter
+       * where the matrix is allocated. If the matrix is allocated on GPU this method
+       * can be called even from device kernels. If the matrix is allocated in GPU device
+       * this method is called from CPU, it transfers values of each matrix element separately and so the
+       * performance is very low. For higher performance see. \ref SparseMatrix::getRow
+       * or \ref SparseMatrix::forElements and \ref SparseMatrix::forAllElements.
+       * The call may fail if the matrix row capacity is exhausted.
+       *
+       * \param row is row index of the element.
+       * \param column is columns index of the element.
+       * \param value is the value the element will be set to.
+       * \param thisElementMultiplicator is multiplicator the original matrix element
+       *   value is multiplied by before addition of given \e value.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_addElement.cpp
+       * \par Output
+       * \include SparseMatrixExample_addElement.out
+       *
+       */
+      __cuda_callable__
+      void addElement( const IndexType row,
+                       const IndexType column,
+                       const RealType& value,
+                       const RealType& thisElementMultiplicator );
+
+      /**
+       * \brief Returns value of matrix element at position given by its row and column index.
+       *
+       * This method can be called from the host system (CPU) no matter
+       * where the matrix is allocated. If the matrix is allocated on GPU this method
+       * can be called even from device kernels. If the matrix is allocated in GPU device
+       * this method is called from CPU, it transfers values of each matrix element separately and so the
+       * performance is very low. For higher performance see. \ref SparseMatrix::getRow
+       * or \ref SparseMatrix::forElements and \ref SparseMatrix::forAllElements.
+       *
+       * \param row is a row index of the matrix element.
+       * \param column i a column index of the matrix element.
+       *
+       * \return value of given matrix element.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_getElement.cpp
+       * \par Output
+       * \include SparseMatrixExample_getElement.out
+       *
+       */
+      __cuda_callable__
+      RealType getElement( const IndexType row,
+                           const IndexType column ) const;
+
+      /**
+       * \brief Method for performing general reduction on matrix rows.
+       *
+       * \tparam Fetch is a type of lambda function for data fetch declared as
+       *          `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`.
+       *          The return type of this lambda can be any non void.
+       * \tparam Reduce is a type of lambda function for reduction declared as
+       *          `reduce( const FetchValue& v1, const FetchValue& v2 ) -> FetchValue`.
+       * \tparam Keep is a type of lambda function for storing results of reduction in each row.
+       *          It is declared as `keep( const IndexType rowIdx, const double& value )`.
+       * \tparam FetchValue is type returned by the Fetch lambda function.
+       *
+       * \param begin defines beginning of the range [begin,end) of rows to be processed.
+       * \param end defines ending of the range [begin,end) of rows to be processed.
+       * \param fetch is an instance of lambda function for data fetch.
+       * \param reduce is an instance of lambda function for reduction.
+       * \param keep in an instance of lambda function for storing results.
+       * \param zero is zero of given reduction operation also known as idempotent element.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_reduceRows.cpp
+       * \par Output
+       * \include SparseMatrixExample_reduceRows.out
+       */
+      template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >
+      void reduceRows( IndexType begin, IndexType end, Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchReal& zero );
+
+      /**
+       * \brief Method for performing general reduction on matrix rows for constant instances.
+       *
+       * \tparam Fetch is a type of lambda function for data fetch declared as
+       *          `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`.
+       *          The return type of this lambda can be any non void.
+       * \tparam Reduce is a type of lambda function for reduction declared as
+       *          `reduce( const FetchValue& v1, const FetchValue& v2 ) -> FetchValue`.
+       * \tparam Keep is a type of lambda function for storing results of reduction in each row.
+       *          It is declared as `keep( const IndexType rowIdx, const double& value )`.
+       * \tparam FetchValue is type returned by the Fetch lambda function.
+       *
+       * \param begin defines beginning of the range [begin,end) of rows to be processed.
+       * \param end defines ending of the range [begin,end) of rows to be processed.
+       * \param fetch is an instance of lambda function for data fetch.
+       * \param reduce is an instance of lambda function for reduction.
+       * \param keep in an instance of lambda function for storing results.
+       * \param zero is zero of given reduction operation also known as idempotent element.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_reduceRows.cpp
+       * \par Output
+       * \include SparseMatrixExample_reduceRows.out
+       */
+      template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >
+      void reduceRows( IndexType begin, IndexType end, Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchReal& zero ) const;
+
+      /**
+       * \brief Method for performing general reduction on all matrix rows.
+       *
+       * \tparam Fetch is a type of lambda function for data fetch declared as
+       *          `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`.
+       *          The return type of this lambda can be any non void.
+       * \tparam Reduce is a type of lambda function for reduction declared as
+       *          `reduce( const FetchValue& v1, const FetchValue& v2 ) -> FetchValue`.
+       * \tparam Keep is a type of lambda function for storing results of reduction in each row.
+       *          It is declared as `keep( const IndexType rowIdx, const double& value )`.
+       * \tparam FetchValue is type returned by the Fetch lambda function.
+       *
+       * \param fetch is an instance of lambda function for data fetch.
+       * \param reduce is an instance of lambda function for reduction.
+       * \param keep in an instance of lambda function for storing results.
+       * \param zero is zero of given reduction operation also known as idempotent element.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_reduceAllRows.cpp
+       * \par Output
+       * \include SparseMatrixExample_reduceAllRows.out
+       */
+      template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >
+      void reduceAllRows( Fetch&& fetch, const Reduce&& reduce, Keep&& keep, const FetchReal& zero );
+
+      /**
+       * \brief Method for performing general reduction on all matrix rows for constant instances.
+       *
+       * \tparam Fetch is a type of lambda function for data fetch declared as
+       *          `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`.
+       *          The return type of this lambda can be any non void.
+       * \tparam Reduce is a type of lambda function for reduction declared as
+       *          `reduce( const FetchValue& v1, const FetchValue& v2 ) -> FetchValue`.
+       * \tparam Keep is a type of lambda function for storing results of reduction in each row.
+       *          It is declared as `keep( const IndexType rowIdx, const double& value )`.
+       * \tparam FetchValue is type returned by the Fetch lambda function.
+       *
+       * \param fetch is an instance of lambda function for data fetch.
+       * \param reduce is an instance of lambda function for reduction.
+       * \param keep in an instance of lambda function for storing results.
+       * \param zero is zero of given reduction operation also known as idempotent element.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_reduceAllRows.cpp
+       * \par Output
+       * \include SparseMatrixExample_reduceAllRows.out
+       */
+      template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >
+      void reduceAllRows( Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchReal& zero ) const;
+
+      /**
+       * \brief Method for parallel iteration over matrix elements of given rows for constant instances.
+       *
+       * \tparam Function is type of lambda function that will operate on matrix elements.
+       *
+       * \param begin defines beginning of the range [ \e begin,\e end ) of rows to be processed.
+       * \param end defines ending of the range [ \e begin, \e end ) of rows to be processed.
+       * \param function is an instance of the lambda function to be called for element of given rows.
+       *
+       * The lambda function `function` should be declared like follows:
+       *
+       * ```
+       * auto function = [] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value ) { ... };
+       * ```
+       *
+       *  The \e localIdx parameter is a rank of the non-zero element in given row.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_forElements.cpp
+       * \par Output
+       * \include SparseMatrixExample_forElements.out
+       */
+      template< typename Function >
+      void forElements( IndexType begin, IndexType end, Function&& function ) const;
+
+      /**
+       * \brief Method for parallel iteration over all matrix elements of given rows for non-constant instances.
+       *
+       * \tparam Function is type of lambda function that will operate on matrix elements.
+       *
+       * \param begin defines beginning of the range [ \e begin,\e end ) of rows to be processed.
+       * \param end defines ending of the range [ \e begin, \e end ) of rows to be processed.
+       * \param function is an instance of the lambda function to be called for each element of given rows.
+       *
+       * The lambda function `function` should be declared like follows:
+       *
+       * ```
+       * auto function = [] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value ) mutable { ... }
+       * ```
+       *
+       *  The \e localIdx parameter is a rank of the non-zero element in given row.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_forElements.cpp
+       * \par Output
+       * \include SparseMatrixExample_forElements.out
+       */
+      template< typename Function >
+      void forElements( IndexType begin, IndexType end, Function&& function );
+
+      /**
+       * \brief Method for parallel iteration over all matrix elements for constant instances.
+       *
+       * See \ref SparseMatrix::forElements.
+       *
+       * \tparam Function is a type of lambda function that will operate on matrix elements.
+       * \param function  is an instance of the lambda function to be called for each matrix element.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_forElements.cpp
+       * \par Output
+       * \include SparseMatrixExample_forElements.out
+       */
+      template< typename Function >
+      void forAllElements( Function&& function ) const;
+
+      /**
+       * \brief Method for parallel iteration over all matrix elements for non-constant instances.
+       *
+       * See \ref SparseMatrix::forElements.
+       *
+       * \tparam Function is a type of lambda function that will operate on matrix elements.
+       * \param function  is an instance of the lambda function to be called for each matrix element.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_forElements.cpp
+       * \par Output
+       * \include SparseMatrixExample_forElements.out
+       */
+      template< typename Function >
+      void forAllElements( Function&& function );
+
+      /**
+       * \brief Method for parallel iteration over matrix rows from interval [ \e begin, \e end).
+       *
+       * In each row, given lambda function is performed. Each row is processed by at most one thread unlike the method
+       * \ref SparseMatrix::forElements where more than one thread can be mapped to each row.
+       *
+       * \tparam Function is type of the lambda function.
+       *
+       * \param begin defines beginning of the range [ \e begin,\e end ) of rows to be processed.
+       * \param end defines ending of the range [ \e begin, \e end ) of rows to be processed.
+       * \param function is an instance of the lambda function to be called for each row.
+       *
+       * ```
+       * auto function = [] __cuda_callable__ ( RowView& row ) mutable { ... };
+       * ```
+       *
+       * \e RowView represents matrix row - see \ref TNL::Matrices::SparseMatrix::RowView.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_forRows.cpp
+       * \par Output
+       * \include SparseMatrixExample_forRows.out
+       */
+      template< typename Function >
+      void forRows( IndexType begin, IndexType end, Function&& function );
+
+      /**
+       * \brief Method for parallel iteration over matrix rows from interval [ \e begin, \e end) for constant instances.
+       *
+       * In each row, given lambda function is performed. Each row is processed by at most one thread unlike the method
+       * \ref SparseMatrix::forElements where more than one thread can be mapped to each row.
+       *
+       * \tparam Function is type of the lambda function.
+       *
+       * \param begin defines beginning of the range [ \e begin,\e end ) of rows to be processed.
+       * \param end defines ending of the range [ \e begin, \e end ) of rows to be processed.
+       * \param function is an instance of the lambda function to be called for each row.
+       *
+       * ```
+       * auto function = [] __cuda_callable__ ( RowView& row ) { ... };
+       * ```
+       *
+       * \e RowView represents matrix row - see \ref TNL::Matrices::SparseMatrix::RowView.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_forRows.cpp
+       * \par Output
+       * \include SparseMatrixExample_forRows.out
+       */
+      template< typename Function >
+      void forRows( IndexType begin, IndexType end, Function&& function ) const;
+
+      /**
+       * \brief Method for parallel iteration over all matrix rows.
+       *
+       * In each row, given lambda function is performed. Each row is processed by at most one thread unlike the method
+       * \ref SparseMatrix::forAllElements where more than one thread can be mapped to each row.
+       *
+       * \tparam Function is type of the lambda function.
+       *
+       * \param function is an instance of the lambda function to be called for each row.
+       *
+       * ```
+       * auto function = [] __cuda_callable__ ( RowView& row ) mutable { ... };
+       * ```
+       *
+       * \e RowView represents matrix row - see \ref TNL::Matrices::SparseMatrix::RowView.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_forRows.cpp
+       * \par Output
+       * \include SparseMatrixExample_forRows.out
+       */
+      template< typename Function >
+      void forAllRows( Function&& function );
+
+      /**
+       * \brief Method for parallel iteration over all matrix rows for constant instances.
+       *
+       * In each row, given lambda function is performed. Each row is processed by at most one thread unlike the method
+       * \ref SparseMatrix::forAllElements where more than one thread can be mapped to each row.
+       *
+       * \tparam Function is type of the lambda function.
+       *
+       * \param function is an instance of the lambda function to be called for each row.
+       *
+       * ```
+       * auto function = [] __cuda_callable__ ( RowView& row ) { ... };
+       * ```
+       *
+       * \e RowView represents matrix row - see \ref TNL::Matrices::SparseMatrix::RowView.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_forRows.cpp
+       * \par Output
+       * \include SparseMatrixExample_forRows.out
+       */
+      template< typename Function >
+      void forAllRows( Function&& function ) const;
+
+      /**
+       * \brief Method for sequential iteration over all matrix rows for constant instances.
+       *
+       * \tparam Function is type of lambda function that will operate on matrix elements.
+       *    It is should have form like
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value )`.
+       *  The column index repeats twice only for compatibility with sparse matrices.
+       *
+       * \param begin defines beginning of the range [begin,end) of rows to be processed.
+       * \param end defines ending of the range [begin,end) of rows to be processed.
+       * \param function is an instance of the lambda function to be called in each row.
+       */
+      template< typename Function >
+      void sequentialForRows( IndexType begin, IndexType end, Function& function ) const;
+
+      /**
+       * \brief Method for sequential iteration over all matrix rows for non-constant instances.
+       *
+       * \tparam Function is type of lambda function that will operate on matrix elements.
+       *    It is should have form like
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, RealType& value )`.
+       *  The column index repeats twice only for compatibility with sparse matrices.
+       *
+       * \param begin defines beginning of the range [begin,end) of rows to be processed.
+       * \param end defines ending of the range [begin,end) of rows to be processed.
+       * \param function is an instance of the lambda function to be called in each row.
+       */
+      template< typename Function >
+      void sequentialForRows( IndexType begin, IndexType end, Function& function );
+
+      /**
+       * \brief This method calls \e sequentialForRows for all matrix rows (for constant instances).
+       *
+       * See \ref SparseMatrix::sequentialForRows.
+       *
+       * \tparam Function is a type of lambda function that will operate on matrix elements.
+       * \param function  is an instance of the lambda function to be called in each row.
+       */
+      template< typename Function >
+      void sequentialForAllRows( Function& function ) const;
+
+      /**
+       * \brief This method calls \e sequentialForRows for all matrix rows.
+       *
+       * See \ref SparseMatrix::sequentialForAllRows.
+       *
+       * \tparam Function is a type of lambda function that will operate on matrix elements.
+       * \param function  is an instance of the lambda function to be called in each row.
+       */
+      template< typename Function >
+      void sequentialForAllRows( Function& function );
+
+      /**
+       * \brief Computes product of matrix and vector.
+       *
+       * More precisely, it computes:
+       *
+       * `outVector = matrixMultiplicator * ( * this ) * inVector + outVectorMultiplicator * outVector`
+       *
+       * \tparam InVector is type of input vector.  It can be \ref Vector,
+       *     \ref VectorView, \ref Array, \ref ArraView or similar container.
+       * \tparam OutVector is type of output vector. It can be \ref Vector,
+       *     \ref VectorView, \ref Array, \ref ArraView or similar container.
+       *
+       * \param inVector is input vector.
+       * \param outVector is output vector.
+       * \param matrixMultiplicator is a factor by which the matrix is multiplied. It is one by default.
+       * \param outVectorMultiplicator is a factor by which the outVector is multiplied before added
+       *    to the result of matrix-vector product. It is zero by default.
+       * \param begin is the beginning of the rows range for which the vector product
+       *    is computed. It is zero by default.
+       * \param end is the end of the rows range for which the vector product
+       *    is computed. It is number if the matrix rows by default.
+       */
+      template< typename InVector,
+                typename OutVector >
+      void vectorProduct( const InVector& inVector,
+                          OutVector& outVector,
+                          const RealType& matrixMultiplicator = 1.0,
+                          const RealType& outVectorMultiplicator = 0.0,
+                          const IndexType firstRow = 0,
+                          const IndexType lastRow = 0 ) const;
+
+      /*template< typename Real2, typename Index2 >
+      void addMatrix( const SparseMatrix< Real2, Segments, Device, Index2 >& matrix,
+                      const RealType& matrixMultiplicator = 1.0,
+                      const RealType& thisMatrixMultiplicator = 1.0 );
+
+      template< typename Real2, typename Index2 >
+      void getTransposition( const SparseMatrix< Real2, Segments, Device, Index2 >& matrix,
+                             const RealType& matrixMultiplicator = 1.0 );
+       */
+
+      template< typename Vector1, typename Vector2 >
+      bool performSORIteration( const Vector1& b,
+                                const IndexType row,
+                                Vector2& x,
+                                const RealType& omega = 1.0 ) const;
+
+      /**
+       * \brief Assignment of exactly the same matrix type.
+       *
+       * \param matrix is input matrix for the assignment.
+       * \return reference to this matrix.
+       */
+      SparseSandboxMatrix& operator=( const SparseSandboxMatrix& matrix );
+
+      /**
+       * \brief Assignment of exactly the same matrix type but different device.
+       *
+       * \param matrix is input matrix for the assignment.
+       * \return reference to this matrix.
+       */
+      template< typename Device_ >
+      SparseSandboxMatrix& operator=( const SparseSandboxMatrix< RealType, Device_, IndexType, MatrixType, RealAllocator, IndexAllocator >& matrix );
+
+      /**
+       * \brief Assignment of dense matrix
+       *
+       * \param matrix is input matrix for the assignment.
+       * \return reference to this matrix.
+       */
+      template< typename Real_, typename Device_, typename Index_, ElementsOrganization Organization, typename RealAllocator_ >
+      SparseSandboxMatrix& operator=( const DenseMatrix< Real_, Device_, Index_, Organization, RealAllocator_ >& matrix );
+
+
+      /**
+       * \brief Assignment of any matrix type other then this and dense.
+       *
+       * **Warning: Assignment of symmetric sparse matrix to general sparse matrix does not give correct result, currently. Only the diagonal and the lower part of the matrix is assigned.**
+       *
+       * \param matrix is input matrix for the assignment.
+       * \return reference to this matrix.
+       */
+      template< typename RHSMatrix >
+      SparseSandboxMatrix& operator=( const RHSMatrix& matrix );
+
+      /**
+       * \brief Comparison operator with another arbitrary matrix type.
+       *
+       * \param matrix is the right-hand side matrix.
+       * \return \e true if the RHS matrix is equal, \e false otherwise.
+       */
+      template< typename Matrix >
+      bool operator==( const Matrix& m ) const;
+
+      /**
+       * \brief Comparison operator with another arbitrary matrix type.
+       *
+       * \param matrix is the right-hand side matrix.
+       * \return \e true if the RHS matrix is equal, \e false otherwise.
+       */
+      template< typename Matrix >
+      bool operator!=( const Matrix& m ) const;
+
+      /**
+       * \brief Method for saving the matrix to the file with given filename.
+       *
+       * \param fileName is name of the file.
+       */
+      void save( const String& fileName ) const;
+
+      /**
+       * \brief Method for loading the matrix from the file with given filename.
+       *
+       * \param fileName is name of the file.
+       */
+      void load( const String& fileName );
+
+      /**
+       * \brief Method for saving the matrix to a file.
+       *
+       * \param file is the output file.
+       */
+      virtual void save( File& file ) const override;
+
+      /**
+       * \brief Method for loading the matrix from a file.
+       *
+       * \param file is the input file.
+       */
+      virtual void load( File& file ) override;
+
+      /**
+       * \brief Method for printing the matrix to output stream.
+       *
+       * \param str is the output stream.
+       */
+      virtual void print( std::ostream& str ) const override;
+
+      /**
+       * \brief Returns a padding index value.
+       *
+       * Padding index is used for column indexes of padding zeros. Padding zeros
+       * are used in some sparse matrix formats for better data alignment in memory.
+       *
+       * \return value of the padding index.
+       */
+      __cuda_callable__
+      IndexType getPaddingIndex() const;
+
+      /**
+       * \brief Getter of segments for non-constant instances.
+       *
+       * \e Segments are a structure for addressing the matrix elements columns and values.
+       * In fact, \e Segments represent the sparse matrix format.
+       *
+       * \return Non-constant reference to segments.
+       */
+      //SegmentsType& getSegments();
+
+      /**
+       * \brief Getter of segments for constant instances.
+       *
+       * \e Segments are a structure for addressing the matrix elements columns and values.
+       * In fact, \e Segments represent the sparse matrix format.
+       *
+       * \return Constant reference to segments.
+       */
+      //const SegmentsType& getSegments() const;
+
+      /**
+       * \brief Getter of column indexes for constant instances.
+       *
+       * \return Constant reference to a vector with matrix elements column indexes.
+       */
+      const ColumnsIndexesVectorType& getColumnIndexes() const;
+
+      /**
+       * \brief Getter of column indexes for nonconstant instances.
+       *
+       * \return Reference to a vector with matrix elements column indexes.
+       */
+      ColumnsIndexesVectorType& getColumnIndexes();
+
+   protected:
+
+      ColumnsIndexesVectorType columnIndexes;
+
+      IndexAllocator indexAllocator;
+
+      ViewType view;
+
+      /**
+       * \brief Container for CSR row pointers.
+       *
+       * SANDBOX_TODO: You may replace it with containers and metadata required by you format.
+       */
+
+      RowPointers rowPointers;
+};
+
+      } // namespace Sandbox
+   } // namespace Matrices
+} // namespace TNL
+
+#include <TNL/Matrices/Sandbox/SparseSandboxMatrix.hpp>
diff --git a/src/TNL/Matrices/Sandbox/SparseSandboxMatrix.hpp b/src/TNL/Matrices/Sandbox/SparseSandboxMatrix.hpp
new file mode 100644
index 0000000000..63f49e6c84
--- /dev/null
+++ b/src/TNL/Matrices/Sandbox/SparseSandboxMatrix.hpp
@@ -0,0 +1,1197 @@
+/***************************************************************************
+                          SparseSandboxMatrix.hpp -  description
+                             -------------------
+    begin                : Apr 19, 2021
+    copyright            : (C) 2021 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#pragma once
+
+#include <functional>
+#include <sstream>
+#include <TNL/Algorithms/Reduction.h>
+#include <TNL/Matrices/Sandbox/SparseSandboxMatrix.h>
+
+namespace TNL {
+   namespace Matrices {
+      namespace Sandbox {
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+SparseSandboxMatrix( const RealAllocatorType& realAllocator,
+                     const IndexAllocatorType& indexAllocator )
+: BaseType( realAllocator ), columnIndexes( indexAllocator ), rowPointers( ( IndexType ) 1, ( IndexType ) 0, indexAllocator )
+{
+   this->view = this->getView();
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Index_t, std::enable_if_t< std::is_integral< Index_t >::value, int > >
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+SparseSandboxMatrix( const Index_t rows,
+                     const Index_t columns,
+                     const RealAllocatorType& realAllocator,
+                     const IndexAllocatorType& indexAllocator )
+: BaseType( rows, columns, realAllocator ), columnIndexes( indexAllocator ), rowPointers( rows + 1, ( IndexType ) 0, indexAllocator )
+{
+   this->view = this->getView();
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename ListIndex >
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+SparseSandboxMatrix( const std::initializer_list< ListIndex >& rowCapacities,
+                     const IndexType columns,
+                     const RealAllocatorType& realAllocator,
+                     const IndexAllocatorType& indexAllocator )
+: BaseType( rowCapacities.size(), columns, realAllocator ), columnIndexes( indexAllocator ), rowPointers( rowCapacities.size() + 1, ( IndexType ) 0, indexAllocator )
+{
+   this->setRowCapacities( RowsCapacitiesType( rowCapacities ) );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename RowCapacitiesVector, std::enable_if_t< TNL::IsArrayType< RowCapacitiesVector >::value, int > >
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+SparseSandboxMatrix( const RowCapacitiesVector& rowCapacities,
+                     const IndexType columns,
+                     const RealAllocatorType& realAllocator,
+                     const IndexAllocatorType& indexAllocator )
+: BaseType( rowCapacities.getSize(), columns, realAllocator ), columnIndexes( indexAllocator ), rowPointers( rowCapacities.getSize() + 1, ( IndexType ) 0, indexAllocator )
+{
+   this->setRowCapacities( rowCapacities );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+SparseSandboxMatrix( const IndexType rows,
+                     const IndexType columns,
+                     const std::initializer_list< std::tuple< IndexType, IndexType, RealType > >& data,
+                     const RealAllocatorType& realAllocator,
+                     const IndexAllocatorType& indexAllocator )
+: BaseType( rows, columns, realAllocator ), columnIndexes( indexAllocator ), rowPointers( rows + 1, ( IndexType ) 0, indexAllocator )
+{
+   this->setElements( data );
+   this->view = this->getView();
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename MapIndex,
+             typename MapValue >
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+SparseSandboxMatrix( const IndexType rows,
+              const IndexType columns,
+              const std::map< std::pair< MapIndex, MapIndex > , MapValue >& map,
+              const RealAllocatorType& realAllocator,
+              const IndexAllocatorType& indexAllocator )
+: BaseType( rows, columns, realAllocator ), columnIndexes( indexAllocator ), rowPointers( rows + 1, ( IndexType ) 0, indexAllocator )
+{
+   this->setDimensions( rows, columns );
+   this->setElements( map );
+   this->view = this->getView();
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+auto
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+getView() const -> ViewType
+{
+   return ViewType( this->getRows(),
+                    this->getColumns(),
+                    const_cast< SparseSandboxMatrix* >( this )->getValues().getView(),  // TODO: remove const_cast
+                    const_cast< SparseSandboxMatrix* >( this )->columnIndexes.getView(),
+                    const_cast< SparseSandboxMatrix* >( this )->rowPointers.getView() );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+auto
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+getConstView() const -> ConstViewType
+{
+   return ConstViewType( this->getRows(),
+                         this->getColumns(),
+                         this->getValues().getConstView(),
+                         this->columnIndexes.getConstView(),
+                         this->segments.getConstView() );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+String
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+getSerializationType()
+{
+   return ViewType::getSerializationType();
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+String
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+getSerializationTypeVirtual() const
+{
+   return this->getSerializationType();
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+setDimensions( const IndexType rows,
+               const IndexType columns )
+{
+   BaseType::setDimensions( rows, columns );
+   this->view = this->getView();
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Matrix_ >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+setLike( const Matrix_& matrix )
+{
+   BaseType::setLike( matrix );
+   // SANDBOX_TODO: Replace the following line with assignment of metadata required by your format. 
+   //               Do not assign matrix elements here.
+   this->rowPointers = matrix.rowPointers;
+   this->view = this->getView();
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename RowsCapacitiesVector >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+setRowCapacities( const RowsCapacitiesVector& rowsCapacities )
+{
+   TNL_ASSERT_EQ( rowsCapacities.getSize(), this->getRows(), "Number of matrix rows does not fit with rowCapacities vector size." );
+   using RowsCapacitiesVectorDevice = typename RowsCapacitiesVector::DeviceType;
+
+   // SANDBOX_TODO: Replace the following lines with the setup of your sparse matrix format based on
+   //               `rowsCapacities`. This container has the same number of elements as is the number of
+   //               rows of this matrix. Each element says how many nonzero elements the user needs to have
+   //               in each row. This number can be increased if the sparse matrix format uses padding zeros.
+   this->rowPointers.setSize( this->getRows() + 1 );
+   if( std::is_same< DeviceType, RowsCapacitiesVectorDevice >::value )
+   {
+      // GOTCHA: when this->getRows() == 0, getView returns a full view with size == 1
+      if( this->getRows() > 0 ) {
+         auto view = this->rowPointers.getView( 0, this->getRows() );
+         view = rowsCapacities;
+      }
+   }
+   else
+   {
+      RowsCapacitiesType thisRowsCapacities;
+      thisRowsCapacities = rowsCapacities;
+      if( this->getRows() > 0 ) {
+         auto view = this->rowPointers.getView( 0, this->getRows() );
+         view = thisRowsCapacities;
+      }
+   }
+   this->rowPointers.setElement( this->getRows(), 0 );
+   this->rowPointers.template scan< Algorithms::ScanType::Exclusive >();
+   // End of sparse matrix format initiation.
+
+   // SANDBOX_TODO: Compute number of all elements that need to be allocated by your format.
+   const auto storageSize = rowPointers.getElement( this->getRows() );
+
+   // The rest of this methods needs no changes.
+   if( ! isBinary() )
+   {
+      this->values.setSize( storageSize );
+      this->values = ( RealType ) 0;
+   }
+   this->columnIndexes.setSize( storageSize );
+   this->columnIndexes = this->getPaddingIndex();
+   this->view = this->getView();
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Vector >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+getRowCapacities( Vector& rowCapacities ) const
+{
+   this->view.getRowCapacities( rowCapacities );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+setElements( const std::initializer_list< std::tuple< IndexType, IndexType, RealType > >& data )
+{
+   const auto& rows = this->getRows();
+   const auto& columns = this->getColumns();
+   Containers::Vector< IndexType, Devices::Host, IndexType > rowCapacities( rows, 0 );
+   for( const auto& i : data )
+   {
+      if( std::get< 0 >( i ) >= rows )
+      {
+         std::stringstream s;
+         s << "Wrong row index " << std::get< 0 >( i ) << " in an initializer list";
+         throw std::logic_error( s.str() );
+      }
+      rowCapacities[ std::get< 0 >( i ) ]++;
+   }
+   SparseSandboxMatrix< Real, Devices::Host, Index, MatrixType > hostMatrix( rows, columns );
+   hostMatrix.setRowCapacities( rowCapacities );
+   for( const auto& i : data )
+   {
+      if( std::get< 1 >( i ) >= columns )
+      {
+         std::stringstream s;
+         s << "Wrong column index " << std::get< 1 >( i ) << " in an initializer list";
+         throw std::logic_error( s.str() );
+      }
+      hostMatrix.setElement( std::get< 0 >( i ), std::get< 1 >( i ), std::get< 2 >( i ) );
+   }
+   ( *this ) = hostMatrix;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename MapIndex,
+             typename MapValue >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+setElements( const std::map< std::pair< MapIndex, MapIndex > , MapValue >& map )
+{
+   Containers::Vector< IndexType, Devices::Host, IndexType > rowsCapacities( this->getRows(), 0 );
+   for( auto element : map )
+      rowsCapacities[ element.first.first ]++;
+   if( !std::is_same< DeviceType, Devices::Host >::value )
+   {
+      SparseSandboxMatrix< Real, Devices::Host, Index, MatrixType > hostMatrix( this->getRows(), this->getColumns() );
+      hostMatrix.setRowCapacities( rowsCapacities );
+      for( auto element : map )
+         hostMatrix.setElement( element.first.first, element.first.second, element.second );
+      *this = hostMatrix;
+   }
+   else
+   {
+      this->setRowCapacities( rowsCapacities );
+      for( auto element : map )
+         this->setElement( element.first.first, element.first.second, element.second );
+   }
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Vector >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+getCompressedRowLengths( Vector& rowLengths ) const
+{
+   this->view.getCompressedRowLengths( rowLengths );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+__cuda_callable__
+Index
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+getRowCapacity( const IndexType row ) const
+{
+   return this->view.getRowCapacity( row );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+Index
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+getNonzeroElementsCount() const
+{
+   return this->view.getNonzeroElementsCount();
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+reset()
+{
+   BaseType::reset();
+   this->columnIndexes.reset();
+   // SANDBOX_TODO: Reset the metadata required by your format here.
+   this->rowPointers.reset();
+   this->view = this->getView();
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+__cuda_callable__ auto
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+getRow( const IndexType& rowIdx ) const -> const ConstRowView
+{
+   return this->view.getRow( rowIdx );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+__cuda_callable__ auto
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+getRow( const IndexType& rowIdx ) -> RowView
+{
+   return this->view.getRow( rowIdx );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+__cuda_callable__ void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+setElement( const IndexType row,
+            const IndexType column,
+            const RealType& value )
+{
+   this->view.setElement( row, column, value );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+__cuda_callable__ void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+addElement( const IndexType row,
+            const IndexType column,
+            const RealType& value,
+            const RealType& thisElementMultiplicator )
+{
+   this->view.addElement( row, column, value, thisElementMultiplicator );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+__cuda_callable__
+auto
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+getElement( const IndexType row,
+            const IndexType column ) const -> RealType
+{
+   return this->view.getElement( row, column );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+template< typename InVector,
+       typename OutVector >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+vectorProduct( const InVector& inVector,
+               OutVector& outVector,
+               const RealType& matrixMultiplicator,
+               const RealType& outVectorMultiplicator,
+               const IndexType firstRow,
+               const IndexType lastRow ) const
+{
+   this->view.vectorProduct( inVector, outVector, matrixMultiplicator, outVectorMultiplicator, firstRow, lastRow );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Fetch, typename Reduce, typename Keep, typename FetchValue >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+reduceRows( IndexType begin, IndexType end, Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchValue& zero )
+{
+   this->view.reduceRows( begin, end, fetch, reduce, keep, zero );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Fetch, typename Reduce, typename Keep, typename FetchValue >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+reduceRows( IndexType begin, IndexType end, Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchValue& zero ) const
+{
+   this->view.reduceRows( begin, end, fetch, reduce, keep, zero );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+reduceAllRows( Fetch&& fetch, const Reduce&& reduce, Keep&& keep, const FetchReal& zero )
+{
+   this->reduceRows( 0, this->getRows(), fetch, reduce, keep, zero );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+reduceAllRows( Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchReal& zero ) const
+{
+   this->reduceRows( 0, this->getRows(), fetch, reduce, keep, zero );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Function >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+forElements( IndexType begin, IndexType end, Function&& function ) const
+{
+   this->view.forElements( begin, end, function );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Function >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+forElements( IndexType begin, IndexType end, Function&& function )
+{
+   this->view.forElements( begin, end, function );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Function >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+forAllElements( Function&& function ) const
+{
+   this->forElements( 0, this->getRows(), function );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Function >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+forAllElements( Function&& function )
+{
+   this->forElements( 0, this->getRows(), function );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Function >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+forRows( IndexType begin, IndexType end, Function&& function )
+{
+   this->getView().forRows( begin, end, function );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Function >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+forRows( IndexType begin, IndexType end, Function&& function ) const
+{
+   this->getConstView().forRows( begin, end, function );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Function >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+forAllRows( Function&& function )
+{
+   this->getView().forAllRows( function );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Function >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+forAllRows( Function&& function ) const
+{
+   this->getConsView().forAllRows( function );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Function >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+sequentialForRows( IndexType begin, IndexType end, Function& function ) const
+{
+   this->view.sequentialForRows( begin, end, function );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Function >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+sequentialForRows( IndexType first, IndexType last, Function& function )
+{
+   this->view.sequentialForRows( first, last, function );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Function >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+sequentialForAllRows( Function& function ) const
+{
+   this->sequentialForRows( 0, this->getRows(), function );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Function >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+sequentialForAllRows( Function& function )
+{
+   this->sequentialForRows( 0, this->getRows(), function );
+}
+
+
+/*template< typename Real,
+          template< typename, typename, typename > class Segments,
+          typename Device,
+          typename Index,
+          typename RealAllocator,
+          typename IndexAllocator >
+template< typename Real2, template< typename, typename > class Segments2, typename Index2, typename RealAllocator2, typename IndexAllocator2 >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+addMatrix( const SparseSandboxMatrix< Real2, Segments2, Device, Index2, RealAllocator2, IndexAllocator2 >& matrix,
+           const RealType& matrixMultiplicator,
+           const RealType& thisMatrixMultiplicator )
+{
+
+}
+
+template< typename Real,
+          template< typename, typename, typename > class Segments,
+          typename Device,
+          typename Index,
+          typename RealAllocator,
+          typename IndexAllocator >
+template< typename Real2, typename Index2 >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+getTransposition( const SparseSandboxMatrix< Real2, Device, Index2 >& matrix,
+                  const RealType& matrixMultiplicator )
+{
+
+}*/
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+template< typename Vector1, typename Vector2 >
+bool
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+performSORIteration( const Vector1& b,
+                     const IndexType row,
+                     Vector2& x,
+                     const RealType& omega ) const
+{
+   return false;
+}
+
+// copy assignment
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >&
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+operator=( const SparseSandboxMatrix& matrix )
+{
+   Matrix< Real, Device, Index >::operator=( matrix );
+   this->columnIndexes = matrix.columnIndexes;
+   // SANDBOX_TODO: Replace the following line with an assignment of metadata required by you sparse matrix format.
+   this->rowPointers = matrix.rowPointers;
+   this->view = this->getView();
+   return *this;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Device_ >
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >&
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+operator=( const SparseSandboxMatrix< RealType, Device_, IndexType, MatrixType, RealAllocator, IndexAllocator >& matrix )
+{
+   Matrix< Real, Device, Index >::operator=( matrix );
+   this->columnIndexes = matrix.columnIndexes;
+   // SANDBOX_TODO: Replace the following line with an assignment of metadata required by you sparse matrix format.
+   this->rowPointers = matrix.rowPointers;
+   this->view = this->getView();
+   return *this;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Real_, typename Device_, typename Index_, ElementsOrganization Organization, typename RealAllocator_ >
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >&
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+operator=( const DenseMatrix< Real_, Device_, Index_, Organization, RealAllocator_ >& matrix )
+{
+   using RHSMatrix = DenseMatrix< Real_, Device_, Index_, Organization, RealAllocator_ >;
+   using RHSIndexType = typename RHSMatrix::IndexType;
+   using RHSRealType = typename RHSMatrix::RealType;
+   using RHSDeviceType = typename RHSMatrix::DeviceType;
+   using RHSRealAllocatorType = typename RHSMatrix::RealAllocatorType;
+
+   Containers::Vector< RHSIndexType, RHSDeviceType, RHSIndexType > rowLengths;
+   matrix.getCompressedRowLengths( rowLengths );
+   this->setLike( matrix );
+   this->setRowCapacities( rowLengths );
+   Containers::Vector< IndexType, DeviceType, IndexType > rowLocalIndexes( matrix.getRows() );
+   rowLocalIndexes = 0;
+
+   // TODO: use getConstView when it works
+   const auto matrixView = const_cast< RHSMatrix& >( matrix ).getView();
+   const IndexType paddingIndex = this->getPaddingIndex();
+   auto columns_view = this->columnIndexes.getView();
+   auto values_view = this->values.getView();
+   auto rowLocalIndexes_view = rowLocalIndexes.getView();
+   columns_view = paddingIndex;
+
+   if( std::is_same< DeviceType, RHSDeviceType >::value )
+   {
+      const auto segments_view = this->segments.getView();
+      auto f = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx, RHSIndexType columnIdx, const RHSRealType& value ) mutable {
+         if( value != 0.0 )
+         {
+            IndexType thisGlobalIdx = segments_view.getGlobalIndex( rowIdx, rowLocalIndexes_view[ rowIdx ]++ );
+            columns_view[ thisGlobalIdx ] = columnIdx;
+            if( ! isBinary() )
+               values_view[ thisGlobalIdx ] = value;
+         }
+      };
+      matrix.forAllElements( f );
+   }
+   else
+   {
+      const IndexType maxRowLength = matrix.getColumns();
+      const IndexType bufferRowsCount( 128 );
+      const size_t bufferSize = bufferRowsCount * maxRowLength;
+      Containers::Vector< RHSRealType, RHSDeviceType, RHSIndexType, RHSRealAllocatorType > matrixValuesBuffer( bufferSize );
+      Containers::Vector< RealType, DeviceType, IndexType, RealAllocatorType > thisValuesBuffer( bufferSize );
+      Containers::Vector< IndexType, DeviceType, IndexType, IndexAllocatorType > thisColumnsBuffer( bufferSize );
+      auto matrixValuesBuffer_view = matrixValuesBuffer.getView();
+      auto thisValuesBuffer_view = thisValuesBuffer.getView();
+
+      IndexType baseRow( 0 );
+      const IndexType rowsCount = this->getRows();
+      while( baseRow < rowsCount )
+      {
+         const IndexType lastRow = min( baseRow + bufferRowsCount, rowsCount );
+         thisColumnsBuffer = paddingIndex;
+
+         ////
+         // Copy matrix elements into buffer
+         auto f1 = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx, RHSIndexType columnIndex, const RHSRealType& value ) mutable {
+            const IndexType bufferIdx = ( rowIdx - baseRow ) * maxRowLength + localIdx;
+            matrixValuesBuffer_view[ bufferIdx ] = value;
+         };
+         matrix.forElements( baseRow, lastRow, f1 );
+
+         ////
+         // Copy the source matrix buffer to this matrix buffer
+         thisValuesBuffer_view = matrixValuesBuffer_view;
+
+         ////
+         // Copy matrix elements from the buffer to the matrix and ignoring
+         // zero matrix elements.
+         const IndexType matrix_columns = this->getColumns();
+         auto f2 = [=] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, IndexType& columnIndex, RealType& value ) mutable {
+            RealType inValue( 0.0 );
+            IndexType bufferIdx, column( rowLocalIndexes_view[ rowIdx ] );
+            while( inValue == 0.0 && column < matrix_columns )
+            {
+               bufferIdx = ( rowIdx - baseRow ) * maxRowLength + column++;
+               inValue = thisValuesBuffer_view[ bufferIdx ];
+            }
+            rowLocalIndexes_view[ rowIdx ] = column;
+            if( inValue == 0.0 )
+            {
+               columnIndex = paddingIndex;
+               value = 0.0;
+            }
+            else
+            {
+               columnIndex = column - 1;
+               value = inValue;
+            }
+         };
+         this->forElements( baseRow, lastRow, f2 );
+         baseRow += bufferRowsCount;
+      }
+      //std::cerr << "This matrix = " << std::endl << *this << std::endl;
+   }
+   this->view = this->getView();
+   return *this;
+
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename RHSMatrix >
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >&
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+operator=( const RHSMatrix& matrix )
+{
+   using RHSIndexType = typename RHSMatrix::IndexType;
+   using RHSRealType = typename RHSMatrix::RealType;
+   using RHSDeviceType = typename RHSMatrix::DeviceType;
+   using RHSRealAllocatorType = typename RHSMatrix::RealAllocatorType;
+
+   Containers::Vector< RHSIndexType, RHSDeviceType, RHSIndexType > rowCapacities;
+   matrix.getRowCapacities( rowCapacities );
+   this->setDimensions( matrix.getRows(), matrix.getColumns() );
+   this->setRowCapacities( rowCapacities );
+   Containers::Vector< IndexType, DeviceType, IndexType > rowLocalIndexes( matrix.getRows() );
+   rowLocalIndexes = 0;
+
+   // TODO: use getConstView when it works
+   const auto matrixView = const_cast< RHSMatrix& >( matrix ).getView();
+   const IndexType paddingIndex = this->getPaddingIndex();
+   auto columns_view = this->columnIndexes.getView();
+   auto values_view = this->values.getView();
+   auto rowLocalIndexes_view = rowLocalIndexes.getView();
+   columns_view = paddingIndex;
+
+   // SANDBOX_TODO: Modify the follwoing accoring to your format
+   auto row_pointers_view = this->rowPointers.getView();
+   if( std::is_same< DeviceType, RHSDeviceType >::value )
+   {
+      auto f = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx_, RHSIndexType columnIndex, const RHSRealType& value ) mutable {
+         IndexType localIdx( rowLocalIndexes_view[ rowIdx ] );
+         IndexType thisRowBegin = row_pointers_view[ rowIdx ];
+         if( value != 0.0 && columnIndex != paddingIndex )
+         {
+            IndexType thisGlobalIdx = thisRowBegin + localIdx++;
+            columns_view[ thisGlobalIdx ] = columnIndex;
+            if( ! isBinary() )
+               values_view[ thisGlobalIdx ] = value;
+            rowLocalIndexes_view[ rowIdx ] = localIdx;
+         }
+      };
+      matrix.forAllElements( f );
+   }
+   else
+   {
+      const IndexType maxRowLength = max( rowCapacities );
+      const IndexType bufferRowsCount( 128 );
+      const size_t bufferSize = bufferRowsCount * maxRowLength;
+      Containers::Vector< RHSRealType, RHSDeviceType, RHSIndexType, RHSRealAllocatorType > matrixValuesBuffer( bufferSize );
+      Containers::Vector< RHSIndexType, RHSDeviceType, RHSIndexType > matrixColumnsBuffer( bufferSize );
+      Containers::Vector< RealType, DeviceType, IndexType, RealAllocatorType > thisValuesBuffer( bufferSize );
+      Containers::Vector< IndexType, DeviceType, IndexType > thisColumnsBuffer( bufferSize );
+      Containers::Vector< IndexType, DeviceType, IndexType > thisRowLengths;
+      Containers::Vector< RHSIndexType, RHSDeviceType, RHSIndexType > rhsRowLengths;
+      matrix.getCompressedRowLengths( rhsRowLengths );
+      thisRowLengths= rhsRowLengths;
+      auto matrixValuesBuffer_view = matrixValuesBuffer.getView();
+      auto matrixColumnsBuffer_view = matrixColumnsBuffer.getView();
+      auto thisValuesBuffer_view = thisValuesBuffer.getView();
+      auto thisColumnsBuffer_view = thisColumnsBuffer.getView();
+      matrixValuesBuffer_view = 0.0;
+
+      IndexType baseRow( 0 );
+      const IndexType rowsCount = this->getRows();
+      while( baseRow < rowsCount )
+      {
+         const IndexType lastRow = min( baseRow + bufferRowsCount, rowsCount );
+         thisColumnsBuffer = paddingIndex;
+         matrixColumnsBuffer_view = paddingIndex;
+
+         ////
+         // Copy matrix elements into buffer
+         auto f1 = [=] __cuda_callable__ ( RHSIndexType rowIdx, RHSIndexType localIdx, RHSIndexType columnIndex, const RHSRealType& value ) mutable {
+            if( columnIndex != paddingIndex )
+            {
+               TNL_ASSERT_LT( rowIdx - baseRow, bufferRowsCount, "" );
+               TNL_ASSERT_LT( localIdx, maxRowLength, "" );
+               const IndexType bufferIdx = ( rowIdx - baseRow ) * maxRowLength + localIdx;
+               TNL_ASSERT_LT( bufferIdx, ( IndexType ) bufferSize, "" );
+               matrixColumnsBuffer_view[ bufferIdx ] = columnIndex;
+               matrixValuesBuffer_view[ bufferIdx ] = value;
+            }
+         };
+         matrix.forElements( baseRow, lastRow, f1 );
+
+         ////
+         // Copy the source matrix buffer to this matrix buffer
+         thisValuesBuffer_view = matrixValuesBuffer_view;
+         thisColumnsBuffer_view = matrixColumnsBuffer_view;
+
+         ////
+         // Copy matrix elements from the buffer to the matrix and ignoring
+         // zero matrix elements
+         //const IndexType matrix_columns = this->getColumns();
+         const auto thisRowLengths_view = thisRowLengths.getConstView();
+         auto f2 = [=] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, IndexType& columnIndex, RealType& value ) mutable {
+            RealType inValue( 0.0 );
+            size_t bufferIdx;
+            IndexType bufferLocalIdx( rowLocalIndexes_view[ rowIdx ] );
+            while( inValue == 0.0 && localIdx < thisRowLengths_view[ rowIdx ] )
+            {
+               bufferIdx = ( rowIdx - baseRow ) * maxRowLength + bufferLocalIdx++;
+               TNL_ASSERT_LT( bufferIdx, bufferSize, "" );
+               inValue = thisValuesBuffer_view[ bufferIdx ];
+            }
+            rowLocalIndexes_view[ rowIdx ] = bufferLocalIdx;
+            if( inValue == 0.0 )
+            {
+               columnIndex = paddingIndex;
+               value = 0.0;
+            }
+            else
+            {
+               columnIndex = thisColumnsBuffer_view[ bufferIdx ];//column - 1;
+               value = inValue;
+            }
+         };
+         this->forElements( baseRow, lastRow, f2 );
+         baseRow += bufferRowsCount;
+      }
+   }
+   this->view = this->getView();
+   return *this;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Matrix >
+bool
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+operator==( const Matrix& m ) const
+{
+   return view == m;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Matrix >
+bool
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+operator!=( const Matrix& m ) const
+{
+   return view != m;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+save( File& file ) const
+{
+   this->view.save( file );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+load( File& file )
+{
+   Matrix< RealType, DeviceType, IndexType >::load( file );
+   file >> this->columnIndexes;
+   // SANDBOX_TODO: Replace the following line with loading of metadata required by your sparse matrix format.
+   file >> rowPointers;
+   this->view = this->getView();
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+save( const String& fileName ) const
+{
+   Object::save( fileName );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+load( const String& fileName )
+{
+   Object::load( fileName );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+void
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+print( std::ostream& str ) const
+{
+   this->view.print( str );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+__cuda_callable__
+Index
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+getPaddingIndex() const
+{
+   return -1;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+auto
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+getColumnIndexes() const -> const ColumnsIndexesVectorType&
+{
+   return this->columnIndexes;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType,
+          typename RealAllocator,
+          typename IndexAllocator >
+auto
+SparseSandboxMatrix< Real, Device, Index, MatrixType, RealAllocator, IndexAllocator >::
+getColumnIndexes() -> ColumnsIndexesVectorType&
+{
+   return this->columnIndexes;
+}
+
+      } // namespace Sandbox
+   } // namespace Matrices
+} // namespace TNL
diff --git a/src/TNL/Matrices/Sandbox/SparseSandboxMatrixRowView.h b/src/TNL/Matrices/Sandbox/SparseSandboxMatrixRowView.h
new file mode 100644
index 0000000000..cabf7b7fd5
--- /dev/null
+++ b/src/TNL/Matrices/Sandbox/SparseSandboxMatrixRowView.h
@@ -0,0 +1,282 @@
+ /***************************************************************************
+                          SparseSandboxMatrixRowView.h -  description
+                             -------------------
+    begin                : Apr 20, 2021
+    copyright            : (C) 2021 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#pragma once
+
+#include <ostream>
+
+#include <TNL/Cuda/CudaCallable.h>
+#include <TNL/Matrices/MatrixRowViewIterator.h>
+
+namespace TNL {
+   namespace Matrices {
+      namespace Sandbox {
+
+/**
+ * \brief RowView is a simple structure for accessing rows of sparse matrix.
+ *
+ * \tparam ValuesView is a vector view storing the matrix elements values.
+ * \tparam ColumnsIndexesView is a vector view storing the column indexes of the matrix element.
+ * \tparam isBinary tells if the the parent matrix is a binary matrix.
+ *
+ * See \ref SparseSandboxMatrix and \ref SparseSandboxMatrixView.
+ *
+ * \par Example
+ * \include Matrices/SparseMatrix/SparseMatrixExample_getRow.cpp
+ * \par Output
+ * \include SparseMatrixExample_getRow.out
+ *
+ * \par Example
+ * \include Matrices/SparseMatrix/SparseMatrixViewExample_getRow.cpp
+ * \par Output
+ * \include SparseMatrixViewExample_getRow.out
+ */
+template< typename ValuesView,
+          typename ColumnsIndexesView,
+          bool isBinary_ >
+class SparseSandboxMatrixRowView
+{
+   public:
+
+      /**
+       * \brief The type of matrix elements.
+       */
+      using RealType = typename ValuesView::RealType;
+
+      /**
+       * \brief The type used for matrix elements indexing.
+       */
+      using IndexType = typename ColumnsIndexesView::IndexType;
+
+      /**
+       * \brief Type of container view used for storing the matrix elements values.
+       */
+      using ValuesViewType = ValuesView;
+
+      /**
+       * \brief Type of container view used for storing the column indexes of the matrix elements.
+       */
+      using ColumnsIndexesViewType = ColumnsIndexesView;
+
+      /**
+       * \brief Type of constant container view used for storing the matrix elements values.
+       */
+      using ConstValuesViewType = typename ValuesViewType::ConstViewType;
+
+      /**
+       * \brief Type of constant container view used for storing the column indexes of the matrix elements.
+       */
+      using ConstColumnsIndexesViewType = typename ColumnsIndexesViewType::ConstViewType;
+
+      /**
+       * \brief Type of sparse matrix row view.
+       */
+      using RowView = SparseSandboxMatrixRowView< ValuesViewType, ColumnsIndexesViewType, isBinary_ >;
+
+      /**
+       * \brief Type of constant sparse matrix row view.
+       */
+      using ConstView = SparseSandboxMatrixRowView< ConstValuesViewType, ConstColumnsIndexesViewType, isBinary_ >;
+
+      /**
+       * \brief The type of related matrix element.
+       */
+      using MatrixElementType = SparseMatrixElement< RealType, IndexType >;
+
+      /**
+       * \brief Type of iterator for the matrix row.
+       */
+      using IteratorType = MatrixRowViewIterator< RowView >;
+
+      /**
+       * \brief Tells whether the parent matrix is a binary matrix.
+       * @return `true` if the matrix is binary.
+       */
+      static constexpr bool isBinary() { return isBinary_; };
+
+      /**
+       * \brief Constructor with \e segmentView, \e values and \e columnIndexes.
+       *
+       * \param rowIdx is row index.
+       * \param offset is the begining of the matrix row in arrays with values and column indexes of matrix elements.
+       * \param size is row size, i.e. number of nonzero matrix elements in the row.
+       * \param values is a container view for storing the matrix elements values.
+       * \param columnIndexes is a container view for storing the column indexes of the matrix elements.
+       */
+      __cuda_callable__
+      SparseSandboxMatrixRowView( IndexType rowIdx,
+                                  IndexType offset,
+                                  IndexType size,
+                                  const ValuesViewType& values,
+                                  const ColumnsIndexesViewType& columnIndexes );
+
+      /**
+       * \brief Returns size of the matrix row, i.e. number of matrix elements in this row.
+       *
+       * \return Size of the matrix row.
+       */
+      __cuda_callable__
+      IndexType getSize() const;
+
+      /**
+       * \brief Returns the matrix row index.
+       *
+       * \return matrix row index.
+       */
+      __cuda_callable__
+      const IndexType& getRowIndex() const;
+
+      /**
+       * \brief Returns constants reference to a column index of an element with given rank in the row.
+       *
+       * \param localIdx is the rank of the non-zero element in given row.
+       *
+       * \return constant reference to the matrix element column index.
+       */
+      __cuda_callable__
+      const IndexType& getColumnIndex( const IndexType localIdx ) const;
+
+      /**
+       * \brief Returns non-constants reference to a column index of an element with given rank in the row.
+       *
+       * \param localIdx is the rank of the non-zero element in given row.
+       *
+       * \return non-constant reference to the matrix element column index.
+       */
+      __cuda_callable__
+      IndexType& getColumnIndex( const IndexType localIdx );
+
+      /**
+       * \brief Returns constants reference to value of an element with given rank in the row.
+       *
+       * \param localIdx is the rank of the non-zero element in given row.
+       *
+       * \return constant reference to the matrix element value.
+       */
+      __cuda_callable__
+      const RealType& getValue( const IndexType localIdx ) const;
+
+      /**
+       * \brief Returns non-constants reference to value of an element with given rank in the row.
+       *
+       * \param localIdx is the rank of the non-zero element in given row.
+       *
+       * \return non-constant reference to the matrix element value.
+       */
+      __cuda_callable__
+      RealType& getValue( const IndexType localIdx );
+
+      /**
+       * \brief Sets a value of matrix element with given rank in the matrix row.
+       *
+       * \param localIdx is the rank of the matrix element in the row.
+       * \param value is the new value of the matrix element.
+       */
+      __cuda_callable__
+      void setValue( const IndexType localIdx,
+                     const RealType& value );
+
+      /**
+       * \brief Sets a column index of matrix element with given rank in the matrix row.
+       *
+       * \param localIdx is the rank of the matrix element in the row.
+       * \param columnIndex is the new column index of the matrix element.
+       */
+      __cuda_callable__
+      void setColumnIndex( const IndexType localIdx,
+                           const IndexType& columnIndex );
+
+      /**
+       * \brief Sets both a value and a column index of matrix element with given rank in the matrix row.
+       *
+       * \param localIdx is the rank of the matrix element in the row.
+       * \param columnIndex is the new column index of the matrix element.
+       * \param value is the new value of the matrix element.
+       */
+      __cuda_callable__
+      void setElement( const IndexType localIdx,
+                       const IndexType columnIndex,
+                       const RealType& value );
+
+      /**
+       * \brief Comparison of two matrix rows.
+       *
+       * The other matrix row can be from any other matrix.
+       *
+       * \param other is another matrix row.
+       * \return \e true if both rows are the same, \e false otherwise.
+       */
+      template< typename _ValuesView,
+                typename _ColumnsIndexesView,
+                bool _isBinary >
+      __cuda_callable__
+      bool operator==( const SparseSandboxMatrixRowView< _ValuesView, _ColumnsIndexesView, _isBinary >& other ) const;
+
+      /**
+       * \brief Returns iterator pointing at the beginning of the matrix row.
+       *
+       * \return iterator pointing at the beginning.
+       */
+      __cuda_callable__
+      IteratorType begin();
+
+      /**
+       * \brief Returns iterator pointing at the end of the matrix row.
+       *
+       * \return iterator pointing at the end.
+       */
+      __cuda_callable__
+      IteratorType end();
+
+      /**
+       * \brief Returns constant iterator pointing at the beginning of the matrix row.
+       *
+       * \return iterator pointing at the beginning.
+       */
+      __cuda_callable__
+      const IteratorType cbegin() const;
+
+      /**
+       * \brief Returns constant iterator pointing at the end of the matrix row.
+       *
+       * \return iterator pointing at the end.
+       */
+      __cuda_callable__
+      const IteratorType cend() const;
+
+   protected:
+
+      IndexType rowIdx, size;
+
+      // SANDBOX_TODO: Replace the following line with data required by your format.
+      IndexType offset;
+
+      ValuesViewType values;
+
+      ColumnsIndexesViewType columnIndexes;
+};
+
+/**
+ * \brief Insertion operator for a sparse matrix row.
+ *
+ * \param str is an output stream.
+ * \param row is an input sparse matrix row.
+ * \return  reference to the output stream.
+ */
+template< typename ValuesView,
+          typename ColumnsIndexesView,
+          bool isBinary_ >
+std::ostream& operator<<( std::ostream& str, const SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >& row );
+
+      } // namespace Sandbox
+   } // namespace Matrices
+} // namespace TNL
+
+#include <TNL/Matrices/Sandbox/SparseSandboxMatrixRowView.hpp>
diff --git a/src/TNL/Matrices/Sandbox/SparseSandboxMatrixRowView.hpp b/src/TNL/Matrices/Sandbox/SparseSandboxMatrixRowView.hpp
new file mode 100644
index 0000000000..09598edd04
--- /dev/null
+++ b/src/TNL/Matrices/Sandbox/SparseSandboxMatrixRowView.hpp
@@ -0,0 +1,240 @@
+/***************************************************************************
+                          SparseSandboxMatrixRowView.hpp -  description
+                             -------------------
+    begin                : Apr 20, 2021
+    copyright            : (C) 2021 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#pragma once
+
+#include <TNL/Matrices/Sandbox/SparseSandboxMatrixRowView.h>
+#include <TNL/Assert.h>
+
+namespace TNL {
+   namespace Matrices {
+      namespace Sandbox {
+
+// SANDBOX_TODO: Modify the follwing constructor by your needs
+template< typename ValuesView,
+          typename ColumnsIndexesView,
+          bool isBinary_ >
+__cuda_callable__
+SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >::
+SparseSandboxMatrixRowView( IndexType rowIdx,
+                            IndexType offset,
+                            IndexType size,
+                            const ValuesViewType& values,
+                            const ColumnsIndexesViewType& columnIndexes )
+ : rowIdx( rowIdx ), offset( offset ), size( size ), values( values ), columnIndexes( columnIndexes )
+{
+}
+
+template< typename ValuesView,
+          typename ColumnsIndexesView,
+          bool isBinary_ >
+__cuda_callable__ auto
+SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >::
+getSize() const -> IndexType
+{
+   return this->size;
+}
+
+template< typename ValuesView,
+          typename ColumnsIndexesView,
+          bool isBinary_ >
+__cuda_callable__
+auto
+SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >::
+getRowIndex() const -> const IndexType&
+{
+   return this->rowIdx;
+}
+
+template< typename ValuesView,
+          typename ColumnsIndexesView,
+          bool isBinary_ >
+__cuda_callable__ auto
+SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >::
+getColumnIndex( const IndexType localIdx ) const -> const IndexType&
+{
+   TNL_ASSERT_LT( localIdx, this->getSize(), "Local index exceeds matrix row capacity." );
+   // SANDBOX_TODO: Modify the following line to match with your sparse format.
+   return columnIndexes[ offset + localIdx ];
+}
+
+template< typename ValuesView,
+          typename ColumnsIndexesView,
+          bool isBinary_ >
+__cuda_callable__ auto
+SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >::
+getColumnIndex( const IndexType localIdx ) -> IndexType&
+{
+   TNL_ASSERT_LT( localIdx, this->getSize(), "Local index exceeds matrix row capacity." );
+   // SANDBOX_TODO: Modify the following line to match with your sparse format.
+   return columnIndexes[ offset + localIdx ];
+}
+
+template< typename ValuesView,
+          typename ColumnsIndexesView,
+          bool isBinary_ >
+__cuda_callable__ auto
+SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >::
+getValue( const IndexType localIdx ) const -> const RealType&
+{
+   TNL_ASSERT_LT( localIdx, this->getSize(), "Local index exceeds matrix row capacity." );
+   TNL_ASSERT_FALSE( isBinary(), "Cannot call this method for binary matrix row." );
+   // SANDBOX_TODO: Modify the following line to match with your sparse format.
+   return values[ offset + localIdx ];
+}
+
+template< typename ValuesView,
+          typename ColumnsIndexesView,
+          bool isBinary_ >
+__cuda_callable__ auto
+SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >::
+getValue( const IndexType localIdx ) -> RealType&
+{
+   TNL_ASSERT_LT( localIdx, this->getSize(), "Local index exceeds matrix row capacity." );
+   TNL_ASSERT_FALSE( isBinary(), "Cannot call this method for binary matrix row." );
+   // SANDBOX_TODO: Modify the following line to match with your sparse format.
+   return values[ offset + localIdx ];
+}
+
+template< typename ValuesView,
+          typename ColumnsIndexesView,
+          bool isBinary_ >
+__cuda_callable__ void
+SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >::
+setValue( const IndexType localIdx,
+          const RealType& value )
+{
+   TNL_ASSERT_LT( localIdx, this->getSize(), "Local index exceeds matrix row capacity." );
+   if( ! isBinary() ) {
+      // SANDBOX_TODO: Modify the following line to match with your sparse format.
+      const IndexType globalIdx = offset + localIdx;
+      values[ globalIdx ] = value;
+   }
+}
+
+template< typename ValuesView,
+          typename ColumnsIndexesView,
+          bool isBinary_ >
+__cuda_callable__ void
+SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >::
+setColumnIndex( const IndexType localIdx,
+                const IndexType& columnIndex )
+{
+   TNL_ASSERT_LT( localIdx, this->getSize(), "Local index exceeds matrix row capacity." );
+   // SANDBOX_TODO: Modify the following line to match with your sparse format.
+   const IndexType globalIdx = offset + localIdx;
+   this->columnIndexes[ globalIdx ] = columnIndex;
+}
+
+template< typename ValuesView,
+          typename ColumnsIndexesView,
+          bool isBinary_ >
+__cuda_callable__ void
+SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >::
+setElement( const IndexType localIdx,
+            const IndexType column,
+            const RealType& value )
+{
+   TNL_ASSERT_LT( localIdx, this->getSize(), "Local index exceeds matrix row capacity." );
+   // SANDBOX_TODO: Modify the following line to match with your sparse format.
+   const IndexType globalIdx = offset + localIdx;
+   columnIndexes[ globalIdx ] = column;
+   if( ! isBinary() )
+      values[ globalIdx ] = value;
+}
+
+template< typename ValuesView,
+          typename ColumnsIndexesView,
+          bool isBinary_ >
+   template< typename _ValuesView,
+             typename _ColumnsIndexesView,
+             bool _isBinary >
+__cuda_callable__
+bool
+SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >::
+operator==( const SparseSandboxMatrixRowView< _ValuesView, _ColumnsIndexesView, _isBinary >& other ) const
+{
+   IndexType i = 0;
+   while( i < getSize() && i < other.getSize() ) {
+      if( getColumnIndex( i ) != other.getColumnIndex( i ) )
+         return false;
+      if( ! _isBinary && getValue( i ) != other.getValue( i ) )
+         return false;
+      ++i;
+   }
+   for( IndexType j = i; j < getSize(); j++ )
+      // TODO: use ... != getPaddingIndex()
+      if( getColumnIndex( j ) >= 0 )
+         return false;
+   for( IndexType j = i; j < other.getSize(); j++ )
+      // TODO: use ... != getPaddingIndex()
+      if( other.getColumnIndex( j ) >= 0 )
+         return false;
+   return true;
+}
+
+template< typename ValuesView,
+          typename ColumnsIndexesView,
+          bool isBinary_ >
+__cuda_callable__ auto
+SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >::
+begin() -> IteratorType
+{
+   return IteratorType( *this, 0 );
+}
+
+template< typename ValuesView,
+          typename ColumnsIndexesView,
+          bool isBinary_ >
+__cuda_callable__ auto
+SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >::
+end() -> IteratorType
+{
+   return IteratorType( *this, this->getSize() );
+}
+
+template< typename ValuesView,
+          typename ColumnsIndexesView,
+          bool isBinary_ >
+__cuda_callable__ auto
+SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >::
+cbegin() const -> const IteratorType
+{
+   return IteratorType( *this, 0 );
+}
+
+template< typename ValuesView,
+          typename ColumnsIndexesView,
+          bool isBinary_ >
+__cuda_callable__ auto
+SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >::
+cend() const -> const IteratorType
+{
+   return IteratorType( *this, this->getSize() );
+}
+
+template< typename ValuesView,
+          typename ColumnsIndexesView,
+          bool isBinary_ >
+std::ostream& operator<<( std::ostream& str, const SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >& row )
+{
+   using NonConstIndex = std::remove_const_t< typename SparseSandboxMatrixRowView< ValuesView, ColumnsIndexesView, isBinary_ >::IndexType >;
+   for( NonConstIndex i = 0; i < row.getSize(); i++ )
+      if( isBinary_ )
+         // TODO: check getPaddingIndex(), print only the column indices of non-zeros but not the values
+         str << " [ " << row.getColumnIndex( i ) << " ] = " << (row.getColumnIndex( i ) >= 0) << ", ";
+      else
+         str << " [ " << row.getColumnIndex( i ) << " ] = " << row.getValue( i ) << ", ";
+   return str;
+}
+
+      } // namespace Sandbox
+   } // namespace Matrices
+} // namespace TNL
diff --git a/src/TNL/Matrices/Sandbox/SparseSandboxMatrixView.h b/src/TNL/Matrices/Sandbox/SparseSandboxMatrixView.h
new file mode 100644
index 0000000000..66247a349d
--- /dev/null
+++ b/src/TNL/Matrices/Sandbox/SparseSandboxMatrixView.h
@@ -0,0 +1,871 @@
+/***************************************************************************
+                          SparseSandboxMatrixView.h -  description
+                             -------------------
+    begin                : Apr 20, 2021
+    copyright            : (C) 2021 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#pragma once
+
+#include <TNL/Matrices/Matrix.h>
+#include <TNL/Matrices/MatrixType.h>
+#include <TNL/Allocators/Default.h>
+#include <TNL/Algorithms/Segments/CSR.h>
+#include <TNL/Matrices/Sandbox/SparseSandboxMatrixRowView.h>
+#include <TNL/TypeTraits.h>
+
+namespace TNL {
+   namespace Matrices {
+      namespace Sandbox {
+
+/**
+ * \brief Implementation of sparse sandbox matrix view.
+ *
+ * It serves as an accessor to \ref SparseSandboxMatrix for example when passing the
+ * matrix to lambda functions. SparseSandboxMatrix view can be also created in CUDA kernels.
+ *
+ * \tparam Real is a type of matrix elements. If \e Real equals \e bool the matrix is treated
+ *    as binary and so the matrix elements values are not stored in the memory since we need
+ *    to remember only coordinates of non-zero elements( which equal one).
+ * \tparam Device is a device where the matrix is allocated.
+ * \tparam Index is a type for indexing of the matrix elements.
+ * \tparam MatrixType specifies a symmetry of matrix. See \ref MatrixType. Symmetric
+ *    matrices store only lower part of the matrix and its diagonal. The upper part is reconstructed on the fly.
+ *    GeneralMatrix with no symmetry is used by default.
+ * \tparam Segments is a structure representing the sparse matrix format. Depending on the pattern of the non-zero elements
+ *    different matrix formats can perform differently especially on GPUs. By default \ref CSR format is used. See also
+ *    \ref Ellpack, \ref SlicedEllpack, \ref ChunkedEllpack or \ref BiEllpack.
+ * \tparam ComputeReal is the same as \e Real mostly but for binary matrices it is set to \e Index type. This can be changed
+ *    bu the user, of course.
+ *
+ */
+template< typename Real,
+          typename Device = Devices::Host,
+          typename Index = int,
+          typename MatrixType = GeneralMatrix >
+class SparseSandboxMatrixView : public MatrixView< Real, Device, Index >
+{
+   static_assert(
+      ! MatrixType::isSymmetric() ||
+      ! std::is_same< Device, Devices::Cuda >::value ||
+      ( std::is_same< Real, float >::value || std::is_same< Real, double >::value || std::is_same< Real, int >::value || std::is_same< Real, long long int >::value ),
+      "Given Real type is not supported by atomic operations on GPU which are necessary for symmetric operations." );
+
+   public:
+
+      // Supporting types - they are not important for the user
+      using BaseType = MatrixView< Real, Device, Index >;
+      using ValuesViewType = typename BaseType::ValuesView;
+      using ConstValuesViewType = typename ValuesViewType::ConstViewType;
+      using ColumnsIndexesViewType = Containers::VectorView< typename TNL::copy_const< Index >::template from< Real >::type, Device, Index >;
+      using ConstColumnsIndexesViewType = typename ColumnsIndexesViewType::ConstViewType;
+      using RowsCapacitiesView = Containers::VectorView< Index, Device, Index >;
+      using ConstRowsCapacitiesView = typename RowsCapacitiesView::ConstViewType;
+
+      /**
+       * \brief Test of symmetric matrix type.
+       *
+       * \return \e true if the matrix is stored as symmetric and \e false otherwise.
+       */
+      static constexpr bool isSymmetric() { return MatrixType::isSymmetric(); };
+
+      /**
+       * \brief Test of binary matrix type.
+       *
+       * \return \e true if the matrix is stored as binary and \e false otherwise.
+       */
+      static constexpr bool isBinary() { return std::is_same< Real, bool >::value; };
+
+      /**
+       * \brief The type of matrix elements.
+       */
+      using RealType = Real;
+
+      //using ComputeRealType = ComputeReal;
+
+      /**
+       * \brief The device where the matrix is allocated.
+       */
+      using DeviceType = Device;
+
+      /**
+       * \brief The type used for matrix elements indexing.
+       */
+      using IndexType = Index;
+
+      /**
+       * \brief Templated type of segments view, i.e. sparse matrix format.
+       */
+      //template< typename Device_, typename Index_ >
+      //using SegmentsViewTemplate = SegmentsView< Device_, Index_ >;
+
+      /**
+       * \brief Type of segments view used by this matrix. It represents the sparse matrix format.
+       */
+      //using SegmentsViewType = SegmentsView< Device, Index >;
+
+      /**
+       * \brief Type of related matrix view.
+       */
+      using ViewType = SparseSandboxMatrixView< Real, Device, Index, MatrixType >;
+
+      /**
+       * \brief Matrix view type for constant instances.
+       */
+      using ConstViewType = SparseSandboxMatrixView< std::add_const_t< Real >, Device, Index, MatrixType >;
+
+      /**
+       * \brief Type for accessing matrix rows.
+       */
+      using RowView = SparseSandboxMatrixRowView< ValuesViewType, ColumnsIndexesViewType, isBinary() >;
+
+      /**
+       * \brief Type for accessing constant matrix rows.
+       */
+      using ConstRowView = SparseSandboxMatrixRowView< ConstValuesViewType, ConstColumnsIndexesViewType, isBinary() >;;
+
+      /**
+       * \brief Helper type for getting self type or its modifications.
+       */
+      template< typename _Real = Real,
+                typename _Device = Device,
+                typename _Index = Index,
+                typename _MatrixType = MatrixType >
+      using Self = SparseSandboxMatrixView< _Real, _Device, _Index, _MatrixType >;
+
+      /**
+       * \brief Type of container view for CSR row pointers.
+       *
+       * SANDBOX_TODO: You may replace it with containers views for metadata of your format.
+       */
+      using RowPointersView = TNL::Containers::VectorView< IndexType, DeviceType, IndexType >;
+
+      /**
+       * \brief Constructor with no parameters.
+       */
+      __cuda_callable__
+      SparseSandboxMatrixView();
+
+      /**
+       * \brief Constructor with all necessary data and views.
+       *
+       * \param rows is a number of matrix rows.
+       * \param columns is a number of matrix columns.
+       * \param values is a vector view with matrix elements values.
+       * \param columnIndexes is a vector view with matrix elements column indexes.
+       * \param rowPointers is a container view with row pointers.
+       *
+       * SANDBOX_TODO: Replace `rowPointers` with metadata by your needs.
+       */
+      __cuda_callable__
+      SparseSandboxMatrixView( const IndexType rows,
+                               const IndexType columns,
+                               const ValuesViewType& values,
+                               const ColumnsIndexesViewType& columnIndexes,
+                               const RowPointersView& rowPointers );
+
+      /**
+       * \brief Copy constructor.
+       *
+       * \param matrix is an input sparse matrix view.
+       */
+      __cuda_callable__
+      SparseSandboxMatrixView( const SparseSandboxMatrixView& matrix ) = default;
+
+      /**
+       * \brief Move constructor.
+       *
+       * \param matrix is an input sparse matrix view.
+       */
+      __cuda_callable__
+      SparseSandboxMatrixView( SparseSandboxMatrixView&& matrix ) = default;
+
+      /**
+       * \brief Returns a modifiable view of the sparse matrix.
+       *
+       * \return sparse matrix view.
+       */
+      __cuda_callable__
+      ViewType getView();
+
+      /**
+       * \brief Returns a non-modifiable view of the sparse matrix.
+       *
+       * \return sparse matrix view.
+       */
+      __cuda_callable__
+      ConstViewType getConstView() const;
+
+      /**
+       * \brief Returns string with serialization type.
+       *
+       * The string has a form `Matrices::SparseSandboxMatrix< RealType,  [any_device], IndexType, General/Symmetric, Format, [any_allocator] >`.
+       *
+       * \return \ref String with the serialization type.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixViewExample_getSerializationType.cpp
+       * \par Output
+       * \include SparseMatrixViewExample_getSerializationType.out
+       */
+      static String getSerializationType();
+
+      /**
+       * \brief Returns string with serialization type.
+       *
+       * See \ref SparseSandboxMatrix::getSerializationType.
+       *
+       * \return \e String with the serialization type.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixExample_getSerializationType.cpp
+       * \par Output
+       * \include SparseMatrixExample_getSerializationType.out
+       */
+      virtual String getSerializationTypeVirtual() const;
+
+      /**
+       * \brief Computes number of non-zeros in each row.
+       *
+       * \param rowLengths is a vector into which the number of non-zeros in each row
+       * will be stored.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixViewExample_getCompressedRowLengths.cpp
+       * \par Output
+       * \include SparseMatrixViewExample_getCompressedRowLengths.out
+       */
+      template< typename Vector >
+      void getCompressedRowLengths( Vector& rowLengths ) const;
+
+      /**
+       * \brief Compute capacities of all rows.
+       *
+       * The row capacities are not stored explicitly and must be computed.
+       *
+       * \param rowCapacities is a vector where the row capacities will be stored.
+       */
+      template< typename Vector >
+      void getRowCapacities( Vector& rowCapacities ) const;
+
+      /**
+       * \brief Returns capacity of given matrix row.
+       *
+       * \param row index of matrix row.
+       * \return number of matrix elements allocated for the row.
+       */
+      __cuda_callable__
+      IndexType getRowCapacity( const IndexType row ) const;
+
+      /**
+       * \brief Returns number of non-zero matrix elements.
+       *
+       * This method really counts the non-zero matrix elements and so
+       * it returns zero for matrix having all allocated elements set to zero.
+       *
+       * \return number of non-zero matrix elements.
+       */
+      IndexType getNonzeroElementsCount() const;
+
+      /**
+       * \brief Constant getter of simple structure for accessing given matrix row.
+       *
+       * \param rowIdx is matrix row index.
+       *
+       * \return RowView for accessing given matrix row.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixViewExample_getConstRow.cpp
+       * \par Output
+       * \include SparseMatrixViewExample_getConstRow.out
+       *
+       * See \ref SparseSandboxMatrixRowView.
+       */
+      __cuda_callable__
+      ConstRowView getRow( const IndexType& rowIdx ) const;
+
+      /**
+       * \brief Non-constant getter of simple structure for accessing given matrix row.
+       *
+       * \param rowIdx is matrix row index.
+       *
+       * \return RowView for accessing given matrix row.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixViewExample_getRow.cpp
+       * \par Output
+       * \include SparseMatrixViewExample_getRow.out
+       *
+       * See \ref SparseSandboxMatrixRowView.
+       */
+      __cuda_callable__
+      RowView getRow( const IndexType& rowIdx );
+
+      /**
+       * \brief Sets element at given \e row and \e column to given \e value.
+       *
+       * This method can be called from the host system (CPU) no matter
+       * where the matrix is allocated. If the matrix is allocated on GPU this method
+       * can be called even from device kernels. If the matrix is allocated in GPU device
+       * this method is called from CPU, it transfers values of each matrix element separately and so the
+       * performance is very low. For higher performance see. \ref SparseSandboxMatrix::getRow
+       * or \ref SparseSandboxMatrix::forElements and \ref SparseSandboxMatrix::forAllElements.
+       * The call may fail if the matrix row capacity is exhausted.
+       *
+       * \param row is row index of the element.
+       * \param column is columns index of the element.
+       * \param value is the value the element will be set to.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixViewExample_setElement.cpp
+       * \par Output
+       * \include SparseMatrixViewExample_setElement.out
+       */
+      __cuda_callable__
+      void setElement( const IndexType row,
+                       const IndexType column,
+                       const RealType& value );
+
+      /**
+       * \brief Add element at given \e row and \e column to given \e value.
+       *
+       * This method can be called from the host system (CPU) no matter
+       * where the matrix is allocated. If the matrix is allocated on GPU this method
+       * can be called even from device kernels. If the matrix is allocated in GPU device
+       * this method is called from CPU, it transfers values of each matrix element separately and so the
+       * performance is very low. For higher performance see. \ref SparseSandboxMatrix::getRow
+       * or \ref SparseSandboxMatrix::forElements and \ref SparseSandboxMatrix::forAllElements.
+       * The call may fail if the matrix row capacity is exhausted.
+       *
+       * \param row is row index of the element.
+       * \param column is columns index of the element.
+       * \param value is the value the element will be set to.
+       * \param thisElementMultiplicator is multiplicator the original matrix element
+       *   value is multiplied by before addition of given \e value.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixViewExample_addElement.cpp
+       * \par Output
+       * \include SparseMatrixViewExample_addElement.out
+       */
+      __cuda_callable__
+      void addElement( IndexType row,
+                       IndexType column,
+                       const RealType& value,
+                       const RealType& thisElementMultiplicator = 1.0 );
+
+      /**
+       * \brief Returns value of matrix element at position given by its row and column index.
+       *
+       * This method can be called from the host system (CPU) no matter
+       * where the matrix is allocated. If the matrix is allocated on GPU this method
+       * can be called even from device kernels. If the matrix is allocated in GPU device
+       * this method is called from CPU, it transfers values of each matrix element separately and so the
+       * performance is very low. For higher performance see. \ref SparseSandboxMatrix::getRow
+       * or \ref SparseSandboxMatrix::forElements and \ref SparseSandboxMatrix::forAllElements.
+       *
+       * \param row is a row index of the matrix element.
+       * \param column i a column index of the matrix element.
+       *
+       * \return value of given matrix element.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixViewExample_getElement.cpp
+       * \par Output
+       * \include SparseMatrixViewExample_getElement.out
+       *
+       */
+      __cuda_callable__
+      RealType getElement( IndexType row,
+                           IndexType column ) const;
+
+      /**
+       * \brief Method for performing general reduction on matrix rows.
+       *
+       * \tparam Fetch is a type of lambda function for data fetch declared as
+       *          `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`.
+       *          The return type of this lambda can be any non void.
+       * \tparam Reduce is a type of lambda function for reduction declared as
+       *          `reduce( const FetchValue& v1, const FetchValue& v2 ) -> FetchValue`.
+       * \tparam Keep is a type of lambda function for storing results of reduction in each row.
+       *          It is declared as `keep( const IndexType rowIdx, const double& value )`.
+       * \tparam FetchValue is type returned by the Fetch lambda function.
+       *
+       * \param begin defines beginning of the range [begin,end) of rows to be processed.
+       * \param end defines ending of the range [begin,end) of rows to be processed.
+       * \param fetch is an instance of lambda function for data fetch.
+       * \param reduce is an instance of lambda function for reduction.
+       * \param keep in an instance of lambda function for storing results.
+       * \param zero is zero of given reduction operation also known as idempotent element.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixViewExample_reduceRows.cpp
+       * \par Output
+       * \include SparseMatrixViewExample_reduceRows.out
+       */
+      template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >
+      void reduceRows( IndexType begin, IndexType end, Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchReal& zero );
+
+      /**
+       * \brief Method for performing general reduction on matrix rows for constant instances.
+       *
+       * \tparam Fetch is a type of lambda function for data fetch declared as
+       *          `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`.
+       *          The return type of this lambda can be any non void.
+       * \tparam Reduce is a type of lambda function for reduction declared as
+       *          `reduce( const FetchValue& v1, const FetchValue& v2 ) -> FetchValue`.
+       * \tparam Keep is a type of lambda function for storing results of reduction in each row.
+       *          It is declared as `keep( const IndexType rowIdx, const double& value )`.
+       * \tparam FetchValue is type returned by the Fetch lambda function.
+       *
+       * \param begin defines beginning of the range [begin,end) of rows to be processed.
+       * \param end defines ending of the range [begin,end) of rows to be processed.
+       * \param fetch is an instance of lambda function for data fetch.
+       * \param reduce is an instance of lambda function for reduction.
+       * \param keep in an instance of lambda function for storing results.
+       * \param zero is zero of given reduction operation also known as idempotent element.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixViewExample_reduceRows.cpp
+       * \par Output
+       * \include SparseMatrixViewExample_reduceRows.out
+       */
+      template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >
+      void reduceRows( IndexType begin, IndexType end, Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchReal& zero ) const;
+
+      /**
+       * \brief Method for performing general reduction on all matrix rows.
+       *
+       * \tparam Fetch is a type of lambda function for data fetch declared as
+       *          `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`.
+       *          The return type of this lambda can be any non void.
+       * \tparam Reduce is a type of lambda function for reduction declared as
+       *          `reduce( const FetchValue& v1, const FetchValue& v2 ) -> FetchValue`.
+       * \tparam Keep is a type of lambda function for storing results of reduction in each row.
+       *          It is declared as `keep( const IndexType rowIdx, const double& value )`.
+       * \tparam FetchValue is type returned by the Fetch lambda function.
+       *
+       * \param fetch is an instance of lambda function for data fetch.
+       * \param reduce is an instance of lambda function for reduction.
+       * \param keep in an instance of lambda function for storing results.
+       * \param zero is zero of given reduction operation also known as idempotent element.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixViewExample_reduceAllRows.cpp
+       * \par Output
+       * \include SparseMatrixViewExample_reduceAllRows.out
+       */
+      template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >
+      void reduceAllRows( Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchReal& zero );
+
+      /**
+       * \brief Method for performing general reduction on all matrix rows for constant instances.
+       *
+       * \tparam Fetch is a type of lambda function for data fetch declared as
+       *          `fetch( IndexType rowIdx, IndexType& columnIdx, RealType& elementValue ) -> FetchValue`.
+       *          The return type of this lambda can be any non void.
+       * \tparam Reduce is a type of lambda function for reduction declared as
+       *          `reduce( const FetchValue& v1, const FetchValue& v2 ) -> FetchValue`.
+       * \tparam Keep is a type of lambda function for storing results of reduction in each row.
+       *          It is declared as `keep( const IndexType rowIdx, const double& value )`.
+       * \tparam FetchValue is type returned by the Fetch lambda function.
+       *
+       * \param fetch is an instance of lambda function for data fetch.
+       * \param reduce is an instance of lambda function for reduction.
+       * \param keep in an instance of lambda function for storing results.
+       * \param zero is zero of given reduction operation also known as idempotent element.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixViewExample_reduceAllRows.cpp
+       * \par Output
+       * \include SparseMatrixViewExample_reduceAllRows.out
+       */
+      template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >
+      void reduceAllRows( Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchReal& zero ) const;
+
+      /**
+       * \brief Method for iteration over all matrix rows for constant instances.
+       *
+       * \tparam Function is type of lambda function that will operate on matrix elements.
+       *    It is should have form like
+       *  `function( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value )`.
+       *  The \e localIdx parameter is a rank of the non-zero element in given row.
+       *
+       * \param begin defines beginning of the range [begin,end) of rows to be processed.
+       * \param end defines ending of the range [begin,end) of rows to be processed.
+       * \param function is an instance of the lambda function to be called in each row.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixViewExample_forRows.cpp
+       * \par Output
+       * \include SparseMatrixViewExample_forRows.out
+       */
+      template< typename Function >
+      void forElements( IndexType begin, IndexType end, Function& function ) const;
+
+      /**
+       * \brief Method for iteration over all matrix rows for non-constant instances.
+       *
+       * \tparam Function is type of lambda function that will operate on matrix elements.
+       *    It is should have form like
+       *  `function( IndexType rowIdx, IndexType localIdx, IndexType columnIdx, const RealType& value )`.
+       *  The \e localIdx parameter is a rank of the non-zero element in given row.
+       *
+       * \param begin defines beginning of the range [begin,end) of rows to be processed.
+       * \param end defines ending of the range [begin,end) of rows to be processed.
+       * \param function is an instance of the lambda function to be called in each row.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixViewExample_forRows.cpp
+       * \par Output
+       * \include SparseMatrixViewExample_forRows.out
+       */
+      template< typename Function >
+      void forElements( IndexType begin, IndexType end, Function& function );
+
+      /**
+       * \brief This method calls \e forElements for all matrix rows (for constant instances).
+       *
+       * See \ref SparseSandboxMatrix::forElements.
+       *
+       * \tparam Function is a type of lambda function that will operate on matrix elements.
+       * \param function  is an instance of the lambda function to be called in each row.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixViewExample_forAllRows.cpp
+       * \par Output
+       * \include SparseMatrixViewExample_forAllRows.out
+       */
+      template< typename Function >
+      void forAllElements( Function& function ) const;
+
+      /**
+       * \brief This method calls \e forElements for all matrix rows.
+       *
+       * See \ref SparseSandboxMatrix::forElements.
+       *
+       * \tparam Function is a type of lambda function that will operate on matrix elements.
+       * \param function  is an instance of the lambda function to be called in each row.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixViewExample_forAllRows.cpp
+       * \par Output
+       * \include SparseMatrixViewExample_forAllRows.out
+       */
+      template< typename Function >
+      void forAllElements( Function& function );
+
+      /**
+       * \brief Method for parallel iteration over matrix rows from interval [ \e begin, \e end).
+       *
+       * In each row, given lambda function is performed. Each row is processed by at most one thread unlike the method
+       * \ref SparseSandboxMatrixView::forElements where more than one thread can be mapped to each row.
+
+       *
+       * \tparam Function is type of the lambda function.
+       *
+       * \param begin defines beginning of the range [ \e begin,\e end ) of rows to be processed.
+       * \param end defines ending of the range [ \e begin, \e end ) of rows to be processed.
+       * \param function is an instance of the lambda function to be called for each row.
+       *
+       * ```
+       * auto function = [] __cuda_callable__ ( RowView& row ) mutable { ... };
+       * ```
+       *
+       * \e RowView represents matrix row - see \ref TNL::Matrices::SparseSandboxMatrixView::RowView.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixViewExample_forRows.cpp
+       * \par Output
+       * \include SparseMatrixViewExample_forRows.out
+       */
+      template< typename Function >
+      void forRows( IndexType begin, IndexType end, Function&& function );
+
+      /**
+       * \brief Method for parallel iteration over matrix rows from interval [ \e begin, \e end) for constant instances.
+       *
+       * In each row, given lambda function is performed. Each row is processed by at most one thread unlike the method
+       * \ref SparseSandboxMatrixView::forElements where more than one thread can be mapped to each row.
+       *
+       * \tparam Function is type of the lambda function.
+       *
+       * \param begin defines beginning of the range [ \e begin,\e end ) of rows to be processed.
+       * \param end defines ending of the range [ \e begin, \e end ) of rows to be processed.
+       * \param function is an instance of the lambda function to be called for each row.
+       *
+       * ```
+       * auto function = [] __cuda_callable__ ( RowView& row ) { ... };
+       * ```
+       *
+       * \e RowView represents matrix row - see \ref TNL::Matrices::SparseSandboxMatrixView::RowView.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixViewExample_forRows.cpp
+       * \par Output
+       * \include SparseMatrixViewExample_forRows.out
+       */
+      template< typename Function >
+      void forRows( IndexType begin, IndexType end, Function&& function ) const;
+
+      /**
+       * \brief Method for parallel iteration over all matrix rows.
+       *
+       * In each row, given lambda function is performed. Each row is processed by at most one thread unlike the method
+       * \ref SparseSandboxMatrixView::forAllElements where more than one thread can be mapped to each row.
+       *
+       * \tparam Function is type of the lambda function.
+       *
+       * \param function is an instance of the lambda function to be called for each row.
+       *
+       * ```
+       * auto function = [] __cuda_callable__ ( RowView& row ) mutable { ... };
+       * ```
+       *
+       * \e RowView represents matrix row - see \ref TNL::Matrices::SparseSandboxMatrixView::RowView.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixViewExample_forRows.cpp
+       * \par Output
+       * \include SparseMatrixViewExample_forRows.out
+       */
+      template< typename Function >
+      void forAllRows( Function&& function );
+
+      /**
+       * \brief Method for parallel iteration over all matrix rows for constant instances.
+       *
+       * In each row, given lambda function is performed. Each row is processed by at most one thread unlike the method
+       * \ref SparseSandboxMatrixView::forAllElements where more than one thread can be mapped to each row.
+       *
+       * \tparam Function is type of the lambda function.
+       *
+       * \param function is an instance of the lambda function to be called for each row.
+       *
+       * ```
+       * auto function = [] __cuda_callable__ ( RowView& row ) { ... };
+       * ```
+       *
+       * \e RowView represents matrix row - see \ref TNL::Matrices::SparseSandboxMatrixView::RowView.
+       *
+       * \par Example
+       * \include Matrices/SparseMatrix/SparseMatrixViewExample_forRows.cpp
+       * \par Output
+       * \include SparseMatrixViewExample_forRows.out
+       */
+      template< typename Function >
+      void forAllRows( Function&& function ) const;
+
+      /**
+       * \brief Method for sequential iteration over all matrix rows for constant instances.
+       *
+       * \tparam Function is type of lambda function that will operate on matrix elements.
+       *    It is should have form like
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, const RealType& value )`.
+       *  The column index repeats twice only for compatibility with sparse matrices.
+       *
+       * \param begin defines beginning of the range [begin,end) of rows to be processed.
+       * \param end defines ending of the range [begin,end) of rows to be processed.
+       * \param function is an instance of the lambda function to be called in each row.
+       */
+      template< typename Function >
+      void sequentialForRows( IndexType begin, IndexType end, Function& function ) const;
+
+      /**
+       * \brief Method for sequential iteration over all matrix rows for non-constant instances.
+       *
+       * \tparam Function is type of lambda function that will operate on matrix elements.
+       *    It is should have form like
+       *  `function( IndexType rowIdx, IndexType columnIdx, IndexType columnIdx_, RealType& value )`.
+       *  The column index repeats twice only for compatibility with sparse matrices.
+       *
+       * \param begin defines beginning of the range [begin,end) of rows to be processed.
+       * \param end defines ending of the range [begin,end) of rows to be processed.
+       * \param function is an instance of the lambda function to be called in each row.
+       */
+      template< typename Function >
+      void sequentialForRows( IndexType begin, IndexType end, Function& function );
+
+      /**
+       * \brief This method calls \e sequentialForRows for all matrix rows (for constant instances).
+       *
+       * See \ref SparseSandboxMatrixView::sequentialForRows.
+       *
+       * \tparam Function is a type of lambda function that will operate on matrix elements.
+       * \param function  is an instance of the lambda function to be called in each row.
+       */
+      template< typename Function >
+      void sequentialForAllRows( Function& function ) const;
+
+      /**
+       * \brief This method calls \e sequentialForRows for all matrix rows.
+       *
+       * See \ref SparseSandboxMatrixView::sequentialForAllRows.
+       *
+       * \tparam Function is a type of lambda function that will operate on matrix elements.
+       * \param function  is an instance of the lambda function to be called in each row.
+       */
+      template< typename Function >
+      void sequentialForAllRows( Function& function );
+
+      /**
+       * \brief Computes product of matrix and vector.
+       *
+       * More precisely, it computes:
+       *
+       * `outVector = matrixMultiplicator * ( * this ) * inVector + outVectorMultiplicator * outVector`
+       *
+       * \tparam InVector is type of input vector.  It can be \ref Vector,
+       *     \ref VectorView, \ref Array, \ref ArraView or similar container.
+       * \tparam OutVector is type of output vector. It can be \ref Vector,
+       *     \ref VectorView, \ref Array, \ref ArraView or similar container.
+       *
+       * \param inVector is input vector.
+       * \param outVector is output vector.
+       * \param matrixMultiplicator is a factor by which the matrix is multiplied. It is one by default.
+       * \param outVectorMultiplicator is a factor by which the outVector is multiplied before added
+       *    to the result of matrix-vector product. It is zero by default.
+       * \param begin is the beginning of the rows range for which the vector product
+       *    is computed. It is zero by default.
+       * \param end is the end of the rows range for which the vector product
+       *    is computed. It is number if the matrix rows by default.
+       */
+      template< typename InVector,
+                typename OutVector >
+      void vectorProduct( const InVector& inVector,
+                          OutVector& outVector,
+                          const RealType matrixMultiplicator = 1.0,
+                          const RealType outVectorMultiplicator = 0.0,
+                          const IndexType begin = 0,
+                          IndexType end = 0 ) const;
+
+      template< typename Vector1, typename Vector2 >
+      bool performSORIteration( const Vector1& b,
+                                const IndexType row,
+                                Vector2& x,
+                                const RealType& omega = 1.0 ) const;
+
+      /**
+       * \brief Assignment of any matrix type.
+       * .
+       * \param matrix is input matrix for the assignment.
+       * \return reference to this matrix.
+       */
+      SparseSandboxMatrixView& operator=( const SparseSandboxMatrixView& matrix );
+
+      /**
+       * \brief Comparison operator with another arbitrary matrix type.
+       *
+       * \param matrix is the right-hand side matrix.
+       * \return \e true if the RHS matrix is equal, \e false otherwise.
+       */
+      template< typename Matrix >
+      bool operator==( const Matrix& m ) const;
+
+      /**
+       * \brief Comparison operator with another arbitrary matrix type.
+       *
+       * \param matrix is the right-hand side matrix.
+       * \return \e true if the RHS matrix is equal, \e false otherwise.
+       */
+      template< typename Matrix >
+      bool operator!=( const Matrix& m ) const;
+
+      /**
+       * \brief Method for saving the matrix to the file with given filename.
+       *
+       * \param fileName is name of the file.
+       */
+      void save( const String& fileName ) const;
+
+      /**
+       * \brief Method for saving the matrix to a file.
+       *
+       * \param file is the output file.
+       */
+      void save( File& file ) const;
+
+      /**
+       * \brief Method for printing the matrix to output stream.
+       *
+       * \param str is the output stream.
+       */
+      void print( std::ostream& str ) const;
+
+      /**
+       * \brief Getter of segments for non-constant instances.
+       *
+       * \e Segments are a structure for addressing the matrix elements columns and values.
+       * In fact, \e Segments represent the sparse matrix format.
+       *
+       * \return Non-constant reference to segments.
+       */
+      //SegmentsViewType& getSegments();
+
+      /**
+       * \brief Getter of segments for constant instances.
+       *
+       * \e Segments are a structure for addressing the matrix elements columns and values.
+       * In fact, \e Segments represent the sparse matrix format.
+       *
+       * \return Constant reference to segments.
+       */
+      //const SegmentsViewType& getSegments() const;
+
+      /**
+       * \brief Getter of column indexes for constant instances.
+       *
+       * \return Constant reference to a vector with matrix elements column indexes.
+       */
+      const ColumnsIndexesViewType& getColumnIndexes() const;
+
+      /**
+       * \brief Getter of column indexes for nonconstant instances.
+       *
+       * \return Reference to a vector with matrix elements column indexes.
+       */
+      ColumnsIndexesViewType& getColumnIndexes();
+
+      /**
+       * \brief Returns a padding index value.
+       *
+       * Padding index is used for column indexes of padding zeros. Padding zeros
+       * are used in some sparse matrix formats for better data alignment in memory.
+       *
+       * \return value of the padding index.
+       */
+      __cuda_callable__
+      IndexType getPaddingIndex() const;
+
+   protected:
+
+      ColumnsIndexesViewType columnIndexes;
+
+      RowPointersView rowPointers;
+      //SegmentsViewType segments;
+
+   private:
+      // TODO: this should be probably moved into a detail namespace
+      template< typename VectorOrView,
+                std::enable_if_t< HasSetSizeMethod< VectorOrView >::value, bool > = true >
+      static void set_size_if_resizable( VectorOrView& v, IndexType size )
+      {
+         v.setSize( size );
+      }
+
+      template< typename VectorOrView,
+                std::enable_if_t< ! HasSetSizeMethod< VectorOrView >::value, bool > = true >
+      static void set_size_if_resizable( VectorOrView& v, IndexType size )
+      {
+         TNL_ASSERT_EQ( v.getSize(), size, "view has wrong size" );
+      }
+};
+
+      } // namespace Sandbox
+   } // namespace Matrices
+} // namespace TNL
+
+#include <TNL/Matrices/Sandbox/SparseSandboxMatrixView.hpp>
diff --git a/src/TNL/Matrices/Sandbox/SparseSandboxMatrixView.hpp b/src/TNL/Matrices/Sandbox/SparseSandboxMatrixView.hpp
new file mode 100644
index 0000000000..fb50a7acf1
--- /dev/null
+++ b/src/TNL/Matrices/Sandbox/SparseSandboxMatrixView.hpp
@@ -0,0 +1,1026 @@
+/***************************************************************************
+                          SparseSandboxMatrixView.hpp -  description
+                             -------------------
+    begin                : Apr 20, 2021
+    copyright            : (C) 2021 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#pragma once
+
+#include <functional>
+#include <TNL/Matrices/Sandbox/SparseSandboxMatrixView.h>
+#include <TNL/Algorithms/Reduction.h>
+#include <TNL/Algorithms/AtomicOperations.h>
+#include <TNL/Matrices/details/SparseMatrix.h>
+
+namespace TNL {
+   namespace Matrices {
+      namespace Sandbox {
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+__cuda_callable__
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+SparseSandboxMatrixView()
+{
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+__cuda_callable__
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+SparseSandboxMatrixView( const IndexType rows,
+                         const IndexType columns,
+                         const ValuesViewType& values,
+                         const ColumnsIndexesViewType& columnIndexes,
+                         const RowPointersView& rowPointers )
+: MatrixView< Real, Device, Index >( rows, columns, values ), columnIndexes( columnIndexes ), rowPointers( rowPointers )
+{
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+__cuda_callable__
+auto
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+getView() -> ViewType
+{
+   return ViewType( this->getRows(),
+                    this->getColumns(),
+                    this->getValues().getView(),
+                    this->columnIndexes.getView(),
+                    this->segments.getView() );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+__cuda_callable__
+auto
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+getConstView() const -> ConstViewType
+{
+   return ConstViewType( this->getRows(),
+                         this->getColumns(),
+                         this->getValues().getConstView(),
+                         this->getColumnsIndexes().getConstView(),
+                         this->segments.getConstView() );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+String
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+getSerializationType()
+{
+   return String( "Matrices::Sandbox::SparseMatrix< " ) +
+             TNL::getSerializationType< RealType >() + ", " +
+             TNL::getSerializationType< IndexType >() + ", " +
+             MatrixType::getSerializationType() + ", [any_allocator], [any_allocator] >";
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+String
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+getSerializationTypeVirtual() const
+{
+   return this->getSerializationType();
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+   template< typename Vector >
+void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+getCompressedRowLengths( Vector& rowLengths ) const
+{
+   details::set_size_if_resizable( rowLengths, this->getRows() );
+   rowLengths = 0;
+   auto rowLengths_view = rowLengths.getView();
+   auto fetch = [] __cuda_callable__ ( IndexType row, IndexType column, const RealType& value ) -> IndexType {
+      return ( value != 0.0 );
+   };
+   auto keep = [=] __cuda_callable__ ( const IndexType rowIdx, const IndexType value ) mutable {
+      rowLengths_view[ rowIdx ] = value;
+   };
+   this->reduceAllRows( fetch, std::plus<>{}, keep, 0 );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+   template< typename Vector >
+void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+getRowCapacities( Vector& rowLengths ) const
+{
+   details::set_size_if_resizable( rowLengths, this->getRows() );
+   rowLengths = 0;
+   auto rowLengths_view = rowLengths.getView();
+   auto fetch = [] __cuda_callable__ ( IndexType row, IndexType column, const RealType& value ) -> IndexType {
+      return 1;
+   };
+   auto keep = [=] __cuda_callable__ ( const IndexType rowIdx, const IndexType value ) mutable {
+      rowLengths_view[ rowIdx ] = value;
+   };
+   this->reduceAllRows( fetch, std::plus<>{}, keep, 0 );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+__cuda_callable__
+Index
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+getRowCapacity( const IndexType row ) const
+{
+   return this->segments.getSegmentSize( row );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+Index
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+getNonzeroElementsCount() const
+{
+   const auto columns_view = this->columnIndexes.getConstView();
+   const IndexType paddingIndex = this->getPaddingIndex();
+   if( ! isSymmetric() )
+   {
+      auto fetch = [=] __cuda_callable__ ( const IndexType i ) -> IndexType {
+         return ( columns_view[ i ] != paddingIndex );
+      };
+      return Algorithms::Reduction< DeviceType >::reduce( ( IndexType ) 0, this->columnIndexes.getSize(), fetch, std::plus<>{}, 0 );
+   }
+   else
+   {
+      const auto rows = this->getRows();
+      const auto columns = this->getColumns();
+      Containers::Vector< IndexType, DeviceType, IndexType > row_sums( this->getRows(), 0 );
+      auto row_sums_view = row_sums.getView();
+      auto row_pointers_view = this->rowPointers.getConstView();
+      const auto columnIndexesView = this->columnIndexes.getConstView();
+      // SANDBOX_TODO: Replace the following lambda function (or more) with code compute number of nonzero matrix elements
+      //               of symmetric matrix. Note, that this is required only by symmetric matrices and that the highest performance
+      //               is not a priority here.
+      auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
+         auto begin = row_pointers_view[ rowIdx ];
+         auto end = row_pointers_view[ rowIdx + 1 ];
+         IndexType sum( 0 );
+         for( IndexType globalIdx = begin; globalIdx < end; globalIdx++ )
+         {
+            const IndexType column = columnIndexesView[ globalIdx ];
+            if( column != paddingIndex )
+               sum += 1 + ( column != rowIdx && column < rows && rowIdx < columns );
+         }
+         row_sums_view[ rowIdx ] = sum;
+      };
+      TNL::Algorithms::ParallelFor< DeviceType >::exec( ( IndexType ) 0, this->getRows(), f );
+      return sum( row_sums );
+   }
+   return 0;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+__cuda_callable__ auto
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+getRow( const IndexType& rowIdx ) const -> ConstRowView
+{
+   TNL_ASSERT_LT( rowIdx, this->getRows(), "Row index is larger than number of matrix rows." );
+   // SANDBOX_TODO: Replace the following with creation of RowView corresponding with your sparse matrix format.
+   return ConstRowView( rowIdx,                                                         // row index
+                        this->rowPointers[ rowIdx ],                                    // row begining
+                        this->rowPointers[ rowIdx + 1 ] - this->rowPointers[ rowIdx ],  // number of elemnts allocated for given row
+                        this->values, this->columnIndexes );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+__cuda_callable__ auto
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+getRow( const IndexType& rowIdx ) -> RowView
+{
+   TNL_ASSERT_LT( rowIdx, this->getRows(), "Row index is larger than number of matrix rows." );
+   // SANDBOX_TODO: Replace this with RowView constructor by your needs.
+   return RowView( rowIdx,                                               // row index
+                   rowPointers[ rowIdx ],                                // index of the first nonzero element in the row
+                   rowPointers[ rowIdx + 1 ] - rowPointers[ rowIdx ],    // number of nonzero elements in the row
+                   this->values,
+                   this->columnIndexes );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+__cuda_callable__ void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+setElement( const IndexType row,
+            const IndexType column,
+            const RealType& value )
+{
+   this->addElement( row, column, value, 0.0 );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+__cuda_callable__ void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+addElement( IndexType row,
+            IndexType column,
+            const RealType& value,
+            const RealType& thisElementMultiplicator )
+{
+   TNL_ASSERT_GE( row, 0, "Sparse matrix row index cannot be negative." );
+   TNL_ASSERT_LT( row, this->getRows(), "Sparse matrix row index is larger than number of matrix rows." );
+   TNL_ASSERT_GE( column, 0, "Sparse matrix column index cannot be negative." );
+   TNL_ASSERT_LT( column, this->getColumns(), "Sparse matrix column index is larger than number of matrix columns." );
+
+   if( isSymmetric() && row < column )
+   {
+      swap( row, column );
+      TNL_ASSERT_LT( row, this->getRows(), "Column index is out of the symmetric part of the matrix after transposition." );
+      TNL_ASSERT_LT( column,this->getColumns(), "Row index is out of the symmetric part of the matrix after transposition." );
+   }
+
+   // SANDBOX_TODO: Replace the following line with a code that computes number of matrix elements allocated for 
+   //               matrix row with indedx `row`. Note that the code works on both host and GPU kernel. To achieve
+   //               the same effect, you may use macro __CUDA_ARCH__ as can be seen bellow in this method.
+   const IndexType rowSize = this->rowPointers.getElement( row + 1 ) - this->rowPointers.getElement( row );
+   IndexType col( this->getPaddingIndex() );
+   IndexType i;
+   IndexType globalIdx;
+   for( i = 0; i < rowSize; i++ )
+   {
+      // SANDBOX_TODO: Replace the following line with a code that computes a global index of `i`-th nonzero matrix element
+      //               in the `row`-th matrix row. The global index is a pointer to arrays `values` and `columnIndexes` storing
+      //               the matrix elements values and column indexes respectively.
+      globalIdx = this->rowPointers.getElement( row ) + i;
+      TNL_ASSERT_LT( globalIdx, this->columnIndexes.getSize(), "" );
+      col = this->columnIndexes.getElement( globalIdx );
+      if( col == column )
+      {
+         if( ! isBinary() )
+            this->values.setElement( globalIdx, thisElementMultiplicator * this->values.getElement( globalIdx ) + value );
+         return;
+      }
+      if( col == this->getPaddingIndex() || col > column )
+         break;
+   }
+   if( i == rowSize )
+   {
+#ifndef __CUDA_ARCH__
+      std::stringstream msg;
+      msg << "The capacity of the sparse matrix row number "  << row << " was exceeded.";
+      throw std::logic_error( msg.str() );
+#else
+      TNL_ASSERT_TRUE( false, "");
+      return;
+#endif
+   }
+   if( col == this->getPaddingIndex() )
+   {
+      this->columnIndexes.setElement( globalIdx, column );
+      if( ! isBinary() )
+         this->values.setElement( globalIdx, value );
+      return;
+   }
+   else
+   {
+      IndexType j = rowSize - 1;
+      while( j > i )
+      {
+         // SANDBOX_TODO: Replace the following two lines with a code that computes a global indexes of `j`-th and `j-1`-th nonzero matrix elements
+         //               in the `row`-th matrix row. The global index is a pointer to arrays `values` and `columnIndexes` storing
+         //               the matrix elements values and column indexes respectively.
+         const IndexType globalIdx1 = this->rowPointers.getElement( row ) + j;
+         const IndexType globalIdx2 = globalIdx1 - 1;
+         // End of code replacement.
+         TNL_ASSERT_LT( globalIdx1, this->columnIndexes.getSize(), "" );
+         TNL_ASSERT_LT( globalIdx2, this->columnIndexes.getSize(), "" );
+         this->columnIndexes.setElement( globalIdx1, this->columnIndexes.getElement( globalIdx2 ) );
+         if( ! isBinary() )
+            this->values.setElement( globalIdx1, this->values.getElement( globalIdx2 ) );
+         j--;
+      }
+
+      this->columnIndexes.setElement( globalIdx, column );
+      if( ! isBinary() )
+         this->values.setElement( globalIdx, value );
+      return;
+   }
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+__cuda_callable__
+auto
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+getElement( IndexType row,
+            IndexType column ) const -> RealType
+{
+   TNL_ASSERT_GE( row, 0, "Sparse matrix row index cannot be negative." );
+   TNL_ASSERT_LT( row, this->getRows(), "Sparse matrix row index is larger than number of matrix rows." );
+   TNL_ASSERT_GE( column, 0, "Sparse matrix column index cannot be negative." );
+   TNL_ASSERT_LT( column, this->getColumns(), "Sparse matrix column index is larger than number of matrix columns." );
+
+   if( isSymmetric() && row < column )
+   {
+      swap( row, column );
+      if( row >= this->getRows() || column >= this->getColumns() )
+         return 0.0;
+   }
+
+   // SANDBOX_TODO: Replace the following lines with a code for getting number of elements allocated for given row.
+   const IndexType rowSize = this->rowPointers.getElement( row + 1 ) - this->rowPointers.getElement( row );
+   for( IndexType i = 0; i < rowSize; i++ )
+   {
+      // SANDBOX_TODO: Replace the following line with a code for getting index of the matrix element in arrays `values` and `columnIdexes`.
+      const IndexType globalIdx = this->rowPointers.getElement( row ) + i;
+      TNL_ASSERT_LT( globalIdx, this->columnIndexes.getSize(), "" );
+      const IndexType col = this->columnIndexes.getElement( globalIdx );
+      if( col == column )
+      {
+         if( isBinary() )
+            return 1;
+         else
+            return this->values.getElement( globalIdx );
+      }
+   }
+   return 0.0;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+template< typename InVector,
+       typename OutVector >
+void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+vectorProduct( const InVector& inVector,
+               OutVector& outVector,
+               const RealType matrixMultiplicator,
+               const RealType outVectorMultiplicator,
+               const IndexType firstRow,
+               IndexType lastRow ) const
+{
+   TNL_ASSERT_EQ( this->getColumns(), inVector.getSize(), "Matrix columns do not fit with input vector." );
+   TNL_ASSERT_EQ( this->getRows(), outVector.getSize(), "Matrix rows do not fit with output vector." );
+
+   using OutVectorReal = typename OutVector::RealType;
+   static_assert(
+         ! MatrixType::isSymmetric() ||
+         ! std::is_same< Device, Devices::Cuda >::value ||
+         ( std::is_same< OutVectorReal, float >::value ||
+           std::is_same< OutVectorReal, double >::value ||
+           std::is_same< OutVectorReal, int >::value ||
+           std::is_same< OutVectorReal, long long int >::value ),
+         "Given Real type is not supported by atomic operations on GPU which are necessary for symmetric operations." );
+
+   const auto inVectorView = inVector.getConstView();
+   auto outVectorView = outVector.getView();
+   const auto valuesView = this->values.getConstView();
+   const auto columnIndexesView = this->columnIndexes.getConstView();
+   const auto rowPointersView = this->rowPointers.getConstView();
+   const IndexType paddingIndex = this->getPaddingIndex();
+#define HAVE_SANDBOX_SIMPLE_SPMV
+   // SANDBOX_TODO: The following is simple direct implementation of SpMV operation with CSR format. We recommend to start by
+   //               replacing this part with SpMV based on your sparse format.
+   if( std::is_same< DeviceType, TNL::Devices::Host >::value )          // this way you may easily specialize for different device types
+   {
+      // SANDBOX_TODO: This simple and naive implementation for CPU.
+      for( IndexType rowIdx = firstRow; rowIdx < lastRow; rowIdx++ )
+      {
+         const auto begin = rowPointers[ rowIdx ];
+         const auto end = rowPointers[ rowIdx + 1 ];
+         RealType sum( 0.0 );
+         for( IndexType globalIdx = begin; globalIdx < end; globalIdx++ )
+               sum += this->values[ globalIdx ] * inVector[ this->columnIndexes[ globalIdx ] ];
+         // SANDBOX_TODO:The following is quite inefficient, its better to specialized the code for cases when
+         // `outVectorMultiplicator` is zero or `matrixMultiplicator` is one - see. the full implementation bellow.
+         outVector[ rowIdx ] = outVector[ rowIdx ] * outVectorMultiplicator + matrixMultiplicator * sum;
+      }
+   }
+   else
+   {
+      //SANDBOX_TODO: The following is general implementation based on ParallelFor and lambda function. It would work even on CPU.
+      auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
+         const auto begin = rowPointersView[ rowIdx ];
+         const auto end = rowPointersView[ rowIdx + 1 ];
+         RealType sum( 0.0 );
+         for( IndexType globalIdx = begin; globalIdx < end; globalIdx++ )
+            sum + valuesView[ globalIdx ] * inVectorView[ columnIndexesView[ globalIdx ] ];
+         outVectorView[ rowIdx ] = outVectorView[ rowIdx ] * outVectorMultiplicator + matrixMultiplicator * sum;
+      };
+      TNL::Algorithms::ParallelFor< DeviceType >::exec( firstRow, lastRow, f );
+   }
+#ifdef HAVE_SANDBOX_SIMPLE_SPMV
+#else
+   // SANDBOX_TODO: The following is fully functional implementation based on method `reduceRows`.
+   if( isSymmetric() )
+      outVector *= outVectorMultiplicator;
+   auto symmetricFetch = [=] __cuda_callable__ ( IndexType row, IndexType localIdx, IndexType globalIdx, bool& compute ) mutable -> RealType {
+      const IndexType column = columnIndexesView[ globalIdx ];
+      compute = ( column != paddingIndex );
+      if( ! compute )
+         return 0.0;
+      if( isSymmetric() && column < row )
+      {
+         if( isBinary() )
+            Algorithms::AtomicOperations< DeviceType >::add( outVectorView[ column ], ( OutVectorReal ) matrixMultiplicator * inVectorView[ row ] );
+         else
+            Algorithms::AtomicOperations< DeviceType >::add( outVectorView[ column ], ( OutVectorReal ) matrixMultiplicator * valuesView[ globalIdx ] * inVectorView[ row ] );
+      }
+      if( isBinary() )
+         return inVectorView[ column ];
+      return valuesView[ globalIdx ] * inVectorView[ column ];
+   };
+   auto fetch = [=] __cuda_callable__ ( IndexType globalIdx, bool& compute ) mutable -> RealType {
+      const IndexType column = columnIndexesView[ globalIdx ];
+      if( isBinary() )
+         return inVectorView[ column ];
+      return valuesView[ globalIdx ] * inVectorView[ column ];
+   };
+
+   auto keeperGeneral = [=] __cuda_callable__ ( IndexType row, const RealType& value ) mutable {
+      if( isSymmetric() )
+      {
+         typename OutVector::RealType aux = matrixMultiplicator * value;
+         Algorithms::AtomicOperations< DeviceType >::add( outVectorView[ row ], aux );
+      }
+      else
+      {
+         if( outVectorMultiplicator == 0.0 )
+            outVectorView[ row ] = matrixMultiplicator * value;
+         else
+            outVectorView[ row ] = outVectorMultiplicator * outVectorView[ row ] + matrixMultiplicator * value;
+      }
+   };
+   auto keeperDirect = [=] __cuda_callable__ ( IndexType row, const RealType& value ) mutable {
+      outVectorView[ row ] = value;
+   };
+   auto keeperMatrixMult = [=] __cuda_callable__ ( IndexType row, const RealType& value ) mutable {
+      outVectorView[ row ] = matrixMultiplicator * value;
+   };
+   auto keeperVectorMult = [=] __cuda_callable__ ( IndexType row, const RealType& value ) mutable {
+      outVectorView[ row ] = outVectorMultiplicator * outVectorView[ row ] + value;
+   };
+
+   if( lastRow == 0 )
+      lastRow = this->getRows();
+   if( isSymmetric() )
+      this->reduceRows( firstRow, lastRow, symmetricFetch, std::plus<>{}, keeperGeneral, ( RealType ) 0.0 );
+   else
+   {
+      if( outVectorMultiplicator == 0.0 )
+      {
+         if( matrixMultiplicator == 1.0 )
+            this->reduceRows( firstRow, lastRow, fetch, std::plus<>{}, keeperDirect, ( RealType ) 0.0 );
+         else
+            this->reduceRows( firstRow, lastRow, fetch, std::plus<>{}, keeperMatrixMult, ( RealType ) 0.0 );
+      }
+      else
+      {
+         if( matrixMultiplicator == 1.0 )
+            this->reduceRows( firstRow, lastRow, fetch, std::plus<>{}, keeperVectorMult, ( RealType ) 0.0 );
+         else
+            this->reduceRows( firstRow, lastRow, fetch, std::plus<>{}, keeperGeneral, ( RealType ) 0.0 );
+      }
+   }
+#endif
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+   template< typename Fetch, typename Reduce, typename Keep, typename FetchValue >
+void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+reduceRows( IndexType begin, IndexType end, Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchValue& zero )
+{
+   auto columns_view = this->columnIndexes.getView();
+   auto values_view = this->values.getView();
+   auto row_pointers_view = this->rowPointers.getConstView();
+   const IndexType paddingIndex_ = this->getPaddingIndex();
+   // SANDBOX_TODO: Replace the following code with the one for computing reduction in rows by your format.
+   //               Note, that this method can be used for implementation of SpMV.
+   auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
+      const auto begin = row_pointers_view[ rowIdx ];
+      const auto end = row_pointers_view[ rowIdx + 1 ];
+      FetchValue sum = zero;
+      for( IndexType globalIdx = begin; globalIdx < end; globalIdx++ )
+      {
+         IndexType& columnIdx = columns_view[ globalIdx ];
+         if( columnIdx != paddingIndex_ )
+         {
+            if( isBinary() )
+               sum = reduce( sum, fetch( rowIdx, columnIdx, 1 ) );
+            else
+               sum = reduce( sum, fetch( rowIdx, columnIdx, values_view[ globalIdx ] ) );
+         }
+      }
+      keep( rowIdx, sum );
+   };
+   TNL::Algorithms::ParallelFor< DeviceType >::exec( begin, end, f );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+   template< typename Fetch, typename Reduce, typename Keep, typename FetchValue >
+void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+reduceRows( IndexType begin, IndexType end, Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchValue& zero ) const
+{
+   auto columns_view = this->columnIndexes.getConstView();
+   auto values_view = this->values.getConstView();
+   const IndexType paddingIndex_ = this->getPaddingIndex();
+   // SANDBOX_TODO: Replace the following code with the one for computing reduction in rows by your format.
+   //               Note, that this method can be used for implementation of SpMV.
+   auto row_pointers_view = this->rowPointers.getConstView();
+   auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
+      const auto begin = row_pointers_view[ rowIdx ];
+      const auto end = row_pointers_view[ rowIdx + 1 ];
+      FetchValue sum = zero;
+      for( IndexType globalIdx = begin; globalIdx < end; globalIdx++ )
+      {
+         const IndexType& columnIdx = columns_view[ globalIdx ];
+         if( columnIdx != paddingIndex_ )
+         {
+            if( isBinary() )
+               sum = reduce( sum, fetch( rowIdx, columnIdx, 1 ) );
+            else
+               sum = reduce( sum, fetch( rowIdx, columnIdx, values_view[ globalIdx ] ) );
+         }
+      }
+      keep( rowIdx, sum );
+   };
+   TNL::Algorithms::ParallelFor< DeviceType >::exec( begin, end, f );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+   template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >
+void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+reduceAllRows( Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchReal& zero )
+{
+   this->reduceRows( 0, this->getRows(), fetch, reduce, keep, zero );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+   template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >
+void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+reduceAllRows( Fetch& fetch, const Reduce& reduce, Keep& keep, const FetchReal& zero ) const
+{
+   this->reduceRows( 0, this->getRows(), fetch, reduce, keep, zero );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+   template< typename Function >
+void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+forElements( IndexType begin, IndexType end, Function& function ) const
+{
+   const auto columns_view = this->columnIndexes.getConstView();
+   const auto values_view = this->values.getConstView();
+   // SANDBOX_TODO: Replace the following code with the one for iterating over all allocated matrix elements.
+   auto row_pointers_view = this->rowPointers.getConstView();
+   auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
+      const auto begin = row_pointers_view[ rowIdx ];
+      const auto end = row_pointers_view[ rowIdx + 1 ];
+      IndexType localIdx( 0 );
+      for( IndexType globalIdx = begin; globalIdx < end; globalIdx++ )
+      {
+         if( isBinary() )
+            function( rowIdx, localIdx, columns_view[ globalIdx ], ( RealType ) 1 );
+         else
+            function( rowIdx, localIdx, columns_view[ globalIdx ], values_view[ globalIdx ] );
+         localIdx++;
+      }
+   };
+   TNL::Algorithms::ParallelFor< DeviceType >::exec( begin, end, f );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+   template< typename Function >
+void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+forElements( IndexType begin, IndexType end, Function& function )
+{
+   auto columns_view = this->columnIndexes.getView();
+   auto values_view = this->values.getView();
+   // SANDBOX_TODO: Replace the following code with the one for iterating over all allocated matrix elements.
+   auto row_pointers_view = this->rowPointers.getConstView();
+   auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
+      const auto begin = row_pointers_view[ rowIdx ];
+      const auto end = row_pointers_view[ rowIdx + 1 ];
+      IndexType localIdx( 0 );
+      RealType one( 1.0 );
+      for( IndexType globalIdx = begin; globalIdx < end; globalIdx++ )
+      {
+         if( isBinary() )
+            function( rowIdx, localIdx, columns_view[ globalIdx ], one ); // TODO: Fix this without using `one`.
+         else
+            function( rowIdx, localIdx, columns_view[ globalIdx ], values_view[ globalIdx ] );
+         localIdx++;
+      }
+   };
+   TNL::Algorithms::ParallelFor< DeviceType >::exec( begin, end, f );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+   template< typename Function >
+void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+forAllElements( Function& function ) const
+{
+   this->forElements( 0, this->getRows(), function );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+   template< typename Function >
+void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+forAllElements( Function& function )
+{
+   this->forElements( 0, this->getRows(), function );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+   template< typename Function >
+void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+forRows( IndexType begin, IndexType end, Function&& function )
+{
+   auto columns_view = this->columnIndexes.getView();
+   auto values_view = this->values.getView();
+   // SANDBOX_TODO: Replace the following code with the one for iteration over matrix rows.
+   auto row_pointers_view = this->rowPointers.getConstView();
+   auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
+      auto rowView = RowView( rowIdx,                                                       // row index
+                              row_pointers_view[ rowIdx ],                                  // row begining
+                              row_pointers_view[ rowIdx + 1 ] -row_pointers_view[ rowIdx ], // number of elemnts allocated for given matrix row
+                              values_view, columns_view );
+      function( rowView );
+   };
+   TNL::Algorithms::ParallelFor< DeviceType >::exec( begin, end, f );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+   template< typename Function >
+void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+forRows( IndexType begin, IndexType end, Function&& function ) const
+{
+   const auto columns_view = this->columnIndexes.getConstView();
+   const auto values_view = this->values.getConstView();
+   // SANDBOX_TODO: Replace the following code with the one for iteration over matrix rows.
+   auto row_pointers_view = this->rowPointers.getConstView();
+   auto f = [=] __cuda_callable__ ( IndexType rowIdx ) {
+      auto rowView = ConstRowView( rowIdx,                                                       // row index
+                                   row_pointers_view[ rowIdx ],                                  // row begining
+                                   row_pointers_view[ rowIdx + 1 ] -row_pointers_view[ rowIdx ], // number of elemnts allocated for given matrix row
+                                   values_view, columns_view );
+      function( rowView );
+   };
+   TNL::Algorithms::ParallelFor< DeviceType >::exec( begin, end, f );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+   template< typename Function >
+void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+forAllRows( Function&& function )
+{
+   this->forRows( 0, this->getRows(), function );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+   template< typename Function >
+void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+forAllRows( Function&& function ) const
+{
+   this->forRows( 0, this->getRows(), function );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+   template< typename Function >
+void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+sequentialForRows( IndexType begin, IndexType end, Function& function ) const
+{
+   for( IndexType row = begin; row < end; row ++ )
+      this->forRows( row, row + 1, function );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+   template< typename Function >
+void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+sequentialForRows( IndexType begin, IndexType end, Function& function )
+{
+   for( IndexType row = begin; row < end; row ++ )
+      this->forRows( row, row + 1, function );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+   template< typename Function >
+void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+sequentialForAllRows( Function& function ) const
+{
+   this->sequentialForRows( 0, this->getRows(), function );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+   template< typename Function >
+void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+sequentialForAllRows( Function& function )
+{
+   this->sequentialForRows( 0, this->getRows(), function );
+}
+
+/*template< typename Real,
+          template< typename, typename > class SegmentsView,
+          typename Device,
+          typename Index,
+          typename RealAllocator,
+          typename IndexAllocator >
+template< typename Real2, template< typename, typename > class Segments2, typename Index2, typename RealAllocator2, typename IndexAllocator2 >
+void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+addMatrix( const SparseSandboxMatrixView< Real2, Segments2, Device, Index2, RealAllocator2, IndexAllocator2 >& matrix,
+           const RealType& matrixMultiplicator,
+           const RealType& thisMatrixMultiplicator )
+{
+
+}
+
+template< typename Real,
+          template< typename, typename > class SegmentsView,
+          typename Device,
+          typename Index,
+          typename RealAllocator,
+          typename IndexAllocator >
+template< typename Real2, typename Index2 >
+void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+getTransposition( const SparseSandboxMatrixView< Real2, Device, Index2 >& matrix,
+                  const RealType& matrixMultiplicator )
+{
+
+}*/
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+template< typename Vector1, typename Vector2 >
+bool
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+performSORIteration( const Vector1& b,
+                     const IndexType row,
+                     Vector2& x,
+                     const RealType& omega ) const
+{
+   return false;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >&
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+operator=( const SparseSandboxMatrixView< Real, Device, Index, MatrixType >& matrix )
+{
+   MatrixView< Real, Device, Index >::operator=( matrix );
+   this->columnIndexes.bind( matrix.columnIndexes );
+   // SANDBOX_TODO: Replace the following line with assignment of metadata required by your
+   //               sparse format.
+   this->rowPointers.bind( matrix.rowPointers );
+   return *this;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+   template< typename Matrix >
+bool
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+operator==( const Matrix& m ) const
+{
+   const auto& view1 = *this;
+   // FIXME: getConstView does not work
+   //const auto view2 = m.getConstView();
+   const auto view2 = m.getView();
+   auto fetch = [=] __cuda_callable__ ( const IndexType i ) -> bool
+   {
+      return view1.getRow( i ) == view2.getRow( i );
+   };
+   return Algorithms::Reduction< DeviceType >::reduce( 0, this->getRows(), fetch, std::logical_and<>{}, true );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+   template< typename Matrix >
+bool
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+operator!=( const Matrix& m ) const
+{
+   return ! operator==( m );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+save( File& file ) const
+{
+   MatrixView< RealType, DeviceType, IndexType >::save( file );
+   file << this->columnIndexes
+        << this->rowPointers;  // SANDBOX_TODO: Replace this with medata required by your format
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+save( const String& fileName ) const
+{
+   Object::save( fileName );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+void
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+print( std::ostream& str ) const
+{
+   if( isSymmetric() )
+   {
+      for( IndexType row = 0; row < this->getRows(); row++ )
+      {
+         str <<"Row: " << row << " -> ";
+         for( IndexType column = 0; column < this->getColumns(); column++ )
+         {
+            auto value = this->getElement( row, column );
+            if( value != ( RealType ) 0 )
+               str << " Col:" << column << "->" << value << "\t";
+         }
+         str << std::endl;
+      }
+   }
+   else
+      for( IndexType row = 0; row < this->getRows(); row++ )
+      {
+         str <<"Row: " << row << " -> ";
+         // SANDBOX_TODO: Replace the followinf line with a code for computing number of elements allocated for given matrix row.
+         const auto rowLength = this->rowPointers.getElement( row + 1 ) - this->rowPointers.getElement( row );
+         for( IndexType i = 0; i < rowLength; i++ )
+         {
+            // SANDBOX_TODO: Replace the following line with a code for getting index of the matrix element in arrays `values` and `columnIdexes`.
+            const IndexType globalIdx = this->rowPointers.getElement( row ) + i;
+            const IndexType column = this->columnIndexes.getElement( globalIdx );
+            if( column == this->getPaddingIndex() )
+               break;
+            RealType value;
+            if( isBinary() )
+               value = ( RealType ) 1.0;
+            else
+               value = this->values.getElement( globalIdx );
+            if( value )
+            {
+               std::stringstream str_;
+               str_ << std::setw( 4 ) << std::right << column << ":" << std::setw( 4 ) << std::left << value;
+               str << std::setw( 10 ) << str_.str();
+            }
+         }
+         str << std::endl;
+      }
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+__cuda_callable__
+Index
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+getPaddingIndex() const
+{
+   return -1;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+auto
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+getColumnIndexes() const -> const ColumnsIndexesViewType&
+{
+   return this->columnIndexes;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index,
+          typename MatrixType >
+auto
+SparseSandboxMatrixView< Real, Device, Index, MatrixType >::
+getColumnIndexes() -> ColumnsIndexesViewType&
+{
+   return this->columnIndexes;
+}
+
+      } // namespace Sandbox
+   } //namespace Matrices
+} // namespace  TNL
diff --git a/src/TNL/Matrices/SparseMatrix.hpp b/src/TNL/Matrices/SparseMatrix.hpp
index 40c568b4df..b58620c2d0 100644
--- a/src/TNL/Matrices/SparseMatrix.hpp
+++ b/src/TNL/Matrices/SparseMatrix.hpp
@@ -417,6 +417,7 @@ SparseMatrix< Real, Device, Index, MatrixType, Segments, ComputeReal, RealAlloca
 reset()
 {
    BaseType::reset();
+   this->columnIndexes.reset();
    this->segments.reset();
    this->view = this->getView();
    TNL_ASSERT_EQ( this->getRows(), segments.getSegmentsCount(), "mismatched segments count" );
diff --git a/src/UnitTests/Matrices/CMakeLists.txt b/src/UnitTests/Matrices/CMakeLists.txt
index a4b06708e6..3210920f58 100644
--- a/src/UnitTests/Matrices/CMakeLists.txt
+++ b/src/UnitTests/Matrices/CMakeLists.txt
@@ -28,6 +28,8 @@ set( COMMON_TESTS
             BinarySparseMatrixCopyTest
             SymmetricSparseMatrixTest_CSR
             LambdaMatrixTest
+            SparseMatrixTest_SandboxMatrix
+            SparseMatrixVectorProductTest_SandboxMatrix
 )
 
 set( CPP_TESTS
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_SandboxMatrix.cpp b/src/UnitTests/Matrices/SparseMatrixTest_SandboxMatrix.cpp
new file mode 100644
index 0000000000..dc856310e5
--- /dev/null
+++ b/src/UnitTests/Matrices/SparseMatrixTest_SandboxMatrix.cpp
@@ -0,0 +1,11 @@
+/***************************************************************************
+                          SparseMatrixTest_SandboxMatrix.cpp -  description
+                             -------------------
+    begin                : Apr 19, 2021
+    copyright            : (C) 2021 by Tomas Oberhuber et al.
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#include "SparseMatrixTest_SandboxMatrix.h"
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_SandboxMatrix.cu b/src/UnitTests/Matrices/SparseMatrixTest_SandboxMatrix.cu
new file mode 120000
index 0000000000..27787fdf24
--- /dev/null
+++ b/src/UnitTests/Matrices/SparseMatrixTest_SandboxMatrix.cu
@@ -0,0 +1 @@
+SparseMatrixTest_SandboxMatrix.cpp
\ No newline at end of file
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_SandboxMatrix.h b/src/UnitTests/Matrices/SparseMatrixTest_SandboxMatrix.h
new file mode 100644
index 0000000000..ad1a0c74d6
--- /dev/null
+++ b/src/UnitTests/Matrices/SparseMatrixTest_SandboxMatrix.h
@@ -0,0 +1,45 @@
+/***************************************************************************
+                          SandboxMatrixTest_SandboxMatrix.h -  description
+                             -------------------
+    begin                : Apr 19, 2021
+    copyright            : (C) 2021 by Tomas Oberhuber et al.
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#include <iostream>
+#include <TNL/Matrices/Sandbox/SparseSandboxMatrix.h>
+
+#ifdef HAVE_GTEST
+#include <gtest/gtest.h>
+
+const char* saveAndLoadFileName = "test_SandboxMatrixTest_SandboxMatrix";
+
+// types for which MatrixTest is instantiated
+using MatrixTypes = ::testing::Types
+<
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< int,     TNL::Devices::Host, int,   TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< long,    TNL::Devices::Host, int,   TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< float,   TNL::Devices::Host, int,   TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< double,  TNL::Devices::Host, int,   TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< int,     TNL::Devices::Host, long,  TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< long,    TNL::Devices::Host, long,  TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< float,   TNL::Devices::Host, long,  TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< double,  TNL::Devices::Host, long,  TNL::Matrices::GeneralMatrix >
+#ifdef HAVE_CUDA
+   ,TNL::Matrices::Sandbox::SparseSandboxMatrix< int,     TNL::Devices::Cuda, int,   TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< long,    TNL::Devices::Cuda, int,   TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< float,   TNL::Devices::Cuda, int,   TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< double,  TNL::Devices::Cuda, int,   TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< int,     TNL::Devices::Cuda, long,  TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< long,    TNL::Devices::Cuda, long,  TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< float,   TNL::Devices::Cuda, long,  TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< double,  TNL::Devices::Cuda, long,  TNL::Matrices::GeneralMatrix >
+#endif
+>;
+
+#endif
+
+#include "SparseMatrixTest.h"
+#include "../main.h"
diff --git a/src/UnitTests/Matrices/SparseMatrixVectorProductTest_SandboxMatrix.cpp b/src/UnitTests/Matrices/SparseMatrixVectorProductTest_SandboxMatrix.cpp
new file mode 100644
index 0000000000..bfa16c02b9
--- /dev/null
+++ b/src/UnitTests/Matrices/SparseMatrixVectorProductTest_SandboxMatrix.cpp
@@ -0,0 +1,11 @@
+/***************************************************************************
+                          SparseMatrixVectorProductTest_CSRScalar.cpp -  description
+                             -------------------
+    begin                : Mar 30, 2021
+    copyright            : (C) 2021 by Tomas Oberhuber et al.
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#include "SparseMatrixVectorProductTest_CSRScalar.h"
diff --git a/src/UnitTests/Matrices/SparseMatrixVectorProductTest_SandboxMatrix.cu b/src/UnitTests/Matrices/SparseMatrixVectorProductTest_SandboxMatrix.cu
new file mode 120000
index 0000000000..bd87e1ad0c
--- /dev/null
+++ b/src/UnitTests/Matrices/SparseMatrixVectorProductTest_SandboxMatrix.cu
@@ -0,0 +1 @@
+SparseMatrixVectorProductTest_SandboxMatrix.cpp
\ No newline at end of file
diff --git a/src/UnitTests/Matrices/SparseMatrixVectorProductTest_SandboxMatrix.h b/src/UnitTests/Matrices/SparseMatrixVectorProductTest_SandboxMatrix.h
new file mode 100644
index 0000000000..7b06af0f3a
--- /dev/null
+++ b/src/UnitTests/Matrices/SparseMatrixVectorProductTest_SandboxMatrix.h
@@ -0,0 +1,45 @@
+/***************************************************************************
+                          SparseMatrixVectorProductTest_SandbxMatrix.h -  description
+                             -------------------
+    begin                : Apr 22, 2021
+    copyright            : (C) 2021 by Tomas Oberhuber et al.
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#include <iostream>
+#include <TNL/Matrices/Sandbox/SparseSandboxMatrix.h>
+
+#ifdef HAVE_GTEST
+#include <gtest/gtest.h>
+
+const char* saveAndLoadFileName = "test_SparseMatrixTest_CSRScalar_segments";
+
+// types for which MatrixTest is instantiated
+using MatrixTypes = ::testing::Types
+<
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< int,     TNL::Devices::Host, int,   TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< long,    TNL::Devices::Host, int,   TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< float,   TNL::Devices::Host, int,   TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< double,  TNL::Devices::Host, int,   TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< int,     TNL::Devices::Host, long,  TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< long,    TNL::Devices::Host, long,  TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< float,   TNL::Devices::Host, long,  TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< double,  TNL::Devices::Host, long,  TNL::Matrices::GeneralMatrix >
+#ifdef HAVE_CUDA
+   ,TNL::Matrices::Sandbox::SparseSandboxMatrix< int,     TNL::Devices::Cuda, int,   TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< long,    TNL::Devices::Cuda, int,   TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< float,   TNL::Devices::Cuda, int,   TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< double,  TNL::Devices::Cuda, int,   TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< int,     TNL::Devices::Cuda, long,  TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< long,    TNL::Devices::Cuda, long,  TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< float,   TNL::Devices::Cuda, long,  TNL::Matrices::GeneralMatrix >,
+    TNL::Matrices::Sandbox::SparseSandboxMatrix< double,  TNL::Devices::Cuda, long,  TNL::Matrices::GeneralMatrix >
+#endif
+>;
+
+#endif
+
+#include "SparseMatrixVectorProductTest.h"
+#include "../main.h"
-- 
GitLab