diff --git a/src/TNL/Containers/Segments/SlicedEllpack.h b/src/TNL/Containers/Segments/SlicedEllpack.h
index 1c110b1f11b18457bf9eeb8a70746b256739f6af..76185bcace2661e85881a0b9dc80fab599f34b2d 100644
--- a/src/TNL/Containers/Segments/SlicedEllpack.h
+++ b/src/TNL/Containers/Segments/SlicedEllpack.h
@@ -10,6 +10,7 @@
 
 #pragma once
 
+#include <TNL/Allocators/Default.h>
 #include <TNL/Containers/Vector.h>
 #include <TNL/Containers/Segments/SlicedEllpackView.h>
 #include <TNL/Containers/Segments/EllpackSegmentView.h>
diff --git a/src/UnitTests/Matrices/CMakeLists.txt b/src/UnitTests/Matrices/CMakeLists.txt
index 9b168bd563781c12dab0c148fc9edea1776d33d3..668e272df62e5f31552854bf41000192a10410ff 100644
--- a/src/UnitTests/Matrices/CMakeLists.txt
+++ b/src/UnitTests/Matrices/CMakeLists.txt
@@ -1,3 +1,5 @@
+ADD_SUBDIRECTORY( Legacy )
+
 IF( BUILD_CUDA )
    CUDA_ADD_EXECUTABLE( SparseMatrixCopyTest SparseMatrixCopyTest.cu OPTIONS ${CXX_TESTS_FLAGS} )
    TARGET_LINK_LIBRARIES( SparseMatrixCopyTest ${GTEST_BOTH_LIBRARIES} )
@@ -5,24 +7,6 @@ IF( BUILD_CUDA )
    CUDA_ADD_EXECUTABLE( SparseMatrixTest SparseMatrixTest.cu OPTIONS ${CXX_TESTS_FLAGS} )
    TARGET_LINK_LIBRARIES( SparseMatrixTest ${GTEST_BOTH_LIBRARIES} )
 
-   CUDA_ADD_EXECUTABLE( SparseMatrixTest_AdEllpack SparseMatrixTest_AdEllpack.cu OPTIONS ${CXX_TESTS_FLAGS} )
-   TARGET_LINK_LIBRARIES( SparseMatrixTest_AdEllpack ${GTEST_BOTH_LIBRARIES} )
-
-   CUDA_ADD_EXECUTABLE( SparseMatrixTest_BiEllpack SparseMatrixTest_BiEllpack.cu OPTIONS ${CXX_TESTS_FLAGS} )
-   TARGET_LINK_LIBRARIES( SparseMatrixTest_BiEllpack ${GTEST_BOTH_LIBRARIES} )
-
-   CUDA_ADD_EXECUTABLE( SparseMatrixTest_ChunkedEllpack SparseMatrixTest_ChunkedEllpack.cu OPTIONS ${CXX_TESTS_FLAGS} )
-   TARGET_LINK_LIBRARIES( SparseMatrixTest_ChunkedEllpack ${GTEST_BOTH_LIBRARIES} )
-
-   CUDA_ADD_EXECUTABLE( SparseMatrixTest_CSR SparseMatrixTest_CSR.cu OPTIONS ${CXX_TESTS_FLAGS} )
-   TARGET_LINK_LIBRARIES( SparseMatrixTest_CSR ${GTEST_BOTH_LIBRARIES} )
-
-   CUDA_ADD_EXECUTABLE( SparseMatrixTest_Ellpack SparseMatrixTest_Ellpack.cu OPTIONS ${CXX_TESTS_FLAGS} )
-   TARGET_LINK_LIBRARIES( SparseMatrixTest_Ellpack ${GTEST_BOTH_LIBRARIES} )
-
-   CUDA_ADD_EXECUTABLE( SparseMatrixTest_SlicedEllpack SparseMatrixTest_SlicedEllpack.cu OPTIONS ${CXX_TESTS_FLAGS} )
-   TARGET_LINK_LIBRARIES( SparseMatrixTest_SlicedEllpack ${GTEST_BOTH_LIBRARIES} )
-
    CUDA_ADD_EXECUTABLE( DenseMatrixTest DenseMatrixTest.cu OPTIONS ${CXX_TESTS_FLAGS} )
    TARGET_LINK_LIBRARIES( DenseMatrixTest ${GTEST_BOTH_LIBRARIES} )
 
@@ -46,30 +30,6 @@ ELSE(  BUILD_CUDA )
    TARGET_COMPILE_OPTIONS( SparseMatrixTest PRIVATE ${CXX_TESTS_FLAGS} )
    TARGET_LINK_LIBRARIES( SparseMatrixTest ${GTEST_BOTH_LIBRARIES} )
 
-   ADD_EXECUTABLE( SparseMatrixTest_AdEllpack SparseMatrixTest_AdEllpack.cpp )
-   TARGET_COMPILE_OPTIONS( SparseMatrixTest_AdEllpack PRIVATE ${CXX_TESTS_FLAGS} )
-   TARGET_LINK_LIBRARIES( SparseMatrixTest_AdEllpack ${GTEST_BOTH_LIBRARIES} )
-
-   ADD_EXECUTABLE( SparseMatrixTest_BiEllpack SparseMatrixTest_BiEllpack.cpp )
-   TARGET_COMPILE_OPTIONS( SparseMatrixTest_BiEllpack PRIVATE ${CXX_TESTS_FLAGS} )
-   TARGET_LINK_LIBRARIES( SparseMatrixTest_BiEllpack ${GTEST_BOTH_LIBRARIES} )
-
-   ADD_EXECUTABLE( SparseMatrixTest_ChunkedEllpack SparseMatrixTest_ChunkedEllpack.cpp )
-   TARGET_COMPILE_OPTIONS( SparseMatrixTest_ChunkedEllpack PRIVATE ${CXX_TESTS_FLAGS} )
-   TARGET_LINK_LIBRARIES( SparseMatrixTest_ChunkedEllpack ${GTEST_BOTH_LIBRARIES} )
-
-   ADD_EXECUTABLE( SparseMatrixTest_CSR SparseMatrixTest_CSR.cpp )
-   TARGET_COMPILE_OPTIONS( SparseMatrixTest_CSR PRIVATE ${CXX_TESTS_FLAGS} )
-   TARGET_LINK_LIBRARIES( SparseMatrixTest_CSR ${GTEST_BOTH_LIBRARIES} )
-
-   ADD_EXECUTABLE( SparseMatrixTest_Ellpack SparseMatrixTest_Ellpack.cpp )
-   TARGET_COMPILE_OPTIONS( SparseMatrixTest_Ellpack PRIVATE ${CXX_TESTS_FLAGS} )
-   TARGET_LINK_LIBRARIES( SparseMatrixTest_Ellpack ${GTEST_BOTH_LIBRARIES} )
-
-   ADD_EXECUTABLE( SparseMatrixTest_SlicedEllpack SparseMatrixTest_SlicedEllpack.cpp )
-   TARGET_COMPILE_OPTIONS( SparseMatrixTest_SlicedEllpack PRIVATE ${CXX_TESTS_FLAGS} )
-   TARGET_LINK_LIBRARIES( SparseMatrixTest_SlicedEllpack ${GTEST_BOTH_LIBRARIES} )
-
    ADD_EXECUTABLE( DenseMatrixTest DenseMatrixTest.cpp )
    TARGET_COMPILE_OPTIONS( DenseMatrixTest PRIVATE ${CXX_TESTS_FLAGS} )
    TARGET_LINK_LIBRARIES( DenseMatrixTest ${GTEST_BOTH_LIBRARIES} )
@@ -95,11 +55,6 @@ ADD_TEST( SparseMatrixCopyTest ${EXECUTABLE_OUTPUT_PATH}/SparseMatrixCopyTest${C
 ADD_TEST( SparseMatrixTest ${EXECUTABLE_OUTPUT_PATH}/SparseMatrixTest${CMAKE_EXECUTABLE_SUFFIX} )
 # TODO: Uncomment the following when AdEllpack works
 #ADD_TEST( SparseMatrixTest_AdEllpack ${EXECUTABLE_OUTPUT_PATH}/SparseMatrixTest_AdEllpack${CMAKE_EXECUTABLE_SUFFIX} )
-ADD_TEST( SparseMatrixTest_BiEllpack ${EXECUTABLE_OUTPUT_PATH}/SparseMatrixTest_BiEllpack${CMAKE_EXECUTABLE_SUFFIX} )
-ADD_TEST( SparseMatrixTest_ChunkedEllpack ${EXECUTABLE_OUTPUT_PATH}/SparseMatrixTest_ChunkedEllpack${CMAKE_EXECUTABLE_SUFFIX} )
-ADD_TEST( SparseMatrixTest_CSR ${EXECUTABLE_OUTPUT_PATH}/SparseMatrixTest_CSR${CMAKE_EXECUTABLE_SUFFIX} )
-ADD_TEST( SparseMatrixTest_Ellpack ${EXECUTABLE_OUTPUT_PATH}/SparseMatrixTest_Ellpack${CMAKE_EXECUTABLE_SUFFIX} )
-ADD_TEST( SparseMatrixTest_SlicedEllpack ${EXECUTABLE_OUTPUT_PATH}/SparseMatrixTest_SlicedEllpack${CMAKE_EXECUTABLE_SUFFIX} )
 # TODO: DenseMatrixTest is not finished
 #ADD_TEST( DenseMatrixTest ${EXECUTABLE_OUTPUT_PATH}/DenseMatrixTest${CMAKE_EXECUTABLE_SUFFIX} )
 
diff --git a/src/UnitTests/Matrices/Legacy/CMakeLists.txt b/src/UnitTests/Matrices/Legacy/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..9cdfe27845a4ef7d75f27b644b18d5b17946080f
--- /dev/null
+++ b/src/UnitTests/Matrices/Legacy/CMakeLists.txt
@@ -0,0 +1,72 @@
+IF( BUILD_CUDA )
+   CUDA_ADD_EXECUTABLE( Legacy_SparseMatrixCopyTest SparseMatrixCopyTest.cu OPTIONS ${CXX_TESTS_FLAGS} )
+   TARGET_LINK_LIBRARIES( Legacy_SparseMatrixCopyTest ${GTEST_BOTH_LIBRARIES} )
+
+   CUDA_ADD_EXECUTABLE( Legacy_SparseMatrixTest SparseMatrixTest.cu OPTIONS ${CXX_TESTS_FLAGS} )
+   TARGET_LINK_LIBRARIES( Legacy_SparseMatrixTest ${GTEST_BOTH_LIBRARIES} )
+
+   CUDA_ADD_EXECUTABLE( Legacy_SparseMatrixTest_AdEllpack SparseMatrixTest_AdEllpack.cu OPTIONS ${CXX_TESTS_FLAGS} )
+   TARGET_LINK_LIBRARIES( Legacy_SparseMatrixTest_AdEllpack ${GTEST_BOTH_LIBRARIES} )
+
+   CUDA_ADD_EXECUTABLE( Legacy_SparseMatrixTest_BiEllpack SparseMatrixTest_BiEllpack.cu OPTIONS ${CXX_TESTS_FLAGS} )
+   TARGET_LINK_LIBRARIES( Legacy_SparseMatrixTest_BiEllpack ${GTEST_BOTH_LIBRARIES} )
+
+   CUDA_ADD_EXECUTABLE( Legacy_SparseMatrixTest_ChunkedEllpack SparseMatrixTest_ChunkedEllpack.cu OPTIONS ${CXX_TESTS_FLAGS} )
+   TARGET_LINK_LIBRARIES( Legacy_SparseMatrixTest_ChunkedEllpack ${GTEST_BOTH_LIBRARIES} )
+
+   CUDA_ADD_EXECUTABLE( Legacy_SparseMatrixTest_CSR SparseMatrixTest_CSR.cu OPTIONS ${CXX_TESTS_FLAGS} )
+   TARGET_LINK_LIBRARIES( Legacy_SparseMatrixTest_CSR ${GTEST_BOTH_LIBRARIES} )
+
+   CUDA_ADD_EXECUTABLE( Legacy_SparseMatrixTest_Ellpack SparseMatrixTest_Ellpack.cu OPTIONS ${CXX_TESTS_FLAGS} )
+   TARGET_LINK_LIBRARIES( Legacy_SparseMatrixTest_Ellpack ${GTEST_BOTH_LIBRARIES} )
+
+   CUDA_ADD_EXECUTABLE( Legacy_SparseMatrixTest_SlicedEllpack SparseMatrixTest_SlicedEllpack.cu OPTIONS ${CXX_TESTS_FLAGS} )
+   TARGET_LINK_LIBRARIES( Legacy_SparseMatrixTest_SlicedEllpack ${GTEST_BOTH_LIBRARIES} )
+
+ELSE(  BUILD_CUDA )
+   ADD_EXECUTABLE( Legacy_SparseMatrixCopyTest SparseMatrixCopyTest.cpp )
+   TARGET_COMPILE_OPTIONS( Legacy_SparseMatrixCopyTest PRIVATE ${CXX_TESTS_FLAGS} )
+   TARGET_LINK_LIBRARIES( Legacy_SparseMatrixCopyTest ${GTEST_BOTH_LIBRARIES} )
+
+   ADD_EXECUTABLE( Legacy_SparseMatrixTest SparseMatrixTest.cpp )
+   TARGET_COMPILE_OPTIONS( Legacy_SparseMatrixTest PRIVATE ${CXX_TESTS_FLAGS} )
+   TARGET_LINK_LIBRARIES( Legacy_SparseMatrixTest ${GTEST_BOTH_LIBRARIES} )
+
+   ADD_EXECUTABLE( Legacy_SparseMatrixTest_AdEllpack SparseMatrixTest_AdEllpack.cpp )
+   TARGET_COMPILE_OPTIONS( Legacy_SparseMatrixTest_AdEllpack PRIVATE ${CXX_TESTS_FLAGS} )
+   TARGET_LINK_LIBRARIES( Legacy_SparseMatrixTest_AdEllpack ${GTEST_BOTH_LIBRARIES} )
+
+   ADD_EXECUTABLE( Legacy_SparseMatrixTest_BiEllpack SparseMatrixTest_BiEllpack.cpp )
+   TARGET_COMPILE_OPTIONS( Legacy_SparseMatrixTest_BiEllpack PRIVATE ${CXX_TESTS_FLAGS} )
+   TARGET_LINK_LIBRARIES( Legacy_SparseMatrixTest_BiEllpack ${GTEST_BOTH_LIBRARIES} )
+
+   ADD_EXECUTABLE( Legacy_SparseMatrixTest_ChunkedEllpack SparseMatrixTest_ChunkedEllpack.cpp )
+   TARGET_COMPILE_OPTIONS( Legacy_SparseMatrixTest_ChunkedEllpack PRIVATE ${CXX_TESTS_FLAGS} )
+   TARGET_LINK_LIBRARIES( Legacy_SparseMatrixTest_ChunkedEllpack ${GTEST_BOTH_LIBRARIES} )
+
+   ADD_EXECUTABLE( Legacy_SparseMatrixTest_CSR SparseMatrixTest_CSR.cpp )
+   TARGET_COMPILE_OPTIONS( Legacy_SparseMatrixTest_CSR PRIVATE ${CXX_TESTS_FLAGS} )
+   TARGET_LINK_LIBRARIES( Legacy_SparseMatrixTest_CSR ${GTEST_BOTH_LIBRARIES} )
+
+   ADD_EXECUTABLE( Legacy_SparseMatrixTest_Ellpack SparseMatrixTest_Ellpack.cpp )
+   TARGET_COMPILE_OPTIONS( Legacy_SparseMatrixTest_Ellpack PRIVATE ${CXX_TESTS_FLAGS} )
+   TARGET_LINK_LIBRARIES( Legacy_SparseMatrixTest_Ellpack ${GTEST_BOTH_LIBRARIES} )
+
+   ADD_EXECUTABLE( Legacy_SparseMatrixTest_SlicedEllpack SparseMatrixTest_SlicedEllpack.cpp )
+   TARGET_COMPILE_OPTIONS( Legacy_SparseMatrixTest_SlicedEllpack PRIVATE ${CXX_TESTS_FLAGS} )
+   TARGET_LINK_LIBRARIES( Legacy_SparseMatrixTest_SlicedEllpack ${GTEST_BOTH_LIBRARIES} )
+
+ENDIF( BUILD_CUDA )
+
+
+ADD_TEST( Legacy_SparseMatrixCopyTest ${EXECUTABLE_OUTPUT_PATH}/SparseMatrixCopyTest${CMAKE_EXECUTABLE_SUFFIX} )
+ADD_TEST( Legacy_SparseMatrixTest ${EXECUTABLE_OUTPUT_PATH}/SparseMatrixTest${CMAKE_EXECUTABLE_SUFFIX} )
+# TODO: Uncomment the following when AdEllpack works
+#ADD_TEST( SparseMatrixTest_AdEllpack ${EXECUTABLE_OUTPUT_PATH}/SparseMatrixTest_AdEllpack${CMAKE_EXECUTABLE_SUFFIX} )
+ADD_TEST( Legacy_SparseMatrixTest_BiEllpack ${EXECUTABLE_OUTPUT_PATH}/SparseMatrixTest_BiEllpack${CMAKE_EXECUTABLE_SUFFIX} )
+ADD_TEST( Legacy_SparseMatrixTest_ChunkedEllpack ${EXECUTABLE_OUTPUT_PATH}/SparseMatrixTest_ChunkedEllpack${CMAKE_EXECUTABLE_SUFFIX} )
+ADD_TEST( Legacy_SparseMatrixTest_CSR ${EXECUTABLE_OUTPUT_PATH}/SparseMatrixTest_CSR${CMAKE_EXECUTABLE_SUFFIX} )
+ADD_TEST( Legacy_SparseMatrixTest_Ellpack ${EXECUTABLE_OUTPUT_PATH}/SparseMatrixTest_Ellpack${CMAKE_EXECUTABLE_SUFFIX} )
+ADD_TEST( Legacy_SparseMatrixTest_SlicedEllpack ${EXECUTABLE_OUTPUT_PATH}/SparseMatrixTest_SlicedEllpack${CMAKE_EXECUTABLE_SUFFIX} )
+# TODO: DenseMatrixTest is not finished
+#ADD_TEST( DenseMatrixTest ${EXECUTABLE_OUTPUT_PATH}/DenseMatrixTest${CMAKE_EXECUTABLE_SUFFIX} )
diff --git a/src/UnitTests/Matrices/Legacy/SparseMatrixCopyTest.cpp b/src/UnitTests/Matrices/Legacy/SparseMatrixCopyTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..30b8f64ecfdbf228856d272a71d3de08980f3987
--- /dev/null
+++ b/src/UnitTests/Matrices/Legacy/SparseMatrixCopyTest.cpp
@@ -0,0 +1,11 @@
+/***************************************************************************
+                          SparseMatrixCopyTest.cpp  -  description
+                             -------------------
+    begin                : Jun 25, 2017
+    copyright            : (C) 2017 by Tomas Oberhuber et al.
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#include "SparseMatrixCopyTest.h"
diff --git a/src/UnitTests/Matrices/Legacy/SparseMatrixCopyTest.cu b/src/UnitTests/Matrices/Legacy/SparseMatrixCopyTest.cu
new file mode 100644
index 0000000000000000000000000000000000000000..431fe481c2db1d5b18cfa849e882c0ed836463c1
--- /dev/null
+++ b/src/UnitTests/Matrices/Legacy/SparseMatrixCopyTest.cu
@@ -0,0 +1,11 @@
+/***************************************************************************
+                          SparseMatrixCopyTest.cu  -  description
+                             -------------------
+    begin                : Jun 25, 2017
+    copyright            : (C) 2017 by Tomas Oberhuber et al.
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#include "SparseMatrixCopyTest.h"
diff --git a/src/UnitTests/Matrices/Legacy/SparseMatrixCopyTest.h b/src/UnitTests/Matrices/Legacy/SparseMatrixCopyTest.h
new file mode 100644
index 0000000000000000000000000000000000000000..7069fd77740473247cf6f04a26bf1f5ec5ff4670
--- /dev/null
+++ b/src/UnitTests/Matrices/Legacy/SparseMatrixCopyTest.h
@@ -0,0 +1,573 @@
+/***************************************************************************
+                          SparseMatrixCopyTest.h -  description
+                             -------------------
+    begin                : Jun 25, 2017
+    copyright            : (C) 2017 by Tomas Oberhuber et al.
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#include <TNL/Matrices/Legacy/CSR.h>
+#include <TNL/Matrices/Legacy/Ellpack.h>
+#include <TNL/Matrices/Legacy/SlicedEllpack.h>
+
+#include <TNL/Matrices/SparseMatrix.h>
+#include <TNL/Matrices/MatrixType.h>
+#include <TNL/Containers/Segments/CSR.h>
+#include <TNL/Containers/Segments/Ellpack.h>
+#include <TNL/Containers/Segments/SlicedEllpack.h>
+
+/*using CSR_host = TNL::Matrices::CSR< int, TNL::Devices::Host, int >;
+using CSR_cuda = TNL::Matrices::CSR< int, TNL::Devices::Cuda, int >;
+using E_host = TNL::Matrices::Ellpack< int, TNL::Devices::Host, int >;
+using E_cuda = TNL::Matrices::Ellpack< int, TNL::Devices::Cuda, int >;
+using SE_host = TNL::Matrices::SlicedEllpack< int, TNL::Devices::Host, int, 2 >;
+using SE_cuda = TNL::Matrices::SlicedEllpack< int, TNL::Devices::Cuda, int, 2 >;*/
+
+template< typename Device, typename Index, typename IndexAllocator >
+using EllpackSegments = TNL::Containers::Segments::Ellpack< Device, Index, IndexAllocator >;
+
+template< typename Device, typename Index, typename IndexAllocator >
+using SlicedEllpackSegments = TNL::Containers::Segments::SlicedEllpack< Device, Index, IndexAllocator >;
+
+using CSR_host = TNL::Matrices::SparseMatrix< int, TNL::Devices::Host, int, TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >;
+using CSR_cuda = TNL::Matrices::SparseMatrix< int, TNL::Devices::Cuda, int, TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >;
+using E_host   = TNL::Matrices::SparseMatrix< int, TNL::Devices::Host, int, TNL::Matrices::GeneralMatrix, EllpackSegments >;
+using E_cuda   = TNL::Matrices::SparseMatrix< int, TNL::Devices::Cuda, int, TNL::Matrices::GeneralMatrix, EllpackSegments >;
+using SE_host  = TNL::Matrices::SparseMatrix< int, TNL::Devices::Host, int, TNL::Matrices::GeneralMatrix, SlicedEllpackSegments >;
+using SE_cuda  = TNL::Matrices::SparseMatrix< int, TNL::Devices::Cuda, int, TNL::Matrices::GeneralMatrix, SlicedEllpackSegments >;
+
+
+#ifdef HAVE_GTEST 
+#include <gtest/gtest.h>
+
+/*
+ * Sets up the following 10x6 sparse matrix:
+ *
+ *    /  1  2             \
+ *    |           3  4  5 |
+ *    |  6  7  8          |
+ *    |     9 10 11 12 13 |
+ *    | 14 15 16 17 18    |
+ *    | 19 20             |
+ *    | 21                |
+ *    | 22                |
+ *    | 23 24 25 26 27    |
+ *    \                28 /
+ */
+template< typename Matrix >
+void setupUnevenRowSizeMatrix( Matrix& m )
+{
+    const int rows = 10;
+    const int cols = 6;
+    m.reset();
+    m.setDimensions( rows, cols );
+    typename Matrix::CompressedRowLengthsVector rowLengths;
+    rowLengths.setSize( rows );
+    rowLengths.setValue( 5 );
+    rowLengths.setElement( 0, 2 );
+    rowLengths.setElement( 1,  3 );
+    rowLengths.setElement( 2,  3 );
+    rowLengths.setElement( 5,  2 );
+    rowLengths.setElement( 6,  1 );
+    rowLengths.setElement( 7,  1 );
+    rowLengths.setElement( 9,  1 );
+    m.setCompressedRowLengths( rowLengths );
+
+    int value = 1;
+    for( int i = 0; i < cols - 4; i++ )  // 0th row
+        m.setElement( 0, i, value++ );
+
+    for( int i = 3; i < cols; i++ )      // 1st row
+        m.setElement( 1, i, value++ );
+
+    for( int i = 0; i < cols - 3; i++ )  // 2nd row
+        m.setElement( 2, i, value++ );
+
+    for( int i = 1; i < cols; i++ )      // 3rd row
+        m.setElement( 3, i, value++ );
+
+    for( int i = 0; i < cols - 1; i++ )  // 4th row
+        m.setElement( 4, i, value++ );
+
+    for( int i = 0; i < cols - 4; i++ )  // 5th row
+        m.setElement( 5, i, value++ );
+
+    m.setElement( 6, 0, value++ );   // 6th row
+
+    m.setElement( 7, 0, value++ );   // 7th row
+
+    for( int i = 0; i < cols - 1; i++ )  // 8th row 
+        m.setElement( 8, i, value++ );
+
+    m.setElement( 9, 5, value++ );   // 9th row
+}
+
+template< typename Matrix >
+void checkUnevenRowSizeMatrix( Matrix& m )
+{
+   ASSERT_EQ( m.getRows(), 10 );
+   ASSERT_EQ( m.getColumns(), 6 );
+
+   EXPECT_EQ( m.getElement( 0, 0 ),  1 );
+   EXPECT_EQ( m.getElement( 0, 1 ),  2 );
+   EXPECT_EQ( m.getElement( 0, 2 ),  0 );
+   EXPECT_EQ( m.getElement( 0, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 0, 4 ),  0 );
+   EXPECT_EQ( m.getElement( 0, 5 ),  0);
+
+   EXPECT_EQ( m.getElement( 1, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 1, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 1, 2 ),  0 );
+   EXPECT_EQ( m.getElement( 1, 3 ),  3 );
+   EXPECT_EQ( m.getElement( 1, 4 ),  4 );
+   EXPECT_EQ( m.getElement( 1, 5 ),  5 );
+
+   EXPECT_EQ( m.getElement( 2, 0 ),  6 );
+   EXPECT_EQ( m.getElement( 2, 1 ),  7 );
+   EXPECT_EQ( m.getElement( 2, 2 ),  8 );
+   EXPECT_EQ( m.getElement( 2, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 2, 4 ),  0 );
+   EXPECT_EQ( m.getElement( 2, 5 ),  0 );
+
+   EXPECT_EQ( m.getElement( 3, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 3, 1 ),  9 );
+   EXPECT_EQ( m.getElement( 3, 2 ), 10 );
+   EXPECT_EQ( m.getElement( 3, 3 ), 11 );
+   EXPECT_EQ( m.getElement( 3, 4 ), 12 );
+   EXPECT_EQ( m.getElement( 3, 5 ), 13 );
+
+   EXPECT_EQ( m.getElement( 4, 0 ), 14 );
+   EXPECT_EQ( m.getElement( 4, 1 ), 15 );
+   EXPECT_EQ( m.getElement( 4, 2 ), 16 );
+   EXPECT_EQ( m.getElement( 4, 3 ), 17 );
+   EXPECT_EQ( m.getElement( 4, 4 ), 18 );
+   EXPECT_EQ( m.getElement( 4, 5 ),  0 );
+
+   EXPECT_EQ( m.getElement( 5, 0 ), 19 );
+   EXPECT_EQ( m.getElement( 5, 1 ), 20 );
+   EXPECT_EQ( m.getElement( 5, 2 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 4 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 5 ),  0 );
+
+   EXPECT_EQ( m.getElement( 6, 0 ), 21 );
+   EXPECT_EQ( m.getElement( 6, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 6, 2 ),  0 );
+   EXPECT_EQ( m.getElement( 6, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 6, 4 ),  0 );
+   EXPECT_EQ( m.getElement( 6, 5 ),  0 );
+   
+   EXPECT_EQ( m.getElement( 7, 0 ), 22 );
+   EXPECT_EQ( m.getElement( 7, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 7, 2 ),  0 );
+   EXPECT_EQ( m.getElement( 7, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 7, 4 ),  0 );
+   EXPECT_EQ( m.getElement( 7, 5 ),  0 );
+   
+   EXPECT_EQ( m.getElement( 8, 0 ), 23 );
+   EXPECT_EQ( m.getElement( 8, 1 ), 24 );
+   EXPECT_EQ( m.getElement( 8, 2 ), 25 );
+   EXPECT_EQ( m.getElement( 8, 3 ), 26 );
+   EXPECT_EQ( m.getElement( 8, 4 ), 27 );
+   EXPECT_EQ( m.getElement( 8, 5 ),  0 );
+   
+   EXPECT_EQ( m.getElement( 9, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 9, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 9, 2 ),  0 );
+   EXPECT_EQ( m.getElement( 9, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 9, 4 ),  0 );
+   EXPECT_EQ( m.getElement( 9, 5 ), 28 );
+}
+
+/*
+ * Sets up the following 7x6 sparse matrix:
+ *
+ *    /              2  1 \
+ *    |           5  4  3 |
+ *    |        8  7  6    |
+ *    |    11 10  9       |
+ *    | 14 13 12          |
+ *    | 16 15             |
+ *    \ 17                /
+ */
+template< typename Matrix >
+void setupAntiTriDiagMatrix( Matrix& m )
+{
+    const int rows = 7;
+    const int cols = 6;
+    m.reset();
+    m.setDimensions( rows, cols );
+    typename Matrix::CompressedRowLengthsVector rowLengths;
+    rowLengths.setSize( rows );
+    rowLengths.setValue( 3 );
+    rowLengths.setElement( 0, 4);
+    rowLengths.setElement( 1,  4 );
+    m.setCompressedRowLengths( rowLengths );
+    
+    int value = 1;
+    for( int i = 0; i < rows; i++ )
+        for( int j = cols - 1; j > 2; j-- )
+            if( j - i + 1 < cols && j - i + 1 >= 0 )
+                m.setElement( i, j - i + 1, value++ );
+}
+
+template< typename Matrix >
+void checkAntiTriDiagMatrix( Matrix& m )
+{
+   ASSERT_EQ( m.getRows(), 7 );
+   ASSERT_EQ( m.getColumns(), 6 );
+
+   EXPECT_EQ( m.getElement( 0, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 0, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 0, 2 ),  0 );
+   EXPECT_EQ( m.getElement( 0, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 0, 4 ),  2 );
+   EXPECT_EQ( m.getElement( 0, 5 ),  1);
+
+   EXPECT_EQ( m.getElement( 1, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 1, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 1, 2 ),  0 );
+   EXPECT_EQ( m.getElement( 1, 3 ),  5 );
+   EXPECT_EQ( m.getElement( 1, 4 ),  4 );
+   EXPECT_EQ( m.getElement( 1, 5 ),  3 );
+
+   EXPECT_EQ( m.getElement( 2, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 2, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 2, 2 ),  8 );
+   EXPECT_EQ( m.getElement( 2, 3 ),  7 );
+   EXPECT_EQ( m.getElement( 2, 4 ),  6 );
+   EXPECT_EQ( m.getElement( 2, 5 ),  0 );
+
+   EXPECT_EQ( m.getElement( 3, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 3, 1 ), 11 );
+   EXPECT_EQ( m.getElement( 3, 2 ), 10 );
+   EXPECT_EQ( m.getElement( 3, 3 ),  9 );
+   EXPECT_EQ( m.getElement( 3, 4 ),  0 );
+   EXPECT_EQ( m.getElement( 3, 5 ),  0 );
+
+   EXPECT_EQ( m.getElement( 4, 0 ), 14 );
+   EXPECT_EQ( m.getElement( 4, 1 ), 13 );
+   EXPECT_EQ( m.getElement( 4, 2 ), 12 );
+   EXPECT_EQ( m.getElement( 4, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 4, 4 ),  0 );
+   EXPECT_EQ( m.getElement( 4, 5 ),  0 );
+
+   EXPECT_EQ( m.getElement( 5, 0 ), 16 );
+   EXPECT_EQ( m.getElement( 5, 1 ), 15 );
+   EXPECT_EQ( m.getElement( 5, 2 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 4 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 5 ),  0 );
+
+   EXPECT_EQ( m.getElement( 6, 0 ), 17 );
+   EXPECT_EQ( m.getElement( 6, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 6, 2 ),  0 );
+   EXPECT_EQ( m.getElement( 6, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 6, 4 ),  0 );
+   EXPECT_EQ( m.getElement( 6, 5 ),  0 );
+}
+
+/*
+ * Sets up the following 7x6 sparse matrix:
+ *
+ *    / 1  2             \
+ *    | 3  4  5          |
+ *    |    6  7  8       |
+ *    |       9 10 11    |
+ *    |         12 13 14 |
+ *    |            15 16 |
+ *    \               17 /
+ */
+template< typename Matrix >
+void setupTriDiagMatrix( Matrix& m )
+{
+   const int rows = 7;
+   const int cols = 6;
+   m.reset();
+   m.setDimensions( rows, cols );
+   typename Matrix::CompressedRowLengthsVector rowLengths;
+   rowLengths.setSize( rows );
+   rowLengths.setValue( 3 );
+   rowLengths.setElement( 0 , 4 );
+   rowLengths.setElement( 1,  4 );
+   m.setCompressedRowLengths( rowLengths );
+
+   int value = 1;
+   for( int i = 0; i < rows; i++ )
+      for( int j = 0; j < 3; j++ )
+         if( i + j - 1 >= 0 && i + j - 1 < cols )
+            m.setElement( i, i + j - 1, value++ );
+}
+
+template< typename Matrix >
+void checkTriDiagMatrix( Matrix& m )
+{
+   ASSERT_EQ( m.getRows(), 7 );
+   ASSERT_EQ( m.getColumns(), 6 );
+
+   EXPECT_EQ( m.getElement( 0, 0 ),  1 );
+   EXPECT_EQ( m.getElement( 0, 1 ),  2 );
+   EXPECT_EQ( m.getElement( 0, 2 ),  0 );
+   EXPECT_EQ( m.getElement( 0, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 0, 4 ),  0 );
+   EXPECT_EQ( m.getElement( 0, 5 ),  0 );
+
+   EXPECT_EQ( m.getElement( 1, 0 ),  3 );
+   EXPECT_EQ( m.getElement( 1, 1 ),  4 );
+   EXPECT_EQ( m.getElement( 1, 2 ),  5 );
+   EXPECT_EQ( m.getElement( 1, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 1, 4 ),  0 );
+   EXPECT_EQ( m.getElement( 1, 5 ),  0 );
+
+   EXPECT_EQ( m.getElement( 2, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 2, 1 ),  6 );
+   EXPECT_EQ( m.getElement( 2, 2 ),  7 );
+   EXPECT_EQ( m.getElement( 2, 3 ),  8 );
+   EXPECT_EQ( m.getElement( 2, 4 ),  0 );
+   EXPECT_EQ( m.getElement( 2, 5 ),  0 );
+
+   EXPECT_EQ( m.getElement( 3, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 3, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 3, 2 ),  9 );
+   EXPECT_EQ( m.getElement( 3, 3 ), 10 );
+   EXPECT_EQ( m.getElement( 3, 4 ), 11 );
+   EXPECT_EQ( m.getElement( 3, 5 ),  0 );
+
+   EXPECT_EQ( m.getElement( 4, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 4, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 4, 2 ),  0 );
+   EXPECT_EQ( m.getElement( 4, 3 ), 12 );
+   EXPECT_EQ( m.getElement( 4, 4 ), 13 );
+   EXPECT_EQ( m.getElement( 4, 5 ), 14 );
+
+   EXPECT_EQ( m.getElement( 5, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 2 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 4 ), 15 );
+   EXPECT_EQ( m.getElement( 5, 5 ), 16 );
+
+   EXPECT_EQ( m.getElement( 6, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 6, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 6, 2 ),  0 );
+   EXPECT_EQ( m.getElement( 6, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 6, 4 ),  0 );
+   EXPECT_EQ( m.getElement( 6, 5 ), 17 );
+}
+
+template< typename Matrix1, typename Matrix2 >
+void testCopyAssignment()
+{
+   {
+      SCOPED_TRACE("Tri Diagonal Matrix");
+
+      Matrix1 triDiag1;
+      setupTriDiagMatrix( triDiag1 );
+      checkTriDiagMatrix( triDiag1 );
+
+      Matrix2 triDiag2;
+      triDiag2 = triDiag1;
+      checkTriDiagMatrix( triDiag2 );
+   }
+   {
+      SCOPED_TRACE("Anti Tri Diagonal Matrix");
+      Matrix1 antiTriDiag1;
+      setupAntiTriDiagMatrix( antiTriDiag1 );
+      checkAntiTriDiagMatrix( antiTriDiag1 );
+
+      Matrix2 antiTriDiag2;
+      antiTriDiag2 = antiTriDiag1;
+      checkAntiTriDiagMatrix( antiTriDiag2 );
+   }
+   {
+      SCOPED_TRACE("Uneven Row Size Matrix");
+      Matrix1 unevenRowSize1;
+      setupUnevenRowSizeMatrix( unevenRowSize1 );
+      checkUnevenRowSizeMatrix( unevenRowSize1 );
+
+      Matrix2 unevenRowSize2;
+      unevenRowSize2 = unevenRowSize1;
+      checkUnevenRowSizeMatrix( unevenRowSize2 );
+   }
+}
+
+template< typename Matrix1, typename Matrix2 >
+void testConversion()
+{
+    
+   {
+        SCOPED_TRACE("Tri Diagonal Matrix");
+        
+        Matrix1 triDiag1;
+        setupTriDiagMatrix( triDiag1 );
+        checkTriDiagMatrix( triDiag1 );
+        
+        Matrix2 triDiag2;
+        //TNL::Matrices::copySparseMatrix( triDiag2, triDiag1 );
+        triDiag2 = triDiag1;
+        checkTriDiagMatrix( triDiag2 );
+   }
+   
+   {
+        SCOPED_TRACE("Anti Tri Diagonal Matrix");
+                
+        Matrix1 antiTriDiag1;
+        setupAntiTriDiagMatrix( antiTriDiag1 );
+        checkAntiTriDiagMatrix( antiTriDiag1 );
+        
+        Matrix2 antiTriDiag2;
+        //TNL::Matrices::copySparseMatrix( antiTriDiag2, antiTriDiag1 );
+        antiTriDiag2 = antiTriDiag1;
+        checkAntiTriDiagMatrix( antiTriDiag2 );
+   }
+   
+   {
+        SCOPED_TRACE("Uneven Row Size Matrix");
+        Matrix1 unevenRowSize1;
+        setupUnevenRowSizeMatrix( unevenRowSize1 );
+        checkUnevenRowSizeMatrix( unevenRowSize1 );
+        
+        Matrix2 unevenRowSize2;
+        //TNL::Matrices::copySparseMatrix( unevenRowSize2, unevenRowSize1 );
+        unevenRowSize2 = unevenRowSize1;
+        checkUnevenRowSizeMatrix( unevenRowSize2 );
+   }
+}
+
+TEST( SparseMatrixCopyTest, CSR_HostToHost )
+{
+   testCopyAssignment< CSR_host, CSR_host >();
+}
+
+#ifdef HAVE_CUDA
+TEST( SparseMatrixCopyTest, CSR_HostToCuda )
+{
+   testCopyAssignment< CSR_host, CSR_cuda >();
+}
+
+TEST( SparseMatrixCopyTest, CSR_CudaToHost )
+{
+   testCopyAssignment< CSR_cuda, CSR_host >();
+}
+
+TEST( SparseMatrixCopyTest, CSR_CudaToCuda )
+{
+   testCopyAssignment< CSR_cuda, CSR_cuda >();
+}
+#endif
+
+
+TEST( SparseMatrixCopyTest, Ellpack_HostToHost )
+{
+   testCopyAssignment< E_host, E_host >();
+}
+
+#ifdef HAVE_CUDA
+TEST( SparseMatrixCopyTest, Ellpack_HostToCuda )
+{
+   testCopyAssignment< E_host, E_cuda >();
+}
+
+TEST( SparseMatrixCopyTest, Ellpack_CudaToHost )
+{
+   testCopyAssignment< E_cuda, E_host >();
+}
+
+TEST( SparseMatrixCopyTest, Ellpack_CudaToCuda )
+{
+   testCopyAssignment< E_cuda, E_cuda >();
+}
+#endif
+
+
+TEST( SparseMatrixCopyTest, SlicedEllpack_HostToHost )
+{
+   testCopyAssignment< SE_host, SE_host >();
+}
+
+#ifdef HAVE_CUDA
+TEST( SparseMatrixCopyTest, SlicedEllpack_HostToCuda )
+{
+   testCopyAssignment< SE_host, SE_cuda >();
+}
+
+TEST( SparseMatrixCopyTest, SlicedEllpack_CudaToHost )
+{
+   testCopyAssignment< SE_cuda, SE_host >();
+}
+
+TEST( SparseMatrixCopyTest, SlicedEllpack_CudaToCuda )
+{
+   testCopyAssignment< SE_cuda, SE_cuda >();
+}
+#endif
+
+
+// test conversion between formats
+TEST( SparseMatrixCopyTest, CSR_to_Ellpack_host )
+{
+   testConversion< CSR_host, E_host >();
+}
+
+TEST( SparseMatrixCopyTest, Ellpack_to_CSR_host )
+{
+   testConversion< E_host, CSR_host >();
+}
+
+TEST( SparseMatrixCopyTest, CSR_to_SlicedEllpack_host )
+{
+   testConversion< CSR_host, SE_host >();
+}
+
+TEST( SparseMatrixCopyTest, SlicedEllpack_to_CSR_host )
+{
+   testConversion< SE_host, CSR_host >();
+}
+
+TEST( SparseMatrixCopyTest, Ellpack_to_SlicedEllpack_host )
+{
+   testConversion< E_host, SE_host >();
+}
+
+TEST( SparseMatrixCopyTest, SlicedEllpack_to_Ellpack_host )
+{
+   testConversion< SE_host, E_host >();
+}
+
+#ifdef HAVE_CUDA
+TEST( SparseMatrixCopyTest, CSR_to_Ellpack_cuda )
+{
+   testConversion< CSR_cuda, E_cuda >();
+}
+
+TEST( SparseMatrixCopyTest, Ellpack_to_CSR_cuda )
+{
+   testConversion< E_cuda, CSR_cuda >();
+}
+
+TEST( SparseMatrixCopyTest, CSR_to_SlicedEllpack_cuda )
+{
+   testConversion< CSR_cuda, SE_cuda >();
+}
+
+TEST( SparseMatrixCopyTest, SlicedEllpack_to_CSR_cuda )
+{
+   testConversion< SE_cuda, CSR_cuda >();
+}
+
+TEST( SparseMatrixCopyTest, Ellpack_to_SlicedEllpack_cuda )
+{
+   testConversion< E_cuda, SE_cuda >();
+}
+
+TEST( SparseMatrixCopyTest, SlicedEllpack_to_Ellpack_cuda )
+{
+   testConversion< SE_cuda, E_cuda >();
+}
+#endif
+
+#endif
+
+#include "../../main.h"
diff --git a/src/UnitTests/Matrices/Legacy/SparseMatrixTest.cpp b/src/UnitTests/Matrices/Legacy/SparseMatrixTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..46f6b9bd3531652feb15f10d60e5736ddfb627fa
--- /dev/null
+++ b/src/UnitTests/Matrices/Legacy/SparseMatrixTest.cpp
@@ -0,0 +1,11 @@
+/***************************************************************************
+                          SparseMatrixTest.cpp -  description
+                             -------------------
+    begin                : Nov 2, 2018
+    copyright            : (C) 2018 by Tomas Oberhuber et al.
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#include "SparseMatrixTest.h"
\ No newline at end of file
diff --git a/src/UnitTests/Matrices/Legacy/SparseMatrixTest.cu b/src/UnitTests/Matrices/Legacy/SparseMatrixTest.cu
new file mode 100644
index 0000000000000000000000000000000000000000..01c23c1937b6ea39c2c99647207d74297ea588c9
--- /dev/null
+++ b/src/UnitTests/Matrices/Legacy/SparseMatrixTest.cu
@@ -0,0 +1,11 @@
+/***************************************************************************
+                          SparseMatrixTest.cu -  description
+                             -------------------
+    begin                : Nov 2, 2018
+    copyright            : (C) 2018 by Tomas Oberhuber et al.
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#include "SparseMatrixTest.h"
\ No newline at end of file
diff --git a/src/UnitTests/Matrices/Legacy/SparseMatrixTest.h b/src/UnitTests/Matrices/Legacy/SparseMatrixTest.h
new file mode 100644
index 0000000000000000000000000000000000000000..ed8bec79603b41e291246428cb59b9a040a56744
--- /dev/null
+++ b/src/UnitTests/Matrices/Legacy/SparseMatrixTest.h
@@ -0,0 +1,39 @@
+/***************************************************************************
+                          SparseMatrixTest.h -  description
+                             -------------------
+    begin                : Nov 2, 2018
+    copyright            : (C) 2018 by Tomas Oberhuber et al.
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#include <TNL/Matrices/Legacy/CSR.h>
+
+#include "SparseMatrixTest.hpp"
+#include <iostream>
+
+#ifdef HAVE_GTEST 
+#include <gtest/gtest.h>
+
+using CSR_host_float = TNL::Matrices::CSR< float, TNL::Devices::Host, int >;
+using CSR_host_int = TNL::Matrices::CSR< int, TNL::Devices::Host, int >;
+
+using CSR_cuda_float = TNL::Matrices::CSR< float, TNL::Devices::Cuda, int >;
+using CSR_cuda_int = TNL::Matrices::CSR< int, TNL::Devices::Cuda, int >;
+
+TEST( SparseMatrixTest, CSR_perforSORIterationTest_Host )
+{
+    test_PerformSORIteration< CSR_host_float >();
+}
+
+#ifdef HAVE_CUDA
+TEST( SparseMatrixTest, CSR_perforSORIterationTest_Cuda )
+{
+   //    test_PerformSORIteration< CSR_cuda_float >();
+}
+#endif
+
+#endif
+
+#include "../../main.h"
diff --git a/src/UnitTests/Matrices/Legacy/SparseMatrixTest.hpp b/src/UnitTests/Matrices/Legacy/SparseMatrixTest.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..c6ff5cbd7349ed52e65d794b3a4df0c7915ba8e6
--- /dev/null
+++ b/src/UnitTests/Matrices/Legacy/SparseMatrixTest.hpp
@@ -0,0 +1,1831 @@
+/***************************************************************************
+                          SparseMatrixTest_impl.h -  description
+                             -------------------
+    begin                : Nov 22, 2018
+    copyright            : (C) 2018 by Tomas Oberhuber et al.
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#include <TNL/Containers/Vector.h>
+#include <TNL/Containers/VectorView.h>
+#include <TNL/Math.h>
+#include <TNL/Algorithms/ParallelFor.h>
+#include <iostream>
+
+// Temporary, until test_OperatorEquals doesn't work for all formats.
+#include <TNL/Matrices/Legacy/ChunkedEllpack.h>
+#include <TNL/Matrices/Legacy/AdEllpack.h>
+#include <TNL/Matrices/Legacy/BiEllpack.h>
+
+#ifdef HAVE_GTEST
+#include <gtest/gtest.h>
+
+template< typename MatrixHostFloat, typename MatrixHostInt >
+void host_test_GetType()
+{
+    bool testRan = false;
+    EXPECT_TRUE( testRan );
+    std::cout << "\nTEST DID NOT RUN. NOT WORKING.\n\n";
+    std::cerr << "This test has not been implemented properly yet.\n" << std::endl;
+}
+
+template< typename MatrixCudaFloat, typename MatrixCudaInt >
+void cuda_test_GetType()
+{
+    bool testRan = false;
+    EXPECT_TRUE( testRan );
+    std::cout << "\nTEST DID NOT RUN. NOT WORKING.\n\n";
+    std::cerr << "This test has not been implemented properly yet.\n" << std::endl;
+}
+
+template< typename Matrix >
+void test_SetDimensions()
+{
+    using RealType = typename Matrix::RealType;
+    using DeviceType = typename Matrix::DeviceType;
+    using IndexType = typename Matrix::IndexType;
+
+    const IndexType rows = 9;
+    const IndexType cols = 8;
+
+    Matrix m;
+    m.setDimensions( rows, cols );
+
+    EXPECT_EQ( m.getRows(), 9 );
+    EXPECT_EQ( m.getColumns(), 8 );
+}
+
+template< typename Matrix >
+void test_SetCompressedRowLengths()
+{
+    using RealType = typename Matrix::RealType;
+    using DeviceType = typename Matrix::DeviceType;
+    using IndexType = typename Matrix::IndexType;
+
+    const IndexType rows = 10;
+    const IndexType cols = 11;
+
+    Matrix m;
+    m.reset();
+    m.setDimensions( rows, cols );
+    typename Matrix::CompressedRowLengthsVector rowLengths;
+    rowLengths.setSize( rows );
+    rowLengths.setValue( 3 );
+
+    IndexType rowLength = 1;
+    for( IndexType i = 2; i < rows; i++ )
+        rowLengths.setElement( i, rowLength++ );
+
+    m.setCompressedRowLengths( rowLengths );
+
+    // Insert values into the rows.
+    RealType value = 1;
+
+    for( IndexType i = 0; i < 3; i++ )      // 0th row
+        m.setElement( 0, i, value++ );
+
+    for( IndexType i = 0; i < 3; i++ )      // 1st row
+        m.setElement( 1, i, value++ );
+
+    for( IndexType i = 0; i < 1; i++ )      // 2nd row
+        m.setElement( 2, i, value++ );
+
+    for( IndexType i = 0; i < 2; i++ )      // 3rd row
+        m.setElement( 3, i, value++ );
+
+    for( IndexType i = 0; i < 3; i++ )      // 4th row
+        m.setElement( 4, i, value++ );
+
+    for( IndexType i = 0; i < 4; i++ )      // 5th row
+        m.setElement( 5, i, value++ );
+
+    for( IndexType i = 0; i < 5; i++ )      // 6th row
+        m.setElement( 6, i, value++ );
+
+    for( IndexType i = 0; i < 6; i++ )      // 7th row
+        m.setElement( 7, i, value++ );
+
+    for( IndexType i = 0; i < 7; i++ )      // 8th row
+        m.setElement( 8, i, value++ );
+
+    for( IndexType i = 0; i < 8; i++ )      // 9th row
+        m.setElement( 9, i, value++ );
+
+
+    EXPECT_EQ( m.getNonZeroRowLength( 0 ), 3 );
+    EXPECT_EQ( m.getNonZeroRowLength( 1 ), 3 );
+    EXPECT_EQ( m.getNonZeroRowLength( 2 ), 1 );
+    EXPECT_EQ( m.getNonZeroRowLength( 3 ), 2 );
+    EXPECT_EQ( m.getNonZeroRowLength( 4 ), 3 );
+    EXPECT_EQ( m.getNonZeroRowLength( 5 ), 4 );
+    EXPECT_EQ( m.getNonZeroRowLength( 6 ), 5 );
+    EXPECT_EQ( m.getNonZeroRowLength( 7 ), 6 );
+    EXPECT_EQ( m.getNonZeroRowLength( 8 ), 7 );
+    EXPECT_EQ( m.getNonZeroRowLength( 9 ), 8 );
+}
+
+template< typename Matrix1, typename Matrix2 >
+void test_SetLike()
+{
+    using RealType = typename Matrix1::RealType;
+    using DeviceType = typename Matrix1::DeviceType;
+    using IndexType = typename Matrix1::IndexType;
+
+    const IndexType rows = 8;
+    const IndexType cols = 7;
+
+    Matrix1 m1;
+    m1.reset();
+    m1.setDimensions( rows + 1, cols + 2 );
+
+    Matrix2 m2;
+    m2.reset();
+    m2.setDimensions( rows, cols );
+
+    m1.setLike( m2 );
+
+
+    EXPECT_EQ( m1.getRows(), m2.getRows() );
+    EXPECT_EQ( m1.getColumns(), m2.getColumns() );
+}
+
+template< typename Matrix >
+void test_GetNumberOfNonzeroMatrixElements()
+{
+   using RealType = typename Matrix::RealType;
+   using DeviceType = typename Matrix::DeviceType;
+   using IndexType = typename Matrix::IndexType;
+
+   /*
+    * Sets up the following 10x10 sparse matrix:
+    *
+    *    /  1  0  2  0  3  0  4  0  0  0  \
+    *    |  5  6  7  0  0  0  0  0  0  0  |
+    *    |  8  9 10 11 12 13 14 15  0  0  |
+    *    | 16 17  0  0  0  0  0  0  0  0  |
+    *    | 18  0  0  0  0  0  0  0  0  0  |
+    *    | 19  0  0  0  0  0  0  0  0  0  |
+    *    | 20  0  0  0  0  0  0  0  0  0  |
+    *    | 21  0  0  0  0  0  0  0  0  0  |
+    *    | 22 23 24 25 26 27 28 29 30 31  |
+    *    \ 32 33 34 35 36 37 38 39 40 41 /
+    */
+
+   const IndexType rows = 10;
+   const IndexType cols = 10;
+
+   Matrix m;
+   m.reset();
+
+   m.setDimensions( rows, cols );
+
+   typename Matrix::CompressedRowLengthsVector rowLengths;
+   rowLengths.setSize( rows );
+   rowLengths.setElement( 0, 4 );
+   rowLengths.setElement( 1, 3 );
+   rowLengths.setElement( 2, 8 );
+   rowLengths.setElement( 3, 2 );
+   for( IndexType i = 4; i < rows - 2; i++ )
+   {
+      rowLengths.setElement( i, 1 );
+   }
+   rowLengths.setElement( 8, 10 );
+   rowLengths.setElement( 9, 10 );
+   m.setCompressedRowLengths( rowLengths );
+
+   RealType value = 1;
+   for( IndexType i = 0; i < 4; i++ )
+      m.setElement( 0, 2 * i, value++ );
+
+   for( IndexType i = 0; i < 3; i++ )
+      m.setElement( 1, i, value++ );
+
+   for( IndexType i = 0; i < 8; i++ )
+      m.setElement( 2, i, value++ );
+
+   for( IndexType i = 0; i < 2; i++ )
+      m.setElement( 3, i, value++ );
+
+   for( IndexType i = 4; i < 8; i++ )
+      m.setElement( i, 0, value++ );
+
+   for( IndexType j = 8; j < rows; j++)
+   {
+      for( IndexType i = 0; i < cols; i++ )
+         m.setElement( j, i, value++ );
+   }
+
+   EXPECT_EQ( m.getNumberOfNonzeroMatrixElements(), 41 );
+}
+
+template< typename Matrix >
+void test_Reset()
+{
+    using RealType = typename Matrix::RealType;
+    using DeviceType = typename Matrix::DeviceType;
+    using IndexType = typename Matrix::IndexType;
+
+/*
+ * Sets up the following 5x4 sparse matrix:
+ *
+ *    /  0  0  0  0 \
+ *    |  0  0  0  0 |
+ *    |  0  0  0  0 |
+ *    |  0  0  0  0 |
+ *    \  0  0  0  0 /
+ */
+
+    const IndexType rows = 5;
+    const IndexType cols = 4;
+
+    Matrix m;
+    m.setDimensions( rows, cols );
+
+    m.reset();
+
+
+    EXPECT_EQ( m.getRows(), 0 );
+    EXPECT_EQ( m.getColumns(), 0 );
+}
+
+template< typename Matrix >
+void test_GetRow()
+{
+    using RealType = typename Matrix::RealType;
+    using DeviceType = typename Matrix::DeviceType;
+    using IndexType = typename Matrix::IndexType;
+
+/*
+ * Sets up the following 10x10 sparse matrix:
+ *
+ *    /  1  0  2  0  3  0  4  0  0  0  \
+ *    |  5  6  7  0  0  0  0  0  0  0  |
+ *    |  8  9 10 11 12 13 14 15  0  0  |
+ *    | 16 17  0  0  0  0  0  0  0  0  |
+ *    | 18  0  0  0  0  0  0  0  0  0  |
+ *    | 19  0  0  0  0  0  0  0  0  0  |
+ *    | 20  0  0  0  0  0  0  0  0  0  |
+ *    | 21  0  0  0  0  0  0  0  0  0  |
+ *    | 22 23 24 25 26 27 28 29 30 31  |
+ *    \ 32 33 34 35 36 37 38 39 40 41 /
+ */
+
+    const IndexType rows = 10;
+    const IndexType cols = 10;
+
+    Matrix m( rows, cols );
+
+    typename Matrix::CompressedRowLengthsVector rowLengths;
+    rowLengths.setSize( rows );
+    rowLengths.setElement( 0, 4 );
+    rowLengths.setElement( 1, 3 );
+    rowLengths.setElement( 2, 8 );
+    rowLengths.setElement( 3, 2 );
+    for( IndexType i = 4; i < rows - 2; i++ )
+    {
+        rowLengths.setElement( i, 1 );
+    }
+    rowLengths.setElement( 8, 10 );
+    rowLengths.setElement( 9, 10 );
+    m.setCompressedRowLengths( rowLengths );
+
+    /*RealType value = 1;
+    for( IndexType i = 0; i < 4; i++ )
+        m.setElement( 0, 2 * i, value++ );
+
+    for( IndexType i = 0; i < 3; i++ )
+        m.setElement( 1, i, value++ );
+
+    for( IndexType i = 0; i < 8; i++ )
+        m.setElement( 2, i, value++ );
+
+    for( IndexType i = 0; i < 2; i++ )
+        m.setElement( 3, i, value++ );
+
+    for( IndexType i = 4; i < 8; i++ )
+        m.setElement( i, 0, value++ );
+
+    for( IndexType j = 8; j < rows; j++)
+    {
+        for( IndexType i = 0; i < cols; i++ )
+            m.setElement( j, i, value++ );
+    }*/
+    auto matrixView = m.getView();
+    auto f = [=] __cuda_callable__ ( const IndexType rowIdx ) mutable {
+       auto row = matrixView.getRow( rowIdx );
+       RealType val;
+       switch( rowIdx )
+       {
+          case 0:
+            val = 1;
+            for( IndexType i = 0; i < 4; i++ )
+               row.setElement( i, 2 * i, val++ );
+            break;
+         case 1:
+            val = 5;
+            for( IndexType i = 0; i < 3; i++ )
+               row.setElement( i, i, val++ );
+            break;
+         case 2:
+            val = 8;
+            for( IndexType i = 0; i < 8; i++ )
+               row.setElement( i, i, val++ );
+            break;
+         case 3:
+            val = 16;
+            for( IndexType i = 0; i < 2; i++ )
+               row.setElement( i, i, val++ );
+            break;
+         case 4:
+            row.setElement( 0, 0, 18 );
+            break;
+         case 5:
+            row.setElement( 0, 0, 19 );
+            break;
+         case 6:
+            row.setElement( 0, 0, 20 );
+            break;
+         case 7:
+            row.setElement( 0, 0, 21 );
+            break;
+         case 8:
+             val = 22;
+             for( IndexType i = 0; i < rows; i++ )
+                row.setElement( i, i, val++ );
+             break;
+         case 9:
+             val = 32;
+             for( IndexType i = 0; i < rows; i++ )
+                row.setElement( i, i, val++ );
+             break;
+       }
+    };
+    TNL::Algorithms::ParallelFor< DeviceType >::exec( ( IndexType ) 0, rows, f );
+
+    EXPECT_EQ( m.getElement( 0, 0 ),  1 );
+    EXPECT_EQ( m.getElement( 0, 1 ),  0 );
+    EXPECT_EQ( m.getElement( 0, 2 ),  2 );
+    EXPECT_EQ( m.getElement( 0, 3 ),  0 );
+    EXPECT_EQ( m.getElement( 0, 4 ),  3 );
+    EXPECT_EQ( m.getElement( 0, 5 ),  0 );
+    EXPECT_EQ( m.getElement( 0, 6 ),  4 );
+    EXPECT_EQ( m.getElement( 0, 7 ),  0 );
+    EXPECT_EQ( m.getElement( 0, 8 ),  0 );
+    EXPECT_EQ( m.getElement( 0, 9 ),  0 );
+
+    EXPECT_EQ( m.getElement( 1, 0 ),  5 );
+    EXPECT_EQ( m.getElement( 1, 1 ),  6 );
+    EXPECT_EQ( m.getElement( 1, 2 ),  7 );
+    EXPECT_EQ( m.getElement( 1, 3 ),  0 );
+    EXPECT_EQ( m.getElement( 1, 4 ),  0 );
+    EXPECT_EQ( m.getElement( 1, 5 ),  0 );
+    EXPECT_EQ( m.getElement( 1, 6 ),  0 );
+    EXPECT_EQ( m.getElement( 1, 7 ),  0 );
+    EXPECT_EQ( m.getElement( 1, 8 ),  0 );
+    EXPECT_EQ( m.getElement( 1, 9 ),  0 );
+
+    EXPECT_EQ( m.getElement( 2, 0 ),  8 );
+    EXPECT_EQ( m.getElement( 2, 1 ),  9 );
+    EXPECT_EQ( m.getElement( 2, 2 ), 10 );
+    EXPECT_EQ( m.getElement( 2, 3 ), 11 );
+    EXPECT_EQ( m.getElement( 2, 4 ), 12 );
+    EXPECT_EQ( m.getElement( 2, 5 ), 13 );
+    EXPECT_EQ( m.getElement( 2, 6 ), 14 );
+    EXPECT_EQ( m.getElement( 2, 7 ), 15 );
+    EXPECT_EQ( m.getElement( 2, 8 ),  0 );
+    EXPECT_EQ( m.getElement( 2, 9 ),  0 );
+
+    EXPECT_EQ( m.getElement( 3, 0 ), 16 );
+    EXPECT_EQ( m.getElement( 3, 1 ), 17 );
+    EXPECT_EQ( m.getElement( 3, 2 ),  0 );
+    EXPECT_EQ( m.getElement( 3, 3 ),  0 );
+    EXPECT_EQ( m.getElement( 3, 4 ),  0 );
+    EXPECT_EQ( m.getElement( 3, 5 ),  0 );
+    EXPECT_EQ( m.getElement( 3, 6 ),  0 );
+    EXPECT_EQ( m.getElement( 3, 7 ),  0 );
+    EXPECT_EQ( m.getElement( 3, 8 ),  0 );
+    EXPECT_EQ( m.getElement( 3, 9 ),  0 );
+
+    EXPECT_EQ( m.getElement( 4, 0 ), 18 );
+    EXPECT_EQ( m.getElement( 4, 1 ),  0 );
+    EXPECT_EQ( m.getElement( 4, 2 ),  0 );
+    EXPECT_EQ( m.getElement( 4, 3 ),  0 );
+    EXPECT_EQ( m.getElement( 4, 4 ),  0 );
+    EXPECT_EQ( m.getElement( 4, 5 ),  0 );
+    EXPECT_EQ( m.getElement( 4, 6 ),  0 );
+    EXPECT_EQ( m.getElement( 4, 7 ),  0 );
+    EXPECT_EQ( m.getElement( 4, 8 ),  0 );
+    EXPECT_EQ( m.getElement( 4, 9 ),  0 );
+
+    EXPECT_EQ( m.getElement( 5, 0 ), 19 );
+    EXPECT_EQ( m.getElement( 5, 1 ),  0 );
+    EXPECT_EQ( m.getElement( 5, 2 ),  0 );
+    EXPECT_EQ( m.getElement( 5, 3 ),  0 );
+    EXPECT_EQ( m.getElement( 5, 4 ),  0 );
+    EXPECT_EQ( m.getElement( 5, 5 ),  0 );
+    EXPECT_EQ( m.getElement( 5, 6 ),  0 );
+    EXPECT_EQ( m.getElement( 5, 7 ),  0 );
+    EXPECT_EQ( m.getElement( 5, 8 ),  0 );
+    EXPECT_EQ( m.getElement( 5, 9 ),  0 );
+
+    EXPECT_EQ( m.getElement( 6, 0 ), 20 );
+    EXPECT_EQ( m.getElement( 6, 1 ),  0 );
+    EXPECT_EQ( m.getElement( 6, 2 ),  0 );
+    EXPECT_EQ( m.getElement( 6, 3 ),  0 );
+    EXPECT_EQ( m.getElement( 6, 4 ),  0 );
+    EXPECT_EQ( m.getElement( 6, 5 ),  0 );
+    EXPECT_EQ( m.getElement( 6, 6 ),  0 );
+    EXPECT_EQ( m.getElement( 6, 7 ),  0 );
+    EXPECT_EQ( m.getElement( 6, 8 ),  0 );
+    EXPECT_EQ( m.getElement( 6, 9 ),  0 );
+
+    EXPECT_EQ( m.getElement( 7, 0 ), 21 );
+    EXPECT_EQ( m.getElement( 7, 1 ),  0 );
+    EXPECT_EQ( m.getElement( 7, 2 ),  0 );
+    EXPECT_EQ( m.getElement( 7, 3 ),  0 );
+    EXPECT_EQ( m.getElement( 7, 4 ),  0 );
+    EXPECT_EQ( m.getElement( 7, 5 ),  0 );
+    EXPECT_EQ( m.getElement( 7, 6 ),  0 );
+    EXPECT_EQ( m.getElement( 7, 7 ),  0 );
+    EXPECT_EQ( m.getElement( 7, 8 ),  0 );
+    EXPECT_EQ( m.getElement( 7, 9 ),  0 );
+
+    EXPECT_EQ( m.getElement( 8, 0 ), 22 );
+    EXPECT_EQ( m.getElement( 8, 1 ), 23 );
+    EXPECT_EQ( m.getElement( 8, 2 ), 24 );
+    EXPECT_EQ( m.getElement( 8, 3 ), 25 );
+    EXPECT_EQ( m.getElement( 8, 4 ), 26 );
+    EXPECT_EQ( m.getElement( 8, 5 ), 27 );
+    EXPECT_EQ( m.getElement( 8, 6 ), 28 );
+    EXPECT_EQ( m.getElement( 8, 7 ), 29 );
+    EXPECT_EQ( m.getElement( 8, 8 ), 30 );
+    EXPECT_EQ( m.getElement( 8, 9 ), 31 );
+
+    EXPECT_EQ( m.getElement( 9, 0 ), 32 );
+    EXPECT_EQ( m.getElement( 9, 1 ), 33 );
+    EXPECT_EQ( m.getElement( 9, 2 ), 34 );
+    EXPECT_EQ( m.getElement( 9, 3 ), 35 );
+    EXPECT_EQ( m.getElement( 9, 4 ), 36 );
+    EXPECT_EQ( m.getElement( 9, 5 ), 37 );
+    EXPECT_EQ( m.getElement( 9, 6 ), 38 );
+    EXPECT_EQ( m.getElement( 9, 7 ), 39 );
+    EXPECT_EQ( m.getElement( 9, 8 ), 40 );
+    EXPECT_EQ( m.getElement( 9, 9 ), 41 );
+}
+
+
+template< typename Matrix >
+void test_SetElement()
+{
+    using RealType = typename Matrix::RealType;
+    using DeviceType = typename Matrix::DeviceType;
+    using IndexType = typename Matrix::IndexType;
+
+/*
+ * Sets up the following 10x10 sparse matrix:
+ *
+ *    /  1  0  2  0  3  0  4  0  0  0  \
+ *    |  5  6  7  0  0  0  0  0  0  0  |
+ *    |  8  9 10 11 12 13 14 15  0  0  |
+ *    | 16 17  0  0  0  0  0  0  0  0  |
+ *    | 18  0  0  0  0  0  0  0  0  0  |
+ *    | 19  0  0  0  0  0  0  0  0  0  |
+ *    | 20  0  0  0  0  0  0  0  0  0  |
+ *    | 21  0  0  0  0  0  0  0  0  0  |
+ *    | 22 23 24 25 26 27 28 29 30 31  |
+ *    \ 32 33 34 35 36 37 38 39 40 41 /
+ */
+
+    const IndexType rows = 10;
+    const IndexType cols = 10;
+
+    Matrix m;
+    m.reset();
+
+    m.setDimensions( rows, cols );
+
+    typename Matrix::CompressedRowLengthsVector rowLengths;
+    rowLengths.setSize( rows );
+    rowLengths.setElement( 0, 4 );
+    rowLengths.setElement( 1, 3 );
+    rowLengths.setElement( 2, 8 );
+    rowLengths.setElement( 3, 2 );
+    for( IndexType i = 4; i < rows - 2; i++ )
+    {
+        rowLengths.setElement( i, 1 );
+    }
+    rowLengths.setElement( 8, 10 );
+    rowLengths.setElement( 9, 10 );
+    m.setCompressedRowLengths( rowLengths );
+
+    RealType value = 1;
+    for( IndexType i = 0; i < 4; i++ )
+        m.setElement( 0, 2 * i, value++ );
+
+    for( IndexType i = 0; i < 3; i++ )
+        m.setElement( 1, i, value++ );
+
+    for( IndexType i = 0; i < 8; i++ )
+        m.setElement( 2, i, value++ );
+
+    for( IndexType i = 0; i < 2; i++ )
+        m.setElement( 3, i, value++ );
+
+    for( IndexType i = 4; i < 8; i++ )
+        m.setElement( i, 0, value++ );
+
+    for( IndexType j = 8; j < rows; j++)
+    {
+        for( IndexType i = 0; i < cols; i++ )
+            m.setElement( j, i, value++ );
+    }
+
+    EXPECT_EQ( m.getElement( 0, 0 ),  1 );
+    EXPECT_EQ( m.getElement( 0, 1 ),  0 );
+    EXPECT_EQ( m.getElement( 0, 2 ),  2 );
+    EXPECT_EQ( m.getElement( 0, 3 ),  0 );
+    EXPECT_EQ( m.getElement( 0, 4 ),  3 );
+    EXPECT_EQ( m.getElement( 0, 5 ),  0 );
+    EXPECT_EQ( m.getElement( 0, 6 ),  4 );
+    EXPECT_EQ( m.getElement( 0, 7 ),  0 );
+    EXPECT_EQ( m.getElement( 0, 8 ),  0 );
+    EXPECT_EQ( m.getElement( 0, 9 ),  0 );
+
+    EXPECT_EQ( m.getElement( 1, 0 ),  5 );
+    EXPECT_EQ( m.getElement( 1, 1 ),  6 );
+    EXPECT_EQ( m.getElement( 1, 2 ),  7 );
+    EXPECT_EQ( m.getElement( 1, 3 ),  0 );
+    EXPECT_EQ( m.getElement( 1, 4 ),  0 );
+    EXPECT_EQ( m.getElement( 1, 5 ),  0 );
+    EXPECT_EQ( m.getElement( 1, 6 ),  0 );
+    EXPECT_EQ( m.getElement( 1, 7 ),  0 );
+    EXPECT_EQ( m.getElement( 1, 8 ),  0 );
+    EXPECT_EQ( m.getElement( 1, 9 ),  0 );
+
+    EXPECT_EQ( m.getElement( 2, 0 ),  8 );
+    EXPECT_EQ( m.getElement( 2, 1 ),  9 );
+    EXPECT_EQ( m.getElement( 2, 2 ), 10 );
+    EXPECT_EQ( m.getElement( 2, 3 ), 11 );
+    EXPECT_EQ( m.getElement( 2, 4 ), 12 );
+    EXPECT_EQ( m.getElement( 2, 5 ), 13 );
+    EXPECT_EQ( m.getElement( 2, 6 ), 14 );
+    EXPECT_EQ( m.getElement( 2, 7 ), 15 );
+    EXPECT_EQ( m.getElement( 2, 8 ),  0 );
+    EXPECT_EQ( m.getElement( 2, 9 ),  0 );
+
+    EXPECT_EQ( m.getElement( 3, 0 ), 16 );
+    EXPECT_EQ( m.getElement( 3, 1 ), 17 );
+    EXPECT_EQ( m.getElement( 3, 2 ),  0 );
+    EXPECT_EQ( m.getElement( 3, 3 ),  0 );
+    EXPECT_EQ( m.getElement( 3, 4 ),  0 );
+    EXPECT_EQ( m.getElement( 3, 5 ),  0 );
+    EXPECT_EQ( m.getElement( 3, 6 ),  0 );
+    EXPECT_EQ( m.getElement( 3, 7 ),  0 );
+    EXPECT_EQ( m.getElement( 3, 8 ),  0 );
+    EXPECT_EQ( m.getElement( 3, 9 ),  0 );
+
+    EXPECT_EQ( m.getElement( 4, 0 ), 18 );
+    EXPECT_EQ( m.getElement( 4, 1 ),  0 );
+    EXPECT_EQ( m.getElement( 4, 2 ),  0 );
+    EXPECT_EQ( m.getElement( 4, 3 ),  0 );
+    EXPECT_EQ( m.getElement( 4, 4 ),  0 );
+    EXPECT_EQ( m.getElement( 4, 5 ),  0 );
+    EXPECT_EQ( m.getElement( 4, 6 ),  0 );
+    EXPECT_EQ( m.getElement( 4, 7 ),  0 );
+    EXPECT_EQ( m.getElement( 4, 8 ),  0 );
+    EXPECT_EQ( m.getElement( 4, 9 ),  0 );
+
+    EXPECT_EQ( m.getElement( 5, 0 ), 19 );
+    EXPECT_EQ( m.getElement( 5, 1 ),  0 );
+    EXPECT_EQ( m.getElement( 5, 2 ),  0 );
+    EXPECT_EQ( m.getElement( 5, 3 ),  0 );
+    EXPECT_EQ( m.getElement( 5, 4 ),  0 );
+    EXPECT_EQ( m.getElement( 5, 5 ),  0 );
+    EXPECT_EQ( m.getElement( 5, 6 ),  0 );
+    EXPECT_EQ( m.getElement( 5, 7 ),  0 );
+    EXPECT_EQ( m.getElement( 5, 8 ),  0 );
+    EXPECT_EQ( m.getElement( 5, 9 ),  0 );
+
+    EXPECT_EQ( m.getElement( 6, 0 ), 20 );
+    EXPECT_EQ( m.getElement( 6, 1 ),  0 );
+    EXPECT_EQ( m.getElement( 6, 2 ),  0 );
+    EXPECT_EQ( m.getElement( 6, 3 ),  0 );
+    EXPECT_EQ( m.getElement( 6, 4 ),  0 );
+    EXPECT_EQ( m.getElement( 6, 5 ),  0 );
+    EXPECT_EQ( m.getElement( 6, 6 ),  0 );
+    EXPECT_EQ( m.getElement( 6, 7 ),  0 );
+    EXPECT_EQ( m.getElement( 6, 8 ),  0 );
+    EXPECT_EQ( m.getElement( 6, 9 ),  0 );
+
+    EXPECT_EQ( m.getElement( 7, 0 ), 21 );
+    EXPECT_EQ( m.getElement( 7, 1 ),  0 );
+    EXPECT_EQ( m.getElement( 7, 2 ),  0 );
+    EXPECT_EQ( m.getElement( 7, 3 ),  0 );
+    EXPECT_EQ( m.getElement( 7, 4 ),  0 );
+    EXPECT_EQ( m.getElement( 7, 5 ),  0 );
+    EXPECT_EQ( m.getElement( 7, 6 ),  0 );
+    EXPECT_EQ( m.getElement( 7, 7 ),  0 );
+    EXPECT_EQ( m.getElement( 7, 8 ),  0 );
+    EXPECT_EQ( m.getElement( 7, 9 ),  0 );
+
+    EXPECT_EQ( m.getElement( 8, 0 ), 22 );
+    EXPECT_EQ( m.getElement( 8, 1 ), 23 );
+    EXPECT_EQ( m.getElement( 8, 2 ), 24 );
+    EXPECT_EQ( m.getElement( 8, 3 ), 25 );
+    EXPECT_EQ( m.getElement( 8, 4 ), 26 );
+    EXPECT_EQ( m.getElement( 8, 5 ), 27 );
+    EXPECT_EQ( m.getElement( 8, 6 ), 28 );
+    EXPECT_EQ( m.getElement( 8, 7 ), 29 );
+    EXPECT_EQ( m.getElement( 8, 8 ), 30 );
+    EXPECT_EQ( m.getElement( 8, 9 ), 31 );
+
+    EXPECT_EQ( m.getElement( 9, 0 ), 32 );
+    EXPECT_EQ( m.getElement( 9, 1 ), 33 );
+    EXPECT_EQ( m.getElement( 9, 2 ), 34 );
+    EXPECT_EQ( m.getElement( 9, 3 ), 35 );
+    EXPECT_EQ( m.getElement( 9, 4 ), 36 );
+    EXPECT_EQ( m.getElement( 9, 5 ), 37 );
+    EXPECT_EQ( m.getElement( 9, 6 ), 38 );
+    EXPECT_EQ( m.getElement( 9, 7 ), 39 );
+    EXPECT_EQ( m.getElement( 9, 8 ), 40 );
+    EXPECT_EQ( m.getElement( 9, 9 ), 41 );
+}
+
+template< typename Matrix >
+void test_AddElement()
+{
+    using RealType = typename Matrix::RealType;
+    using DeviceType = typename Matrix::DeviceType;
+    using IndexType = typename Matrix::IndexType;
+
+/*
+ * Sets up the following 6x5 sparse matrix:
+ *
+ *    /  1  2  3  0  0 \
+ *    |  0  4  5  6  0 |
+ *    |  0  0  7  8  9 |
+ *    | 10  0  0  0  0 |
+ *    |  0 11  0  0  0 |
+ *    \  0  0  0 12  0 /
+ */
+
+    const IndexType rows = 6;
+    const IndexType cols = 5;
+
+    Matrix m;
+    m.reset();
+    m.setDimensions( rows, cols );
+    typename Matrix::CompressedRowLengthsVector rowLengths;
+    rowLengths.setSize( rows );
+    rowLengths.setValue( 3 );
+    m.setCompressedRowLengths( rowLengths );
+
+    RealType value = 1;
+    for( IndexType i = 0; i < cols - 2; i++ )     // 0th row
+        m.setElement( 0, i, value++ );
+
+    for( IndexType i = 1; i < cols - 1; i++ )     // 1st row
+        m.setElement( 1, i, value++ );
+
+    for( IndexType i = 2; i < cols; i++ )         // 2nd row
+        m.setElement( 2, i, value++ );
+
+    m.setElement( 3, 0, value++ );      // 3rd row
+
+    m.setElement( 4, 1, value++ );      // 4th row
+
+    m.setElement( 5, 3, value++ );      // 5th row
+
+
+    // Check the set elements
+    EXPECT_EQ( m.getElement( 0, 0 ),  1 );
+    EXPECT_EQ( m.getElement( 0, 1 ),  2 );
+    EXPECT_EQ( m.getElement( 0, 2 ),  3 );
+    EXPECT_EQ( m.getElement( 0, 3 ),  0 );
+    EXPECT_EQ( m.getElement( 0, 4 ),  0 );
+
+    EXPECT_EQ( m.getElement( 1, 0 ),  0 );
+    EXPECT_EQ( m.getElement( 1, 1 ),  4 );
+    EXPECT_EQ( m.getElement( 1, 2 ),  5 );
+    EXPECT_EQ( m.getElement( 1, 3 ),  6 );
+    EXPECT_EQ( m.getElement( 1, 4 ),  0 );
+
+    EXPECT_EQ( m.getElement( 2, 0 ),  0 );
+    EXPECT_EQ( m.getElement( 2, 1 ),  0 );
+    EXPECT_EQ( m.getElement( 2, 2 ),  7 );
+    EXPECT_EQ( m.getElement( 2, 3 ),  8 );
+    EXPECT_EQ( m.getElement( 2, 4 ),  9 );
+
+    EXPECT_EQ( m.getElement( 3, 0 ), 10 );
+    EXPECT_EQ( m.getElement( 3, 1 ),  0 );
+    EXPECT_EQ( m.getElement( 3, 2 ),  0 );
+    EXPECT_EQ( m.getElement( 3, 3 ),  0 );
+    EXPECT_EQ( m.getElement( 3, 4 ),  0 );
+
+    EXPECT_EQ( m.getElement( 4, 0 ),  0 );
+    EXPECT_EQ( m.getElement( 4, 1 ), 11 );
+    EXPECT_EQ( m.getElement( 4, 2 ),  0 );
+    EXPECT_EQ( m.getElement( 4, 3 ),  0 );
+    EXPECT_EQ( m.getElement( 4, 4 ),  0 );
+
+    EXPECT_EQ( m.getElement( 5, 0 ),  0 );
+    EXPECT_EQ( m.getElement( 5, 1 ),  0 );
+    EXPECT_EQ( m.getElement( 5, 2 ),  0 );
+    EXPECT_EQ( m.getElement( 5, 3 ), 12 );
+    EXPECT_EQ( m.getElement( 5, 4 ),  0 );
+
+    // Add new elements to the old elements with a multiplying factor applied to the old elements.
+
+/*
+ * Sets up the following 6x5 sparse matrix:
+ *
+ *    /  1  2  3  0  0 \
+ *    |  0  4  5  6  0 |
+ *    |  0  0  7  8  9 |
+ *    | 10  0  0  0  0 |
+ *    |  0 11  0  0  0 |
+ *    \  0  0  0 12  0 /
+ */
+
+/*
+ * The following setup results in the following 6x5 sparse matrix:
+ *
+ *    /  3  6  9  0  0 \
+ *    |  0 12 15 18  0 |
+ *    |  0  0 21 24 27 |
+ *    | 30 11 12  0  0 |
+ *    |  0 35 14 15  0 |
+ *    \  0  0 16 41 18 /
+ */
+
+    RealType newValue = 1;
+    for( IndexType i = 0; i < cols - 2; i++ )         // 0th row
+        m.addElement( 0, i, newValue++, 2.0 );
+
+    for( IndexType i = 1; i < cols - 1; i++ )         // 1st row
+        m.addElement( 1, i, newValue++, 2.0 );
+
+    for( IndexType i = 2; i < cols; i++ )             // 2nd row
+        m.addElement( 2, i, newValue++, 2.0 );
+
+    for( IndexType i = 0; i < cols - 2; i++ )         // 3rd row
+        m.addElement( 3, i, newValue++, 2.0 );
+
+    for( IndexType i = 1; i < cols - 1; i++ )         // 4th row
+        m.addElement( 4, i, newValue++, 2.0 );
+
+    for( IndexType i = 2; i < cols; i++ )             // 5th row
+        m.addElement( 5, i, newValue++, 2.0 );
+
+
+    EXPECT_EQ( m.getElement( 0, 0 ),  3 );
+    EXPECT_EQ( m.getElement( 0, 1 ),  6 );
+    EXPECT_EQ( m.getElement( 0, 2 ),  9 );
+    EXPECT_EQ( m.getElement( 0, 3 ),  0 );
+    EXPECT_EQ( m.getElement( 0, 4 ),  0 );
+
+    EXPECT_EQ( m.getElement( 1, 0 ),  0 );
+    EXPECT_EQ( m.getElement( 1, 1 ), 12 );
+    EXPECT_EQ( m.getElement( 1, 2 ), 15 );
+    EXPECT_EQ( m.getElement( 1, 3 ), 18 );
+    EXPECT_EQ( m.getElement( 1, 4 ),  0 );
+
+    EXPECT_EQ( m.getElement( 2, 0 ),  0 );
+    EXPECT_EQ( m.getElement( 2, 1 ),  0 );
+    EXPECT_EQ( m.getElement( 2, 2 ), 21 );
+    EXPECT_EQ( m.getElement( 2, 3 ), 24 );
+    EXPECT_EQ( m.getElement( 2, 4 ), 27 );
+
+    EXPECT_EQ( m.getElement( 3, 0 ), 30 );
+    EXPECT_EQ( m.getElement( 3, 1 ), 11 );
+    EXPECT_EQ( m.getElement( 3, 2 ), 12 );
+    EXPECT_EQ( m.getElement( 3, 3 ),  0 );
+    EXPECT_EQ( m.getElement( 3, 4 ),  0 );
+
+    EXPECT_EQ( m.getElement( 4, 0 ),  0 );
+    EXPECT_EQ( m.getElement( 4, 1 ), 35 );
+    EXPECT_EQ( m.getElement( 4, 2 ), 14 );
+    EXPECT_EQ( m.getElement( 4, 3 ), 15 );
+    EXPECT_EQ( m.getElement( 4, 4 ),  0 );
+
+    EXPECT_EQ( m.getElement( 5, 0 ),  0 );
+    EXPECT_EQ( m.getElement( 5, 1 ),  0 );
+    EXPECT_EQ( m.getElement( 5, 2 ), 16 );
+    EXPECT_EQ( m.getElement( 5, 3 ), 41 );
+    EXPECT_EQ( m.getElement( 5, 4 ), 18 );
+}
+
+template< typename Matrix >
+void test_SetRow()
+{
+    using RealType = typename Matrix::RealType;
+    using DeviceType = typename Matrix::DeviceType;
+    using IndexType = typename Matrix::IndexType;
+
+/*
+ * Sets up the following 3x7 sparse matrix:
+ *
+ *    /  0  0  0  1  1  1  0 \
+ *    |  2  2  2  0  0  0  0 |
+ *    \  3  3  3  0  0  0  0 /
+ */
+
+    const IndexType rows = 3;
+    const IndexType cols = 7;
+
+    Matrix m;
+    m.reset();
+    m.setDimensions( rows, cols );
+    typename Matrix::CompressedRowLengthsVector rowLengths;
+    rowLengths.setSize( rows );
+    rowLengths.setValue( 6 );
+    rowLengths.setElement( 1, 3 );
+    m.setCompressedRowLengths( rowLengths );
+
+    RealType value = 1;
+    for( IndexType i = 0; i < 3; i++ )
+    {
+        m.setElement( 0, i + 3, value );
+        m.setElement( 1, i, value + 1 );
+        m.setElement( 2, i, value + 2 );
+    }
+
+    RealType row1 [ 3 ] = { 11, 11, 11 }; IndexType colIndexes1 [ 3 ] = { 0, 1, 2 };
+    RealType row2 [ 3 ] = { 22, 22, 22 }; IndexType colIndexes2 [ 3 ] = { 0, 1, 2 };
+    RealType row3 [ 3 ] = { 33, 33, 33 }; IndexType colIndexes3 [ 3 ] = { 3, 4, 5 };
+
+    RealType row = 0;
+    IndexType elements = 3;
+
+    m.setRow( row++, colIndexes1, row1, elements );
+    m.setRow( row++, colIndexes2, row2, elements );
+    m.setRow( row++, colIndexes3, row3, elements );
+
+
+    EXPECT_EQ( m.getElement( 0, 0 ), 11 );
+    EXPECT_EQ( m.getElement( 0, 1 ), 11 );
+    EXPECT_EQ( m.getElement( 0, 2 ), 11 );
+    EXPECT_EQ( m.getElement( 0, 3 ),  0 );
+    EXPECT_EQ( m.getElement( 0, 4 ),  0 );
+    EXPECT_EQ( m.getElement( 0, 5 ),  0 );
+    EXPECT_EQ( m.getElement( 0, 6 ),  0 );
+
+    EXPECT_EQ( m.getElement( 1, 0 ), 22 );
+    EXPECT_EQ( m.getElement( 1, 1 ), 22 );
+    EXPECT_EQ( m.getElement( 1, 2 ), 22 );
+    EXPECT_EQ( m.getElement( 1, 3 ),  0 );
+    EXPECT_EQ( m.getElement( 1, 4 ),  0 );
+    EXPECT_EQ( m.getElement( 1, 5 ),  0 );
+    EXPECT_EQ( m.getElement( 1, 6 ),  0 );
+
+    EXPECT_EQ( m.getElement( 2, 0 ),  0 );
+    EXPECT_EQ( m.getElement( 2, 1 ),  0 );
+    EXPECT_EQ( m.getElement( 2, 2 ),  0 );
+    EXPECT_EQ( m.getElement( 2, 3 ), 33 );
+    EXPECT_EQ( m.getElement( 2, 4 ), 33 );
+    EXPECT_EQ( m.getElement( 2, 5 ), 33 );
+    EXPECT_EQ( m.getElement( 2, 6 ),  0 );
+}
+
+template< typename Matrix >
+void test_VectorProduct()
+{
+    using RealType = typename Matrix::RealType;
+    using DeviceType = typename Matrix::DeviceType;
+    using IndexType = typename Matrix::IndexType;
+    using VectorType = TNL::Containers::Vector< RealType, DeviceType, IndexType >;
+
+/*
+ * Sets up the following 4x4 sparse matrix:
+ *
+ *    /  1  0  0  0 \
+ *    |  0  2  0  3 |
+ *    |  0  4  0  0 |
+ *    \  0  0  5  0 /
+ */
+
+    const IndexType m_rows_1 = 4;
+    const IndexType m_cols_1 = 4;
+
+    Matrix m_1;
+    m_1.reset();
+    m_1.setDimensions( m_rows_1, m_cols_1 );
+    typename Matrix::CompressedRowLengthsVector rowLengths_1;
+    rowLengths_1.setSize( m_rows_1 );
+    rowLengths_1.setElement( 0, 1 );
+    rowLengths_1.setElement( 1, 2 );
+    rowLengths_1.setElement( 2, 1 );
+    rowLengths_1.setElement( 3, 1 );
+    m_1.setCompressedRowLengths( rowLengths_1 );
+
+    RealType value_1 = 1;
+    m_1.setElement( 0, 0, value_1++ );      // 0th row
+
+    m_1.setElement( 1, 1, value_1++ );      // 1st row
+    m_1.setElement( 1, 3, value_1++ );
+
+    m_1.setElement( 2, 1, value_1++ );      // 2nd row
+
+    m_1.setElement( 3, 2, value_1++ );      // 3rd row
+
+    VectorType inVector_1;
+    inVector_1.setSize( m_cols_1 );
+    for( IndexType i = 0; i < inVector_1.getSize(); i++ )
+        inVector_1.setElement( i, 2 );
+
+    VectorType outVector_1;
+    outVector_1.setSize( m_rows_1 );
+    for( IndexType j = 0; j < outVector_1.getSize(); j++ )
+        outVector_1.setElement( j, 0 );
+
+
+    m_1.vectorProduct( inVector_1, outVector_1 );
+
+
+    EXPECT_EQ( outVector_1.getElement( 0 ),  2 );
+    EXPECT_EQ( outVector_1.getElement( 1 ), 10 );
+    EXPECT_EQ( outVector_1.getElement( 2 ),  8 );
+    EXPECT_EQ( outVector_1.getElement( 3 ), 10 );
+
+
+/*
+ * Sets up the following 4x4 sparse matrix:
+ *
+ *    /  1  2  3  0 \
+ *    |  0  0  0  4 |
+ *    |  5  6  7  0 |
+ *    \  0  8  0  0 /
+ */
+
+    const IndexType m_rows_2 = 4;
+    const IndexType m_cols_2 = 4;
+
+    Matrix m_2;
+    m_2.reset();
+    m_2.setDimensions( m_rows_2, m_cols_2 );
+    typename Matrix::CompressedRowLengthsVector rowLengths_2;
+    rowLengths_2.setSize( m_rows_2 );
+    rowLengths_2.setValue( 3 );
+    rowLengths_2.setElement( 1, 1 );
+    rowLengths_2.setElement( 3, 1 );
+    m_2.setCompressedRowLengths( rowLengths_2 );
+
+    RealType value_2 = 1;
+    for( IndexType i = 0; i < 3; i++ )   // 0th row
+        m_2.setElement( 0, i, value_2++ );
+
+    m_2.setElement( 1, 3, value_2++ );      // 1st row
+
+    for( IndexType i = 0; i < 3; i++ )   // 2nd row
+        m_2.setElement( 2, i, value_2++ );
+
+    for( IndexType i = 1; i < 2; i++ )       // 3rd row
+        m_2.setElement( 3, i, value_2++ );
+
+    VectorType inVector_2;
+    inVector_2.setSize( m_cols_2 );
+    for( IndexType i = 0; i < inVector_2.getSize(); i++ )
+        inVector_2.setElement( i, 2 );
+
+    VectorType outVector_2;
+    outVector_2.setSize( m_rows_2 );
+    for( IndexType j = 0; j < outVector_2.getSize(); j++ )
+        outVector_2.setElement( j, 0 );
+
+
+    m_2.vectorProduct( inVector_2, outVector_2 );
+
+
+    EXPECT_EQ( outVector_2.getElement( 0 ), 12 );
+    EXPECT_EQ( outVector_2.getElement( 1 ),  8 );
+    EXPECT_EQ( outVector_2.getElement( 2 ), 36 );
+    EXPECT_EQ( outVector_2.getElement( 3 ), 16 );
+
+
+/*
+ * Sets up the following 4x4 sparse matrix:
+ *
+ *    /  1  2  3  0 \
+ *    |  0  4  5  6 |
+ *    |  7  8  9  0 |
+ *    \  0 10 11 12 /
+ */
+
+    const IndexType m_rows_3 = 4;
+    const IndexType m_cols_3 = 4;
+
+    Matrix m_3;
+    m_3.reset();
+    m_3.setDimensions( m_rows_3, m_cols_3 );
+    typename Matrix::CompressedRowLengthsVector rowLengths_3;
+    rowLengths_3.setSize( m_rows_3 );
+    rowLengths_3.setValue( 3 );
+    m_3.setCompressedRowLengths( rowLengths_3 );
+
+    RealType value_3 = 1;
+    for( IndexType i = 0; i < 3; i++ )          // 0th row
+        m_3.setElement( 0, i, value_3++ );
+
+    for( IndexType i = 1; i < 4; i++ )
+        m_3.setElement( 1, i, value_3++ );      // 1st row
+
+    for( IndexType i = 0; i < 3; i++ )          // 2nd row
+        m_3.setElement( 2, i, value_3++ );
+
+    for( IndexType i = 1; i < 4; i++ )          // 3rd row
+        m_3.setElement( 3, i, value_3++ );
+
+    VectorType inVector_3;
+    inVector_3.setSize( m_cols_3 );
+    for( IndexType i = 0; i < inVector_3.getSize(); i++ )
+        inVector_3.setElement( i, 2 );
+
+    VectorType outVector_3;
+    outVector_3.setSize( m_rows_3 );
+    for( IndexType j = 0; j < outVector_3.getSize(); j++ )
+        outVector_3.setElement( j, 0 );
+
+
+    m_3.vectorProduct( inVector_3, outVector_3 );
+
+
+    EXPECT_EQ( outVector_3.getElement( 0 ), 12 );
+    EXPECT_EQ( outVector_3.getElement( 1 ), 30 );
+    EXPECT_EQ( outVector_3.getElement( 2 ), 48 );
+    EXPECT_EQ( outVector_3.getElement( 3 ), 66 );
+
+
+/*
+ * Sets up the following 8x8 sparse matrix:
+ *
+ *    /  1  2  3  0  0  4  0  0 \
+ *    |  0  5  6  7  8  0  0  0 |
+ *    |  9 10 11 12 13  0  0  0 |
+ *    |  0 14 15 16 17  0  0  0 |
+ *    |  0  0 18 19 20 21  0  0 |
+ *    |  0  0  0 22 23 24 25  0 |
+ *    | 26 27 28 29 30  0  0  0 |
+ *    \ 31 32 33 34 35  0  0  0 /
+ */
+
+    const IndexType m_rows_4 = 8;
+    const IndexType m_cols_4 = 8;
+
+    Matrix m_4;
+    m_4.reset();
+    m_4.setDimensions( m_rows_4, m_cols_4 );
+    typename Matrix::CompressedRowLengthsVector rowLengths_4;
+    rowLengths_4.setSize( m_rows_4 );
+    rowLengths_4.setValue( 4 );
+    rowLengths_4.setElement( 2, 5 );
+    rowLengths_4.setElement( 6, 5 );
+    rowLengths_4.setElement( 7, 5 );
+    m_4.setCompressedRowLengths( rowLengths_4 );
+
+    RealType value_4 = 1;
+    for( IndexType i = 0; i < 3; i++ )       // 0th row
+        m_4.setElement( 0, i, value_4++ );
+
+    m_4.setElement( 0, 5, value_4++ );
+
+    for( IndexType i = 1; i < 5; i++ )       // 1st row
+        m_4.setElement( 1, i, value_4++ );
+
+    for( IndexType i = 0; i < 5; i++ )       // 2nd row
+        m_4.setElement( 2, i, value_4++ );
+
+    for( IndexType i = 1; i < 5; i++ )       // 3rd row
+        m_4.setElement( 3, i, value_4++ );
+
+    for( IndexType i = 2; i < 6; i++ )       // 4th row
+        m_4.setElement( 4, i, value_4++ );
+
+    for( IndexType i = 3; i < 7; i++ )       // 5th row
+        m_4.setElement( 5, i, value_4++ );
+
+    for( IndexType i = 0; i < 5; i++ )       // 6th row
+        m_4.setElement( 6, i, value_4++ );
+
+    for( IndexType i = 0; i < 5; i++ )       // 7th row
+        m_4.setElement( 7, i, value_4++ );
+
+    VectorType inVector_4;
+    inVector_4.setSize( m_cols_4 );
+    for( IndexType i = 0; i < inVector_4.getSize(); i++ )
+        inVector_4.setElement( i, 2 );
+
+    VectorType outVector_4;
+    outVector_4.setSize( m_rows_4 );
+    for( IndexType j = 0; j < outVector_4.getSize(); j++ )
+        outVector_4.setElement( j, 0 );
+
+
+    m_4.vectorProduct( inVector_4, outVector_4 );
+
+
+    EXPECT_EQ( outVector_4.getElement( 0 ),  20 );
+    EXPECT_EQ( outVector_4.getElement( 1 ),  52 );
+    EXPECT_EQ( outVector_4.getElement( 2 ), 110 );
+    EXPECT_EQ( outVector_4.getElement( 3 ), 124 );
+    EXPECT_EQ( outVector_4.getElement( 4 ), 156 );
+    EXPECT_EQ( outVector_4.getElement( 5 ), 188 );
+    EXPECT_EQ( outVector_4.getElement( 6 ), 280 );
+    EXPECT_EQ( outVector_4.getElement( 7 ), 330 );
+
+
+   /*
+    * Sets up the following 8x8 sparse matrix:
+    *
+    *    /  1  2  3  0  4  5  0  1 \   6
+    *    |  0  6  0  7  0  0  0  1 |   3
+    *    |  0  8  9  0 10  0  0  1 |   4
+    *    |  0 11 12 13 14  0  0  1 |   5
+    *    |  0 15  0  0  0  0  0  1 |   2
+    *    |  0 16 17 18 19 20 21  1 |   7
+    *    | 22 23 24 25 26 27 28  1 |   8
+    *    \ 29 30 31 32 33 34 35 36 /   8
+    */
+
+    const IndexType m_rows_5 = 8;
+    const IndexType m_cols_5 = 8;
+
+    Matrix m_5;
+    m_5.reset();
+    m_5.setDimensions( m_rows_5, m_cols_5 );
+    typename Matrix::CompressedRowLengthsVector rowLengths_5;
+    rowLengths_5.setSize( m_rows_5 );
+    rowLengths_5.setElement(0, 6);
+    rowLengths_5.setElement(1, 3);
+    rowLengths_5.setElement(2, 4);
+    rowLengths_5.setElement(3, 5);
+    rowLengths_5.setElement(4, 2);
+    rowLengths_5.setElement(5, 7);
+    rowLengths_5.setElement(6, 8);
+    rowLengths_5.setElement(7, 8);
+    m_5.setCompressedRowLengths( rowLengths_5 );
+
+    RealType value_5 = 1;
+    for( IndexType i = 0; i < 3; i++ )   // 0th row
+        m_5.setElement( 0, i, value_5++ );
+
+    m_5.setElement( 0, 4, value_5++ );           // 0th row
+    m_5.setElement( 0, 5, value_5++ );
+
+    m_5.setElement( 1, 1, value_5++ );           // 1st row
+    m_5.setElement( 1, 3, value_5++ );
+
+    for( IndexType i = 1; i < 3; i++ )            // 2nd row
+        m_5.setElement( 2, i, value_5++ );
+
+    m_5.setElement( 2, 4, value_5++ );           // 2nd row
+
+    for( IndexType i = 1; i < 5; i++ )            // 3rd row
+        m_5.setElement( 3, i, value_5++ );
+
+    m_5.setElement( 4, 1, value_5++ );           // 4th row
+
+    for( IndexType i = 1; i < 7; i++ )            // 5th row
+        m_5.setElement( 5, i, value_5++ );
+
+    for( IndexType i = 0; i < 7; i++ )            // 6th row
+        m_5.setElement( 6, i, value_5++ );
+
+    for( IndexType i = 0; i < 8; i++ )            // 7th row
+        m_5.setElement( 7, i, value_5++ );
+
+    for( IndexType i = 0; i < 7; i++ )            // 1s at the end of rows
+        m_5.setElement( i, 7, 1);
+
+    VectorType inVector_5;
+    inVector_5.setSize( m_cols_5 );
+    for( IndexType i = 0; i < inVector_5.getSize(); i++ )
+        inVector_5.setElement( i, 2 );
+
+    VectorType outVector_5;
+    outVector_5.setSize( m_rows_5 );
+    for( IndexType j = 0; j < outVector_5.getSize(); j++ )
+        outVector_5.setElement( j, 0 );
+
+    m_5.vectorProduct( inVector_5, outVector_5 );
+
+    EXPECT_EQ( outVector_5.getElement( 0 ),  32 );
+    EXPECT_EQ( outVector_5.getElement( 1 ),  28 );
+    EXPECT_EQ( outVector_5.getElement( 2 ),  56 );
+    EXPECT_EQ( outVector_5.getElement( 3 ), 102 );
+    EXPECT_EQ( outVector_5.getElement( 4 ),  32 );
+    EXPECT_EQ( outVector_5.getElement( 5 ), 224 );
+    EXPECT_EQ( outVector_5.getElement( 6 ), 352 );
+    EXPECT_EQ( outVector_5.getElement( 7 ), 520 );
+}
+
+template< typename Matrix >
+void test_RowsReduction()
+{
+   using RealType = typename Matrix::RealType;
+   using DeviceType = typename Matrix::DeviceType;
+   using IndexType = typename Matrix::IndexType;
+
+   /*
+    * Sets up the following 8x8 sparse matrix:
+    *
+    *    /  1  2  3  0  4  5  0  1 \   6
+    *    |  0  6  0  7  0  0  0  1 |   3
+    *    |  0  8  9  0 10  0  0  1 |   4
+    *    |  0 11 12 13 14  0  0  1 |   5
+    *    |  0 15  0  0  0  0  0  1 |   2
+    *    |  0 16 17 18 19 20 21  1 |   7
+    *    | 22 23 24 25 26 27 28  1 |   8
+    *    \ 29 30 31 32 33 34 35 36 /   8
+    */
+
+   const IndexType rows = 8;
+   const IndexType cols = 8;
+
+   Matrix m;
+   m.setDimensions( rows, cols );
+   typename Matrix::RowsCapacitiesType rowsCapacities( rows );
+   //rowLengths.setSize( rows );
+   rowsCapacities.setElement(0, 6);
+   rowsCapacities.setElement(1, 3);
+   rowsCapacities.setElement(2, 4);
+   rowsCapacities.setElement(3, 5);
+   rowsCapacities.setElement(4, 2);
+   rowsCapacities.setElement(5, 7);
+   rowsCapacities.setElement(6, 8);
+   rowsCapacities.setElement(7, 8);
+   m.setCompressedRowLengths( rowsCapacities );
+
+   RealType value = 1;
+   for( IndexType i = 0; i < 3; i++ )   // 0th row
+      m.setElement( 0, i, value++ );
+
+   m.setElement( 0, 4, value++ );           // 0th row
+   m.setElement( 0, 5, value++ );
+
+   m.setElement( 1, 1, value++ );           // 1st row
+   m.setElement( 1, 3, value++ );
+
+   for( IndexType i = 1; i < 3; i++ )            // 2nd row
+      m.setElement( 2, i, value++ );
+
+   m.setElement( 2, 4, value++ );           // 2nd row
+
+   for( IndexType i = 1; i < 5; i++ )            // 3rd row
+      m.setElement( 3, i, value++ );
+
+   m.setElement( 4, 1, value++ );           // 4th row
+
+   for( IndexType i = 1; i < 7; i++ )            // 5th row
+      m.setElement( 5, i, value++ );
+
+   for( IndexType i = 0; i < 7; i++ )            // 6th row
+      m.setElement( 6, i, value++ );
+
+   for( IndexType i = 0; i < 8; i++ )            // 7th row
+       m.setElement( 7, i, value++ );
+
+   for( IndexType i = 0; i < 7; i++ )            // 1s at the end of rows
+      m.setElement( i, 7, 1);
+
+   ////
+   // Compute number of non-zero elements in rows.
+   typename Matrix::RowsCapacitiesType rowLengths( rows );
+   auto rowLengths_view = rowLengths.getView();
+   auto fetch = [] __cuda_callable__ ( IndexType row, IndexType column, const RealType& value ) -> IndexType {
+      return ( value != 0.0 );
+   };
+   auto reduce = [] __cuda_callable__ ( IndexType& aux, const IndexType a ) {
+      aux += a;
+   };
+   auto keep = [=] __cuda_callable__ ( const IndexType rowIdx, const IndexType value ) mutable {
+      rowLengths_view[ rowIdx ] = value;
+   };
+   m.allRowsReduction( fetch, reduce, keep, 0 );
+   EXPECT_EQ( rowsCapacities, rowLengths );
+   m.getCompressedRowLengths( rowLengths );
+   EXPECT_EQ( rowsCapacities, rowLengths );
+
+   ////
+   // Compute max norm
+   TNL::Containers::Vector< RealType, DeviceType, IndexType > rowSums( rows );
+   auto rowSums_view = rowSums.getView();
+   auto max_fetch = [] __cuda_callable__ ( IndexType row, IndexType column, const RealType& value ) -> IndexType {
+      return abs( value );
+   };
+   auto max_reduce = [] __cuda_callable__ ( IndexType& aux, const IndexType a ) {
+      aux += a;
+   };
+   auto max_keep = [=] __cuda_callable__ ( const IndexType rowIdx, const IndexType value ) mutable {
+      rowSums_view[ rowIdx ] = value;
+   };
+   m.allRowsReduction( max_fetch, max_reduce, max_keep, 0 );
+   const RealType maxNorm = TNL::max( rowSums );
+   EXPECT_EQ( maxNorm, 260 ) ; // 29+30+31+32+33+34+35+36
+}
+
+template< typename Matrix >
+void test_PerformSORIteration()
+{
+    using RealType = typename Matrix::RealType;
+    using DeviceType = typename Matrix::DeviceType;
+    using IndexType = typename Matrix::IndexType;
+
+/*
+ * Sets up the following 4x4 sparse matrix:
+ *
+ *    /  4  1  0  0 \
+ *    |  1  4  1  0 |
+ *    |  0  1  4  1 |
+ *    \  0  0  1  4 /
+ */
+
+    const IndexType m_rows = 4;
+    const IndexType m_cols = 4;
+
+    Matrix m;
+    m.reset();
+    m.setDimensions( m_rows, m_cols );
+    typename Matrix::CompressedRowLengthsVector rowLengths;
+    rowLengths.setSize( m_rows );
+    rowLengths.setValue( 3 );
+    m.setCompressedRowLengths( rowLengths );
+
+    m.setElement( 0, 0, 4.0 );        // 0th row
+    m.setElement( 0, 1, 1.0);
+
+    m.setElement( 1, 0, 1.0 );        // 1st row
+    m.setElement( 1, 1, 4.0 );
+    m.setElement( 1, 2, 1.0 );
+
+    m.setElement( 2, 1, 1.0 );        // 2nd row
+    m.setElement( 2, 2, 4.0 );
+    m.setElement( 2, 3, 1.0 );
+
+    m.setElement( 3, 2, 1.0 );        // 3rd row
+    m.setElement( 3, 3, 4.0 );
+
+    RealType bVector [ 4 ] = { 1, 1, 1, 1 };
+    RealType xVector [ 4 ] = { 1, 1, 1, 1 };
+
+    IndexType row = 0;
+    RealType omega = 1;
+
+
+    m.performSORIteration( bVector, row++, xVector, omega);
+
+    EXPECT_EQ( xVector[ 0 ], 0.0 );
+    EXPECT_EQ( xVector[ 1 ], 1.0 );
+    EXPECT_EQ( xVector[ 2 ], 1.0 );
+    EXPECT_EQ( xVector[ 3 ], 1.0 );
+
+
+    m.performSORIteration( bVector, row++, xVector, omega);
+
+    EXPECT_EQ( xVector[ 0 ], 0.0 );
+    EXPECT_EQ( xVector[ 1 ], 0.0 );
+    EXPECT_EQ( xVector[ 2 ], 1.0 );
+    EXPECT_EQ( xVector[ 3 ], 1.0 );
+
+
+    m.performSORIteration( bVector, row++, xVector, omega);
+
+    EXPECT_EQ( xVector[ 0 ], 0.0 );
+    EXPECT_EQ( xVector[ 1 ], 0.0 );
+    EXPECT_EQ( xVector[ 2 ], 0.0 );
+    EXPECT_EQ( xVector[ 3 ], 1.0 );
+
+
+    m.performSORIteration( bVector, row++, xVector, omega);
+
+    EXPECT_EQ( xVector[ 0 ], 0.0 );
+    EXPECT_EQ( xVector[ 1 ], 0.0 );
+    EXPECT_EQ( xVector[ 2 ], 0.0 );
+    EXPECT_EQ( xVector[ 3 ], 0.25 );
+}
+
+// This test is only for AdEllpack
+template< typename Matrix >
+void test_OperatorEquals()
+{
+   using RealType = typename Matrix::RealType;
+   using DeviceType = typename Matrix::DeviceType;
+   using IndexType = typename Matrix::IndexType;
+
+   if( std::is_same< DeviceType, TNL::Devices::Cuda >::value )
+       return;
+   else
+   {
+       using AdELL_host = TNL::Matrices::AdEllpack< RealType, TNL::Devices::Host, IndexType >;
+       using AdELL_cuda = TNL::Matrices::AdEllpack< RealType, TNL::Devices::Cuda, IndexType >;
+
+       /*
+        * Sets up the following 8x8 sparse matrix:
+        *
+        *    /  1  2  3  0  4  5  0  1 \   6
+        *    |  0  6  0  7  0  0  0  1 |   3
+        *    |  0  8  9  0 10  0  0  1 |   4
+        *    |  0 11 12 13 14  0  0  1 |   5
+        *    |  0 15  0  0  0  0  0  1 |   2
+        *    |  0 16 17 18 19 20 21  1 |   7
+        *    | 22 23 24 25 26 27 28  1 |   8
+        *    \ 29 30 31 32 33 34 35 36 /   8
+        */
+
+        const IndexType m_rows = 8;
+        const IndexType m_cols = 8;
+
+        AdELL_host m_host;
+
+        m_host.reset();
+        m_host.setDimensions( m_rows, m_cols );
+        typename AdELL_host::CompressedRowLengthsVector rowLengths;
+        rowLengths.setSize( m_rows );
+        rowLengths.setElement(0, 6);
+        rowLengths.setElement(1, 3);
+        rowLengths.setElement(2, 4);
+        rowLengths.setElement(3, 5);
+        rowLengths.setElement(4, 2);
+        rowLengths.setElement(5, 7);
+        rowLengths.setElement(6, 8);
+        rowLengths.setElement(7, 8);
+        m_host.setCompressedRowLengths( rowLengths );
+
+        RealType value = 1;
+        for( IndexType i = 0; i < 3; i++ )   // 0th row
+            m_host.setElement( 0, i, value++ );
+
+        m_host.setElement( 0, 4, value++ );           // 0th row
+        m_host.setElement( 0, 5, value++ );
+
+        m_host.setElement( 1, 1, value++ );           // 1st row
+        m_host.setElement( 1, 3, value++ );
+
+        for( IndexType i = 1; i < 3; i++ )            // 2nd row
+            m_host.setElement( 2, i, value++ );
+
+        m_host.setElement( 2, 4, value++ );           // 2nd row
+
+
+        for( IndexType i = 1; i < 5; i++ )            // 3rd row
+            m_host.setElement( 3, i, value++ );
+
+        m_host.setElement( 4, 1, value++ );           // 4th row
+
+        for( IndexType i = 1; i < 7; i++ )            // 5th row
+            m_host.setElement( 5, i, value++ );
+
+        for( IndexType i = 0; i < 7; i++ )            // 6th row
+            m_host.setElement( 6, i, value++ );
+
+        for( IndexType i = 0; i < 8; i++ )            // 7th row
+            m_host.setElement( 7, i, value++ );
+
+        for( IndexType i = 0; i < 7; i++ )            // 1s at the end or rows: 5, 6
+            m_host.setElement( i, 7, 1);
+
+        EXPECT_EQ( m_host.getElement( 0, 0 ),  1 );
+        EXPECT_EQ( m_host.getElement( 0, 1 ),  2 );
+        EXPECT_EQ( m_host.getElement( 0, 2 ),  3 );
+        EXPECT_EQ( m_host.getElement( 0, 3 ),  0 );
+        EXPECT_EQ( m_host.getElement( 0, 4 ),  4 );
+        EXPECT_EQ( m_host.getElement( 0, 5 ),  5 );
+        EXPECT_EQ( m_host.getElement( 0, 6 ),  0 );
+        EXPECT_EQ( m_host.getElement( 0, 7 ),  1 );
+
+        EXPECT_EQ( m_host.getElement( 1, 0 ),  0 );
+        EXPECT_EQ( m_host.getElement( 1, 1 ),  6 );
+        EXPECT_EQ( m_host.getElement( 1, 2 ),  0 );
+        EXPECT_EQ( m_host.getElement( 1, 3 ),  7 );
+        EXPECT_EQ( m_host.getElement( 1, 4 ),  0 );
+        EXPECT_EQ( m_host.getElement( 1, 5 ),  0 );
+        EXPECT_EQ( m_host.getElement( 1, 6 ),  0 );
+        EXPECT_EQ( m_host.getElement( 1, 7 ),  1 );
+
+        EXPECT_EQ( m_host.getElement( 2, 0 ),  0 );
+        EXPECT_EQ( m_host.getElement( 2, 1 ),  8 );
+        EXPECT_EQ( m_host.getElement( 2, 2 ),  9 );
+        EXPECT_EQ( m_host.getElement( 2, 3 ),  0 );
+        EXPECT_EQ( m_host.getElement( 2, 4 ), 10 );
+        EXPECT_EQ( m_host.getElement( 2, 5 ),  0 );
+        EXPECT_EQ( m_host.getElement( 2, 6 ),  0 );
+        EXPECT_EQ( m_host.getElement( 2, 7 ),  1 );
+
+        EXPECT_EQ( m_host.getElement( 3, 0 ),  0 );
+        EXPECT_EQ( m_host.getElement( 3, 1 ), 11 );
+        EXPECT_EQ( m_host.getElement( 3, 2 ), 12 );
+        EXPECT_EQ( m_host.getElement( 3, 3 ), 13 );
+        EXPECT_EQ( m_host.getElement( 3, 4 ), 14 );
+        EXPECT_EQ( m_host.getElement( 3, 5 ),  0 );
+        EXPECT_EQ( m_host.getElement( 3, 6 ),  0 );
+        EXPECT_EQ( m_host.getElement( 3, 7 ),  1 );
+
+        EXPECT_EQ( m_host.getElement( 4, 0 ),  0 );
+        EXPECT_EQ( m_host.getElement( 4, 1 ), 15 );
+        EXPECT_EQ( m_host.getElement( 4, 2 ),  0 );
+        EXPECT_EQ( m_host.getElement( 4, 3 ),  0 );
+        EXPECT_EQ( m_host.getElement( 4, 4 ),  0 );
+        EXPECT_EQ( m_host.getElement( 4, 5 ),  0 );
+        EXPECT_EQ( m_host.getElement( 4, 6 ),  0 );
+        EXPECT_EQ( m_host.getElement( 4, 7 ),  1 );
+
+        EXPECT_EQ( m_host.getElement( 5, 0 ),  0 );
+        EXPECT_EQ( m_host.getElement( 5, 1 ), 16 );
+        EXPECT_EQ( m_host.getElement( 5, 2 ), 17 );
+        EXPECT_EQ( m_host.getElement( 5, 3 ), 18 );
+        EXPECT_EQ( m_host.getElement( 5, 4 ), 19 );
+        EXPECT_EQ( m_host.getElement( 5, 5 ), 20 );
+        EXPECT_EQ( m_host.getElement( 5, 6 ), 21 );
+        EXPECT_EQ( m_host.getElement( 5, 7 ),  1 );
+
+        EXPECT_EQ( m_host.getElement( 6, 0 ), 22 );
+        EXPECT_EQ( m_host.getElement( 6, 1 ), 23 );
+        EXPECT_EQ( m_host.getElement( 6, 2 ), 24 );
+        EXPECT_EQ( m_host.getElement( 6, 3 ), 25 );
+        EXPECT_EQ( m_host.getElement( 6, 4 ), 26 );
+        EXPECT_EQ( m_host.getElement( 6, 5 ), 27 );
+        EXPECT_EQ( m_host.getElement( 6, 6 ), 28 );
+        EXPECT_EQ( m_host.getElement( 6, 7 ),  1 );
+
+        EXPECT_EQ( m_host.getElement( 7, 0 ), 29 );
+        EXPECT_EQ( m_host.getElement( 7, 1 ), 30 );
+        EXPECT_EQ( m_host.getElement( 7, 2 ), 31 );
+        EXPECT_EQ( m_host.getElement( 7, 3 ), 32 );
+        EXPECT_EQ( m_host.getElement( 7, 4 ), 33 );
+        EXPECT_EQ( m_host.getElement( 7, 5 ), 34 );
+        EXPECT_EQ( m_host.getElement( 7, 6 ), 35 );
+        EXPECT_EQ( m_host.getElement( 7, 7 ), 36 );
+
+        AdELL_cuda m_cuda;
+
+        // Copy the host matrix into the cuda matrix
+        m_cuda = m_host;
+
+        // Reset the host matrix
+        m_host.reset();
+
+        // Copy the cuda matrix back into the host matrix
+        m_host = m_cuda;
+
+        // Check the newly created double-copy host matrix
+        EXPECT_EQ( m_host.getElement( 0, 0 ),  1 );
+        EXPECT_EQ( m_host.getElement( 0, 1 ),  2 );
+        EXPECT_EQ( m_host.getElement( 0, 2 ),  3 );
+        EXPECT_EQ( m_host.getElement( 0, 3 ),  0 );
+        EXPECT_EQ( m_host.getElement( 0, 4 ),  4 );
+        EXPECT_EQ( m_host.getElement( 0, 5 ),  5 );
+        EXPECT_EQ( m_host.getElement( 0, 6 ),  0 );
+        EXPECT_EQ( m_host.getElement( 0, 7 ),  1 );
+
+        EXPECT_EQ( m_host.getElement( 1, 0 ),  0 );
+        EXPECT_EQ( m_host.getElement( 1, 1 ),  6 );
+        EXPECT_EQ( m_host.getElement( 1, 2 ),  0 );
+        EXPECT_EQ( m_host.getElement( 1, 3 ),  7 );
+        EXPECT_EQ( m_host.getElement( 1, 4 ),  0 );
+        EXPECT_EQ( m_host.getElement( 1, 5 ),  0 );
+        EXPECT_EQ( m_host.getElement( 1, 6 ),  0 );
+        EXPECT_EQ( m_host.getElement( 1, 7 ),  1 );
+
+        EXPECT_EQ( m_host.getElement( 2, 0 ),  0 );
+        EXPECT_EQ( m_host.getElement( 2, 1 ),  8 );
+        EXPECT_EQ( m_host.getElement( 2, 2 ),  9 );
+        EXPECT_EQ( m_host.getElement( 2, 3 ),  0 );
+        EXPECT_EQ( m_host.getElement( 2, 4 ), 10 );
+        EXPECT_EQ( m_host.getElement( 2, 5 ),  0 );
+        EXPECT_EQ( m_host.getElement( 2, 6 ),  0 );
+        EXPECT_EQ( m_host.getElement( 2, 7 ),  1 );
+
+        EXPECT_EQ( m_host.getElement( 3, 0 ),  0 );
+        EXPECT_EQ( m_host.getElement( 3, 1 ), 11 );
+        EXPECT_EQ( m_host.getElement( 3, 2 ), 12 );
+        EXPECT_EQ( m_host.getElement( 3, 3 ), 13 );
+        EXPECT_EQ( m_host.getElement( 3, 4 ), 14 );
+        EXPECT_EQ( m_host.getElement( 3, 5 ),  0 );
+        EXPECT_EQ( m_host.getElement( 3, 6 ),  0 );
+        EXPECT_EQ( m_host.getElement( 3, 7 ),  1 );
+
+        EXPECT_EQ( m_host.getElement( 4, 0 ),  0 );
+        EXPECT_EQ( m_host.getElement( 4, 1 ), 15 );
+        EXPECT_EQ( m_host.getElement( 4, 2 ),  0 );
+        EXPECT_EQ( m_host.getElement( 4, 3 ),  0 );
+        EXPECT_EQ( m_host.getElement( 4, 4 ),  0 );
+        EXPECT_EQ( m_host.getElement( 4, 5 ),  0 );
+        EXPECT_EQ( m_host.getElement( 4, 6 ),  0 );
+        EXPECT_EQ( m_host.getElement( 4, 7 ),  1 );
+
+        EXPECT_EQ( m_host.getElement( 5, 0 ),  0 );
+        EXPECT_EQ( m_host.getElement( 5, 1 ), 16 );
+        EXPECT_EQ( m_host.getElement( 5, 2 ), 17 );
+        EXPECT_EQ( m_host.getElement( 5, 3 ), 18 );
+        EXPECT_EQ( m_host.getElement( 5, 4 ), 19 );
+        EXPECT_EQ( m_host.getElement( 5, 5 ), 20 );
+        EXPECT_EQ( m_host.getElement( 5, 6 ), 21 );
+        EXPECT_EQ( m_host.getElement( 5, 7 ),  1 );
+
+        EXPECT_EQ( m_host.getElement( 6, 0 ), 22 );
+        EXPECT_EQ( m_host.getElement( 6, 1 ), 23 );
+        EXPECT_EQ( m_host.getElement( 6, 2 ), 24 );
+        EXPECT_EQ( m_host.getElement( 6, 3 ), 25 );
+        EXPECT_EQ( m_host.getElement( 6, 4 ), 26 );
+        EXPECT_EQ( m_host.getElement( 6, 5 ), 27 );
+        EXPECT_EQ( m_host.getElement( 6, 6 ), 28 );
+        EXPECT_EQ( m_host.getElement( 6, 7 ),  1 );
+
+        EXPECT_EQ( m_host.getElement( 7, 0 ), 29 );
+        EXPECT_EQ( m_host.getElement( 7, 1 ), 30 );
+        EXPECT_EQ( m_host.getElement( 7, 2 ), 31 );
+        EXPECT_EQ( m_host.getElement( 7, 3 ), 32 );
+        EXPECT_EQ( m_host.getElement( 7, 4 ), 33 );
+        EXPECT_EQ( m_host.getElement( 7, 5 ), 34 );
+        EXPECT_EQ( m_host.getElement( 7, 6 ), 35 );
+        EXPECT_EQ( m_host.getElement( 7, 7 ), 36 );
+
+        // Try vectorProduct with copied cuda matrix to see if it works correctly.
+        using VectorType = TNL::Containers::Vector< RealType, TNL::Devices::Cuda, IndexType >;
+
+        VectorType inVector;
+        inVector.setSize( m_cols );
+        for( IndexType i = 0; i < inVector.getSize(); i++ )
+            inVector.setElement( i, 2 );
+
+        VectorType outVector;
+        outVector.setSize( m_rows );
+        for( IndexType j = 0; j < outVector.getSize(); j++ )
+            outVector.setElement( j, 0 );
+
+        m_cuda.vectorProduct( inVector, outVector );
+
+        EXPECT_EQ( outVector.getElement( 0 ),  32 );
+        EXPECT_EQ( outVector.getElement( 1 ),  28 );
+        EXPECT_EQ( outVector.getElement( 2 ),  56 );
+        EXPECT_EQ( outVector.getElement( 3 ), 102 );
+        EXPECT_EQ( outVector.getElement( 4 ),  32 );
+        EXPECT_EQ( outVector.getElement( 5 ), 224 );
+        EXPECT_EQ( outVector.getElement( 6 ), 352 );
+        EXPECT_EQ( outVector.getElement( 7 ), 520 );
+   }
+}
+
+template< typename Matrix >
+void test_SaveAndLoad( const char* filename )
+{
+   using RealType = typename Matrix::RealType;
+   using DeviceType = typename Matrix::DeviceType;
+   using IndexType = typename Matrix::IndexType;
+
+   /*
+    * Sets up the following 4x4 sparse matrix:
+    *
+    *    /  1  2  3  0 \
+    *    |  0  4  0  5 |
+    *    |  6  7  8  0 |
+    *    \  0  9 10 11 /
+    */
+
+    const IndexType m_rows = 4;
+    const IndexType m_cols = 4;
+
+    Matrix savedMatrix;
+    savedMatrix.reset();
+    savedMatrix.setDimensions( m_rows, m_cols );
+    typename Matrix::CompressedRowLengthsVector rowLengths;
+    rowLengths.setSize( m_rows );
+    rowLengths.setValue( 3 );
+    savedMatrix.setCompressedRowLengths( rowLengths );
+
+    RealType value = 1;
+    for( IndexType i = 0; i < m_cols - 1; i++ )   // 0th row
+        savedMatrix.setElement( 0, i, value++ );
+
+    savedMatrix.setElement( 1, 1, value++ );
+    savedMatrix.setElement( 1, 3, value++ );      // 1st row
+
+    for( IndexType i = 0; i < m_cols - 1; i++ )   // 2nd row
+        savedMatrix.setElement( 2, i, value++ );
+
+    for( IndexType i = 1; i < m_cols; i++ )       // 3rd row
+        savedMatrix.setElement( 3, i, value++ );
+
+    ASSERT_NO_THROW( savedMatrix.save( filename ) );
+
+    Matrix loadedMatrix;
+    loadedMatrix.reset();
+    loadedMatrix.setDimensions( m_rows, m_cols );
+    typename Matrix::CompressedRowLengthsVector rowLengths2;
+    rowLengths2.setSize( m_rows );
+    rowLengths2.setValue( 3 );
+    loadedMatrix.setCompressedRowLengths( rowLengths2 );
+
+
+    ASSERT_NO_THROW( loadedMatrix.load( filename ) );
+
+
+    EXPECT_EQ( savedMatrix.getElement( 0, 0 ), loadedMatrix.getElement( 0, 0 ) );
+    EXPECT_EQ( savedMatrix.getElement( 0, 1 ), loadedMatrix.getElement( 0, 1 ) );
+    EXPECT_EQ( savedMatrix.getElement( 0, 2 ), loadedMatrix.getElement( 0, 2 ) );
+    EXPECT_EQ( savedMatrix.getElement( 0, 3 ), loadedMatrix.getElement( 0, 3 ) );
+
+    EXPECT_EQ( savedMatrix.getElement( 1, 0 ), loadedMatrix.getElement( 1, 0 ) );
+    EXPECT_EQ( savedMatrix.getElement( 1, 1 ), loadedMatrix.getElement( 1, 1 ) );
+    EXPECT_EQ( savedMatrix.getElement( 1, 2 ), loadedMatrix.getElement( 1, 2 ) );
+    EXPECT_EQ( savedMatrix.getElement( 1, 3 ), loadedMatrix.getElement( 1, 3 ) );
+
+    EXPECT_EQ( savedMatrix.getElement( 2, 0 ), loadedMatrix.getElement( 2, 0 ) );
+    EXPECT_EQ( savedMatrix.getElement( 2, 1 ), loadedMatrix.getElement( 2, 1 ) );
+    EXPECT_EQ( savedMatrix.getElement( 2, 2 ), loadedMatrix.getElement( 2, 2 ) );
+    EXPECT_EQ( savedMatrix.getElement( 2, 3 ), loadedMatrix.getElement( 2, 3 ) );
+
+    EXPECT_EQ( savedMatrix.getElement( 3, 0 ), loadedMatrix.getElement( 3, 0 ) );
+    EXPECT_EQ( savedMatrix.getElement( 3, 1 ), loadedMatrix.getElement( 3, 1 ) );
+    EXPECT_EQ( savedMatrix.getElement( 3, 2 ), loadedMatrix.getElement( 3, 2 ) );
+    EXPECT_EQ( savedMatrix.getElement( 3, 3 ), loadedMatrix.getElement( 3, 3 ) );
+
+    EXPECT_EQ( savedMatrix.getElement( 0, 0 ),  1 );
+    EXPECT_EQ( savedMatrix.getElement( 0, 1 ),  2 );
+    EXPECT_EQ( savedMatrix.getElement( 0, 2 ),  3 );
+    EXPECT_EQ( savedMatrix.getElement( 0, 3 ),  0 );
+
+    EXPECT_EQ( savedMatrix.getElement( 1, 0 ),  0 );
+    EXPECT_EQ( savedMatrix.getElement( 1, 1 ),  4 );
+    EXPECT_EQ( savedMatrix.getElement( 1, 2 ),  0 );
+    EXPECT_EQ( savedMatrix.getElement( 1, 3 ),  5 );
+
+    EXPECT_EQ( savedMatrix.getElement( 2, 0 ),  6 );
+    EXPECT_EQ( savedMatrix.getElement( 2, 1 ),  7 );
+    EXPECT_EQ( savedMatrix.getElement( 2, 2 ),  8 );
+    EXPECT_EQ( savedMatrix.getElement( 2, 3 ),  0 );
+
+    EXPECT_EQ( savedMatrix.getElement( 3, 0 ),  0 );
+    EXPECT_EQ( savedMatrix.getElement( 3, 1 ),  9 );
+    EXPECT_EQ( savedMatrix.getElement( 3, 2 ), 10 );
+    EXPECT_EQ( savedMatrix.getElement( 3, 3 ), 11 );
+
+    EXPECT_EQ( std::remove( filename ), 0 );
+}
+
+template< typename Matrix >
+void test_Print()
+{
+    using RealType = typename Matrix::RealType;
+    using DeviceType = typename Matrix::DeviceType;
+    using IndexType = typename Matrix::IndexType;
+
+/*
+ * Sets up the following 5x4 sparse matrix:
+ *
+ *    /  1  2  3  0 \
+ *    |  0  0  0  4 |
+ *    |  5  6  7  0 |
+ *    |  0  8  9 10 |
+ *    \  0  0 11 12 /
+ */
+
+    const IndexType m_rows = 5;
+    const IndexType m_cols = 4;
+
+    Matrix m;
+    m.reset();
+    m.setDimensions( m_rows, m_cols );
+    typename Matrix::CompressedRowLengthsVector rowLengths;
+    rowLengths.setSize( m_rows );
+    rowLengths.setValue( 3 );
+    m.setCompressedRowLengths( rowLengths );
+
+    RealType value = 1;
+    for( IndexType i = 0; i < m_cols - 1; i++ )   // 0th row
+        m.setElement( 0, i, value++ );
+
+    m.setElement( 1, 3, value++ );      // 1st row
+
+    for( IndexType i = 0; i < m_cols - 1; i++ )   // 2nd row
+        m.setElement( 2, i, value++ );
+
+    for( IndexType i = 1; i < m_cols; i++ )       // 3rd row
+        m.setElement( 3, i, value++ );
+
+    for( IndexType i = 2; i < m_cols; i++ )       // 4th row
+        m.setElement( 4, i, value++ );
+
+    #include <sstream>
+    std::stringstream printed;
+    std::stringstream couted;
+
+    //change the underlying buffer and save the old buffer
+    auto old_buf = std::cout.rdbuf(printed.rdbuf());
+
+    m.print( std::cout ); //all the std::cout goes to ss
+
+    std::cout.rdbuf(old_buf); //reset
+
+    couted << "Row: 0 ->  Col:0->1	 Col:1->2	 Col:2->3\t\n"
+               "Row: 1 ->  Col:3->4\t\n"
+               "Row: 2 ->  Col:0->5	 Col:1->6	 Col:2->7\t\n"
+               "Row: 3 ->  Col:1->8	 Col:2->9	 Col:3->10\t\n"
+               "Row: 4 ->  Col:2->11	 Col:3->12\t\n";
+
+
+    EXPECT_EQ( printed.str(), couted.str() );
+}
+
+#endif
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_AdEllpack.cpp b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_AdEllpack.cpp
similarity index 100%
rename from src/UnitTests/Matrices/SparseMatrixTest_AdEllpack.cpp
rename to src/UnitTests/Matrices/Legacy/SparseMatrixTest_AdEllpack.cpp
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_AdEllpack.cu b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_AdEllpack.cu
similarity index 100%
rename from src/UnitTests/Matrices/SparseMatrixTest_AdEllpack.cu
rename to src/UnitTests/Matrices/Legacy/SparseMatrixTest_AdEllpack.cu
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_AdEllpack.h b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_AdEllpack.h
similarity index 99%
rename from src/UnitTests/Matrices/SparseMatrixTest_AdEllpack.h
rename to src/UnitTests/Matrices/Legacy/SparseMatrixTest_AdEllpack.h
index 2169b96df3015ba1d5a501cedfe755165c5ea2d0..d2d268dacef5b602d3ac7208bdd6f0148ed6068a 100644
--- a/src/UnitTests/Matrices/SparseMatrixTest_AdEllpack.h
+++ b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_AdEllpack.h
@@ -133,4 +133,4 @@ TYPED_TEST( AdEllpackMatrixTest, printTest )
 #endif
 
 
-#include "../main.h"
+#include "../../main.h"
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_BiEllpack.cpp b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_BiEllpack.cpp
similarity index 100%
rename from src/UnitTests/Matrices/SparseMatrixTest_BiEllpack.cpp
rename to src/UnitTests/Matrices/Legacy/SparseMatrixTest_BiEllpack.cpp
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_BiEllpack.cu b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_BiEllpack.cu
similarity index 100%
rename from src/UnitTests/Matrices/SparseMatrixTest_BiEllpack.cu
rename to src/UnitTests/Matrices/Legacy/SparseMatrixTest_BiEllpack.cu
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_BiEllpack.h b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_BiEllpack.h
similarity index 99%
rename from src/UnitTests/Matrices/SparseMatrixTest_BiEllpack.h
rename to src/UnitTests/Matrices/Legacy/SparseMatrixTest_BiEllpack.h
index c74fa635f38af1a883b80969194a64eaafa73984..9dab63c1a983988fdf5b5722c7df73b1331e5ca9 100644
--- a/src/UnitTests/Matrices/SparseMatrixTest_BiEllpack.h
+++ b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_BiEllpack.h
@@ -142,4 +142,4 @@ TYPED_TEST( BiEllpackMatrixTest, printTest )
 }
 #endif // HAVE_GTEST
 
-#include "../main.h"
+#include "../../main.h"
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_CSR.cpp b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_CSR.cpp
similarity index 100%
rename from src/UnitTests/Matrices/SparseMatrixTest_CSR.cpp
rename to src/UnitTests/Matrices/Legacy/SparseMatrixTest_CSR.cpp
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_CSR.cu b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_CSR.cu
similarity index 100%
rename from src/UnitTests/Matrices/SparseMatrixTest_CSR.cu
rename to src/UnitTests/Matrices/Legacy/SparseMatrixTest_CSR.cu
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_CSR.h b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_CSR.h
similarity index 99%
rename from src/UnitTests/Matrices/SparseMatrixTest_CSR.h
rename to src/UnitTests/Matrices/Legacy/SparseMatrixTest_CSR.h
index c9dfc770f86dc071ef1ecc8ccd7e9c76111308ff..3cae12e3a5be3da5547e5153ed8c32649b7dbe07 100644
--- a/src/UnitTests/Matrices/SparseMatrixTest_CSR.h
+++ b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_CSR.h
@@ -136,4 +136,4 @@ TYPED_TEST( CSRMatrixTest, printTest )
 
 #endif
 
-#include "../main.h"
+#include "../../main.h"
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_ChunkedEllpack.cpp b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_ChunkedEllpack.cpp
similarity index 100%
rename from src/UnitTests/Matrices/SparseMatrixTest_ChunkedEllpack.cpp
rename to src/UnitTests/Matrices/Legacy/SparseMatrixTest_ChunkedEllpack.cpp
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_ChunkedEllpack.cu b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_ChunkedEllpack.cu
similarity index 100%
rename from src/UnitTests/Matrices/SparseMatrixTest_ChunkedEllpack.cu
rename to src/UnitTests/Matrices/Legacy/SparseMatrixTest_ChunkedEllpack.cu
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_ChunkedEllpack.h b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_ChunkedEllpack.h
similarity index 99%
rename from src/UnitTests/Matrices/SparseMatrixTest_ChunkedEllpack.h
rename to src/UnitTests/Matrices/Legacy/SparseMatrixTest_ChunkedEllpack.h
index 45801fa3a9adfd5353fcb6dcec551855b698f52e..a3c049910b5752d8fa6f7e32f4090e20364794fb 100644
--- a/src/UnitTests/Matrices/SparseMatrixTest_ChunkedEllpack.h
+++ b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_ChunkedEllpack.h
@@ -144,4 +144,4 @@ TYPED_TEST( ChunkedEllpackMatrixTest, printTest )
 
 #endif
 
-#include "../main.h"
+#include "../../main.h"
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_Ellpack.cpp b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_Ellpack.cpp
similarity index 100%
rename from src/UnitTests/Matrices/SparseMatrixTest_Ellpack.cpp
rename to src/UnitTests/Matrices/Legacy/SparseMatrixTest_Ellpack.cpp
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_Ellpack.cu b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_Ellpack.cu
similarity index 100%
rename from src/UnitTests/Matrices/SparseMatrixTest_Ellpack.cu
rename to src/UnitTests/Matrices/Legacy/SparseMatrixTest_Ellpack.cu
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_Ellpack.h b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_Ellpack.h
similarity index 99%
rename from src/UnitTests/Matrices/SparseMatrixTest_Ellpack.h
rename to src/UnitTests/Matrices/Legacy/SparseMatrixTest_Ellpack.h
index 26d270a3d3e62061f7cfcbb7c4bf4015c36c01d6..fa6b2027c63380c08b5c7a8e59cf9cb27f7cdd69 100644
--- a/src/UnitTests/Matrices/SparseMatrixTest_Ellpack.h
+++ b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_Ellpack.h
@@ -136,4 +136,4 @@ TYPED_TEST( EllpackMatrixTest, printTest )
 
 #endif
 
-#include "../main.h"
+#include "../../main.h"
diff --git a/src/UnitTests/Matrices/Legacy/SparseMatrixTest_SlicedEllpack.cpp b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_SlicedEllpack.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..40e2e94b81ca64051ddceee82f46dd2d20e66e42
--- /dev/null
+++ b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_SlicedEllpack.cpp
@@ -0,0 +1 @@
+#include "SparseMatrixTest_SlicedEllpack.h"
diff --git a/src/UnitTests/Matrices/Legacy/SparseMatrixTest_SlicedEllpack.cu b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_SlicedEllpack.cu
new file mode 100644
index 0000000000000000000000000000000000000000..40e2e94b81ca64051ddceee82f46dd2d20e66e42
--- /dev/null
+++ b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_SlicedEllpack.cu
@@ -0,0 +1 @@
+#include "SparseMatrixTest_SlicedEllpack.h"
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_SlicedEllpack.h b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_SlicedEllpack.h
similarity index 53%
rename from src/UnitTests/Matrices/SparseMatrixTest_SlicedEllpack.h
rename to src/UnitTests/Matrices/Legacy/SparseMatrixTest_SlicedEllpack.h
index 00184754c6c26632b3e2ffbccd3b244be58a4771..7f5ad546f6fd1f0a76551db46dfaa71928a4e5b6 100644
--- a/src/UnitTests/Matrices/SparseMatrixTest_SlicedEllpack.h
+++ b/src/UnitTests/Matrices/Legacy/SparseMatrixTest_SlicedEllpack.h
@@ -1,5 +1,5 @@
 /***************************************************************************
-                          SparseMatrixTest_SlicedEllpack_segments.h -  description
+                          SparseMatrixTest_SlicedEllpack.h -  description
                              -------------------
     begin                : Dec 9, 2019
     copyright            : (C) 2019 by Tomas Oberhuber et al.
@@ -8,8 +8,7 @@
 
 /* See Copyright Notice in tnl/Copyright */
 
-#include <TNL/Containers/Segments/SlicedEllpack.h>
-#include <TNL/Matrices/SparseMatrix.h>
+#include <TNL/Matrices/Legacy/SlicedEllpack.h>
 
 
 #include "SparseMatrixTest.hpp"
@@ -26,45 +25,38 @@ protected:
    using SlicedEllpackMatrixType = Matrix;
 };
 
-////
-// Row-major format is used for the host system
-template< typename Device, typename Index >
-using RowMajorSlicedEllpack = TNL::Containers::Segments::SlicedEllpack< Device, Index, true, 32 >;
+template< typename Real, typename Device, typename Index >
+using SlicedEllpackType = TNL::Matrices::SlicedEllpack< Real, Device, Index, 32 >;
 
 
-////
-// Column-major format is used for GPUs
-template< typename Device, typename Index >
-using ColumnMajorSlicedEllpack = TNL::Containers::Segments::SlicedEllpack< Device, Index, false, 32 >;
-
 // types for which MatrixTest is instantiated
 using SlicedEllpackMatrixTypes = ::testing::Types
 <
-    TNL::Matrices::SparseMatrix< int,     RowMajorSlicedEllpack, TNL::Devices::Host, short >,
-    TNL::Matrices::SparseMatrix< long,    RowMajorSlicedEllpack, TNL::Devices::Host, short >,
-    TNL::Matrices::SparseMatrix< float,   RowMajorSlicedEllpack, TNL::Devices::Host, short >,
-    TNL::Matrices::SparseMatrix< double,  RowMajorSlicedEllpack, TNL::Devices::Host, short >,
-    TNL::Matrices::SparseMatrix< int,     RowMajorSlicedEllpack, TNL::Devices::Host, int   >,
-    TNL::Matrices::SparseMatrix< long,    RowMajorSlicedEllpack, TNL::Devices::Host, int   >,
-    TNL::Matrices::SparseMatrix< float,   RowMajorSlicedEllpack, TNL::Devices::Host, int   >,
-    TNL::Matrices::SparseMatrix< double,  RowMajorSlicedEllpack, TNL::Devices::Host, int   >,
-    TNL::Matrices::SparseMatrix< int,     RowMajorSlicedEllpack, TNL::Devices::Host, long  >,
-    TNL::Matrices::SparseMatrix< long,    RowMajorSlicedEllpack, TNL::Devices::Host, long  >,
-    TNL::Matrices::SparseMatrix< float,   RowMajorSlicedEllpack, TNL::Devices::Host, long  >,
-    TNL::Matrices::SparseMatrix< double,  RowMajorSlicedEllpack, TNL::Devices::Host, long  >
+    SlicedEllpackType< int,     TNL::Devices::Host, short >,
+    SlicedEllpackType< long,    TNL::Devices::Host, short >,
+    SlicedEllpackType< float,   TNL::Devices::Host, short >,
+    SlicedEllpackType< double,  TNL::Devices::Host, short >,
+    SlicedEllpackType< int,     TNL::Devices::Host, int   >,
+    SlicedEllpackType< long,    TNL::Devices::Host, int   >,
+    SlicedEllpackType< float,   TNL::Devices::Host, int   >,
+    SlicedEllpackType< double,  TNL::Devices::Host, int   >,
+    SlicedEllpackType< int,     TNL::Devices::Host, long  >,
+    SlicedEllpackType< long,    TNL::Devices::Host, long  >,
+    SlicedEllpackType< float,   TNL::Devices::Host, long  >,
+    SlicedEllpackType< double,  TNL::Devices::Host, long  >
 #ifdef HAVE_CUDA
-   ,TNL::Matrices::SparseMatrix< int,     ColumnMajorSlicedEllpack, TNL::Devices::Cuda, short >,
-    TNL::Matrices::SparseMatrix< long,    ColumnMajorSlicedEllpack, TNL::Devices::Cuda, short >,
-    TNL::Matrices::SparseMatrix< float,   ColumnMajorSlicedEllpack, TNL::Devices::Cuda, short >,
-    TNL::Matrices::SparseMatrix< double,  ColumnMajorSlicedEllpack, TNL::Devices::Cuda, short >,
-    TNL::Matrices::SparseMatrix< int,     ColumnMajorSlicedEllpack, TNL::Devices::Cuda, int   >,
-    TNL::Matrices::SparseMatrix< long,    ColumnMajorSlicedEllpack, TNL::Devices::Cuda, int   >,
-    TNL::Matrices::SparseMatrix< float,   ColumnMajorSlicedEllpack, TNL::Devices::Cuda, int   >,
-    TNL::Matrices::SparseMatrix< double,  ColumnMajorSlicedEllpack, TNL::Devices::Cuda, int   >,
-    TNL::Matrices::SparseMatrix< int,     ColumnMajorSlicedEllpack, TNL::Devices::Cuda, long  >,
-    TNL::Matrices::SparseMatrix< long,    ColumnMajorSlicedEllpack, TNL::Devices::Cuda, long  >,
-    TNL::Matrices::SparseMatrix< float,   ColumnMajorSlicedEllpack, TNL::Devices::Cuda, long  >,
-    TNL::Matrices::SparseMatrix< double,  ColumnMajorSlicedEllpack, TNL::Devices::Cuda, long  >
+   ,SlicedEllpackType< int,     TNL::Devices::Cuda, short >,
+    SlicedEllpackType< long,    TNL::Devices::Cuda, short >,
+    SlicedEllpackType< float,   TNL::Devices::Cuda, short >,
+    SlicedEllpackType< double,  TNL::Devices::Cuda, short >,
+    SlicedEllpackType< int,     TNL::Devices::Cuda, int   >,
+    SlicedEllpackType< long,    TNL::Devices::Cuda, int   >,
+    SlicedEllpackType< float,   TNL::Devices::Cuda, int   >,
+    SlicedEllpackType< double,  TNL::Devices::Cuda, int   >,
+    SlicedEllpackType< int,     TNL::Devices::Cuda, long  >,
+    SlicedEllpackType< long,    TNL::Devices::Cuda, long  >,
+    SlicedEllpackType< float,   TNL::Devices::Cuda, long  >,
+    SlicedEllpackType< double,  TNL::Devices::Cuda, long  >
 #endif
 >;
 
@@ -149,4 +141,4 @@ TYPED_TEST( SlicedEllpackMatrixTest, printTest )
 
 #endif
 
-#include "../main.h"
+#include "../../main.h"
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_SlicedEllpack.cpp b/src/UnitTests/Matrices/SparseMatrixTest_SlicedEllpack.cpp
deleted file mode 100644
index a88301100dd94be8d859e545ea9a03cc98f4da73..0000000000000000000000000000000000000000
--- a/src/UnitTests/Matrices/SparseMatrixTest_SlicedEllpack.cpp
+++ /dev/null
@@ -1 +0,0 @@
-#include "SparseMatrixTest_SlicedEllpack_segments.h"
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_SlicedEllpack.cu b/src/UnitTests/Matrices/SparseMatrixTest_SlicedEllpack.cu
deleted file mode 100644
index a88301100dd94be8d859e545ea9a03cc98f4da73..0000000000000000000000000000000000000000
--- a/src/UnitTests/Matrices/SparseMatrixTest_SlicedEllpack.cu
+++ /dev/null
@@ -1 +0,0 @@
-#include "SparseMatrixTest_SlicedEllpack_segments.h"