diff --git a/src/TNL/Matrices/CSR.h b/src/TNL/Matrices/CSR.h
index ef7ba5d6f925d2e56e0df64b2951fe1752a7f84f..e4564562522765f1806846c4a619a049ec3f768b 100644
--- a/src/TNL/Matrices/CSR.h
+++ b/src/TNL/Matrices/CSR.h
@@ -13,6 +13,9 @@
 #include <TNL/Matrices/Sparse.h>
 #include <TNL/Containers/Vector.h>
 
+#include <TNL/Devices/Cuda.h>
+#include <TNL/Exceptions/CudaBadAlloc.h>
+
 namespace TNL {
 namespace Matrices {
    
@@ -41,17 +44,21 @@ private:
 
 public:
 
-   typedef Real RealType;
-   typedef Device DeviceType;
-   typedef Index IndexType;
+   using RealType = Real;
+   using DeviceType = Device;
+   using IndexType = Index;
    typedef typename Sparse< RealType, DeviceType, IndexType >:: CompressedRowLengthsVector CompressedRowLengthsVector;
    typedef typename Sparse< RealType, DeviceType, IndexType >::ConstCompressedRowLengthsVectorView ConstCompressedRowLengthsVectorView;
    typedef CSR< Real, Device, Index > ThisType;
    typedef CSR< Real, Devices::Host, Index > HostType;
    typedef CSR< Real, Devices::Cuda, Index > CudaType;
    typedef Sparse< Real, Device, Index > BaseType;
-   typedef typename BaseType::MatrixRow MatrixRow;
-   typedef SparseRow< const RealType, const IndexType > ConstMatrixRow;
+   //typedef typename BaseType::MatrixRow MatrixRow;
+   
+   using MatrixRow = typename BaseType::MatrixRow;
+   using ConstMatrixRow = typename BaseType::ConstMatrixRow;
+   //using typename BaseType::ConstMatrixRow;
+   //typedef SparseRow< const RealType, const IndexType > ConstMatrixRow;
 
 
    enum SPMVCudaKernel { scalar, vector, hybrid };
@@ -75,7 +82,11 @@ public:
 
    __cuda_callable__
    IndexType getRowLengthFast( const IndexType row ) const;
-
+   
+   IndexType getNonZeroRowLength( const IndexType row ) const;
+   
+   IndexType getNonZeroRowLengthFast( const IndexType row ) const;
+   
    template< typename Real2, typename Device2, typename Index2 >
    void setLike( const CSR< Real2, Device2, Index2 >& matrix );
 
diff --git a/src/TNL/Matrices/CSR_impl.h b/src/TNL/Matrices/CSR_impl.h
index b4dff85470bf86021c69478138eb3be86f74d593..0a682a9dc25cbea51617bb8c5ab89b0a76b0d846 100644
--- a/src/TNL/Matrices/CSR_impl.h
+++ b/src/TNL/Matrices/CSR_impl.h
@@ -131,6 +131,60 @@ Index CSR< Real, Device, Index >::getRowLengthFast( const IndexType row ) const
    return this->rowPointers[ row + 1 ] - this->rowPointers[ row ];
 }
 
+template< typename Real,
+          typename Device,
+          typename Index >
+Index CSR< Real, Device, Index >::getNonZeroRowLength( const IndexType row ) const
+{
+    // TODO: Fix/Implement
+    TNL_ASSERT( false, std::cerr << "TODO: Fix/Implement" );
+    return 0;
+//    if( std::is_same< DeviceType, Devices::Host >::value )
+//    {
+//       ConstMatrixRow matrixRow = this->getRow( row );
+//       return matrixRow.getNonZeroElementsCount();
+//    }
+//    if( std::is_same< DeviceType, Devices::Cuda >::value )
+//    {
+//       IndexType *cols = new IndexType[4];
+//       RealType *vals = new RealType[4];
+//       for( int i = 0; i < 4; i++ )
+//       {
+//           cols[i] = i;
+//           vals[i] = 1.0;
+//       }
+//       ConstMatrixRow matrixRow(cols, vals, 4, 1);
+// //      ConstMatrixRow matrixRow = this->getRow( row );// If the program even compiles, this line fails because a segfault is thrown on the first line of getRow()
+//       // WHEN debugging with GDB:
+//       //  (gdb) p this->rowPointers[0]
+//       //    Could not find operator[].
+//       //  (gdb) p rowPointers.getElement(0)
+//       //    Attempt to take address of value not located in memory.
+//       IndexType resultHost ( 0 );
+//       IndexType* resultCuda = Devices::Cuda::passToDevice( resultHost );
+//       // PROBLEM: If the second parameter of getNonZeroRowLengthCudaKernel is '&resultCuda', the following issue is thrown:
+//       //          'error: no instance of function template "TNL::Matrices::getNonZeroRowLengthCudaKernel" matches the argument list'
+//       TNL::Matrices::getNonZeroRowLengthCudaKernel< ConstMatrixRow, IndexType ><<< 1, 1 >>>( matrixRow, resultCuda ); // matrixRow works fine, tested them both separately
+//       delete []cols;
+//       delete []vals;
+//       std::cout << "Checkpoint BEFORE passFromDevice" << std::endl;
+//       resultHost = Devices::Cuda::passFromDevice( resultCuda ); // This causes a crash: Illegal memory address in Cuda_impl.h at TNL_CHECK_CUDA_DEVICE
+//       std::cout << "Checkpoint AFTER passFromDevice" << std::endl;
+//       Devices::Cuda::freeFromDevice( resultCuda );
+//       return resultHost;
+//   }
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+__cuda_callable__
+Index CSR< Real, Device, Index >::getNonZeroRowLengthFast( const IndexType row ) const
+{  
+   ConstMatrixRow matrixRow = this->getRow( row );
+   return matrixRow.getNonZeroElementsCount();
+}
+
 template< typename Real,
           typename Device,
           typename Index >
@@ -430,12 +484,12 @@ typename CSR< Real, Device, Index >::ConstMatrixRow
 CSR< Real, Device, Index >::
 getRow( const IndexType rowIndex ) const
 {
-   const IndexType rowOffset = this->rowPointers[ rowIndex ];
-   const IndexType rowLength = this->rowPointers[ rowIndex + 1 ] - rowOffset;
-   return ConstMatrixRow( &this->columnIndexes[ rowOffset ],
-                          &this->values[ rowOffset ],
-                          rowLength,
-                          1 );
+    const IndexType rowOffset = this->rowPointers[ rowIndex ];
+    const IndexType rowLength = this->rowPointers[ rowIndex + 1 ] - rowOffset;
+    return ConstMatrixRow( &this->columnIndexes[ rowOffset ],
+                           &this->values[ rowOffset ],
+                           rowLength,
+                           1 );
 }
 
 template< typename Real,
diff --git a/src/TNL/Matrices/ChunkedEllpack.h b/src/TNL/Matrices/ChunkedEllpack.h
index 35bbfa89799eff2b248283cda4ef141bcf7eb039..ff889a49fe0d7b71c4e1fb1e6e1165e7d8970b86 100644
--- a/src/TNL/Matrices/ChunkedEllpack.h
+++ b/src/TNL/Matrices/ChunkedEllpack.h
@@ -104,6 +104,8 @@ public:
 
    __cuda_callable__
    IndexType getRowLengthFast( const IndexType row ) const;
+   
+   IndexType getNonZeroRowLength( const IndexType row ) const;
 
    template< typename Real2, typename Device2, typename Index2 >
    void setLike( const ChunkedEllpack< Real2, Device2, Index2 >& matrix );
diff --git a/src/TNL/Matrices/ChunkedEllpack_impl.h b/src/TNL/Matrices/ChunkedEllpack_impl.h
index 20dbfa68349b7f23ad66e36c300d1894ea7e40be..ff1dd0742f1f2010b441206395e5cf278deaf37e 100644
--- a/src/TNL/Matrices/ChunkedEllpack_impl.h
+++ b/src/TNL/Matrices/ChunkedEllpack_impl.h
@@ -179,10 +179,28 @@ bool ChunkedEllpack< Real, Device, Index >::setSlice( ConstCompressedRowLengthsV
     */
    IndexType maxChunkInSlice( 0 );
    for( IndexType i = sliceBegin; i < sliceEnd; i++ )
-      maxChunkInSlice = max( maxChunkInSlice,
-                          ceil( ( RealType ) rowLengths[ i ] /
-                                ( RealType ) this->rowToChunkMapping[ i ] ) );
-   TNL_ASSERT( maxChunkInSlice > 0,
+   {       
+//       ALL OF THE FOLLOWING std::couts are for troubleshooting purposes, can be deleted.
+//       std::cout << "Troubleshooting invalid ceil operation: " << std::endl;
+//       std::cout << "maxChunkInSlice = " << maxChunkInSlice << std::endl;
+//       std::cout << "( RealType ) rowLengths[ i ] = " <<  ( RealType ) rowLengths[ i ] << std::endl;
+//       std::cout << "( RealType ) this->rowToChunkMapping[ i ] = " <<  ( RealType ) this->rowToChunkMapping[ i ] << std::endl;
+//       std::cout << " ceil( RealType / RealType ) = " << ceil( ( RealType ) rowLengths[ i ] / ( RealType ) this->rowToChunkMapping[ i ] ) << std::endl;
+//       std::cout << "( int ) rowLengths[ i ] = " <<  ( int ) rowLengths[ i ] << std::endl;
+//       std::cout << "( int ) this->rowToChunkMapping[ i ] = " <<  ( int ) this->rowToChunkMapping[ i ] << std::endl;
+//       std::cout << " ceil( int / int ) = " << ceil( ( int ) rowLengths[ i ] / ( int ) this->rowToChunkMapping[ i ] ) << std::endl;
+//       std::cout << "( float ) rowLengths[ i ] = " <<  ( float ) rowLengths[ i ] << std::endl;
+//       std::cout << "( float ) this->rowToChunkMapping[ i ] = " <<  ( float ) this->rowToChunkMapping[ i ] << std::endl;
+//       std::cout << " ceil( float / float ) = " << ceil( ( float ) rowLengths[ i ] / ( float ) this->rowToChunkMapping[ i ] ) << std::endl;
+//       The ceil function doesn't work when rowLengths and the other this.->... is 
+//       typecasted into ( RealType ), because when RealType is int, it will perform 
+//       an integer division and return the int as a double, which in this case 
+//       will be zero and make the assertion fail ( https://stackoverflow.com/questions/33273359/in-c-using-the-ceil-a-division-is-not-working ).
+//       To fix this, typecast them to ( float ), instead of ( RealType )
+       maxChunkInSlice = max( maxChunkInSlice,
+                          roundUpDivision( rowLengths[ i ], this->rowToChunkMapping[ i ] ) );
+   }
+      TNL_ASSERT( maxChunkInSlice > 0,
               std::cerr << " maxChunkInSlice = " << maxChunkInSlice << std::endl );
 
    /****
@@ -231,6 +249,13 @@ void ChunkedEllpack< Real, Device, Index >::setCompressedRowLengths( ConstCompre
          this->setSlice( rowLengths, sliceIndex, elementsToAllocation );
       this->rowPointers.computePrefixSum();
    }
+   
+//   std::cout << "\ngetRowLength after first if: " << std::endl;
+//   for( IndexType i = 0; i < rowLengths.getSize(); i++ )
+//   {
+//       std::cout << getRowLength( i ) << std::endl;
+//   }
+//   std::cout << "\n";
 
    if( std::is_same< Device, Devices::Cuda >::value )
    {
@@ -255,6 +280,7 @@ void ChunkedEllpack< Real, Device, Index >::setCompressedRowLengths( ConstCompre
       elementsToAllocation = hostMatrix.values.getSize();
    }
    this->maxRowLength = rowLengths.max();
+//   std::cout << "\nrowLengths.max() = " << rowLengths.max() << std::endl;
    Sparse< Real, Device, Index >::allocateMatrixElements( elementsToAllocation );
 }
 
@@ -281,6 +307,30 @@ Index ChunkedEllpack< Real, Device, Index >::getRowLengthFast( const IndexType r
    return rowPointers[ row + 1 ] - rowPointers[ row ];
 }
 
+template< typename Real,
+          typename Device,
+          typename Index >
+Index ChunkedEllpack< Real, Device, Index >::getNonZeroRowLength( const IndexType row ) const
+{
+    ConstMatrixRow matrixRow = getRow( row );
+    return matrixRow.getNonZeroElementsCount( Device::getDeviceType() );
+    
+//    IndexType elementCount ( 0 );
+//    ConstMatrixRow matrixRow = this->getRow( row );
+//    
+//    auto computeNonZeros = [&] /*__cuda_callable__*/ ( IndexType i ) mutable
+//    {
+//        std::cout << "matrixRow.getElementValue( i ) = " << matrixRow.getElementValue( i ) << " != 0.0" << std::endl;
+//        if( matrixRow.getElementValue( i ) !=  0.0 )
+//            elementCount++;
+//        
+//        std::cout << "End of lambda elementCount = " << elementCount << std::endl;
+//    };
+//   
+//    ParallelFor< DeviceType >::exec( ( IndexType ) 0, matrixRow.getLength(), computeNonZeros );
+//    return elementCount;
+}
+
 template< typename Real,
           typename Device,
           typename Index >
@@ -952,10 +1002,10 @@ getRow( const IndexType rowIndex ) const
 {
    const IndexType rowOffset = this->rowPointers[ rowIndex ];
    const IndexType rowLength = this->rowPointers[ rowIndex + 1 ] - rowOffset;
-   return MatrixRow( &this->columnIndexes[ rowOffset ],
-                     &this->values[ rowOffset ],
-                     rowLength,
-                     1 );
+   return ConstMatrixRow( &this->columnIndexes[ rowOffset ],
+                          &this->values[ rowOffset ],
+                          rowLength,
+                          1 );
 }
 
 
diff --git a/src/TNL/Matrices/Ellpack.h b/src/TNL/Matrices/Ellpack.h
index 1646db1c5c8b37bb635af0ee3501afd3fce6e431..7d17ff07e8649b7c3db872f08672cf0f36fb92ce 100644
--- a/src/TNL/Matrices/Ellpack.h
+++ b/src/TNL/Matrices/Ellpack.h
@@ -67,6 +67,8 @@ public:
 
    __cuda_callable__
    IndexType getRowLengthFast( const IndexType row ) const;
+   
+   IndexType getNonZeroRowLength( const IndexType row ) const;
 
    template< typename Real2, typename Device2, typename Index2 >
    void setLike( const Ellpack< Real2, Device2, Index2 >& matrix );
diff --git a/src/TNL/Matrices/Ellpack_impl.h b/src/TNL/Matrices/Ellpack_impl.h
index 6186206439474d97d8edced12b4671b257b6f0ed..f3c05c492d8e269834b523436e79fd25edab2d0f 100644
--- a/src/TNL/Matrices/Ellpack_impl.h
+++ b/src/TNL/Matrices/Ellpack_impl.h
@@ -123,6 +123,15 @@ Index Ellpack< Real, Device, Index >::getRowLengthFast( const IndexType row ) co
    return this->rowLengths;
 }
 
+template< typename Real,
+          typename Device,
+          typename Index >
+Index Ellpack< Real, Device, Index >::getNonZeroRowLength( const IndexType row ) const
+{
+    ConstMatrixRow matrixRow = getRow( row );
+    return matrixRow.getNonZeroElementsCount( Device::getDeviceType() );
+}
+
 template< typename Real,
           typename Device,
           typename Index >
diff --git a/src/TNL/Matrices/SlicedEllpack.h b/src/TNL/Matrices/SlicedEllpack.h
index 6f68f2fa8aea4979b8f4685d2ee25d3039653ea7..0fc9ccb0b2ca42bd2ec6441f9fe76cd661679f8e 100644
--- a/src/TNL/Matrices/SlicedEllpack.h
+++ b/src/TNL/Matrices/SlicedEllpack.h
@@ -95,6 +95,8 @@ public:
 
    __cuda_callable__
    IndexType getRowLengthFast( const IndexType row ) const;
+   
+   IndexType getNonZeroRowLength( const IndexType row ) const;
 
    template< typename Real2, typename Device2, typename Index2 >
    void setLike( const SlicedEllpack< Real2, Device2, Index2, SliceSize >& matrix );
diff --git a/src/TNL/Matrices/SlicedEllpack_impl.h b/src/TNL/Matrices/SlicedEllpack_impl.h
index 95a601a00a01ead3f11d0cd0ca0f96a0373b9606..4d7593d3f33a1e1140f87fd8eb876172e70e42a5 100644
--- a/src/TNL/Matrices/SlicedEllpack_impl.h
+++ b/src/TNL/Matrices/SlicedEllpack_impl.h
@@ -121,6 +121,16 @@ Index SlicedEllpack< Real, Device, Index, SliceSize >::getRowLengthFast( const I
    return this->sliceCompressedRowLengths[ slice ];
 }
 
+template< typename Real,
+          typename Device,
+          typename Index ,
+          int SliceSize >
+Index SlicedEllpack< Real, Device, Index, SliceSize >::getNonZeroRowLength( const IndexType row ) const
+{
+    ConstMatrixRow matrixRow = getRow( row );
+    return matrixRow.getNonZeroElementsCount( Device::getDeviceType() );
+}
+
 template< typename Real,
           typename Device,
           typename Index,
@@ -212,7 +222,7 @@ bool SlicedEllpack< Real, Device, Index, SliceSize >::addElementFast( const Inde
                                                                                const RealType& thisElementMultiplicator )
 {
    TNL_ASSERT( row >= 0 && row < this->rows &&
-              column >= 0 && column <= this->rows,
+              column >= 0 && column <= this->columns,
               std::cerr << " row = " << row
                    << " column = " << column
                    << " this->rows = " << this->rows
@@ -260,7 +270,7 @@ bool SlicedEllpack< Real, Device, Index, SliceSize >::addElement( const IndexTyp
                                                                            const RealType& thisElementMultiplicator )
 {
    TNL_ASSERT( row >= 0 && row < this->rows &&
-              column >= 0 && column <= this->rows,
+              column >= 0 && column <= this->columns,
               std::cerr << " row = " << row
                    << " column = " << column
                    << " this->rows = " << this->rows
diff --git a/src/TNL/Matrices/Sparse.h b/src/TNL/Matrices/Sparse.h
index 2ee49219ee2d8fa8be6662c417dd68c3c3a6c690..069ade36cb989e18e1b6e1cc9821af5df50de8c1 100644
--- a/src/TNL/Matrices/Sparse.h
+++ b/src/TNL/Matrices/Sparse.h
@@ -30,6 +30,7 @@ class Sparse : public Matrix< Real, Device, Index >
    typedef Containers::Vector< IndexType, DeviceType, IndexType > ColumnIndexesVector;
    typedef Matrix< Real, Device, Index > BaseType;
    typedef SparseRow< RealType, IndexType > MatrixRow;
+   typedef SparseRow< const RealType, const IndexType > ConstMatrixRow;
 
    Sparse();
 
diff --git a/src/TNL/Matrices/SparseRow.h b/src/TNL/Matrices/SparseRow.h
index e7547ee679e89cef5e2c607720d13d538ba3c31a..fac855eae71a26cdb4dbf62f927f12d0f24b5af1 100644
--- a/src/TNL/Matrices/SparseRow.h
+++ b/src/TNL/Matrices/SparseRow.h
@@ -21,6 +21,9 @@ namespace Matrices {
 template< typename Real, typename Index >
 class SparseRow
 {
+   using RealType = Real;
+   using IndexType = Index;
+   
    public:
 
       __cuda_callable__
@@ -51,6 +54,9 @@ class SparseRow
 
       __cuda_callable__
       Index getLength() const;
+      
+      __cuda_callable__
+      Index getNonZeroElementsCount() const;
 
       void print( std::ostream& str ) const;
 
diff --git a/src/TNL/Matrices/SparseRow_impl.h b/src/TNL/Matrices/SparseRow_impl.h
index f6921b15bbf8f282f0dee1a2f6c842a230cf0dd3..6c86b9d5176e6b9a52a8a4177907e134b39a3394 100644
--- a/src/TNL/Matrices/SparseRow_impl.h
+++ b/src/TNL/Matrices/SparseRow_impl.h
@@ -11,6 +11,11 @@
 #pragma once
 
 #include <TNL/Matrices/SparseRow.h>
+#include <TNL/ParallelFor.h>
+
+// Following includes are here to enable usage of std::vector and std::cout. To avoid having to include Device type (HOW would this be done anyway)
+#include <iostream>
+#include <vector>
 
 namespace TNL {
 namespace Matrices {   
@@ -107,6 +112,46 @@ getLength() const
    return length;
 }
 
+#ifdef HAVE_CUDA
+template< typename MatrixRow, typename Index >
+__global__
+void getNonZeroRowLengthCudaKernel( const MatrixRow row, Index* result )
+{
+//    TODO: Fix/Implement
+    TNL_ASSERT( false, std::cerr << "TODO: Fix/Implement" );
+//    int threadId = blockIdx.x * blockDim.x + threadIdx.x;
+//    if( threadId == 0 )
+//    {
+//       *result = row.getNonZeroElementsCount();
+//    }
+}
+#endif
+
+template< typename Real, typename Index >
+__cuda_callable__
+Index
+SparseRow< Real, Index >::
+getNonZeroElementsCount() const
+{
+//    TODO: Fix/Implement
+    TNL_ASSERT( false, std::cerr << "TODO: Fix/Implement" );
+    return 0;
+//    using NonConstIndex = typename std::remove_const< Index >::type;
+//    
+//    NonConstIndex elementCount ( 0 );
+//   
+//    for( NonConstIndex i = 0; i < length; i++ )
+//    {
+////        std::cout << "this->values[ i * step ] = " << this->values[ i * step ] << " != 0.0" << std::endl;
+//        if( this->values[ i * step ] != 0.0 ) // Returns the same amount of elements in a row as does getRowLength() in ChunkedEllpack. WHY?
+//            elementCount++;
+//    }
+//    
+////    std::cout << "Element Count = " << elementCount << "\n";
+//    
+//    return elementCount;
+}
+
 template< typename Real, typename Index >
 void
 SparseRow< Real, Index >::
diff --git a/src/UnitTests/Matrices/CMakeLists.txt b/src/UnitTests/Matrices/CMakeLists.txt
index 8fbf4e4c8ba5d36e63b5d80819853dccfd7b0a68..adb189ac62164b1ea46874d2732d3e1ed0b66b16 100644
--- a/src/UnitTests/Matrices/CMakeLists.txt
+++ b/src/UnitTests/Matrices/CMakeLists.txt
@@ -2,12 +2,32 @@ IF( BUILD_CUDA )
    CUDA_ADD_EXECUTABLE( SparseMatrixCopyTest SparseMatrixCopyTest.h SparseMatrixCopyTest.cu OPTIONS ${CXX_TESTS_FLAGS} )
    TARGET_LINK_LIBRARIES( SparseMatrixCopyTest ${GTEST_BOTH_LIBRARIES}
                                                            tnl )
+
+   CUDA_ADD_EXECUTABLE( SparseMatrixTest SparseMatrixTest.h SparseMatrixTest.cu OPTIONS ${CXX_TESTS_FLAGS} )
+   TARGET_LINK_LIBRARIES( SparseMatrixTest ${GTEST_BOTH_LIBRARIES}
+                                                           tnl )
+
+   CUDA_ADD_EXECUTABLE( DenseMatrixTest DenseMatrixTest.h DenseMatrixTest.cu OPTIONS ${CXX_TESTS_FLAGS} )
+   TARGET_LINK_LIBRARIES( DenseMatrixTest ${GTEST_BOTH_LIBRARIES}
+                                                           tnl )
 ELSE(  BUILD_CUDA )
    ADD_EXECUTABLE( SparseMatrixCopyTest SparseMatrixCopyTest.h SparseMatrixCopyTest.cpp )
    TARGET_COMPILE_OPTIONS( SparseMatrixCopyTest PRIVATE ${CXX_TESTS_FLAGS} )
    TARGET_LINK_LIBRARIES( SparseMatrixCopyTest ${GTEST_BOTH_LIBRARIES}
                                                            tnl )
+
+   ADD_EXECUTABLE( SparseMatrixTest SparseMatrixTest.h SparseMatrixTest.cpp )
+   TARGET_COMPILE_OPTIONS( SparseMatrixTest PRIVATE ${CXX_TESTS_FLAGS} )
+   TARGET_LINK_LIBRARIES( SparseMatrixTest ${GTEST_BOTH_LIBRARIES}
+                                                           tnl )
+
+   ADD_EXECUTABLE( DenseMatrixTest DenseMatrixTest.h DenseMatrixTest.cpp )
+   TARGET_COMPILE_OPTIONS( DenseMatrixTest PRIVATE ${CXX_TESTS_FLAGS} )
+   TARGET_LINK_LIBRARIES( DenseMatrixTest ${GTEST_BOTH_LIBRARIES}
+                                                           tnl )
 ENDIF( BUILD_CUDA )
 
 
 ADD_TEST( SparseMatrixCopyTest ${EXECUTABLE_OUTPUT_PATH}/SparseMatrixCopyTest${CMAKE_EXECUTABLE_SUFFIX} )
+ADD_TEST( SparseMatrixTest ${EXECUTABLE_OUTPUT_PATH}/SparseMatrixTest${CMAKE_EXECUTABLE_SUFFIX} )
+ADD_TEST( DenseMatrixTest ${EXECUTABLE_OUTPUT_PATH}/SparseMatrixTest${CMAKE_EXECUTABLE_SUFFIX} )
diff --git a/src/UnitTests/Matrices/DenseMatrixTest.cpp b/src/UnitTests/Matrices/DenseMatrixTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a56349360a8ba6667dfb567c64a6f4fa0d7ff48f
--- /dev/null
+++ b/src/UnitTests/Matrices/DenseMatrixTest.cpp
@@ -0,0 +1,11 @@
+/***************************************************************************
+                          DenseMatrixTest.cpp -  description
+                             -------------------
+    begin                : Nov 10, 2018
+    copyright            : (C) 2018 by Tomas Oberhuber et al.
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#include "DenseMatrixTest.h"
\ No newline at end of file
diff --git a/src/UnitTests/Matrices/DenseMatrixTest.cu b/src/UnitTests/Matrices/DenseMatrixTest.cu
new file mode 100644
index 0000000000000000000000000000000000000000..11d45efdb5f30ecc236e8086d9149f8f35b8c0d3
--- /dev/null
+++ b/src/UnitTests/Matrices/DenseMatrixTest.cu
@@ -0,0 +1,11 @@
+/***************************************************************************
+                          DenseMatrixTest.cu -  description
+                             -------------------
+    begin                : Nov 10, 2018
+    copyright            : (C) 2018 by Tomas Oberhuber et al.
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+#include "DenseMatrixTest.h"
\ No newline at end of file
diff --git a/src/UnitTests/Matrices/DenseMatrixTest.h b/src/UnitTests/Matrices/DenseMatrixTest.h
new file mode 100644
index 0000000000000000000000000000000000000000..c67386e8199fb486f7625279902a45736d67ebdd
--- /dev/null
+++ b/src/UnitTests/Matrices/DenseMatrixTest.h
@@ -0,0 +1,1664 @@
+/***************************************************************************
+                          DenseMatrixTest.h -  description
+                             -------------------
+    begin                : Nov 10, 2018
+    copyright            : (C) 2018 by Tomas Oberhuber et al.
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/* See Copyright Notice in tnl/Copyright */
+
+// TODO
+/*
+  * getType()                           ::HOW?  How to test this for each format? edit string how?
+ *      MISTAKE! found it for Cuda instead of Devices::Cuda. Incorrect String in src/TNL/Devices/Cuda.cpp
+ * getTypeVirtual()                     ::TEST? This just calls getType().
+ * getSerializationType()               ::TEST? This just calls getType().
+ * getSerializationTypeVirtual()        ::TEST? This just calls getSerializationType().
+ * setDimensions()                          ::DONE
+ * setLike()                                ::DONE
+ * setCompressedRowLengths()            ::NOT IMPLEMENTED! The function body is empty.
+ * getRowLength()                           ::DONE
+ * getRowLengthFast()                   ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * getMaxRowLength()                    ::TEST? This function is identical to getRowLength().
+ * getNumberOfMatrixElements()              ::DONE
+ * getNumberOfNonZeroMatrixElements()       ::DONE
+ * reset()                                  ::DONE
+ * setValue()                               ::DONE
+ * operator()                           ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * const operator()                     ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * setElementFast()                     ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * setElement()                             ::DONE ; USED! in any test with individual value assignment.
+ * addElementFast()                     ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * addElement()                             ::DONE
+ * setRowFast()                         ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * setRow()                                 ::DONE
+ *      MISTAKE! This function unlike the setRow() for CSR, doesn't replace all the elements of a row, it only replaces the elements it has values for in its arrays.
+ * addRowFast()                         ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * addRow()                                 ::DONE
+ * getElementFast()                     ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * getElement()                             ::USED! in any test with individual value reading.
+ * getRowFast()                         ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * getRow()                             ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * const getRow()                       ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * MatrixRow getRow()                   ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * ConstMatrixRow getRow()              ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * rowVectorProduct()                   ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * vectorProduct()                          ::DONE
+ *      This used to throw illegal memory access, but instead of using ints for vectors, using Types, helped.
+ * addMatrix()                              ::DONE
+ * DenseMatrixProductKernel()           ::HOW? How to test __global__?
+ * getMatrixProdut()                    ::HOW? It won't build: When testing CPU: no parameters match function DenseMatrixProductKernel(); when testing GPU: identifier tnlCudaMin is undefined. 
+ * DenseTranspositionAlignedKernel()    ::HOW? How to test __global__?
+ * DenseTranspositionNonAlignedKernel() ::HOW? How to test __global__?
+ * getTransposition()                   ::HOW? It won't build when testing CPU: no parameters match functions DenseTranspositionAlignedKernel() and DenseTranspositionNonAlignedKernel(). On GPU if will throw terminate and (core dumped).
+ *      MISTAKE! For GPU it works completely fine, when rows == cols. Otherwise it throws assertion failed.
+ * performSORIteration()                ::HOW? Throws segmentation fault CUDA.
+ * operator=()                          ::HOW? What is this supposed to enable? Overloading operators?
+ * save( String& fileName )                 ::DONE
+ * load( String& fileName )                 ::DONE
+ * save( File& file)                    ::USED! In save( String& fileName )
+ * load( File& file )                   ::USED! In load( String& fileName )
+ * print()                                  ::DONE
+ * getElementIndex()                    ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ */
+
+// GENERAL TODO
+/*
+ * Template tests for all formats.
+ * Figure out __cuda_callable_. When trying to call __cuda_callable__ functions
+ *          a segmentation fault (core dumped) is thrown.
+ *      ==>__cuda_callable__ works only for CPU at the moment. (for loops vs thread kernel assignment)
+ */
+
+
+#include <TNL/Devices/Host.h>
+#include <TNL/Matrices/Matrix.h>
+#include <TNL/Matrices/Dense.h>
+#include <TNL/Containers/Array.h>
+
+#include <TNL/Containers/Vector.h>
+#include <TNL/Containers/VectorView.h>
+#include <TNL/Math.h>
+#include <iostream>
+
+using Dense_host_float = TNL::Matrices::Dense< float, TNL::Devices::Host, int >;
+using Dense_host_int = TNL::Matrices::Dense< int, TNL::Devices::Host, int >;
+
+using Dense_cuda_float = TNL::Matrices::Dense< float, TNL::Devices::Cuda, int >;
+using Dense_cuda_int = TNL::Matrices::Dense< int, TNL::Devices::Cuda, int >;
+
+#ifdef HAVE_GTEST 
+#include <type_traits>
+
+#include <gtest/gtest.h>
+
+template< typename MatrixHostFloat, typename MatrixHostInt >
+void host_test_GetType()
+{
+    MatrixHostFloat mtrxHostFloat;
+    MatrixHostInt mtrxHostInt;
+    
+    EXPECT_EQ( mtrxHostFloat.getType(), TNL::String( "Matrices::Dense< float, Devices::Host, int >" ) );
+    EXPECT_EQ( mtrxHostInt.getType(), TNL::String( "Matrices::Dense< int, Devices::Host, int >" ) );
+}
+
+// QUESITON: Cant these two functions be combined into one? Because if no CUDA is present and we were to call
+//           CUDA into the function in the TEST, to be tested, then we could have a problem.
+
+template< typename MatrixCudaFloat, typename MatrixCudaInt >
+void cuda_test_GetType()
+{
+    MatrixCudaFloat mtrxCudaFloat;
+    MatrixCudaInt mtrxCudaInt;
+
+    EXPECT_EQ( mtrxCudaFloat.getType(), TNL::String( "Matrices::Dense< float, Cuda, int >" ) );    // This is mistakenly labeled in /src/TNL/Devices/Cuda.cpp
+    EXPECT_EQ( mtrxCudaInt.getType(), TNL::String( "Matrices::Dense< int, Cuda, int >" ) );        // Should be Devices::Cuda
+}
+
+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 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_GetRowLength()
+{
+    using RealType = typename Matrix::RealType;
+    using DeviceType = typename Matrix::DeviceType;
+    using IndexType = typename Matrix::IndexType;
+    
+    const IndexType rows = 8;
+    const IndexType cols = 7;
+    
+    Matrix m;
+    m.reset();
+    m.setDimensions( rows, cols );
+    
+    EXPECT_EQ( m.getRowLength( 0 ), 7 );
+    EXPECT_EQ( m.getRowLength( 1 ), 7 );
+    EXPECT_EQ( m.getRowLength( 2 ), 7 );
+    EXPECT_EQ( m.getRowLength( 3 ), 7 );
+    EXPECT_EQ( m.getRowLength( 4 ), 7 );
+    EXPECT_EQ( m.getRowLength( 5 ), 7 );
+    EXPECT_EQ( m.getRowLength( 6 ), 7 );
+    EXPECT_EQ( m.getRowLength( 7 ), 7 );
+}
+
+template< typename Matrix >
+void test_GetNumberOfMatrixElements()
+{
+    using RealType = typename Matrix::RealType;
+    using DeviceType = typename Matrix::DeviceType;
+    using IndexType = typename Matrix::IndexType;
+    
+    const IndexType rows = 7;
+    const IndexType cols = 6;
+    
+    Matrix m;
+    m.reset();
+    m.setDimensions( rows, cols );
+    
+    EXPECT_EQ( m.getNumberOfMatrixElements(), 42 );
+}
+
+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 7x6 dense matrix:
+ *
+ *    /  0  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 29 30 |
+ *    | 31 32 33 34 35 36 |
+ *    \ 37 38 39 40 41  0 /
+ */
+    const IndexType rows = 7;
+    const IndexType cols = 6;
+    
+    Matrix m;
+    m.reset();
+    m.setDimensions( rows, cols );
+    
+    RealType value = 1;
+    for( IndexType i = 0; i < rows; i++ )
+        for( IndexType j = 0; j < cols; j++ )
+            m.setElement( i, j, value++ );
+    
+    m.setElement( 0, 0, 0); // Set the first element of the diagonal to 0.
+    m.setElement( 6, 5, 0); // Set the last element of the diagonal to 0.
+    
+    EXPECT_EQ( m.getNumberOfNonzeroMatrixElements(), 40 );
+}
+
+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 dense 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_SetValue()
+{
+    using RealType = typename Matrix::RealType;
+    using DeviceType = typename Matrix::DeviceType;
+    using IndexType = typename Matrix::IndexType;  
+/*
+ * Sets up the following 7x6 dense 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 29 30 |
+ *    | 31 32 33 34 35 36 |
+ *    \ 37 38 39 40 41 42 /
+ */
+    const IndexType rows = 7;
+    const IndexType cols = 6;
+    
+    Matrix m;
+    m.reset();
+    m.setDimensions( rows, cols );
+    
+    RealType value = 1;
+    for( IndexType i = 0; i < rows; i++ )
+        for( IndexType j = 0; j < cols; j++ )
+            m.setElement( i, j, value++ );
+    
+    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 ),  4 );
+    EXPECT_EQ( m.getElement( 0, 4 ),  5 );
+    EXPECT_EQ( m.getElement( 0, 5 ),  6 );
+    
+    EXPECT_EQ( m.getElement( 1, 0 ),  7 );
+    EXPECT_EQ( m.getElement( 1, 1 ),  8 );
+    EXPECT_EQ( m.getElement( 1, 2 ),  9 );
+    EXPECT_EQ( m.getElement( 1, 3 ), 10 );
+    EXPECT_EQ( m.getElement( 1, 4 ), 11 );
+    EXPECT_EQ( m.getElement( 1, 5 ), 12 );
+    
+    EXPECT_EQ( m.getElement( 2, 0 ), 13 );
+    EXPECT_EQ( m.getElement( 2, 1 ), 14 );
+    EXPECT_EQ( m.getElement( 2, 2 ), 15 );
+    EXPECT_EQ( m.getElement( 2, 3 ), 16 );
+    EXPECT_EQ( m.getElement( 2, 4 ), 17 );
+    EXPECT_EQ( m.getElement( 2, 5 ), 18 );
+    
+    EXPECT_EQ( m.getElement( 3, 0 ), 19 );
+    EXPECT_EQ( m.getElement( 3, 1 ), 20 );
+    EXPECT_EQ( m.getElement( 3, 2 ), 21 );
+    EXPECT_EQ( m.getElement( 3, 3 ), 22 );
+    EXPECT_EQ( m.getElement( 3, 4 ), 23 );
+    EXPECT_EQ( m.getElement( 3, 5 ), 24 );
+    
+    EXPECT_EQ( m.getElement( 4, 0 ), 25 );
+    EXPECT_EQ( m.getElement( 4, 1 ), 26 );
+    EXPECT_EQ( m.getElement( 4, 2 ), 27 );
+    EXPECT_EQ( m.getElement( 4, 3 ), 28 );
+    EXPECT_EQ( m.getElement( 4, 4 ), 29 );
+    EXPECT_EQ( m.getElement( 4, 5 ), 30 );
+    
+    EXPECT_EQ( m.getElement( 5, 0 ), 31 );
+    EXPECT_EQ( m.getElement( 5, 1 ), 32 );
+    EXPECT_EQ( m.getElement( 5, 2 ), 33 );
+    EXPECT_EQ( m.getElement( 5, 3 ), 34 );
+    EXPECT_EQ( m.getElement( 5, 4 ), 35 );
+    EXPECT_EQ( m.getElement( 5, 5 ), 36 );
+    
+    EXPECT_EQ( m.getElement( 6, 0 ), 37 );
+    EXPECT_EQ( m.getElement( 6, 1 ), 38 );
+    EXPECT_EQ( m.getElement( 6, 2 ), 39 );
+    EXPECT_EQ( m.getElement( 6, 3 ), 40 );
+    EXPECT_EQ( m.getElement( 6, 4 ), 41 );
+    EXPECT_EQ( m.getElement( 6, 5 ), 42 );
+    
+    // Set the values of all elements to a certain number
+    m.setValue( 42 );
+    
+    EXPECT_EQ( m.getElement( 0, 0 ), 42 );
+    EXPECT_EQ( m.getElement( 0, 1 ), 42 );
+    EXPECT_EQ( m.getElement( 0, 2 ), 42 );
+    EXPECT_EQ( m.getElement( 0, 3 ), 42 );
+    EXPECT_EQ( m.getElement( 0, 4 ), 42 );
+    EXPECT_EQ( m.getElement( 0, 5 ), 42 );
+    
+    EXPECT_EQ( m.getElement( 1, 0 ), 42 );
+    EXPECT_EQ( m.getElement( 1, 1 ), 42 );
+    EXPECT_EQ( m.getElement( 1, 2 ), 42 );
+    EXPECT_EQ( m.getElement( 1, 3 ), 42 );
+    EXPECT_EQ( m.getElement( 1, 4 ), 42 );
+    EXPECT_EQ( m.getElement( 1, 5 ), 42 );
+    
+    EXPECT_EQ( m.getElement( 2, 0 ), 42 );
+    EXPECT_EQ( m.getElement( 2, 1 ), 42 );
+    EXPECT_EQ( m.getElement( 2, 2 ), 42 );
+    EXPECT_EQ( m.getElement( 2, 3 ), 42 );
+    EXPECT_EQ( m.getElement( 2, 4 ), 42 );
+    EXPECT_EQ( m.getElement( 2, 5 ), 42 );
+    
+    EXPECT_EQ( m.getElement( 3, 0 ), 42 );
+    EXPECT_EQ( m.getElement( 3, 1 ), 42 );
+    EXPECT_EQ( m.getElement( 3, 2 ), 42 );
+    EXPECT_EQ( m.getElement( 3, 3 ), 42 );
+    EXPECT_EQ( m.getElement( 3, 4 ), 42 );
+    EXPECT_EQ( m.getElement( 3, 5 ), 42 );
+    
+    EXPECT_EQ( m.getElement( 4, 0 ), 42 );
+    EXPECT_EQ( m.getElement( 4, 1 ), 42 );
+    EXPECT_EQ( m.getElement( 4, 2 ), 42 );
+    EXPECT_EQ( m.getElement( 4, 3 ), 42 );
+    EXPECT_EQ( m.getElement( 4, 4 ), 42 );
+    EXPECT_EQ( m.getElement( 4, 5 ), 42 );
+    
+    EXPECT_EQ( m.getElement( 5, 0 ), 42 );
+    EXPECT_EQ( m.getElement( 5, 1 ), 42 );
+    EXPECT_EQ( m.getElement( 5, 2 ), 42 );
+    EXPECT_EQ( m.getElement( 5, 3 ), 42 );
+    EXPECT_EQ( m.getElement( 5, 4 ), 42 );
+    EXPECT_EQ( m.getElement( 5, 5 ), 42 );
+    
+    EXPECT_EQ( m.getElement( 6, 0 ), 42 );
+    EXPECT_EQ( m.getElement( 6, 1 ), 42 );
+    EXPECT_EQ( m.getElement( 6, 2 ), 42 );
+    EXPECT_EQ( m.getElement( 6, 3 ), 42 );
+    EXPECT_EQ( m.getElement( 6, 4 ), 42 );
+    EXPECT_EQ( m.getElement( 6, 5 ), 42 );
+}
+
+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 5x5 dense 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 /
+ */
+    const IndexType rows = 5;
+    const IndexType cols = 5;
+    
+    Matrix m;
+    m.reset();
+    m.setDimensions( rows, cols );    
+    
+    RealType value = 1;
+    for( IndexType i = 0; i < rows; i++ )
+        for( IndexType j = 0; j < cols; j++ )
+            m.setElement( i, j, value++ );
+    
+    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 ),  4 );
+    EXPECT_EQ( m.getElement( 0, 4 ),  5 );
+    
+    EXPECT_EQ( m.getElement( 1, 0 ),  6 );
+    EXPECT_EQ( m.getElement( 1, 1 ),  7 );
+    EXPECT_EQ( m.getElement( 1, 2 ),  8 );
+    EXPECT_EQ( m.getElement( 1, 3 ),  9 );
+    EXPECT_EQ( m.getElement( 1, 4 ), 10 );
+    
+    EXPECT_EQ( m.getElement( 2, 0 ), 11 );
+    EXPECT_EQ( m.getElement( 2, 1 ), 12 );
+    EXPECT_EQ( m.getElement( 2, 2 ), 13 );
+    EXPECT_EQ( m.getElement( 2, 3 ), 14 );
+    EXPECT_EQ( m.getElement( 2, 4 ), 15 );
+    
+    EXPECT_EQ( m.getElement( 3, 0 ), 16 );
+    EXPECT_EQ( m.getElement( 3, 1 ), 17 );
+    EXPECT_EQ( m.getElement( 3, 2 ), 18 );
+    EXPECT_EQ( m.getElement( 3, 3 ), 19 );
+    EXPECT_EQ( m.getElement( 3, 4 ), 20 );
+    
+    EXPECT_EQ( m.getElement( 4, 0 ), 21 );
+    EXPECT_EQ( m.getElement( 4, 1 ), 22 );
+    EXPECT_EQ( m.getElement( 4, 2 ), 23 );
+    EXPECT_EQ( m.getElement( 4, 3 ), 24 );
+    EXPECT_EQ( m.getElement( 4, 4 ), 25 );
+}
+
+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 dense 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 29 30 /
+ */
+    const IndexType rows = 6;
+    const IndexType cols = 5;
+    
+    Matrix m;
+    m.reset();
+    m.setDimensions( rows, cols );
+    
+    RealType value = 1;
+    for( IndexType i = 0; i < rows; i++ )
+        for( IndexType j = 0; j < cols; j++ )
+            m.setElement( i, j, value++ );
+    
+    // Check the added 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 ),  4 );
+    EXPECT_EQ( m.getElement( 0, 4 ),  5 );
+    
+    EXPECT_EQ( m.getElement( 1, 0 ),  6 );
+    EXPECT_EQ( m.getElement( 1, 1 ),  7 );
+    EXPECT_EQ( m.getElement( 1, 2 ),  8 );
+    EXPECT_EQ( m.getElement( 1, 3 ),  9 );
+    EXPECT_EQ( m.getElement( 1, 4 ), 10 );
+    
+    EXPECT_EQ( m.getElement( 2, 0 ), 11 );
+    EXPECT_EQ( m.getElement( 2, 1 ), 12 );
+    EXPECT_EQ( m.getElement( 2, 2 ), 13 );
+    EXPECT_EQ( m.getElement( 2, 3 ), 14 );
+    EXPECT_EQ( m.getElement( 2, 4 ), 15 );
+    
+    EXPECT_EQ( m.getElement( 3, 0 ), 16 );
+    EXPECT_EQ( m.getElement( 3, 1 ), 17 );
+    EXPECT_EQ( m.getElement( 3, 2 ), 18 );
+    EXPECT_EQ( m.getElement( 3, 3 ), 19 );
+    EXPECT_EQ( m.getElement( 3, 4 ), 20 );
+    
+    EXPECT_EQ( m.getElement( 4, 0 ), 21 );
+    EXPECT_EQ( m.getElement( 4, 1 ), 22 );
+    EXPECT_EQ( m.getElement( 4, 2 ), 23 );
+    EXPECT_EQ( m.getElement( 4, 3 ), 24 );
+    EXPECT_EQ( m.getElement( 4, 4 ), 25 );
+    
+    EXPECT_EQ( m.getElement( 5, 0 ), 26 );
+    EXPECT_EQ( m.getElement( 5, 1 ), 27 );
+    EXPECT_EQ( m.getElement( 5, 2 ), 28 );
+    EXPECT_EQ( m.getElement( 5, 3 ), 29 );
+    EXPECT_EQ( m.getElement( 5, 4 ), 30 );
+    
+    // Add new elements to the old elements with a multiplying factor applied to the old elements.
+/*
+ * The following setup results in the following 6x5 dense matrix:
+ *
+ *    /  3  6  9 12 15 \
+ *    | 18 21 24 27 30 |
+ *    | 33 36 39 42 45 |
+ *    | 48 51 54 57 60 |
+ *    | 63 66 69 72 75 |
+ *    \ 78 81 84 87 90 /
+ */
+    RealType newValue = 1;
+    RealType multiplicator = 2;
+    for( IndexType i = 0; i < rows; i++ )
+        for( IndexType j = 0; j < cols; j++ )
+            m.addElement( i, j, newValue++, multiplicator );    
+    
+    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 ), 12 );
+    EXPECT_EQ( m.getElement( 0, 4 ), 15 );
+    
+    EXPECT_EQ( m.getElement( 1, 0 ), 18 );
+    EXPECT_EQ( m.getElement( 1, 1 ), 21 );
+    EXPECT_EQ( m.getElement( 1, 2 ), 24 );
+    EXPECT_EQ( m.getElement( 1, 3 ), 27 );
+    EXPECT_EQ( m.getElement( 1, 4 ), 30 );
+    
+    EXPECT_EQ( m.getElement( 2, 0 ), 33 );
+    EXPECT_EQ( m.getElement( 2, 1 ), 36 );
+    EXPECT_EQ( m.getElement( 2, 2 ), 39 );
+    EXPECT_EQ( m.getElement( 2, 3 ), 42 );
+    EXPECT_EQ( m.getElement( 2, 4 ), 45 );
+    
+    EXPECT_EQ( m.getElement( 3, 0 ), 48 );
+    EXPECT_EQ( m.getElement( 3, 1 ), 51 );
+    EXPECT_EQ( m.getElement( 3, 2 ), 54 );
+    EXPECT_EQ( m.getElement( 3, 3 ), 57 );
+    EXPECT_EQ( m.getElement( 3, 4 ), 60 );
+    
+    EXPECT_EQ( m.getElement( 4, 0 ), 63 );
+    EXPECT_EQ( m.getElement( 4, 1 ), 66 );
+    EXPECT_EQ( m.getElement( 4, 2 ), 69 );
+    EXPECT_EQ( m.getElement( 4, 3 ), 72 );
+    EXPECT_EQ( m.getElement( 4, 4 ), 75 );
+    
+    EXPECT_EQ( m.getElement( 5, 0 ), 78 );
+    EXPECT_EQ( m.getElement( 5, 1 ), 81 );
+    EXPECT_EQ( m.getElement( 5, 2 ), 84 );
+    EXPECT_EQ( m.getElement( 5, 3 ), 87 );
+    EXPECT_EQ( m.getElement( 5, 4 ), 90 );
+}
+
+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 dense matrix:
+ *
+ *    /  1  2  3  4  5  6  7 \
+ *    |  8  9 10 11 12 13 14 |
+ *    \ 15 16 17 18 19 20 21 /
+ */
+    const IndexType rows = 3;
+    const IndexType cols = 7;
+    
+    Matrix m;
+    m.reset();
+    m.setDimensions( rows, cols );
+    
+    RealType value = 1;
+    for( IndexType i = 0; i < rows; i++ )
+        for( IndexType j = 0; j < cols; j++ )
+            m.setElement( i, j, value++ );       
+    
+    RealType row1 [ 5 ] = { 11, 11, 11, 11, 11 }; IndexType colIndexes1 [ 5 ] = { 0, 1, 2, 3, 4 };
+    RealType row2 [ 5 ] = { 22, 22, 22, 22, 22 }; IndexType colIndexes2 [ 5 ] = { 0, 1, 2, 3, 4 };
+    RealType row3 [ 5 ] = { 33, 33, 33, 33, 33 }; IndexType colIndexes3 [ 5 ] = { 2, 3, 4, 5, 6 };
+    
+    IndexType row = 0;
+    IndexType elements = 5;
+    
+    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 ), 11 );
+    EXPECT_EQ( m.getElement( 0, 4 ), 11 );
+    EXPECT_EQ( m.getElement( 0, 5 ),  6 );
+    EXPECT_EQ( m.getElement( 0, 6 ),  7 );
+    
+    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 ), 22 );
+    EXPECT_EQ( m.getElement( 1, 4 ), 22 );
+    EXPECT_EQ( m.getElement( 1, 5 ), 13 );
+    EXPECT_EQ( m.getElement( 1, 6 ), 14 );
+    
+    EXPECT_EQ( m.getElement( 2, 0 ), 15 );
+    EXPECT_EQ( m.getElement( 2, 1 ), 16 );
+    EXPECT_EQ( m.getElement( 2, 2 ), 33 );
+    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 ), 33 );
+}
+
+template< typename Matrix >
+void test_AddRow()
+{
+    using RealType = typename Matrix::RealType;
+    using DeviceType = typename Matrix::DeviceType;
+    using IndexType = typename Matrix::IndexType;
+/*
+ * Sets up the following 6x5 dense 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 29 30 /
+ */
+    const IndexType rows = 6;
+    const IndexType cols = 5;
+    
+    Matrix m;
+    m.reset();
+    m.setDimensions( rows, cols );
+    
+    RealType value = 1;
+    for( IndexType i = 0; i < rows; i++ )
+        for( IndexType j = 0; j < cols; j++ )
+            m.setElement( i, j, value++ );
+    
+    // Check the added 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 ),  4 );
+    EXPECT_EQ( m.getElement( 0, 4 ),  5 );
+    
+    EXPECT_EQ( m.getElement( 1, 0 ),  6 );
+    EXPECT_EQ( m.getElement( 1, 1 ),  7 );
+    EXPECT_EQ( m.getElement( 1, 2 ),  8 );
+    EXPECT_EQ( m.getElement( 1, 3 ),  9 );
+    EXPECT_EQ( m.getElement( 1, 4 ), 10 );
+    
+    EXPECT_EQ( m.getElement( 2, 0 ), 11 );
+    EXPECT_EQ( m.getElement( 2, 1 ), 12 );
+    EXPECT_EQ( m.getElement( 2, 2 ), 13 );
+    EXPECT_EQ( m.getElement( 2, 3 ), 14 );
+    EXPECT_EQ( m.getElement( 2, 4 ), 15 );
+    
+    EXPECT_EQ( m.getElement( 3, 0 ), 16 );
+    EXPECT_EQ( m.getElement( 3, 1 ), 17 );
+    EXPECT_EQ( m.getElement( 3, 2 ), 18 );
+    EXPECT_EQ( m.getElement( 3, 3 ), 19 );
+    EXPECT_EQ( m.getElement( 3, 4 ), 20 );
+    
+    EXPECT_EQ( m.getElement( 4, 0 ), 21 );
+    EXPECT_EQ( m.getElement( 4, 1 ), 22 );
+    EXPECT_EQ( m.getElement( 4, 2 ), 23 );
+    EXPECT_EQ( m.getElement( 4, 3 ), 24 );
+    EXPECT_EQ( m.getElement( 4, 4 ), 25 );
+    
+    EXPECT_EQ( m.getElement( 5, 0 ), 26 );
+    EXPECT_EQ( m.getElement( 5, 1 ), 27 );
+    EXPECT_EQ( m.getElement( 5, 2 ), 28 );
+    EXPECT_EQ( m.getElement( 5, 3 ), 29 );
+    EXPECT_EQ( m.getElement( 5, 4 ), 30 );
+    
+    // Add new elements to the old elements with a multiplying factor applied to the old elements.
+/*
+ * The following setup results in the following 6x5 sparse matrix:
+ *
+ *    /  3  6  9 12 15 \
+ *    | 18 21 24 27 30 |
+ *    | 33 36 39 42 45 |
+ *    | 48 51 54 57 60 |
+ *    | 63 66 69 72 75 |
+ *    \ 78 81 84 87 90 /
+ */
+    
+    RealType row0 [ 5 ] = { 11, 11, 11, 11, 0 }; IndexType colIndexes0 [ 5 ] = { 0, 1, 2, 3, 4 };
+    RealType row1 [ 5 ] = { 22, 22, 22, 22, 0 }; IndexType colIndexes1 [ 5 ] = { 0, 1, 2, 3, 4 };
+    RealType row2 [ 5 ] = { 33, 33, 33, 33, 0 }; IndexType colIndexes2 [ 5 ] = { 0, 1, 2, 3, 4 };
+    RealType row3 [ 5 ] = { 44, 44, 44, 44, 0 }; IndexType colIndexes3 [ 5 ] = { 0, 1, 2, 3, 4 };
+    RealType row4 [ 5 ] = { 55, 55, 55, 55, 0 }; IndexType colIndexes4 [ 5 ] = { 0, 1, 2, 3, 4 };
+    RealType row5 [ 5 ] = { 66, 66, 66, 66, 0 }; IndexType colIndexes5 [ 5 ] = { 0, 1, 2, 3, 4 };
+    
+    IndexType row = 0;
+    IndexType elements = 5;
+    RealType thisRowMultiplicator = 0;
+    
+    m.addRow( row++, colIndexes0, row0, elements, thisRowMultiplicator++ );
+    m.addRow( row++, colIndexes1, row1, elements, thisRowMultiplicator++ );
+    m.addRow( row++, colIndexes2, row2, elements, thisRowMultiplicator++ );
+    m.addRow( row++, colIndexes3, row3, elements, thisRowMultiplicator++ );
+    m.addRow( row++, colIndexes4, row4, elements, thisRowMultiplicator++ );
+    m.addRow( row++, colIndexes5, row5, elements, thisRowMultiplicator++ );
+    
+    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 ),  11 );
+    EXPECT_EQ( m.getElement( 0, 4 ),   0 );
+    
+    EXPECT_EQ( m.getElement( 1, 0 ),  28 );
+    EXPECT_EQ( m.getElement( 1, 1 ),  29 );
+    EXPECT_EQ( m.getElement( 1, 2 ),  30 );
+    EXPECT_EQ( m.getElement( 1, 3 ),  31 );
+    EXPECT_EQ( m.getElement( 1, 4 ),  10 );
+    
+    EXPECT_EQ( m.getElement( 2, 0 ),  55 );
+    EXPECT_EQ( m.getElement( 2, 1 ),  57 );
+    EXPECT_EQ( m.getElement( 2, 2 ),  59 );
+    EXPECT_EQ( m.getElement( 2, 3 ),  61 );
+    EXPECT_EQ( m.getElement( 2, 4 ),  30 );
+    
+    EXPECT_EQ( m.getElement( 3, 0 ),  92 );
+    EXPECT_EQ( m.getElement( 3, 1 ),  95 );
+    EXPECT_EQ( m.getElement( 3, 2 ),  98 );
+    EXPECT_EQ( m.getElement( 3, 3 ), 101 );
+    EXPECT_EQ( m.getElement( 3, 4 ),  60 );
+    
+    EXPECT_EQ( m.getElement( 4, 0 ), 139 );
+    EXPECT_EQ( m.getElement( 4, 1 ), 143 );
+    EXPECT_EQ( m.getElement( 4, 2 ), 147 );
+    EXPECT_EQ( m.getElement( 4, 3 ), 151 );
+    EXPECT_EQ( m.getElement( 4, 4 ), 100 );
+    
+    EXPECT_EQ( m.getElement( 5, 0 ), 196 );
+    EXPECT_EQ( m.getElement( 5, 1 ), 201 );
+    EXPECT_EQ( m.getElement( 5, 2 ), 206 );
+    EXPECT_EQ( m.getElement( 5, 3 ), 211 );
+    EXPECT_EQ( m.getElement( 5, 4 ), 150 );
+}
+
+template< typename Matrix >
+void test_VectorProduct()
+{
+    using RealType = typename Matrix::RealType;
+    using DeviceType = typename Matrix::DeviceType;
+    using IndexType = typename Matrix::IndexType;
+/*
+ * Sets up the following 5x4 dense matrix:
+ *
+ *    /  1  2  3  4 \
+ *    |  5  6  7  8 |
+ *    |  9 10 11 12 |
+ *    | 13 14 15 16 |
+ *    \ 17 18 19 20 /
+ */
+    const IndexType rows = 5;
+    const IndexType cols = 4;
+    
+    Matrix m;
+    m.reset();
+    m.setDimensions( rows, cols );
+    
+    RealType value = 1;
+    for( IndexType i = 0; i < rows; i++ )
+        for( IndexType j = 0; j < cols; j++)
+            m.setElement( i, j, value++ );
+
+    using VectorType = TNL::Containers::Vector< RealType, DeviceType, IndexType >;
+    
+    VectorType inVector;
+    inVector.setSize( 4 );
+    for( IndexType i = 0; i < inVector.getSize(); i++ )        
+        inVector.setElement( i, 2 );
+
+    VectorType outVector;  
+    outVector.setSize( 5 );
+    for( IndexType j = 0; j < outVector.getSize(); j++ )
+        outVector.setElement( j, 0 );
+ 
+    
+    m.vectorProduct( inVector, outVector);
+   
+    EXPECT_EQ( outVector.getElement( 0 ),  20 );
+    EXPECT_EQ( outVector.getElement( 1 ),  52 );
+    EXPECT_EQ( outVector.getElement( 2 ),  84 );
+    EXPECT_EQ( outVector.getElement( 3 ), 116 );
+    EXPECT_EQ( outVector.getElement( 4 ), 148 );
+}
+
+template< typename Matrix >
+void test_AddMatrix()
+{
+    using RealType = typename Matrix::RealType;
+    using DeviceType = typename Matrix::DeviceType;
+    using IndexType = typename Matrix::IndexType;
+/*
+ * Sets up the following 5x4 dense matrix:
+ *
+ *    /  1  2  3  4 \
+ *    |  5  6  7  8 |
+ *    |  9 10 11 12 |
+ *    | 13 14 15 16 |
+ *    \ 17 18 19 20 /
+ */
+    const IndexType rows = 5;
+    const IndexType cols = 4;
+    
+    Matrix m;                           // We need this matrix to preserve the values for EXPECT_EQ statements comparing the actual operation;
+    m.reset();
+    m.setDimensions( rows, cols );
+    
+    RealType value = 1;
+    for( IndexType i = 0; i < rows; i++ )
+        for( IndexType j = 0; j < cols; j++)
+            m.setElement( i, j, value++ );
+    
+/*
+ * Sets up the following 5x4 dense matrix:
+ *
+ *    /  1  2  3  4 \
+ *    |  5  6  7  8 |
+ *    |  9 10 11 12 |
+ *    | 13 14 15 16 |
+ *    \ 17 18 19 20 /
+ */
+    
+    Matrix m2;
+    m2.reset();
+    m2.setDimensions( rows, cols );
+    
+    RealType newValue = 1;
+    for( IndexType i = 0; i < rows; i++ )
+        for( IndexType j = 0; j < cols; j++)
+            m2.setElement( i, j, newValue++ );
+    
+    /*
+ * Sets up the following 5x4 dense matrix:
+ *
+ *    /  1  2  3  4 \
+ *    |  5  6  7  8 |
+ *    |  9 10 11 12 |
+ *    | 13 14 15 16 |
+ *    \ 17 18 19 20 /
+ */
+    
+    Matrix mResult;
+    mResult.reset();
+    mResult.setDimensions( rows, cols );
+    
+    mResult = m;
+    
+    RealType matrixMultiplicator = 2;
+    RealType thisMatrixMultiplicator = 1;
+    
+    mResult.addMatrix( m2, matrixMultiplicator, thisMatrixMultiplicator );
+    
+    EXPECT_EQ( mResult.getElement( 0, 0 ), matrixMultiplicator * m2.getElement( 0, 0 ) + thisMatrixMultiplicator * m.getElement( 0, 0 ) );
+    EXPECT_EQ( mResult.getElement( 0, 1 ), matrixMultiplicator * m2.getElement( 0, 1 ) + thisMatrixMultiplicator * m.getElement( 0, 1 ) );
+    EXPECT_EQ( mResult.getElement( 0, 2 ), matrixMultiplicator * m2.getElement( 0, 2 ) + thisMatrixMultiplicator * m.getElement( 0, 2 ) );
+    EXPECT_EQ( mResult.getElement( 0, 3 ), matrixMultiplicator * m2.getElement( 0, 3 ) + thisMatrixMultiplicator * m.getElement( 0, 3 ) );
+    
+    EXPECT_EQ( mResult.getElement( 1, 0 ), matrixMultiplicator * m2.getElement( 1, 0 ) + thisMatrixMultiplicator * m.getElement( 1, 0 ) );
+    EXPECT_EQ( mResult.getElement( 1, 1 ), matrixMultiplicator * m2.getElement( 1, 1 ) + thisMatrixMultiplicator * m.getElement( 1, 1 ) );
+    EXPECT_EQ( mResult.getElement( 1, 2 ), matrixMultiplicator * m2.getElement( 1, 2 ) + thisMatrixMultiplicator * m.getElement( 1, 2 ) );
+    EXPECT_EQ( mResult.getElement( 1, 3 ), matrixMultiplicator * m2.getElement( 1, 3 ) + thisMatrixMultiplicator * m.getElement( 1, 3 ) );
+    
+    EXPECT_EQ( mResult.getElement( 2, 0 ), matrixMultiplicator * m2.getElement( 2, 0 ) + thisMatrixMultiplicator * m.getElement( 2, 0 ) );
+    EXPECT_EQ( mResult.getElement( 2, 1 ), matrixMultiplicator * m2.getElement( 2, 1 ) + thisMatrixMultiplicator * m.getElement( 2, 1 ) );
+    EXPECT_EQ( mResult.getElement( 2, 2 ), matrixMultiplicator * m2.getElement( 2, 2 ) + thisMatrixMultiplicator * m.getElement( 2, 2 ) );
+    EXPECT_EQ( mResult.getElement( 2, 3 ), matrixMultiplicator * m2.getElement( 2, 3 ) + thisMatrixMultiplicator * m.getElement( 2, 3 ) );
+    
+    EXPECT_EQ( mResult.getElement( 3, 0 ), matrixMultiplicator * m2.getElement( 3, 0 ) + thisMatrixMultiplicator * m.getElement( 3, 0 ) );
+    EXPECT_EQ( mResult.getElement( 3, 1 ), matrixMultiplicator * m2.getElement( 3, 1 ) + thisMatrixMultiplicator * m.getElement( 3, 1 ) );
+    EXPECT_EQ( mResult.getElement( 3, 2 ), matrixMultiplicator * m2.getElement( 3, 2 ) + thisMatrixMultiplicator * m.getElement( 3, 2 ) );
+    EXPECT_EQ( mResult.getElement( 3, 3 ), matrixMultiplicator * m2.getElement( 3, 3 ) + thisMatrixMultiplicator * m.getElement( 3, 3 ) );
+    
+    EXPECT_EQ( mResult.getElement( 4, 0 ), matrixMultiplicator * m2.getElement( 4, 0 ) + thisMatrixMultiplicator * m.getElement( 4, 0 ) );
+    EXPECT_EQ( mResult.getElement( 4, 1 ), matrixMultiplicator * m2.getElement( 4, 1 ) + thisMatrixMultiplicator * m.getElement( 4, 1 ) );
+    EXPECT_EQ( mResult.getElement( 4, 2 ), matrixMultiplicator * m2.getElement( 4, 2 ) + thisMatrixMultiplicator * m.getElement( 4, 2 ) );
+    EXPECT_EQ( mResult.getElement( 4, 3 ), matrixMultiplicator * m2.getElement( 4, 3 ) + thisMatrixMultiplicator * m.getElement( 4, 3 ) );
+    
+    EXPECT_EQ( mResult.getElement( 0, 0 ),  3 );
+    EXPECT_EQ( mResult.getElement( 0, 1 ),  6 );
+    EXPECT_EQ( mResult.getElement( 0, 2 ),  9 );
+    EXPECT_EQ( mResult.getElement( 0, 3 ), 12 );
+    
+    EXPECT_EQ( mResult.getElement( 1, 0 ), 15 );
+    EXPECT_EQ( mResult.getElement( 1, 1 ), 18 );
+    EXPECT_EQ( mResult.getElement( 1, 2 ), 21 );
+    EXPECT_EQ( mResult.getElement( 1, 3 ), 24 );
+    
+    EXPECT_EQ( mResult.getElement( 2, 0 ), 27 );
+    EXPECT_EQ( mResult.getElement( 2, 1 ), 30 );
+    EXPECT_EQ( mResult.getElement( 2, 2 ), 33 );
+    EXPECT_EQ( mResult.getElement( 2, 3 ), 36 );
+    
+    EXPECT_EQ( mResult.getElement( 3, 0 ), 39 );
+    EXPECT_EQ( mResult.getElement( 3, 1 ), 42 );
+    EXPECT_EQ( mResult.getElement( 3, 2 ), 45 );
+    EXPECT_EQ( mResult.getElement( 3, 3 ), 48 );
+    
+    EXPECT_EQ( mResult.getElement( 4, 0 ), 51 );
+    EXPECT_EQ( mResult.getElement( 4, 1 ), 54 );
+    EXPECT_EQ( mResult.getElement( 4, 2 ), 57 );
+    EXPECT_EQ( mResult.getElement( 4, 3 ), 60 );
+}
+
+template< typename Matrix >
+void test_GetMatrixProduct()
+{
+    using RealType = typename Matrix::RealType;
+    using DeviceType = typename Matrix::DeviceType;
+    using IndexType = typename Matrix::IndexType;
+/*
+ * Sets up the following 5x4 dense matrix:
+ *
+ *    /  1  2  3  4 \
+ *    |  5  6  7  8 |
+ *    |  9 10 11 12 |
+ *    | 13 14 15 16 |
+ *    \ 17 18 19 20 /
+ */
+    const IndexType leftRows = 5;
+    const IndexType leftCols = 4;
+    
+    Matrix leftMatrix;
+    leftMatrix.reset();
+    leftMatrix.setDimensions( leftRows, leftCols );
+    
+    RealType value = 1;
+    for( IndexType i = 0; i < leftRows; i++ )
+        for( IndexType j = 0; j < leftCols; j++)
+            leftMatrix.setElement( i, j, value++ );
+    
+/*
+ * Sets up the following 4x5 dense matrix:
+ *
+ *    /  1  2  3  4  5 \
+ *    |  6  7  8  9 10 |
+ *    | 11 12 13 14 15 |
+ *    \ 16 17 18 19 20 /
+ */
+    const IndexType rightRows = 4;
+    const IndexType rightCols = 5;
+    
+    Matrix rightMatrix;
+    rightMatrix.reset();
+    rightMatrix.setDimensions( rightRows, rightCols );
+    
+    RealType newValue = 1;
+    for( IndexType i = 0; i < rightRows; i++ )
+        for( IndexType j = 0; j < rightCols; j++)
+            rightMatrix.setElement( i, j, newValue++ );
+    
+/*
+ * Sets up the following 5x5 resulting dense matrix:
+ *
+ *    /  0  0  0  0 \
+ *    |  0  0  0  0 |
+ *    |  0  0  0  0 |
+ *    |  0  0  0  0 |
+ *    \  0  0  0  0 /
+ */
+    
+    Matrix mResult;
+    mResult.reset();
+    mResult.setDimensions( leftRows, rightCols );
+    mResult.setValue( 0 );
+    
+    RealType leftMatrixMultiplicator = 1;
+    RealType rightMatrixMultiplicator = 2;
+/*   
+ *      /  1  2  3  4 \                            /  220  240  260  280  300 \
+ *      |  5  6  7  8 |       /  1  2  3  4  5 \   |  492  544  596  648  700 |
+ *  1 * |  9 10 11 12 | * 2 * |  6  7  8  9 10 | = |  764  848  932 1016 1100 |
+ *      | 13 14 15 16 |       | 11 12 13 14 15 |   | 1036 1152 1268 1384 1500 |
+ *      \ 17 18 19 20 /       \ 16 17 18 19 20 /   \ 1308 1456 1604 1752 1900 /
+ */
+    
+    mResult.getMatrixProduct( leftMatrix, rightMatrix, leftMatrixMultiplicator, rightMatrixMultiplicator );
+    
+    EXPECT_EQ( mResult.getElement( 0, 0 ),  220 );
+    EXPECT_EQ( mResult.getElement( 0, 1 ),  240 );
+    EXPECT_EQ( mResult.getElement( 0, 2 ),  260 );
+    EXPECT_EQ( mResult.getElement( 0, 3 ),  280 );
+    EXPECT_EQ( mResult.getElement( 0, 4 ),  300 );
+    
+    EXPECT_EQ( mResult.getElement( 1, 0 ),  492 );
+    EXPECT_EQ( mResult.getElement( 1, 1 ),  544 );
+    EXPECT_EQ( mResult.getElement( 1, 2 ),  596 );
+    EXPECT_EQ( mResult.getElement( 1, 3 ),  648 );
+    EXPECT_EQ( mResult.getElement( 1, 4 ),  700 );
+    
+    EXPECT_EQ( mResult.getElement( 2, 0 ),  764 );
+    EXPECT_EQ( mResult.getElement( 2, 1 ),  848 );
+    EXPECT_EQ( mResult.getElement( 2, 2 ),  932 );
+    EXPECT_EQ( mResult.getElement( 2, 3 ), 1016 );
+    EXPECT_EQ( mResult.getElement( 2, 4 ), 1100 );
+    
+    EXPECT_EQ( mResult.getElement( 3, 0 ), 1036 );
+    EXPECT_EQ( mResult.getElement( 3, 1 ), 1152 );
+    EXPECT_EQ( mResult.getElement( 3, 2 ), 1268 );
+    EXPECT_EQ( mResult.getElement( 3, 3 ), 1384 );
+    EXPECT_EQ( mResult.getElement( 3, 4 ), 1500 );
+    
+    EXPECT_EQ( mResult.getElement( 4, 0 ), 1308 );
+    EXPECT_EQ( mResult.getElement( 4, 1 ), 1456 );
+    EXPECT_EQ( mResult.getElement( 4, 2 ), 1604 );
+    EXPECT_EQ( mResult.getElement( 4, 3 ), 1752 );
+    EXPECT_EQ( mResult.getElement( 4, 4 ), 1900 );
+}
+
+template< typename Matrix >
+void test_GetTransposition()
+{
+    using RealType = typename Matrix::RealType;
+    using DeviceType = typename Matrix::DeviceType;
+    using IndexType = typename Matrix::IndexType;
+/*
+ * Sets up the following 3x2 dense matrix:
+ *
+ *    /  1  2 \
+ *    |  3  4 |
+ *    \  5  6 /
+ */
+    const IndexType rows = 3;
+    const IndexType cols = 2;
+
+    Matrix m;
+    m.reset();
+    m.setDimensions( rows, cols );
+
+    RealType value = 1;
+    for( IndexType i = 0; i < rows; i++ )
+        for( IndexType j = 0; j < cols; j++ )
+            m.setElement( i, j, value++ );
+
+    m.print( std::cout );
+    
+/*
+ * Sets up the following 2x3 dense matrix:
+ *
+ *    /  0  0  0 \
+ *    \  0  0  0 /
+ */ 
+    Matrix mTransposed;
+    mTransposed.reset();
+    mTransposed.setDimensions( cols, rows );
+    
+    mTransposed.print( std::cout );
+    
+    RealType matrixMultiplicator = 1;
+    
+    mTransposed.getTransposition( m, matrixMultiplicator );
+    
+    mTransposed.print( std::cout );
+    
+/*
+ * Should result in the following 2x3 dense matrix:
+ *
+ *    /  1  3  5 \
+ *    \  2  4  6 /
+ */ 
+    
+    EXPECT_EQ( mTransposed.getElement( 0, 0 ), 1 );
+    EXPECT_EQ( mTransposed.getElement( 0, 1 ), 3 );
+    EXPECT_EQ( mTransposed.getElement( 0, 2 ), 5 );
+    
+    EXPECT_EQ( mTransposed.getElement( 1, 0 ), 2 );
+    EXPECT_EQ( mTransposed.getElement( 1, 1 ), 4 );
+    EXPECT_EQ( mTransposed.getElement( 1, 2 ), 6 );
+    
+/*
+ * Sets up the following 5x5 dense 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 /
+ */
+            //    const IndexType rows = 5;
+            //    const IndexType cols = 5;
+            //    
+            //    Matrix m;
+            //    m.reset();
+            //    m.setDimensions( rows, cols );
+            //    
+            //    RealType value = 1;
+            //    for( IndexType i = 0; i < rows; i++ )
+            //        for( IndexType j = 0; j < cols; j++)
+            //            m.setElement( i, j, value++ );
+    
+/*
+ * Sets up the following 5x5 dense matrix:
+ *
+ *    /  2 12 22 32 42 \
+ *    |  4 14 24 34 44 |
+ *    |  6 16 26 36 46 |
+ *    |  8 18 28 38 48 |
+ *    \ 10 20 30 40 50 /
+ */
+            //    const IndexType resultRows = cols;
+            //    const IndexType resultCols = rows;
+            //    
+            //    Matrix mResult;
+            //    mResult.reset();
+            //    mResult.setDimensions( resultRows, resultCols );
+            //    mResult.setValue( 0 );
+            //    
+            //    RealType matrixMultiplicator = 2;
+            //    
+            //    mResult.getTransposition( m, matrixMultiplicator );
+    
+/*
+ * Should result in the following 5x5 resulting dense matrix:
+ *
+ *    /  0  0  0  0  0 \
+ *    |  0  0  0  0  0 |
+ *    |  0  0  0  0  0 |
+ *    |  0  0  0  0  0 |
+ *    \  0  0  0  0  0 /
+ */
+            //    
+            //    EXPECT_EQ( mResult.getElement( 0, 0 ),  2 );
+            //    EXPECT_EQ( mResult.getElement( 0, 1 ), 12 );
+            //    EXPECT_EQ( mResult.getElement( 0, 2 ), 22 );
+            //    EXPECT_EQ( mResult.getElement( 0, 3 ), 32 );
+            //    EXPECT_EQ( mResult.getElement( 0, 4 ), 42 );
+            //    
+            //    EXPECT_EQ( mResult.getElement( 1, 0 ),  4 );
+            //    EXPECT_EQ( mResult.getElement( 1, 1 ), 14 );
+            //    EXPECT_EQ( mResult.getElement( 1, 2 ), 24 );
+            //    EXPECT_EQ( mResult.getElement( 1, 3 ), 34 );
+            //    EXPECT_EQ( mResult.getElement( 1, 4 ), 44 );
+            //    
+            //    EXPECT_EQ( mResult.getElement( 2, 0 ),  6 );
+            //    EXPECT_EQ( mResult.getElement( 2, 1 ), 16 );
+            //    EXPECT_EQ( mResult.getElement( 2, 2 ), 26 );
+            //    EXPECT_EQ( mResult.getElement( 2, 3 ), 36 );
+            //    EXPECT_EQ( mResult.getElement( 2, 4 ), 46 );
+            //    
+            //    EXPECT_EQ( mResult.getElement( 3, 0 ),  8 );
+            //    EXPECT_EQ( mResult.getElement( 3, 1 ), 18 );
+            //    EXPECT_EQ( mResult.getElement( 3, 2 ), 28 );
+            //    EXPECT_EQ( mResult.getElement( 3, 3 ), 38 );
+            //    EXPECT_EQ( mResult.getElement( 3, 4 ), 48 );
+            //    
+            //    EXPECT_EQ( mResult.getElement( 4, 0 ), 10 );
+            //    EXPECT_EQ( mResult.getElement( 4, 1 ), 20 );
+            //    EXPECT_EQ( mResult.getElement( 4, 2 ), 30 );
+            //    EXPECT_EQ( mResult.getElement( 4, 3 ), 40 );
+            //    EXPECT_EQ( mResult.getElement( 4, 4 ), 50 );
+}
+
+
+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 dense matrix:
+ *
+ *    /  4  1  1  1 \
+ *    |  1  4  1  1 |
+ *    |  1  1  4  1 |
+ *    \  1  1  1  4 /
+ */
+    const IndexType rows = 4;
+    const IndexType cols = 4;
+    
+    Matrix m;
+    m.reset();
+    m.setDimensions( rows, cols );
+    
+    m.setElement( 0, 0, 4.0 );        // 0th row
+    m.setElement( 0, 1, 1.0 );
+    m.setElement( 0, 2, 1.0 );
+    m.setElement( 0, 3, 1.0 );
+        
+    m.setElement( 1, 0, 1.0 );        // 1st row
+    m.setElement( 1, 1, 4.0 );
+    m.setElement( 1, 2, 1.0 );
+    m.setElement( 1, 3, 1.0 );
+        
+    m.setElement( 2, 0, 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, 0, 1.0 );        // 3rd row
+    m.setElement( 3, 1, 1.0 );
+    m.setElement( 3, 2, 1.0 );
+    m.setElement( 3, 3, 4.0 );
+    
+    RealType bVector [ 4 ] = { 1.0, 1.0, 1.0, 1.0 };
+    RealType xVector [ 4 ] = { 1.0, 1.0, 1.0, 1.0 };
+    
+    IndexType row = 0;
+    RealType omega = 1;
+    
+    m.performSORIteration( bVector, row++, xVector, omega);
+    
+    EXPECT_EQ( xVector[ 0 ], -0.5 );
+    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.5 );
+    EXPECT_EQ( xVector[ 1 ], -0.125 );
+    EXPECT_EQ( xVector[ 2 ],  1.0 );
+    EXPECT_EQ( xVector[ 3 ],  1.0 );
+    
+    m.performSORIteration( bVector, row++, xVector, omega);
+    
+    EXPECT_EQ( xVector[ 0 ], -0.5 );
+    EXPECT_EQ( xVector[ 1 ], -0.125 );
+    EXPECT_EQ( xVector[ 2 ],  0.15625 );
+    EXPECT_EQ( xVector[ 3 ],  1.0 );
+    
+    m.performSORIteration( bVector, row++, xVector, omega);
+    
+    EXPECT_EQ( xVector[ 0 ], -0.5 );
+    EXPECT_EQ( xVector[ 1 ], -0.125 );
+    EXPECT_EQ( xVector[ 2 ], 0.15625 );
+    EXPECT_EQ( xVector[ 3 ], 0.3671875 );
+}
+
+template< typename Matrix >
+void test_SaveAndLoad()
+{
+    using RealType = typename Matrix::RealType;
+    using DeviceType = typename Matrix::DeviceType;
+    using IndexType = typename Matrix::IndexType;
+/*
+ * Sets up the following 4x4 dense matrix:
+ *
+ *    /  1  2  3  4 \
+ *    |  5  6  7  8 |
+ *    |  9 10 11 12 |
+ *    \ 13 14 15 16 /
+ */
+    const IndexType rows = 4;
+    const IndexType cols = 4;
+    
+    Matrix savedMatrix;
+    savedMatrix.reset();
+    savedMatrix.setDimensions( rows, cols );
+    
+    RealType value = 1;
+    for( IndexType i = 0; i < rows; i++ )
+        for( IndexType j = 0; j < cols; j++ )
+            savedMatrix.setElement( i, j, value++ );
+        
+    savedMatrix.save( "denseMatrixFile" );
+    
+    Matrix loadedMatrix;
+    loadedMatrix.reset();
+    loadedMatrix.setDimensions( rows, cols );
+    
+    loadedMatrix.load( "denseMatrixFile" );
+    
+    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 ),  4 );
+    
+    EXPECT_EQ( savedMatrix.getElement( 1, 0 ),  5 );
+    EXPECT_EQ( savedMatrix.getElement( 1, 1 ),  6 );
+    EXPECT_EQ( savedMatrix.getElement( 1, 2 ),  7 );
+    EXPECT_EQ( savedMatrix.getElement( 1, 3 ),  8 );
+    
+    EXPECT_EQ( savedMatrix.getElement( 2, 0 ),  9 );
+    EXPECT_EQ( savedMatrix.getElement( 2, 1 ), 10 );
+    EXPECT_EQ( savedMatrix.getElement( 2, 2 ), 11 );
+    EXPECT_EQ( savedMatrix.getElement( 2, 3 ), 12 );
+    
+    EXPECT_EQ( savedMatrix.getElement( 3, 0 ), 13 );
+    EXPECT_EQ( savedMatrix.getElement( 3, 1 ), 14 );
+    EXPECT_EQ( savedMatrix.getElement( 3, 2 ), 15 );
+    EXPECT_EQ( savedMatrix.getElement( 3, 3 ), 16 );
+    
+    std::cout << "\nThis will create a file called 'denseMatrixFile' (of the matrix created in the test function), in .../tnl-dev/Debug/bin/\n\n";
+}
+
+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  4 \
+ *    |  5  6  7  8 |
+ *    |  9 10 11 12 |
+ *    | 13 14 15 16 |
+ *    \ 17 18 19 20 /
+ */
+    const IndexType rows = 5;
+    const IndexType cols = 4;
+    
+    Matrix m;
+    m.reset();
+    m.setDimensions( rows, cols );
+    
+    RealType value = 1;
+    for( IndexType i = 0; i < rows; i++)
+        for( IndexType j = 0; j < cols; j++)
+            m.setElement( i, j, value++ );
+    
+    // This is from: https://stackoverflow.com/questions/5193173/getting-cout-output-to-a-stdstring
+    #include <sstream>
+    std::stringstream printed;
+    std::stringstream couted;
+    
+    // This is from: https://stackoverflow.com/questions/19485536/redirect-output-of-an-function-printing-to-console-to-string
+    //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
+    
+    //printed << printed.str() << std::endl;
+    couted << "Row: 0 ->  Col:0->1	 Col:1->2	 Col:2->3	 Col:3->4\t\n"
+              "Row: 1 ->  Col:0->5	 Col:1->6	 Col:2->7	 Col:3->8\t\n"
+              "Row: 2 ->  Col:0->9	 Col:1->10	 Col:2->11	 Col:3->12\t\n"
+              "Row: 3 ->  Col:0->13	 Col:1->14	 Col:2->15	 Col:3->16\t\n"
+              "Row: 4 ->  Col:0->17	 Col:1->18	 Col:2->19	 Col:3->20\t\n";
+    
+    EXPECT_EQ( printed.str(), couted.str() );
+}
+
+// test fixture for typed tests
+template< typename Matrix >
+class MatrixTest : public ::testing::Test
+{
+protected:
+   using MatrixType = Matrix;
+};
+
+// types for which MatrixTest is instantiated
+using MatrixTypes = ::testing::Types
+<
+    TNL::Matrices::Dense< int,    TNL::Devices::Host, short >,
+    TNL::Matrices::Dense< long,   TNL::Devices::Host, short >,
+    TNL::Matrices::Dense< float,  TNL::Devices::Host, short >,
+    TNL::Matrices::Dense< double, TNL::Devices::Host, short >,
+    TNL::Matrices::Dense< int,    TNL::Devices::Host, int >,
+    TNL::Matrices::Dense< long,   TNL::Devices::Host, int >,
+    TNL::Matrices::Dense< float,  TNL::Devices::Host, int >,
+    TNL::Matrices::Dense< double, TNL::Devices::Host, int >,
+    TNL::Matrices::Dense< int,    TNL::Devices::Host, long >,
+    TNL::Matrices::Dense< long,   TNL::Devices::Host, long >,
+    TNL::Matrices::Dense< float,  TNL::Devices::Host, long >,
+    TNL::Matrices::Dense< double, TNL::Devices::Host, long >
+#ifdef HAVE_CUDA
+    ,TNL::Matrices::Dense< int,    TNL::Devices::Cuda, short >,
+    TNL::Matrices::Dense< long,   TNL::Devices::Cuda, short >,
+    TNL::Matrices::Dense< float,  TNL::Devices::Cuda, short >,
+    TNL::Matrices::Dense< double, TNL::Devices::Cuda, short >,
+    TNL::Matrices::Dense< int,    TNL::Devices::Cuda, int >,
+    TNL::Matrices::Dense< long,   TNL::Devices::Cuda, int >,
+    TNL::Matrices::Dense< float,  TNL::Devices::Cuda, int >,
+    TNL::Matrices::Dense< double, TNL::Devices::Cuda, int >,
+    TNL::Matrices::Dense< int,    TNL::Devices::Cuda, long >,
+    TNL::Matrices::Dense< long,   TNL::Devices::Cuda, long >,
+    TNL::Matrices::Dense< float,  TNL::Devices::Cuda, long >,
+    TNL::Matrices::Dense< double, TNL::Devices::Cuda, long >
+#endif
+>;
+
+TYPED_TEST_CASE( MatrixTest, MatrixTypes );
+
+TYPED_TEST( MatrixTest, setDimensionsTest )
+{
+    using MatrixType = typename TestFixture::MatrixType;
+    
+    test_SetDimensions< MatrixType >();
+}
+
+TYPED_TEST( MatrixTest, setLikeTest )
+{
+    using MatrixType = typename TestFixture::MatrixType;
+    
+    test_SetLike< MatrixType, MatrixType >();
+}
+
+TYPED_TEST( MatrixTest, getRowLengthTest )
+{
+    using MatrixType = typename TestFixture::MatrixType;
+    
+    test_GetRowLength< MatrixType >();
+}
+
+TYPED_TEST( MatrixTest, getNumberOfMatrixElementsTest )
+{
+    using MatrixType = typename TestFixture::MatrixType;
+    
+    test_GetNumberOfMatrixElements< MatrixType >();
+}
+
+TYPED_TEST( MatrixTest, getNumberOfNonzeroMatrixElementsTest )
+{
+    using MatrixType = typename TestFixture::MatrixType;
+    
+    test_GetNumberOfNonzeroMatrixElements< MatrixType >();
+}
+
+TYPED_TEST( MatrixTest, resetTest )
+{
+    using MatrixType = typename TestFixture::MatrixType;
+    
+    test_Reset< MatrixType >();
+}
+
+TYPED_TEST( MatrixTest, setValueTest )
+{
+    using MatrixType = typename TestFixture::MatrixType;
+    
+    test_SetValue< MatrixType >();
+}
+
+TYPED_TEST( MatrixTest, setElementTest )
+{
+    using MatrixType = typename TestFixture::MatrixType;
+    
+    test_SetElement< MatrixType >();
+}
+
+TYPED_TEST( MatrixTest, addElementTest )
+{
+    using MatrixType = typename TestFixture::MatrixType;
+    
+    test_AddElement< MatrixType >();
+}
+
+TYPED_TEST( MatrixTest, setRowTest )
+{
+    using MatrixType = typename TestFixture::MatrixType;
+    
+    test_SetRow< MatrixType >();
+}
+
+TYPED_TEST( MatrixTest, addRowTest )
+{
+    using MatrixType = typename TestFixture::MatrixType;
+    
+    test_AddRow< MatrixType >();
+}
+
+TYPED_TEST( MatrixTest, vectorProductTest )
+{
+    using MatrixType = typename TestFixture::MatrixType;
+    
+    test_VectorProduct< MatrixType >();
+}
+
+TYPED_TEST( MatrixTest, addMatrixTest )
+{
+    using MatrixType = typename TestFixture::MatrixType;
+    
+    test_AddMatrix< MatrixType >();
+}
+
+TYPED_TEST( MatrixTest, saveAndLoadTest )
+{
+    using MatrixType = typename TestFixture::MatrixType;
+    
+    test_SaveAndLoad< MatrixType >();
+}
+
+TYPED_TEST( MatrixTest, printTest )
+{
+    using MatrixType = typename TestFixture::MatrixType;
+    
+    test_Print< MatrixType >();
+}
+
+//// test_getType is not general enough yet. DO NOT TEST IT YET.
+
+//TEST( DenseMatrixTest, Dense_GetTypeTest_Host )
+//{
+//    host_test_GetType< Dense_host_float, Dense_host_int >();
+//}
+//
+//#ifdef HAVE_CUDA
+//TEST( DenseMatrixTest, Dense_GetTypeTest_Cuda )
+//{
+//    cuda_test_GetType< Dense_cuda_float, Dense_cuda_int >();
+//}
+//#endif
+
+TEST( DenseMatrixTest, Dense_getMatrixProductTest_Host )
+{
+//    test_GetMatrixProduct< Dense_host_int >();
+    bool testRan = false;
+    EXPECT_TRUE( testRan );
+    std::cout << "\nTEST DID NOT RUN. NOT WORKING.\n\n";
+    std::cout << "If launched on CPU, this test will not build, but will print the following message: \n";
+    std::cout << "      /home/lukas/tnl-dev/src/TNL/Matrices/Dense_impl.h(609): error: no instance of function template \"TNL::Matrices::DenseMatrixProductKernel\" matches the argument list\n";
+    std::cout << "              argument types are: (TNL::Matrices::Dense<int, TNL::Devices::Host, int> *, Dense_host_int *, Dense_host_int *, const int, const int, int, int)\n";
+    std::cout << "          detected during:\n";
+    std::cout << "              instantiation of \"void TNL::Matrices::Dense<Real, Device, Index>::getMatrixProduct(const Matrix1 &, const Matrix2 &, const TNL::Matrices::Dense<Real, Device, Index>::RealType &, const TNL::Matrices::Dense<Real, Device, Index>::RealType &) [with Real=int, Device=TNL::Devices::Host, Index=int, Matrix1=Dense_host_int, Matrix2=Dense_host_int, tileDim=32]\"\n";
+    std::cout << "              /home/lukas/tnl-dev/src/UnitTests/Matrices/DenseMatrixTest.h(901): here\n";
+    std::cout << "                  instantiation of \"void test_GetMatrixProduct<Matrix>() [with Matrix=Dense_host_int]\"\n";
+    std::cout << "              /home/lukas/tnl-dev/src/UnitTests/Matrices/DenseMatrixTest.h(1315): here\n\n";
+}
+
+#ifdef HAVE_CUDA
+TEST( DenseMatrixTest, Dense_getMatrixProductTest_Cuda )
+{
+//    test_GetMatrixProduct< Dense_cuda_int >();
+    bool testRan = false;
+    EXPECT_TRUE( testRan );
+    std::cout << "\nTEST DID NOT RUN. NOT WORKING.\n\n";
+    std::cout << "If launched on GPU, this test will not build, but will print the following message: \n";
+    std::cout << "      /home/lukas/tnl-dev/src/TNL/Matrices/Dense_impl.h(510): error: identifier \"tnlCudaMin\" is undefined\n";
+    std::cout << "          detected during:\n";
+    std::cout << "              instantiation of \"void TNL::Matrices::DenseMatrixProductKernel<Real,Index,Matrix1,Matrix2,tileDim,tileRowBlockSize>(TNL::Matrices::Dense<Real, TNL::Devices::Cuda, Index> *, const Matrix1 *, const Matrix2 *, Real, Real, Index, Index) [with Real=int, Index=int, Matrix1=Dense_cuda_int, Matrix2=Dense_cuda_int, tileDim=32, tileRowBlockSize=8]\"\n";
+    std::cout << "              instantiation of \"void TNL::Matrices::Dense<Real, Device, Index>::getMatrixProduct(const Matrix1 &, const Matrix2 &, const TNL::Matrices::Dense<Real, Device, Index>::RealType &, const TNL::Matrices::Dense<Real, Device, Index>::RealType &) [with Real=int, Device=TNL::Devices::Cuda, Index=int, Matrix1=Dense_cuda_int, Matrix2=Dense_cuda_int, tileDim=32]\"\n";
+    std::cout << "              /home/lukas/tnl-dev/src/UnitTests/Matrices/DenseMatrixTest.h(901): here\n";
+    std::cout << "                  instantiation of \"void test_GetMatrixProduct<Matrix>() [with Matrix=Dense_cuda_int]\"\n";
+    std::cout << "              /home/lukas/tnl-dev/src/UnitTests/Matrices/DenseMatrixTest.h(1332): here\n\n";
+}
+#endif
+
+TEST( DenseMatrixTest, Dense_getTranspositionTest_Host )
+{
+//    test_GetTransposition< Dense_host_int >();
+    bool testRan = false;
+    EXPECT_TRUE( testRan );
+    std::cout << "\nTEST DID NOT RUN. NOT WORKING.\n\n";
+    std::cout << "If launched on CPU, this test will not build, but will print the following message: \n";
+    std::cout << "      /home/lukas/tnl-dev/src/TNL/Matrices/Dense_impl.h(836): error: no instance of function template \"TNL::Matrices::DenseTranspositionAlignedKernel\" matches the argument list\n";
+    std::cout << "              argument types are: (TNL::Matrices::Dense<int, TNL::Devices::Host, int> *, Dense_host_int *, const int, int, int)\n";
+    std::cout << "          detected during:\n";
+    std::cout << "              instantiation of \"void TNL::Matrices::Dense<Real, Device, Index>::getTransposition(const Matrix &, const TNL::Matrices::Dense<Real, Device, Index>::RealType &) [with Real=int, Device=TNL::Devices::Host, Index=int, Matrix=Dense_host_int, tileDim=32]\"\n";
+    std::cout << "              /home/lukas/tnl-dev/src/UnitTests/Matrices/DenseMatrixTest.h(977): here\n";
+    std::cout << "                  instantiation of \"void test_GetTransposition<Matrix>() [with Matrix=Dense_host_int]\"\n";
+    std::cout << "              /home/lukas/tnl-dev/src/UnitTests/Matrices/DenseMatrixTest.h(1420): here\n\n";
+    std::cout << "AND this message: \n";
+    std::cout << "      /home/lukas/tnl-dev/src/TNL/Matrices/Dense_impl.h(852): error: no instance of function template \"TNL::Matrices::DenseTranspositionNonAlignedKernel\" matches the argument list\n";
+    std::cout << "              argument types are: (TNL::Matrices::Dense<int, TNL::Devices::Host, int> *, Dense_host_int *, const int, int, int)\n";
+    std::cout << "          detected during:\n";
+    std::cout << "              instantiation of \"void TNL::Matrices::Dense<Real, Device, Index>::getTransposition(const Matrix &, const TNL::Matrices::Dense<Real, Device, Index>::RealType &) [with Real=int, Device=TNL::Devices::Host, Index=int, Matrix=Dense_host_int, tileDim=32]\"\n";
+    std::cout << "              /home/lukas/tnl-dev/src/UnitTests/Matrices/DenseMatrixTest.h(977): here\n";
+    std::cout << "                  instantiation of \"void test_GetTransposition<Matrix>() [with Matrix=Dense_host_int]\"\n";
+    std::cout << "              /home/lukas/tnl-dev/src/UnitTests/Matrices/DenseMatrixTest.h(1420): here\n\n";
+}
+
+#ifdef HAVE_CUDA
+TEST( DenseMatrixTest, Dense_getTranspositionTest_Cuda )
+{
+//    test_GetTransposition< Dense_cuda_int >();
+    bool testRan = false;
+    EXPECT_TRUE( testRan );
+    std::cout << "\nTEST DID NOT RUN. NOT WORKING.\n\n";
+    std::cout << "If launched on GPU, this test throws the following message: \n";
+    std::cout << "  Assertion 'row >= 0 && row < this->getRows() && column >= 0 && column < this->getColumns()' failed !!!\n";
+    std::cout << "      File: /home/lukas/tnl-dev/src/TNL/Matrices/Dense_impl.h \n";
+    std::cout << "      Line: 329 \n";
+    std::cout << "      Diagnostics: Not supported with CUDA.\n";
+    std::cout << "  Assertion 'row >= 0 && row < this->getRows() && column >= 0 && column < this->getColumns()' failed !!! \n";
+    std::cout << "      File: /home/lukas/tnl-dev/src/TNL/Matrices/Dense_impl.h \n";
+    std::cout << "      Line: 329 \n";
+    std::cout << "      Diagnostics: Not supported with CUDA.\n";
+    std::cout << "  Assertion 'row >= 0 && row < this->getRows() && column >= 0 && column < this->getColumns()' failed !!! \n";
+    std::cout << "      File: /home/lukas/tnl-dev/src/TNL/Matrices/Dense_impl.h \n";
+    std::cout << "      Line: 329 \n";
+    std::cout << "      Diagnostics: Not supported with CUDA.\n";
+    std::cout << "  Assertion 'row >= 0 && row < this->getRows() && column >= 0 && column < this->getColumns()' failed !!! \n";
+    std::cout << "      File: /home/lukas/tnl-dev/src/TNL/Matrices/Dense_impl.h \n";
+    std::cout << "      Line: 329 \n";
+    std::cout << "      Diagnostics: Not supported with CUDA.\n";
+    std::cout << "  terminate called after throwing an instance of 'TNL::Exceptions::CudaRuntimeError'\n";
+    std::cout << "          what():  CUDA ERROR 4 (cudaErrorLaunchFailure): unspecified launch failure.\n";
+    std::cout << "  Source: line 57 in /home/lukas/tnl-dev/src/TNL/Containers/Algorithms/ArrayOperationsCuda_impl.h: unspecified launch failure\n";
+    std::cout << "  [1]    4003 abort (core dumped)  ./DenseMatrixTest-dbg\n";
+}
+#endif
+
+TEST( DenseMatrixTest, Dense_performSORIterationTest_Host )
+{
+    test_PerformSORIteration< Dense_host_float >();
+}
+
+#ifdef HAVE_CUDA
+TEST( DenseMatrixTest, Dense_performSORIterationTest_Cuda )
+{
+//    test_PerformSORIteration< Dense_cuda_float >();
+    bool testRan = false;
+    EXPECT_TRUE( testRan );
+    std::cout << "\nTEST DID NOT RUN. NOT WORKING.\n\n";
+    std::cout << "If launched, this test throws the following message: \n";
+    std::cout << "      [1]    6992 segmentation fault (core dumped)  ./SparseMatrixTest-dbg\n\n";
+    std::cout << "\n THIS IS NOT IMPLEMENTED FOR CUDA YET!!\n\n";
+}
+#endif
+
+#endif
+
+#include "../GtestMissingError.h"
+int main( int argc, char* argv[] )
+{
+#ifdef HAVE_GTEST
+   ::testing::InitGoogleTest( &argc, argv );
+   return RUN_ALL_TESTS();
+#else
+   throw GtestMissingError();
+#endif
+}
+
diff --git a/src/UnitTests/Matrices/SparseMatrixCopyTest.h b/src/UnitTests/Matrices/SparseMatrixCopyTest.h
index a11a8b4442527b371603b8f5d43b70b8e4ff558d..2885bac093d8a3e59144535baf892f27ef0cc24e 100644
--- a/src/UnitTests/Matrices/SparseMatrixCopyTest.h
+++ b/src/UnitTests/Matrices/SparseMatrixCopyTest.h
@@ -22,6 +22,233 @@ using SE_cuda = TNL::Matrices::SlicedEllpack< int, TNL::Devices::Cuda, int, 2 >;
 #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:
  *
@@ -34,7 +261,7 @@ using SE_cuda = TNL::Matrices::SlicedEllpack< int, TNL::Devices::Cuda, int, 2 >;
  *    \               17 /
  */
 template< typename Matrix >
-void setupMatrix( Matrix& m )
+void setupTriDiagMatrix( Matrix& m )
 {
    const int rows = 7;
    const int cols = 6;
@@ -55,7 +282,7 @@ void setupMatrix( Matrix& m )
 }
 
 template< typename Matrix >
-void checkMatrix( Matrix& m )
+void checkTriDiagMatrix( Matrix& m )
 {
    ASSERT_EQ( m.getRows(), 7 );
    ASSERT_EQ( m.getColumns(), 6 );
@@ -113,25 +340,80 @@ void checkMatrix( Matrix& m )
 template< typename Matrix1, typename Matrix2 >
 void testCopyAssignment()
 {
-   Matrix1 m1;
-   setupMatrix( m1 );
-   checkMatrix( m1 );
-
-   Matrix2 m2;
-   m2 = m1;
-   checkMatrix( m2 );
+   {
+        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()
 {
-   Matrix1 m1;
-   setupMatrix( m1 );
-   checkMatrix( m1 );
-
-   Matrix2 m2;
-   TNL::Matrices::copySparseMatrix( m2, m1 );
-   checkMatrix( m2 );
+    
+   {
+        SCOPED_TRACE("Tri Diagonal Matrix");
+        
+        Matrix1 triDiag1;
+        setupTriDiagMatrix( triDiag1 );
+        checkTriDiagMatrix( triDiag1 );
+        
+        Matrix2 triDiag2;
+        TNL::Matrices::copySparseMatrix( triDiag2, triDiag1 );
+        checkTriDiagMatrix( triDiag2 );
+   }
+   
+   {
+        SCOPED_TRACE("Anti Tri Diagonal Matrix");
+                
+        Matrix1 antiTriDiag1;
+        setupAntiTriDiagMatrix( antiTriDiag1 );
+        checkAntiTriDiagMatrix( antiTriDiag1 );
+        
+        Matrix2 antiTriDiag2;
+        TNL::Matrices::copySparseMatrix( antiTriDiag2, antiTriDiag1 );
+        checkAntiTriDiagMatrix( antiTriDiag2 );
+   }
+   
+   {
+        SCOPED_TRACE("Uneven Row Size Matrix");
+        Matrix1 unevenRowSize1;
+        setupUnevenRowSizeMatrix( unevenRowSize1 );
+        checkUnevenRowSizeMatrix( unevenRowSize1 );
+        
+        Matrix2 unevenRowSize2;
+        TNL::Matrices::copySparseMatrix( unevenRowSize2, unevenRowSize1 );
+        checkUnevenRowSizeMatrix( unevenRowSize2 );
+   }
 }
 
 
diff --git a/src/UnitTests/Matrices/SparseMatrixTest.cpp b/src/UnitTests/Matrices/SparseMatrixTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..46f6b9bd3531652feb15f10d60e5736ddfb627fa
--- /dev/null
+++ b/src/UnitTests/Matrices/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/SparseMatrixTest.cu b/src/UnitTests/Matrices/SparseMatrixTest.cu
new file mode 100644
index 0000000000000000000000000000000000000000..01c23c1937b6ea39c2c99647207d74297ea588c9
--- /dev/null
+++ b/src/UnitTests/Matrices/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/SparseMatrixTest.h b/src/UnitTests/Matrices/SparseMatrixTest.h
new file mode 100644
index 0000000000000000000000000000000000000000..1145570b2df29ef8de343f7ca17f161d3a4cf0db
--- /dev/null
+++ b/src/UnitTests/Matrices/SparseMatrixTest.h
@@ -0,0 +1,844 @@
+/***************************************************************************
+                          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 */
+
+// TODO
+/*
+ * getType()                        ::HOW?  How to test this for each format? edit string how?
+ *      Found the mistake for Cuda instead of Devices::Cuda. Incorrect String in src/TNL/Devices/Cuda.cpp
+ *      MISSING: indexType is missing in CSR_impl.h
+ * getTypeVirtual()                 ::TEST? This just calls getType().
+ * getSerializationType()           ::TEST? This just calls HostType::getType().
+ * getSerializationTypeVirtual()    ::TEST? This just calls getSerializationType().
+ * setDimensions()                      ::DONE
+ * setCompressedRowLengths()            ::DONE
+ * getRowLength()                   ::USED! In test_SetCompressedRowLengths() to verify the test itself.
+ * getRowLengthFast()               ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * setLike()                            ::DONE
+ * reset()                              ::DONE
+ * setElementFast()                 ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * setElement()                         ::DONE
+ * addElementFast()                 ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * addElement()                         ::DONE
+ * setRowFast()                     ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * setRow()                             ::DONE
+ *      MISTAKE!!! In SlicedEllpack: addElement(), line 263, "column <= this->rows" shouldn't it be: "column <= this->columns", otherwise test_SetRow causes the assertion to fail.
+ * addRowFast()                     ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * addRow()                         ::NOT IMPLEMENTED! This calls addRowFast() which isn't implemented. Implement? Is it supposed to add an extra row to the matrix or add elements of a row to another row in the matrix?
+ * getElementFast()                 ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * getElement()                     ::USED! In test_SetElement(), test_AddElement() and test_setRow() to verify the test itself.
+ * getRowFast()                     ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * MatrixRow getRow()               ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * ConstMatrixRow getRow()          ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * rowVectorProduct()               ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * vectorProduct()                      ::DONE
+ *      This used to throw illegal memory access, but instead of using ints for vectors, using Types, helped.
+ * addMatrix()                      ::NOT IMPLEMENTED!
+ * getTransposition()               ::NOT IMPLMENETED!
+ * performSORIteration()            ::HOW? Throws segmentation fault CUDA.
+ * operator=()                      ::HOW? What is this supposed to enable? Overloading operators?
+ * save( File& file)                ::USED! In save( String& fileName )
+ * load( File& file )               ::USED! In load( String& fileName )
+ * save( String& fileName )             ::DONE
+ * load( String& fileName )             ::DONE
+ * print()                              ::DONE
+ * setCudaKernelType()              ::NOT SUPPOSED TO TEST! via notes from 1.11.2018 supervisor meeting.
+ * getCudaKernelType()              ::NOT SUPPOSED TO TEST! via notes from 1.11.2018 supervisor meeting.
+ * setCudaWarpSize()                ::NOT SUPPOSED TO TEST! via notes from 1.11.2018 supervisor meeting.
+ * getCudaWarpSize()                ::NOT SUPPOSED TO TEST! via notes from 1.11.2018 supervisor meeting.
+ * setHybridModeSplit()             ::NOT SUPPOSED TO TEST! via notes from 1.11.2018 supervisor meeting.
+ * getHybridModeSplit()             ::NOT SUPPOSED TO TEST! via notes from 1.11.2018 supervisor meeting.
+ * spmvCudaVectorized()             ::TEST? How to test __device__?
+ * vectorProductCuda()              ::TEST? How to test __device__?
+ */
+
+// GENERAL TODO
+/*
+ * For every function, EXPECT_EQ needs to be done, even for zeros in matrices.
+ * Figure out __cuda_callable_. When trying to call __cuda_callable__ functions
+ *      a segmentation fault (core dumped) is thrown.
+ *  ==>__cuda_callable__ works only for CPU at the moment. (for loops vs thread kernel assignment)
+ */
+
+#include <TNL/Matrices/CSR.h>
+#include <TNL/Matrices/Ellpack.h>
+#include <TNL/Matrices/SlicedEllpack.h>
+
+#include <TNL/Matrices/ChunkedEllpack.h>
+#include <TNL/Matrices/AdEllpack.h>
+#include <TNL/Matrices/BiEllpack.h>
+
+#include <UnitTests/Matrices/SparseMatrixTest_impl.h>
+#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 >;
+
+#ifdef NOT_WORKING
+// test fixture for typed tests
+template< typename Matrix >
+class AdEllpackMatrixTest : public ::testing::Test
+{
+protected:
+   using AdEllpackMatrixType = Matrix;
+};
+
+// types for which MatrixTest is instantiated
+using AdEllpackMatrixTypes = ::testing::Types
+<
+    TNL::Matrices::AdEllpack< int,    TNL::Devices::Host, short >,
+    TNL::Matrices::AdEllpack< long,   TNL::Devices::Host, short >,
+    TNL::Matrices::AdEllpack< float,  TNL::Devices::Host, short >,
+    TNL::Matrices::AdEllpack< double, TNL::Devices::Host, short >,
+    TNL::Matrices::AdEllpack< int,    TNL::Devices::Host, int >,
+    TNL::Matrices::AdEllpack< long,   TNL::Devices::Host, int >,
+    TNL::Matrices::AdEllpack< float,  TNL::Devices::Host, int >,
+    TNL::Matrices::AdEllpack< double, TNL::Devices::Host, int >,
+    TNL::Matrices::AdEllpack< int,    TNL::Devices::Host, long >,
+    TNL::Matrices::AdEllpack< long,   TNL::Devices::Host, long >,
+    TNL::Matrices::AdEllpack< float,  TNL::Devices::Host, long >,
+    TNL::Matrices::AdEllpack< double, TNL::Devices::Host, long >,
+#ifdef HAVE_CUDA
+    TNL::Matrices::AdEllpack< int,    TNL::Devices::Cuda, short >,
+    TNL::Matrices::AdEllpack< long,   TNL::Devices::Cuda, short >,
+    TNL::Matrices::AdEllpack< float,  TNL::Devices::Cuda, short >,
+    TNL::Matrices::AdEllpack< double, TNL::Devices::Cuda, short >,
+    TNL::Matrices::AdEllpack< int,    TNL::Devices::Cuda, int >,
+    TNL::Matrices::AdEllpack< long,   TNL::Devices::Cuda, int >,
+    TNL::Matrices::AdEllpack< float,  TNL::Devices::Cuda, int >,
+    TNL::Matrices::AdEllpack< double, TNL::Devices::Cuda, int >,
+    TNL::Matrices::AdEllpack< int,    TNL::Devices::Cuda, long >,
+    TNL::Matrices::AdEllpack< long,   TNL::Devices::Cuda, long >,
+    TNL::Matrices::AdEllpack< float,  TNL::Devices::Cuda, long >,
+    TNL::Matrices::AdEllpack< double, TNL::Devices::Cuda, long >
+#endif
+>;
+
+TYPED_TEST_CASE( AdEllpackMatrixTest, AdEllpackMatrixTypes);
+
+TYPED_TEST( AdEllpackMatrixTest, setDimensionsTest )
+{
+    using AdEllpackMatrixType = typename TestFixture::AdEllpackMatrixType;
+    
+    test_SetDimensions< AdEllpackMatrixType >();
+}
+
+TYPED_TEST( AdEllpackMatrixTest, setCompressedRowLengthsTest )
+{
+//    using AdEllpackMatrixType = typename TestFixture::AdEllpackMatrixType;
+    
+//    test_SetCompressedRowLengths< AdEllpackMatrixType >();
+    
+    bool testRan = false;
+    EXPECT_TRUE( testRan );
+    std::cout << "\nTEST DID NOT RUN. NOT WORKING.\n\n";
+    std::cout << "      This test is dependent on the input format. \n";
+    std::cout << "      Almost every format allocates elements per row differently.\n\n";
+    std::cout << "\n    TODO: Finish implementation of getNonZeroRowLength (Only non-zero elements, not the number of allocated elements.)\n\n";
+}
+
+TYPED_TEST( AdEllpackMatrixTest, setLikeTest )
+{
+    using AdEllpackMatrixType = typename TestFixture::AdEllpackMatrixType;
+    
+    test_SetLike< AdEllpackMatrixType, AdEllpackMatrixType >();
+}
+
+TYPED_TEST( AdEllpackMatrixTest, resetTest )
+{
+    using AdEllpackMatrixType = typename TestFixture::AdEllpackMatrixType;
+    
+    test_Reset< AdEllpackMatrixType >();
+}
+
+TYPED_TEST( AdEllpackMatrixTest, setElementTest )
+{
+    using AdEllpackMatrixType = typename TestFixture::AdEllpackMatrixType;
+    
+    test_SetElement< AdEllpackMatrixType >();
+}
+
+TYPED_TEST( AdEllpackMatrixTest, addElementTest )
+{
+    using AdEllpackMatrixType = typename TestFixture::AdEllpackMatrixType;
+    
+    test_AddElement< AdEllpackMatrixType >();
+}
+
+TYPED_TEST( AdEllpackMatrixTest, setRowTest )
+{
+    using AdEllpackMatrixType = typename TestFixture::AdEllpackMatrixType;
+    
+    test_SetRow< AdEllpackMatrixType >();
+}
+
+TYPED_TEST( AdEllpackMatrixTest, vectorProductTest )
+{
+    using AdEllpackMatrixType = typename TestFixture::AdEllpackMatrixType;
+    
+    test_VectorProduct< AdEllpackMatrixType >();
+}
+
+TYPED_TEST( AdEllpackMatrixTest, saveAndLoadTest )
+{
+    using AdEllpackMatrixType = typename TestFixture::AdEllpackMatrixType;
+    
+    test_SaveAndLoad< AdEllpackMatrixType >();
+}
+
+TYPED_TEST( AdEllpackMatrixTest, printTest )
+{
+    using AdEllpackMatrixType = typename TestFixture::AdEllpackMatrixType;
+    
+    test_Print< AdEllpackMatrixType >();
+}
+
+// test fixture for typed tests
+template< typename Matrix >
+class BiEllpackMatrixTest : public ::testing::Test
+{
+protected:
+   using BiEllpackMatrixType = Matrix;
+};
+
+// types for which MatrixTest is instantiated
+using BiEllpackMatrixTypes = ::testing::Types
+<
+    TNL::Matrices::BiEllpack< int,    TNL::Devices::Host, short >,
+    TNL::Matrices::BiEllpack< long,   TNL::Devices::Host, short >,
+    TNL::Matrices::BiEllpack< float,  TNL::Devices::Host, short >,
+    TNL::Matrices::BiEllpack< double, TNL::Devices::Host, short >,
+    TNL::Matrices::BiEllpack< int,    TNL::Devices::Host, int >,
+    TNL::Matrices::BiEllpack< long,   TNL::Devices::Host, int >,
+    TNL::Matrices::BiEllpack< float,  TNL::Devices::Host, int >,
+    TNL::Matrices::BiEllpack< double, TNL::Devices::Host, int >,
+    TNL::Matrices::BiEllpack< int,    TNL::Devices::Host, long >,
+    TNL::Matrices::BiEllpack< long,   TNL::Devices::Host, long >,
+    TNL::Matrices::BiEllpack< float,  TNL::Devices::Host, long >,
+    TNL::Matrices::BiEllpack< double, TNL::Devices::Host, long >//,
+//#ifdef HAVE_CUDA
+//    TNL::Matrices::BiEllpack< int,    TNL::Devices::Cuda, short >,
+//    TNL::Matrices::BiEllpack< long,   TNL::Devices::Cuda, short >,
+//    TNL::Matrices::BiEllpack< float,  TNL::Devices::Cuda, short >,
+//    TNL::Matrices::BiEllpack< double, TNL::Devices::Cuda, short >,
+//    TNL::Matrices::BiEllpack< int,    TNL::Devices::Cuda, int >,
+//    TNL::Matrices::BiEllpack< long,   TNL::Devices::Cuda, int >,
+//    TNL::Matrices::BiEllpack< float,  TNL::Devices::Cuda, int >,
+//    TNL::Matrices::BiEllpack< double, TNL::Devices::Cuda, int >,
+//    TNL::Matrices::BiEllpack< int,    TNL::Devices::Cuda, long >,
+//    TNL::Matrices::BiEllpack< long,   TNL::Devices::Cuda, long >,
+//    TNL::Matrices::BiEllpack< float,  TNL::Devices::Cuda, long >,
+//    TNL::Matrices::BiEllpack< double, TNL::Devices::Cuda, long >
+//#endif
+>;
+
+TYPED_TEST_CASE( BiEllpackMatrixTest, BiEllpackMatrixTypes);
+
+TYPED_TEST( BiEllpackMatrixTest, setDimensionsTest )
+{
+    using BiEllpackMatrixType = typename TestFixture::BiEllpackMatrixType;
+    
+    test_SetDimensions< BiEllpackMatrixType >();
+}
+
+TYPED_TEST( BiEllpackMatrixTest, setCompressedRowLengthsTest )
+{
+//    using BiEllpackMatrixType = typename TestFixture::BiEllpackMatrixType;
+    
+//    test_SetCompressedRowLengths< BiEllpackMatrixType >();
+    
+    bool testRan = false;
+    EXPECT_TRUE( testRan );
+    std::cout << "\nTEST DID NOT RUN. NOT WORKING.\n\n";
+    std::cout << "      This test is dependent on the input format. \n";
+    std::cout << "      Almost every format allocates elements per row differently.\n\n";
+    std::cout << "\n    TODO: Finish implementation of getNonZeroRowLength (Only non-zero elements, not the number of allocated elements.)\n\n";
+}
+
+TYPED_TEST( BiEllpackMatrixTest, setLikeTest )
+{
+    using BiEllpackMatrixType = typename TestFixture::BiEllpackMatrixType;
+    
+    test_SetLike< BiEllpackMatrixType, BiEllpackMatrixType >();
+}
+
+TYPED_TEST( BiEllpackMatrixTest, resetTest )
+{
+    using BiEllpackMatrixType = typename TestFixture::BiEllpackMatrixType;
+    
+    test_Reset< BiEllpackMatrixType >();
+}
+
+TYPED_TEST( BiEllpackMatrixTest, setElementTest )
+{
+    using BiEllpackMatrixType = typename TestFixture::BiEllpackMatrixType;
+    
+    test_SetElement< BiEllpackMatrixType >();
+}
+
+TYPED_TEST( BiEllpackMatrixTest, addElementTest )
+{
+    using BiEllpackMatrixType = typename TestFixture::BiEllpackMatrixType;
+    
+    test_AddElement< BiEllpackMatrixType >();
+}
+
+TYPED_TEST( BiEllpackMatrixTest, setRowTest )
+{
+    using BiEllpackMatrixType = typename TestFixture::BiEllpackMatrixType;
+    
+    test_SetRow< BiEllpackMatrixType >();
+}
+
+TYPED_TEST( BiEllpackMatrixTest, vectorProductTest )
+{
+    using BiEllpackMatrixType = typename TestFixture::BiEllpackMatrixType;
+    
+    test_VectorProduct< BiEllpackMatrixType >();
+}
+
+TYPED_TEST( BiEllpackMatrixTest, saveAndLoadTest )
+{
+    using BiEllpackMatrixType = typename TestFixture::BiEllpackMatrixType;
+    
+    test_SaveAndLoad< BiEllpackMatrixType >();
+}
+
+TYPED_TEST( BiEllpackMatrixTest, printTest )
+{
+    using BiEllpackMatrixType = typename TestFixture::BiEllpackMatrixType;
+    
+    test_Print< BiEllpackMatrixType >();
+}
+
+#endif
+
+// GTEST ::testing::Types<> has a limit of 38.
+
+// test fixture for typed tests
+template< typename Matrix >
+class ChunkedEllpackMatrixTest : public ::testing::Test
+{
+protected:
+   using ChunkedEllpackMatrixType = Matrix;
+};
+
+// columnIndexes of ChunkedEllpack appear to be broken, when printed, it prints out a bunch of 4s.
+// rowPointers have interesting elements? 0 18 36 42 54 72 96 126 162 204 256 when rows = 10, cols = 11; rowLengths = 3 3 1 2 3 4 5 6 7 8
+// and 0 52 103 154 205 256 when rows = 5, cols = 4; rowLengths = 3 3 3 3 3
+
+
+// types for which MatrixTest is instantiated
+using ChEllpackMatrixTypes = ::testing::Types
+<
+    TNL::Matrices::ChunkedEllpack< int,    TNL::Devices::Host, short >,
+    TNL::Matrices::ChunkedEllpack< long,   TNL::Devices::Host, short >,
+    TNL::Matrices::ChunkedEllpack< float,  TNL::Devices::Host, short >,
+    TNL::Matrices::ChunkedEllpack< double, TNL::Devices::Host, short >,
+    TNL::Matrices::ChunkedEllpack< int,    TNL::Devices::Host, int >,
+    TNL::Matrices::ChunkedEllpack< long,   TNL::Devices::Host, int >,
+    TNL::Matrices::ChunkedEllpack< float,  TNL::Devices::Host, int >,
+    TNL::Matrices::ChunkedEllpack< double, TNL::Devices::Host, int >,
+    TNL::Matrices::ChunkedEllpack< int,    TNL::Devices::Host, long >,
+    TNL::Matrices::ChunkedEllpack< long,   TNL::Devices::Host, long >,
+    TNL::Matrices::ChunkedEllpack< float,  TNL::Devices::Host, long >,
+    TNL::Matrices::ChunkedEllpack< double, TNL::Devices::Host, long >
+#ifdef HAVE_CUDA
+    ,TNL::Matrices::ChunkedEllpack< int,    TNL::Devices::Cuda, short >,
+    TNL::Matrices::ChunkedEllpack< long,   TNL::Devices::Cuda, short >,
+    TNL::Matrices::ChunkedEllpack< float,  TNL::Devices::Cuda, short >,
+    TNL::Matrices::ChunkedEllpack< double, TNL::Devices::Cuda, short >,
+    TNL::Matrices::ChunkedEllpack< int,    TNL::Devices::Cuda, int >,
+    TNL::Matrices::ChunkedEllpack< long,   TNL::Devices::Cuda, int >,
+    TNL::Matrices::ChunkedEllpack< float,  TNL::Devices::Cuda, int >,
+    TNL::Matrices::ChunkedEllpack< double, TNL::Devices::Cuda, int >,
+    TNL::Matrices::ChunkedEllpack< int,    TNL::Devices::Cuda, long >,
+    TNL::Matrices::ChunkedEllpack< long,   TNL::Devices::Cuda, long >,
+    TNL::Matrices::ChunkedEllpack< float,  TNL::Devices::Cuda, long >,
+    TNL::Matrices::ChunkedEllpack< double, TNL::Devices::Cuda, long >
+#endif
+>;
+
+TYPED_TEST_CASE( ChunkedEllpackMatrixTest, ChEllpackMatrixTypes);
+
+TYPED_TEST( ChunkedEllpackMatrixTest, setDimensionsTest )
+{
+    using ChunkedEllpackMatrixType = typename TestFixture::ChunkedEllpackMatrixType;
+    
+    test_SetDimensions< ChunkedEllpackMatrixType >();
+}
+
+//TYPED_TEST( ChunkedEllpackMatrixTest, setCompressedRowLengthsTest )
+//{
+////    using ChunkedEllpackMatrixType = typename TestFixture::ChunkedEllpackMatrixType;
+//    
+////    test_SetCompressedRowLengths< ChunkedEllpackMatrixType >();
+//    
+//    bool testRan = false;
+//    EXPECT_TRUE( testRan );
+//    std::cout << "\nTEST DID NOT RUN. NOT WORKING.\n\n";
+//    std::cout << "      This test is dependent on the input format. \n";
+//    std::cout << "      Almost every format allocates elements per row differently.\n\n";
+//    std::cout << "\n    TODO: Finish implementation of getNonZeroRowLength (Only non-zero elements, not the number of allocated elements.)\n\n";
+//}
+
+TYPED_TEST( ChunkedEllpackMatrixTest, setLikeTest )
+{
+    using ChunkedEllpackMatrixType = typename TestFixture::ChunkedEllpackMatrixType;
+    
+    test_SetLike< ChunkedEllpackMatrixType, ChunkedEllpackMatrixType >();
+}
+
+TYPED_TEST( ChunkedEllpackMatrixTest, resetTest )
+{
+    using ChunkedEllpackMatrixType = typename TestFixture::ChunkedEllpackMatrixType;
+    
+    test_Reset< ChunkedEllpackMatrixType >();
+}
+
+TYPED_TEST( ChunkedEllpackMatrixTest, setElementTest )
+{
+    using ChunkedEllpackMatrixType = typename TestFixture::ChunkedEllpackMatrixType;
+    
+    test_SetElement< ChunkedEllpackMatrixType >();
+}
+
+TYPED_TEST( ChunkedEllpackMatrixTest, addElementTest )
+{
+    using ChunkedEllpackMatrixType = typename TestFixture::ChunkedEllpackMatrixType;
+    
+    test_AddElement< ChunkedEllpackMatrixType >();
+}
+
+TYPED_TEST( ChunkedEllpackMatrixTest, setRowTest )
+{
+    using ChunkedEllpackMatrixType = typename TestFixture::ChunkedEllpackMatrixType;
+    
+    test_SetRow< ChunkedEllpackMatrixType >();
+}
+
+TYPED_TEST( ChunkedEllpackMatrixTest, vectorProductTest )
+{
+    using ChunkedEllpackMatrixType = typename TestFixture::ChunkedEllpackMatrixType;
+    
+    test_VectorProduct< ChunkedEllpackMatrixType >();
+}
+
+TYPED_TEST( ChunkedEllpackMatrixTest, saveAndLoadTest )
+{
+    using ChunkedEllpackMatrixType = typename TestFixture::ChunkedEllpackMatrixType;
+    
+    test_SaveAndLoad< ChunkedEllpackMatrixType >();
+}
+
+TYPED_TEST( ChunkedEllpackMatrixTest, printTest )
+{
+    using ChunkedEllpackMatrixType = typename TestFixture::ChunkedEllpackMatrixType;
+    
+    test_Print< ChunkedEllpackMatrixType >();
+}
+
+// test fixture for typed tests
+template< typename Matrix >
+class CSRMatrixTest : public ::testing::Test
+{
+protected:
+   using CSRMatrixType = Matrix;
+};
+
+// types for which MatrixTest is instantiated
+using CSRMatrixTypes = ::testing::Types
+<
+    TNL::Matrices::CSR< int,    TNL::Devices::Host, short >,
+    TNL::Matrices::CSR< long,   TNL::Devices::Host, short >,
+    TNL::Matrices::CSR< float,  TNL::Devices::Host, short >,
+    TNL::Matrices::CSR< double, TNL::Devices::Host, short >,
+    TNL::Matrices::CSR< int,    TNL::Devices::Host, int >,
+    TNL::Matrices::CSR< long,   TNL::Devices::Host, int >,
+    TNL::Matrices::CSR< float,  TNL::Devices::Host, int >,
+    TNL::Matrices::CSR< double, TNL::Devices::Host, int >,
+    TNL::Matrices::CSR< int,    TNL::Devices::Host, long >,
+    TNL::Matrices::CSR< long,   TNL::Devices::Host, long >,
+    TNL::Matrices::CSR< float,  TNL::Devices::Host, long >,
+    TNL::Matrices::CSR< double, TNL::Devices::Host, long >
+#ifdef HAVE_CUDA
+    ,TNL::Matrices::CSR< int,    TNL::Devices::Cuda, short >,
+    TNL::Matrices::CSR< long,   TNL::Devices::Cuda, short >,
+    TNL::Matrices::CSR< float,  TNL::Devices::Cuda, short >,
+    TNL::Matrices::CSR< double, TNL::Devices::Cuda, short >,
+    TNL::Matrices::CSR< int,    TNL::Devices::Cuda, int >,
+    TNL::Matrices::CSR< long,   TNL::Devices::Cuda, int >,
+    TNL::Matrices::CSR< float,  TNL::Devices::Cuda, int >,
+    TNL::Matrices::CSR< double, TNL::Devices::Cuda, int >,
+    TNL::Matrices::CSR< int,    TNL::Devices::Cuda, long >,
+    TNL::Matrices::CSR< long,   TNL::Devices::Cuda, long >,
+    TNL::Matrices::CSR< float,  TNL::Devices::Cuda, long >,
+    TNL::Matrices::CSR< double, TNL::Devices::Cuda, long >
+#endif
+>;
+
+TYPED_TEST_CASE( CSRMatrixTest, CSRMatrixTypes);
+
+TYPED_TEST( CSRMatrixTest, setDimensionsTest )
+{
+    using CSRMatrixType = typename TestFixture::CSRMatrixType;
+    
+    test_SetDimensions< CSRMatrixType >();
+}
+
+//TYPED_TEST( CSRMatrixTest, setCompressedRowLengthsTest )
+//{
+////    using CSRMatrixType = typename TestFixture::CSRMatrixType;
+//    
+////    test_SetCompressedRowLengths< CSRMatrixType >();
+//    
+//    bool testRan = false;
+//    EXPECT_TRUE( testRan );
+//    std::cout << "\nTEST DID NOT RUN. NOT WORKING.\n\n";
+//    std::cout << "      This test is dependent on the input format. \n";
+//    std::cout << "      Almost every format allocates elements per row differently.\n\n";
+//    std::cout << "\n    TODO: Finish implementation of getNonZeroRowLength (Only non-zero elements, not the number of allocated elements.)\n\n";
+//}
+
+TYPED_TEST( CSRMatrixTest, setLikeTest )
+{
+    using CSRMatrixType = typename TestFixture::CSRMatrixType;
+    
+    test_SetLike< CSRMatrixType, CSRMatrixType >();
+}
+
+TYPED_TEST( CSRMatrixTest, resetTest )
+{
+    using CSRMatrixType = typename TestFixture::CSRMatrixType;
+    
+    test_Reset< CSRMatrixType >();
+}
+
+TYPED_TEST( CSRMatrixTest, setElementTest )
+{
+    using CSRMatrixType = typename TestFixture::CSRMatrixType;
+    
+    test_SetElement< CSRMatrixType >();
+}
+
+TYPED_TEST( CSRMatrixTest, addElementTest )
+{
+    using CSRMatrixType = typename TestFixture::CSRMatrixType;
+    
+    test_AddElement< CSRMatrixType >();
+}
+
+TYPED_TEST( CSRMatrixTest, setRowTest )
+{
+    using CSRMatrixType = typename TestFixture::CSRMatrixType;
+    
+    test_SetRow< CSRMatrixType >();
+}
+
+TYPED_TEST( CSRMatrixTest, vectorProductTest )
+{
+    using CSRMatrixType = typename TestFixture::CSRMatrixType;
+    
+    test_VectorProduct< CSRMatrixType >();
+}
+
+TYPED_TEST( CSRMatrixTest, saveAndLoadTest )
+{
+    using CSRMatrixType = typename TestFixture::CSRMatrixType;
+    
+    test_SaveAndLoad< CSRMatrixType >();
+}
+
+TYPED_TEST( CSRMatrixTest, printTest )
+{
+    using CSRMatrixType = typename TestFixture::CSRMatrixType;
+    
+    test_Print< CSRMatrixType >();
+}
+
+// test fixture for typed tests
+template< typename Matrix >
+class EllpackMatrixTest : public ::testing::Test
+{
+protected:
+   using EllpackMatrixType = Matrix;
+};
+
+// types for which MatrixTest is instantiated
+using EllpackMatrixTypes = ::testing::Types
+<
+    TNL::Matrices::Ellpack< int,    TNL::Devices::Host, short >,
+    TNL::Matrices::Ellpack< long,   TNL::Devices::Host, short >,
+    TNL::Matrices::Ellpack< float,  TNL::Devices::Host, short >,
+    TNL::Matrices::Ellpack< double, TNL::Devices::Host, short >,
+    TNL::Matrices::Ellpack< int,    TNL::Devices::Host, int >,
+    TNL::Matrices::Ellpack< long,   TNL::Devices::Host, int >,
+    TNL::Matrices::Ellpack< float,  TNL::Devices::Host, int >,
+    TNL::Matrices::Ellpack< double, TNL::Devices::Host, int >,
+    TNL::Matrices::Ellpack< int,    TNL::Devices::Host, long >,
+    TNL::Matrices::Ellpack< long,   TNL::Devices::Host, long >,
+    TNL::Matrices::Ellpack< float,  TNL::Devices::Host, long >,
+    TNL::Matrices::Ellpack< double, TNL::Devices::Host, long >
+#ifdef HAVE_CUDA
+    ,TNL::Matrices::Ellpack< int,    TNL::Devices::Cuda, short >,
+    TNL::Matrices::Ellpack< long,   TNL::Devices::Cuda, short >,
+    TNL::Matrices::Ellpack< float,  TNL::Devices::Cuda, short >,
+    TNL::Matrices::Ellpack< double, TNL::Devices::Cuda, short >,
+    TNL::Matrices::Ellpack< int,    TNL::Devices::Cuda, int >,
+    TNL::Matrices::Ellpack< long,   TNL::Devices::Cuda, int >,
+    TNL::Matrices::Ellpack< float,  TNL::Devices::Cuda, int >,
+    TNL::Matrices::Ellpack< double, TNL::Devices::Cuda, int >,
+    TNL::Matrices::Ellpack< int,    TNL::Devices::Cuda, long >,
+    TNL::Matrices::Ellpack< long,   TNL::Devices::Cuda, long >,
+    TNL::Matrices::Ellpack< float,  TNL::Devices::Cuda, long >,
+    TNL::Matrices::Ellpack< double, TNL::Devices::Cuda, long >
+#endif
+>;
+
+TYPED_TEST_CASE( EllpackMatrixTest, EllpackMatrixTypes );
+
+TYPED_TEST( EllpackMatrixTest, setDimensionsTest )
+{
+    using EllpackMatrixType = typename TestFixture::EllpackMatrixType;
+    
+    test_SetDimensions< EllpackMatrixType >();
+}
+
+//TYPED_TEST( EllpackMatrixTest, setCompressedRowLengthsTest )
+//{
+////    using EllpackMatrixType = typename TestFixture::EllpackMatrixType;
+//    
+////    test_SetCompressedRowLengths< EllpackMatrixType >();
+//    
+//    bool testRan = false;
+//    EXPECT_TRUE( testRan );
+//    std::cout << "\nTEST DID NOT RUN. NOT WORKING.\n\n";
+//    std::cout << "      This test is dependent on the input format. \n";
+//    std::cout << "      Almost every format allocates elements per row differently.\n\n";
+//    std::cout << "\n    TODO: Finish implementation of getNonZeroRowLength (Only non-zero elements, not the number of allocated elements.)\n\n";
+//}
+
+TYPED_TEST( EllpackMatrixTest, setLikeTest )
+{
+    using EllpackMatrixType = typename TestFixture::EllpackMatrixType;
+    
+    test_SetLike< EllpackMatrixType, EllpackMatrixType >();
+}
+
+TYPED_TEST( EllpackMatrixTest, resetTest )
+{
+    using EllpackMatrixType = typename TestFixture::EllpackMatrixType;
+    
+    test_Reset< EllpackMatrixType >();
+}
+
+TYPED_TEST( EllpackMatrixTest, setElementTest )
+{
+    using EllpackMatrixType = typename TestFixture::EllpackMatrixType;
+    
+    test_SetElement< EllpackMatrixType >();
+}
+
+TYPED_TEST( EllpackMatrixTest, addElementTest )
+{
+    using EllpackMatrixType = typename TestFixture::EllpackMatrixType;
+    
+    test_AddElement< EllpackMatrixType >();
+}
+
+TYPED_TEST( EllpackMatrixTest, setRowTest )
+{
+    using EllpackMatrixType = typename TestFixture::EllpackMatrixType;
+    
+    test_SetRow< EllpackMatrixType >();
+}
+
+TYPED_TEST( EllpackMatrixTest, vectorProductTest )
+{
+    using EllpackMatrixType = typename TestFixture::EllpackMatrixType;
+    
+    test_VectorProduct< EllpackMatrixType >();
+}
+
+TYPED_TEST( EllpackMatrixTest, saveAndLoadTest )
+{
+    using EllpackMatrixType = typename TestFixture::EllpackMatrixType;
+    
+    test_SaveAndLoad< EllpackMatrixType >();
+}
+
+TYPED_TEST( EllpackMatrixTest, printTest )
+{
+    using EllpackMatrixType = typename TestFixture::EllpackMatrixType;
+    
+    test_Print< EllpackMatrixType >();
+}
+
+// test fixture for typed tests
+template< typename Matrix >
+class SlicedEllpackMatrixTest : public ::testing::Test
+{
+protected:
+   using SlicedEllpackMatrixType = Matrix;
+};
+
+// types for which MatrixTest is instantiated
+using SlicedEllpackMatrixTypes = ::testing::Types
+<
+    TNL::Matrices::SlicedEllpack< int,    TNL::Devices::Host, short >,
+    TNL::Matrices::SlicedEllpack< long,   TNL::Devices::Host, short >,
+    TNL::Matrices::SlicedEllpack< float,  TNL::Devices::Host, short >,
+    TNL::Matrices::SlicedEllpack< double, TNL::Devices::Host, short >,
+    TNL::Matrices::SlicedEllpack< int,    TNL::Devices::Host, int >,
+    TNL::Matrices::SlicedEllpack< long,   TNL::Devices::Host, int >,
+    TNL::Matrices::SlicedEllpack< float,  TNL::Devices::Host, int >,
+    TNL::Matrices::SlicedEllpack< double, TNL::Devices::Host, int >,
+    TNL::Matrices::SlicedEllpack< int,    TNL::Devices::Host, long >,
+    TNL::Matrices::SlicedEllpack< long,   TNL::Devices::Host, long >,
+    TNL::Matrices::SlicedEllpack< float,  TNL::Devices::Host, long >,
+    TNL::Matrices::SlicedEllpack< double, TNL::Devices::Host, long >
+#ifdef HAVE_CUDA
+    ,TNL::Matrices::SlicedEllpack< int,    TNL::Devices::Cuda, short >,
+    TNL::Matrices::SlicedEllpack< long,   TNL::Devices::Cuda, short >,
+    TNL::Matrices::SlicedEllpack< float,  TNL::Devices::Cuda, short >,
+    TNL::Matrices::SlicedEllpack< double, TNL::Devices::Cuda, short >,
+    TNL::Matrices::SlicedEllpack< int,    TNL::Devices::Cuda, int >,
+    TNL::Matrices::SlicedEllpack< long,   TNL::Devices::Cuda, int >,
+    TNL::Matrices::SlicedEllpack< float,  TNL::Devices::Cuda, int >,
+    TNL::Matrices::SlicedEllpack< double, TNL::Devices::Cuda, int >,
+    TNL::Matrices::SlicedEllpack< int,    TNL::Devices::Cuda, long >,
+    TNL::Matrices::SlicedEllpack< long,   TNL::Devices::Cuda, long >,
+    TNL::Matrices::SlicedEllpack< float,  TNL::Devices::Cuda, long >,
+    TNL::Matrices::SlicedEllpack< double, TNL::Devices::Cuda, long >
+#endif
+>;
+
+TYPED_TEST_CASE( SlicedEllpackMatrixTest, SlicedEllpackMatrixTypes );
+
+TYPED_TEST( SlicedEllpackMatrixTest, setDimensionsTest )
+{
+    using SlicedEllpackMatrixType = typename TestFixture::SlicedEllpackMatrixType;
+    
+    test_SetDimensions< SlicedEllpackMatrixType >();
+}
+
+//TYPED_TEST( SlicedEllpackMatrixTest, setCompressedRowLengthsTest )
+//{
+////    using SlicedEllpackMatrixType = typename TestFixture::SlicedEllpackMatrixType;
+//    
+////    test_SetCompressedRowLengths< SlicedEllpackMatrixType >();
+//    
+//    bool testRan = false;
+//    EXPECT_TRUE( testRan );
+//    std::cout << "\nTEST DID NOT RUN. NOT WORKING.\n\n";
+//    std::cout << "      This test is dependent on the input format. \n";
+//    std::cout << "      Almost every format allocates elements per row differently.\n\n";
+//    std::cout << "\n    TODO: Finish implementation of getNonZeroRowLength (Only non-zero elements, not the number of allocated elements.)\n\n";
+//}
+
+TYPED_TEST( SlicedEllpackMatrixTest, setLikeTest )
+{
+    using SlicedEllpackMatrixType = typename TestFixture::SlicedEllpackMatrixType;
+    
+    test_SetLike< SlicedEllpackMatrixType, SlicedEllpackMatrixType >();
+}
+
+TYPED_TEST( SlicedEllpackMatrixTest, resetTest )
+{
+    using SlicedEllpackMatrixType = typename TestFixture::SlicedEllpackMatrixType;
+    
+    test_Reset< SlicedEllpackMatrixType >();
+}
+
+TYPED_TEST( SlicedEllpackMatrixTest, setElementTest )
+{
+    using SlicedEllpackMatrixType = typename TestFixture::SlicedEllpackMatrixType;
+    
+    test_SetElement< SlicedEllpackMatrixType >();
+}
+
+TYPED_TEST( SlicedEllpackMatrixTest, addElementTest )
+{
+    using SlicedEllpackMatrixType = typename TestFixture::SlicedEllpackMatrixType;
+    
+    test_AddElement< SlicedEllpackMatrixType >();
+}
+
+TYPED_TEST( SlicedEllpackMatrixTest, setRowTest )
+{
+    using SlicedEllpackMatrixType = typename TestFixture::SlicedEllpackMatrixType;
+    
+    test_SetRow< SlicedEllpackMatrixType >();
+}
+
+TYPED_TEST( SlicedEllpackMatrixTest, vectorProductTest )
+{
+    using SlicedEllpackMatrixType = typename TestFixture::SlicedEllpackMatrixType;
+    
+    test_VectorProduct< SlicedEllpackMatrixType >();
+}
+
+TYPED_TEST( SlicedEllpackMatrixTest, saveAndLoadTest )
+{
+    using SlicedEllpackMatrixType = typename TestFixture::SlicedEllpackMatrixType;
+    
+    test_SaveAndLoad< SlicedEllpackMatrixType >();
+}
+
+TYPED_TEST( SlicedEllpackMatrixTest, printTest )
+{
+    using SlicedEllpackMatrixType = typename TestFixture::SlicedEllpackMatrixType;
+    
+    test_Print< SlicedEllpackMatrixType >();
+}
+
+//// test_getType is not general enough yet. DO NOT TEST IT YET.
+
+//TEST( SparseMatrixTest, CSR_GetTypeTest_Host )
+//{
+//    host_test_GetType< CSR_host_float, CSR_host_int >();
+//}
+//
+//#ifdef HAVE_CUDA
+//TEST( SparseMatrixTest, CSR_GetTypeTest_Cuda )
+//{
+//    cuda_test_GetType< CSR_cuda_float, CSR_cuda_int >();
+//}
+//#endif
+
+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 "../GtestMissingError.h"
+int main( int argc, char* argv[] )
+{
+#ifdef HAVE_GTEST
+   ::testing::InitGoogleTest( &argc, argv );
+   return RUN_ALL_TESTS();
+#else
+   throw GtestMissingError();
+#endif
+}
\ No newline at end of file
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_impl.h b/src/UnitTests/Matrices/SparseMatrixTest_impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..8c3858641d8622424e42cb182646f2d8c854b10b
--- /dev/null
+++ b/src/UnitTests/Matrices/SparseMatrixTest_impl.h
@@ -0,0 +1,926 @@
+/***************************************************************************
+                          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 */
+
+// TODO
+/*
+ * getType()                        ::HOW?  How to test this for each format? edit string how?
+ *      Found the mistake for Cuda instead of Devices::Cuda. Incorrect String in src/TNL/Devices/Cuda.cpp
+ *      MISSING: indexType is missing in CSR_impl.h
+ * getTypeVirtual()                 ::TEST? This just calls getType().
+ * getSerializationType()           ::TEST? This just calls HostType::getType().
+ * getSerializationTypeVirtual()    ::TEST? This just calls getSerializationType().
+ * setDimensions()                      ::DONE
+ * setCompressedRowLengths()            ::DONE
+ * getRowLength()                   ::USED! In test_SetCompressedRowLengths() to verify the test itself.
+ * getRowLengthFast()               ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * setLike()                            ::DONE
+ * reset()                              ::DONE
+ * setElementFast()                 ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * setElement()                         ::DONE
+ * addElementFast()                 ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * addElement()                         ::DONE
+ * setRowFast()                     ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * setRow()                             ::DONE
+ *      MISTAKE!!! In SlicedEllpack: addElement(), line 263, "column <= this->rows" shouldn't it be: "column <= this->columns", otherwise test_SetRow causes the assertion to fail.
+ * addRowFast()                     ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * addRow()                         ::NOT IMPLEMENTED! This calls addRowFast() which isn't implemented. Implement? Is it supposed to add an extra row to the matrix or add elements of a row to another row in the matrix?
+ * getElementFast()                 ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * getElement()                     ::USED! In test_SetElement(), test_AddElement() and test_setRow() to verify the test itself.
+ * getRowFast()                     ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * MatrixRow getRow()               ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * ConstMatrixRow getRow()          ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * rowVectorProduct()               ::TEST? How to test __cuda_callable__? ONLY TEST ON CPU FOR NOW
+ * vectorProduct()                      ::DONE
+ *      This used to throw illegal memory access, but instead of using ints for vectors, using Types, helped.
+ * addMatrix()                      ::NOT IMPLEMENTED!
+ * getTransposition()               ::NOT IMPLMENETED!
+ * performSORIteration()            ::HOW? Throws segmentation fault CUDA.
+ * operator=()                      ::HOW? What is this supposed to enable? Overloading operators?
+ * save( File& file)                ::USED! In save( String& fileName )
+ * load( File& file )               ::USED! In load( String& fileName )
+ * save( String& fileName )             ::DONE
+ * load( String& fileName )             ::DONE
+ * print()                              ::DONE
+ * setCudaKernelType()              ::NOT SUPPOSED TO TEST! via notes from 1.11.2018 supervisor meeting.
+ * getCudaKernelType()              ::NOT SUPPOSED TO TEST! via notes from 1.11.2018 supervisor meeting.
+ * setCudaWarpSize()                ::NOT SUPPOSED TO TEST! via notes from 1.11.2018 supervisor meeting.
+ * getCudaWarpSize()                ::NOT SUPPOSED TO TEST! via notes from 1.11.2018 supervisor meeting.
+ * setHybridModeSplit()             ::NOT SUPPOSED TO TEST! via notes from 1.11.2018 supervisor meeting.
+ * getHybridModeSplit()             ::NOT SUPPOSED TO TEST! via notes from 1.11.2018 supervisor meeting.
+ * spmvCudaVectorized()             ::TEST? How to test __device__?
+ * vectorProductCuda()              ::TEST? How to test __device__?
+ */
+
+// GENERAL TODO
+/*
+ * For every function, EXPECT_EQ needs to be done, even for zeros in matrices.
+ * Figure out __cuda_callable_. When trying to call __cuda_callable__ functions
+ *      a segmentation fault (core dumped) is thrown.
+ *  ==>__cuda_callable__ works only for CPU at the moment. (for loops vs thread kernel assignment).
+ *                       If we want to use __cuda_callable__ on the GPU, we need to call it as a kernel.
+ */
+
+#include <TNL/Containers/Vector.h>
+#include <TNL/Containers/VectorView.h>
+#include <TNL/Math.h>
+#include <iostream>
+
+#ifdef HAVE_GTEST 
+#include <gtest/gtest.h>
+
+template< typename MatrixHostFloat, typename MatrixHostInt >
+void host_test_GetType()
+{
+    MatrixHostFloat mtrxHostFloat;
+    MatrixHostInt mtrxHostInt;
+    
+    
+    EXPECT_EQ( mtrxHostFloat.getType(), TNL::String( "Matrices::CSR< float, Devices::Host >" ) );
+    EXPECT_EQ( mtrxHostInt.getType(), TNL::String( "Matrices::CSR< int, Devices::Host >" ) ); 
+}
+
+template< typename MatrixCudaFloat, typename MatrixCudaInt >
+void cuda_test_GetType()
+{
+    MatrixCudaFloat mtrxCudaFloat;
+    MatrixCudaInt mtrxCudaInt;
+    
+
+    EXPECT_EQ( mtrxCudaFloat.getType(), TNL::String( "Matrices::CSR< float, Cuda >" ) );    // This is mistakenly labeled in /src/TNL/Devices/Cuda.cpp
+    EXPECT_EQ( mtrxCudaInt.getType(), TNL::String( "Matrices::CSR< int, Cuda >" ) );        // Should be Devices::Cuda
+}
+
+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 );
+    
+//    if( m.getType() == TNL::String( TNL::String( "Matrices::CSR< ") +
+//                       TNL::String( TNL::getType< RealType >() ) +
+//                       TNL::String( ", " ) +
+//                       TNL::String( Matrix::DeviceType::getDeviceType() ) +
+//                       //TNL::String( ", " ) +
+//                       //TNL::String( TNL::getType< IndexType >() ) +
+//                       TNL::String( " >" ) )
+//      )
+//    {
+//        EXPECT_EQ( m.getRowLength( 0 ), 3 );
+//        EXPECT_EQ( m.getRowLength( 1 ), 3 );
+//        EXPECT_EQ( m.getRowLength( 2 ), 1 );
+//        EXPECT_EQ( m.getRowLength( 3 ), 2 );
+//        EXPECT_EQ( m.getRowLength( 4 ), 3 );
+//        EXPECT_EQ( m.getRowLength( 5 ), 4 );
+//        EXPECT_EQ( m.getRowLength( 6 ), 5 );
+//        EXPECT_EQ( m.getRowLength( 7 ), 6 );
+//        EXPECT_EQ( m.getRowLength( 8 ), 7 );
+//        EXPECT_EQ( m.getRowLength( 9 ), 8 );
+//    }
+//    else if( m.getType() == TNL::String( TNL::String( "Matrices::AdEllpack< ") +
+//                            TNL::String( TNL::getType< RealType >() ) +
+//                            TNL::String( ", " ) +
+//                            TNL::String( Matrix::DeviceType::getDeviceType() ) +
+//                            TNL::String( ", " ) +
+//                            TNL::String( TNL::getType< IndexType >() ) +
+//                            TNL::String( " >" ) ) 
+//                            || 
+//             m.getType() == TNL::String( TNL::String( "Matrices::SlicedEllpack< ") +
+//                            TNL::String( TNL::getType< RealType >() ) +
+//                            TNL::String( ", " ) +
+//                            TNL::String( Matrix::DeviceType::getDeviceType() ) +
+//                            TNL::String( " >" ) )
+//           )
+//    {
+//        EXPECT_EQ( m.getRowLength( 0 ), 8 );
+//        EXPECT_EQ( m.getRowLength( 1 ), 8 );
+//        EXPECT_EQ( m.getRowLength( 2 ), 8 );
+//        EXPECT_EQ( m.getRowLength( 3 ), 8 );
+//        EXPECT_EQ( m.getRowLength( 4 ), 8 );
+//        EXPECT_EQ( m.getRowLength( 5 ), 8 );
+//        EXPECT_EQ( m.getRowLength( 6 ), 8 );
+//        EXPECT_EQ( m.getRowLength( 7 ), 8 );
+//        EXPECT_EQ( m.getRowLength( 8 ), 8 );
+//        EXPECT_EQ( m.getRowLength( 9 ), 8 );
+//    }
+//    else if( m.getType() == TNL::String( TNL::String( "Matrices::Ellpack< ") +
+//                            TNL::String( TNL::getType< RealType >() ) +
+//                            TNL::String( ", " ) +
+//                            TNL::String( Matrix::DeviceType::getDeviceType() ) +
+//                            TNL::String( ", " ) +
+//                            TNL::String( TNL::getType< IndexType >() ) +
+//                            TNL::String( " >" ) ) 
+//                            ||
+//             m.getType() == TNL::String( TNL::String( "Matrices::ChunkedEllpack< ") +
+//                            TNL::String( TNL::getType< RealType >() ) +
+//                            TNL::String( ", " ) +
+//                            TNL::String( Matrix::DeviceType::getDeviceType() ) +
+//                            TNL::String( " >" ) )
+//           )
+//    {
+//        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 );
+//    }
+//    else
+//    {
+//        EXPECT_EQ( m.getRowLength( 0 ), 3 );
+//        EXPECT_EQ( m.getRowLength( 1 ), 3 );
+//        EXPECT_EQ( m.getRowLength( 2 ), 1 );
+//        EXPECT_EQ( m.getRowLength( 3 ), 2 );
+//        EXPECT_EQ( m.getRowLength( 4 ), 3 );
+//        EXPECT_EQ( m.getRowLength( 5 ), 4 );
+//        EXPECT_EQ( m.getRowLength( 6 ), 5 );
+//        EXPECT_EQ( m.getRowLength( 7 ), 6 );
+//        EXPECT_EQ( m.getRowLength( 8 ), 7 );
+//        EXPECT_EQ( m.getRowLength( 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_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_SetElement()
+{
+    using RealType = typename Matrix::RealType;
+    using DeviceType = typename Matrix::DeviceType;
+    using IndexType = typename Matrix::IndexType;
+    
+/*
+ * Sets up the following 5x5 sparse matrix:
+ *
+ *    /  1  0  0  0  0 \
+ *    |  0  2  0  0  0 |
+ *    |  0  0  3  0  0 |
+ *    |  0  0  0  4  0 |
+ *    \  0  0  0  0  5 /
+ */
+    
+    const IndexType rows = 5;
+    const IndexType cols = 5;
+    
+    Matrix m;
+    m.reset();
+    m.setDimensions( rows, cols );
+    typename Matrix::CompressedRowLengthsVector rowLengths;
+    rowLengths.setSize( rows );
+    rowLengths.setValue( 1 );
+    m.setCompressedRowLengths( rowLengths );    
+    
+    RealType value = 1;
+    for( IndexType i = 0; i < rows; i++ )
+        m.setElement( i, i, value++ );
+    
+    
+    EXPECT_EQ( m.getElement( 0, 0 ), 1 );
+    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 ), 0 );
+    
+    EXPECT_EQ( m.getElement( 1, 0 ), 0 );
+    EXPECT_EQ( m.getElement( 1, 1 ), 2 );
+    EXPECT_EQ( m.getElement( 1, 2 ), 0 );
+    EXPECT_EQ( m.getElement( 1, 3 ), 0 );
+    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 ), 3 );
+    EXPECT_EQ( m.getElement( 2, 3 ), 0 );
+    EXPECT_EQ( m.getElement( 2, 4 ), 0 );
+    
+    EXPECT_EQ( m.getElement( 3, 0 ), 0 );
+    EXPECT_EQ( m.getElement( 3, 1 ), 0 );
+    EXPECT_EQ( m.getElement( 3, 2 ), 0 );
+    EXPECT_EQ( m.getElement( 3, 3 ), 4 );
+    EXPECT_EQ( m.getElement( 3, 4 ), 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 ), 0 );
+    EXPECT_EQ( m.getElement( 4, 4 ), 5 );
+}
+
+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.
+
+/*
+ * 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;
+    
+/*
+ * 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++ );
+
+    using VectorType = TNL::Containers::Vector< RealType, DeviceType, IndexType >;
+    
+    VectorType inVector;
+    inVector.setSize( 4 );
+    for( IndexType i = 0; i < inVector.getSize(); i++ )        
+        inVector.setElement( i, 2 );
+
+    VectorType outVector;  
+    outVector.setSize( 5 );
+    for( IndexType j = 0; j < outVector.getSize(); j++ )
+        outVector.setElement( j, 0 );
+ 
+    
+    m.vectorProduct( inVector, outVector);
+    
+   
+    EXPECT_EQ( outVector.getElement( 0 ), 12 );
+    EXPECT_EQ( outVector.getElement( 1 ),  8 );
+    EXPECT_EQ( outVector.getElement( 2 ), 36 );
+    EXPECT_EQ( outVector.getElement( 3 ), 54 );
+    EXPECT_EQ( outVector.getElement( 4 ), 46 );
+}
+
+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 );
+}
+
+template< typename Matrix >
+void test_SaveAndLoad()
+{
+   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++ );
+        
+    savedMatrix.save( "sparseMatrixFile" );
+    
+    Matrix loadedMatrix;
+    loadedMatrix.reset();
+    loadedMatrix.setDimensions( m_rows, m_cols );
+    typename Matrix::CompressedRowLengthsVector rowLengths2;
+    rowLengths2.setSize( m_rows );
+    rowLengths2.setValue( 3 );
+    loadedMatrix.setCompressedRowLengths( rowLengths2 );
+    
+    
+    loadedMatrix.load( "sparseMatrixFile" );
+    
+    
+    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 );
+    
+    std::remove( "sparseMatrixFile" );
+}
+
+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++ );
+    
+    // This is from: https://stackoverflow.com/questions/5193173/getting-cout-output-to-a-stdstring
+    #include <sstream>
+    std::stringstream printed;
+    std::stringstream couted;
+    
+    // This is from: https://stackoverflow.com/questions/19485536/redirect-output-of-an-function-printing-to-console-to-string
+    //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
+    
+    //printed << printed.str() << std::endl;
+    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
\ No newline at end of file