From 8573f541dd5dd163e2f47b824787fe9464764c27 Mon Sep 17 00:00:00 2001
From: Tomas Oberhuber <tomas.oberhuber@fjfi.cvut.cz>
Date: Wed, 11 Dec 2019 18:00:42 +0100
Subject: [PATCH] Implementing sparse matrcies assignment.

---
 src/TNL/Containers/Segments/CSR.hpp           |  11 +-
 src/TNL/Containers/Segments/Ellpack.hpp       |  18 +-
 src/TNL/Containers/Segments/SlicedEllpack.hpp |   6 +-
 src/TNL/Matrices/Matrix.h                     |   9 +-
 src/TNL/Matrices/SparseMatrix.h               |  19 +-
 src/TNL/Matrices/SparseMatrix.hpp             |  88 +++-
 src/UnitTests/Matrices/SparseMatrixTest.hpp   | 494 +++++++++---------
 7 files changed, 374 insertions(+), 271 deletions(-)

diff --git a/src/TNL/Containers/Segments/CSR.hpp b/src/TNL/Containers/Segments/CSR.hpp
index ef7431038c..ccb483125d 100644
--- a/src/TNL/Containers/Segments/CSR.hpp
+++ b/src/TNL/Containers/Segments/CSR.hpp
@@ -154,11 +154,12 @@ CSR< Device, Index >::
 forSegments( IndexType first, IndexType last, Function& f, Args... args ) const
 {
    const auto offsetsView = this->offsets.getConstView();
-   auto l = [=] __cuda_callable__ ( const IndexType i, Args... args ) mutable {
-      const IndexType begin = offsetsView[ i ];
-      const IndexType end = offsetsView[ i + 1 ];
-      for( IndexType j = begin; j < end; j++  )
-         if( ! f( i, j, args... ) )
+   auto l = [=] __cuda_callable__ ( const IndexType segmentIdx, Args... args ) mutable {
+      const IndexType begin = offsetsView[ segmentIdx ];
+      const IndexType end = offsetsView[ segmentIdx + 1 ];
+      IndexType localIdx( 0 );
+      for( IndexType globalIdx = begin; globalIdx < end; globalIdx++  )
+         if( ! f( segmentIdx, localIdx++, globalIdx, args... ) )
             break;
    };
    Algorithms::ParallelFor< Device >::exec( first, last, l, args... );
diff --git a/src/TNL/Containers/Segments/Ellpack.hpp b/src/TNL/Containers/Segments/Ellpack.hpp
index d3d90be5e0..337009e99a 100644
--- a/src/TNL/Containers/Segments/Ellpack.hpp
+++ b/src/TNL/Containers/Segments/Ellpack.hpp
@@ -192,11 +192,12 @@ forSegments( IndexType first, IndexType last, Function& f, Args... args ) const
    if( RowMajorOrder )
    {
       const IndexType segmentSize = this->segmentSize;
-      auto l = [=] __cuda_callable__ ( const IndexType i, Args... args ) mutable {
-         const IndexType begin = i * segmentSize;
+      auto l = [=] __cuda_callable__ ( const IndexType segmentIdx, Args... args ) mutable {
+         const IndexType begin = segmentIdx * segmentSize;
          const IndexType end = begin + segmentSize;
-         for( IndexType j = begin; j < end; j++  )
-            if( ! f( i, j, args... ) )
+         IndexType localIdx( 0 );
+         for( IndexType globalIdx = begin; globalIdx < end; globalIdx++  )
+            if( ! f( segmentIdx, localIdx++, globalIdx,  args... ) )
                break;
       };
       Algorithms::ParallelFor< Device >::exec( first, last, l, args... );
@@ -205,11 +206,12 @@ forSegments( IndexType first, IndexType last, Function& f, Args... args ) const
    {
       const IndexType storageSize = this->getStorageSize();
       const IndexType alignedSize = this->alignedSize;
-      auto l = [=] __cuda_callable__ ( const IndexType i, Args... args ) mutable {
-         const IndexType begin = i;
+      auto l = [=] __cuda_callable__ ( const IndexType segmentIdx, Args... args ) mutable {
+         const IndexType begin = segmentIdx;
          const IndexType end = storageSize;
-         for( IndexType j = begin; j < end; j += alignedSize )
-            if( ! f( i, j, args... ) )
+         IndexType localIdx( 0 );
+         for( IndexType globalIdx = begin; globalIdx < end; globalIdx += alignedSize )
+            if( ! f( segmentIdx, localIdx++, globalIdx, args... ) )
                break;
       };
       Algorithms::ParallelFor< Device >::exec( first, last, l, args... );
diff --git a/src/TNL/Containers/Segments/SlicedEllpack.hpp b/src/TNL/Containers/Segments/SlicedEllpack.hpp
index c8e74ec591..d721edb00b 100644
--- a/src/TNL/Containers/Segments/SlicedEllpack.hpp
+++ b/src/TNL/Containers/Segments/SlicedEllpack.hpp
@@ -225,8 +225,9 @@ forSegments( IndexType first, IndexType last, Function& f, Args... args ) const
          const IndexType segmentSize = sliceSegmentSizes_view[ sliceIdx ];
          const IndexType begin = sliceOffsets_view[ sliceIdx ] + segmentInSliceIdx * segmentSize;
          const IndexType end = begin + segmentSize;
+         IndexType localIdx( 0 );
          for( IndexType globalIdx = begin; globalIdx < end; globalIdx++  )
-            if( ! f( segmentIdx, globalIdx, args... ) )
+            if( ! f( segmentIdx, localIdx++, globalIdx, args... ) )
                break;
       };
       Algorithms::ParallelFor< Device >::exec( first, last, l, args... );
@@ -239,8 +240,9 @@ forSegments( IndexType first, IndexType last, Function& f, Args... args ) const
          const IndexType segmentSize = sliceSegmentSizes_view[ sliceIdx ];
          const IndexType begin = sliceOffsets_view[ sliceIdx ] + segmentInSliceIdx;
          const IndexType end = sliceOffsets_view[ sliceIdx + 1 ];
+         IndexType localIdx( 0 );
          for( IndexType globalIdx = begin; globalIdx < end; globalIdx += SliceSize )
-            if( ! f( segmentIdx, globalIdx, args... ) )
+            if( ! f( segmentIdx, localIdx++, globalIdx, args... ) )
                break;
       };
       Algorithms::ParallelFor< Device >::exec( first, last, l, args... );
diff --git a/src/TNL/Matrices/Matrix.h b/src/TNL/Matrices/Matrix.h
index a877fd5c24..4a038eb2e8 100644
--- a/src/TNL/Matrices/Matrix.h
+++ b/src/TNL/Matrices/Matrix.h
@@ -39,7 +39,7 @@ public:
    using RealAllocatorType = RealAllocator;
 
    Matrix( const RealAllocatorType& allocator = RealAllocatorType() );
-   
+
    Matrix( const IndexType rows,
            const IndexType columns,
            const RealAllocatorType& allocator = RealAllocatorType() );
@@ -100,9 +100,9 @@ public:
 
    virtual Real getElement( const IndexType row,
                             const IndexType column ) const = 0;
-   
+
    const ValuesVector& getValues() const;
-   
+
    ValuesVector& getValues();
 
    // TODO: parallelize and optimize for sparse matrices
@@ -137,7 +137,8 @@ public:
    __cuda_callable__
    Index getValuesSize() const;
 
-   protected:
+   // TODO: restore this
+   //protected:
 
    IndexType rows, columns, numberOfColors;
 
diff --git a/src/TNL/Matrices/SparseMatrix.h b/src/TNL/Matrices/SparseMatrix.h
index 268fba6d39..b6a618e105 100644
--- a/src/TNL/Matrices/SparseMatrix.h
+++ b/src/TNL/Matrices/SparseMatrix.h
@@ -39,7 +39,7 @@ class SparseMatrix : public Matrix< Real, Device, Index, RealAllocator >
       using ConstRowsCapacitiesView = typename RowsCapacitiesView::ConstViewType;
       using ValuesVectorType = typename Matrix< Real, Device, Index, RealAllocator >::ValuesVector;
       using ColumnsVectorType = Containers::Vector< IndexType, DeviceType, IndexType, IndexAllocatorType >;
-      
+
       // TODO: remove this - it is here only for compatibility with original matrix implementation
       typedef Containers::Vector< IndexType, DeviceType, IndexType > CompressedRowLengthsVector;
       typedef Containers::VectorView< IndexType, DeviceType, IndexType > CompressedRowLengthsVectorView;
@@ -64,6 +64,9 @@ class SparseMatrix : public Matrix< Real, Device, Index, RealAllocator >
 
       void setCompressedRowLengths( ConstCompressedRowLengthsVectorView rowLengths );
 
+      template< typename Vector >
+      void getCompressedRowLengths( Vector& rowLengths ) const;
+
       IndexType getRowLength( const IndexType row ) const;
 
       __cuda_callable__
@@ -167,10 +170,16 @@ class SparseMatrix : public Matrix< Real, Device, Index, RealAllocator >
 
       template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >
       void rowsReduction( IndexType first, IndexType last, Fetch& fetch, Reduce& reduce, Keep& keep, const FetchReal& zero ) const;
-      
+
       template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >
       void allRowsReduction( Fetch& fetch, Reduce& reduce, Keep& keep, const FetchReal& zero ) const;
 
+      template< typename Function >
+      void forRows( IndexType first, IndexType last, Function& function ) const;
+
+      template< typename Function >
+      void forAllRows( Function& function ) const;
+
       template< typename Vector1, typename Vector2 >
       bool performSORIteration( const Vector1& b,
                                 const IndexType row,
@@ -201,7 +210,9 @@ class SparseMatrix : public Matrix< Real, Device, Index, RealAllocator >
 
       __cuda_callable__
       IndexType getPaddingIndex() const;
-   protected:
+
+// TODO: restore it and also in Matrix
+//   protected:
 
       ColumnsVectorType columnIndexes;
 
@@ -210,6 +221,8 @@ class SparseMatrix : public Matrix< Real, Device, Index, RealAllocator >
       IndexAllocator indexAlloctor;
 
       RealAllocator realAllocator;
+
+
 };
 
 }  // namespace Conatiners
diff --git a/src/TNL/Matrices/SparseMatrix.hpp b/src/TNL/Matrices/SparseMatrix.hpp
index e26693f6a1..3605daaef4 100644
--- a/src/TNL/Matrices/SparseMatrix.hpp
+++ b/src/TNL/Matrices/SparseMatrix.hpp
@@ -116,6 +116,32 @@ setCompressedRowLengths( ConstCompressedRowLengthsVectorView rowLengths )
    this->columnIndexes = this->getPaddingIndex();
 }
 
+template< typename Real,
+          template< typename, typename > class Segments,
+          typename Device,
+          typename Index,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Vector >
+void
+SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+getCompressedRowLengths( Vector& rowLengths ) const
+{
+   rowLengths.setSize( this->getRows() );
+   rowLengths = 0;
+   auto rowLengths_view = rowLengths.getView();
+   auto fetch = [] __cuda_callable__ ( IndexType row, IndexType column, const RealType& value ) -> IndexType {
+      return ( value != 0.0 );
+   };
+   auto reduce = [] __cuda_callable__ ( IndexType& aux, const IndexType a ) {
+      aux += a;
+   };
+   auto keep = [=] __cuda_callable__ ( const IndexType rowIdx, const IndexType value ) mutable {
+      rowLengths_view[ rowIdx ] = value;
+   };
+   this->allRowsReduction( fetch, reduce, keep, 0 );
+}
+
 template< typename Real,
           template< typename, typename > class Segments,
           typename Device,
@@ -554,6 +580,43 @@ allRowsReduction( Fetch& fetch, Reduce& reduce, Keep& keep, const FetchReal& zer
    this->rowsReduction( 0, this->getRows(), fetch, reduce, keep, zero );
 }
 
+template< typename Real,
+          template< typename, typename > class Segments,
+          typename Device,
+          typename Index,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Function >
+void
+SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+forRows( IndexType first, IndexType last, Function& function ) const
+{
+   const auto columns_view = this->columnIndexes.getConstView();
+   const auto values_view = this->values.getConstView();
+   const IndexType paddingIndex_ = this->getPaddingIndex();
+   auto fetch_ = [=] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, IndexType globalIdx ) mutable -> decltype( fetch( IndexType(), IndexType(), RealType() ) ) {
+      IndexType columnIdx = columns_view[ globalIdx ];
+      if( columnIdx != paddingIndex_ )
+         return fetch( rowIdx, columnIdx, values_view[ globalIdx ] );
+      return zero;
+   };
+   this->segments.segmentsReduction( first, last, fetch_, reduce, keep, zero );
+
+}
+
+template< typename Real,
+          template< typename, typename > class Segments,
+          typename Device,
+          typename Index,
+          typename RealAllocator,
+          typename IndexAllocator >
+   template< typename Function >
+void
+SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+forAllRows( Function& function ) const
+{
+   this->forRows( 0, this->getRows(), function );
+}
 
 /*template< typename Real,
           template< typename, typename > class Segments,
@@ -638,12 +701,31 @@ SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >&
 SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
 operator=( const SparseMatrix< Real2, Segments2, Device2, Index2, RealAllocator2, IndexAllocator2 >& matrix )
 {
+   using RHSMatrixType = SparseMatrix< Real2, Segments2, Device2, Index2, RealAllocator2, IndexAllocator2 >;
    if( std::is_same< Device, Device2 >::value )
    {
-      
+      RowsCapacitiesType rowLengths;
+      matrix.getCompressedRowLengths( rowLengths );
+      this->setCompressedRowLengths( rowLengths );
+      // TODO: Replace this with SparseMatrixView
+      const auto matrix_columns_view = matrix.columnIndexes.getConstView();
+      const auto matrix_values_view = matrix.values.getConstView();
+      const auto segments_view = this->segments.getConstView();
+      auto this_columns_view = this->columnIndexes.getView();
+      auto this_values_view = this->values.getView();
+      const IndexType paddingIndex = this->getPaddingIndex();
+      auto f = [=] __cuda_callable__ ( IndexType rowIdx, IndexType localIdx, IndexType globalIdx ) {
+         const IndexType column = matrix_columns_view[ globalIdx ];
+         if( column != paddingIndex )
+         {
+            const RealType value = matrix_values_view[ globalIdx ];
+            IndexType thisGlobalIdx = segments_view.getGlobalIdx( rowIdx, localIdx );
+            this_columns_view[ thisGlobalIdx ] = column;
+            this_values_view[ thisGlobalIdx ] = value;
+         }
+      };
+      matrix.forAllRows( f );
    }
-      
-   
 }
 
 template< typename Real,
diff --git a/src/UnitTests/Matrices/SparseMatrixTest.hpp b/src/UnitTests/Matrices/SparseMatrixTest.hpp
index b366f4e2fe..07a60178f1 100644
--- a/src/UnitTests/Matrices/SparseMatrixTest.hpp
+++ b/src/UnitTests/Matrices/SparseMatrixTest.hpp
@@ -18,7 +18,7 @@
 #include <TNL/Matrices/AdEllpack.h>
 #include <TNL/Matrices/BiEllpack.h>
 
-#ifdef HAVE_GTEST 
+#ifdef HAVE_GTEST
 #include <gtest/gtest.h>
 
 template< typename MatrixHostFloat, typename MatrixHostInt >
@@ -36,7 +36,7 @@ void cuda_test_GetType()
     bool testRan = false;
     EXPECT_TRUE( testRan );
     std::cout << "\nTEST DID NOT RUN. NOT WORKING.\n\n";
-    std::cerr << "This test has not been implemented properly yet.\n" << std::endl;    
+    std::cerr << "This test has not been implemented properly yet.\n" << std::endl;
 }
 
 template< typename Matrix >
@@ -45,13 +45,13 @@ 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 );
 }
@@ -62,41 +62,41 @@ 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++ );
 
@@ -111,8 +111,8 @@ void test_SetCompressedRowLengths()
 
     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 );
@@ -131,21 +131,21 @@ 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() );
 }
@@ -174,7 +174,7 @@ void test_GetNumberOfNonzeroMatrixElements()
 
    const IndexType rows = 10;
    const IndexType cols = 10;
-    
+
    Matrix m;
    m.reset();
 
@@ -225,7 +225,7 @@ 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:
  *
@@ -235,16 +235,16 @@ void test_Reset()
  *    |  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 );
 }
@@ -255,7 +255,7 @@ void test_SetElement()
     using RealType = typename Matrix::RealType;
     using DeviceType = typename Matrix::DeviceType;
     using IndexType = typename Matrix::IndexType;
-    
+
 /*
  * Sets up the following 10x10 sparse matrix:
  *
@@ -270,15 +270,15 @@ void test_SetElement()
  *    | 22 23 24 25 26 27 28 29 30 31  |
  *    \ 32 33 34 35 36 37 38 39 40 41 /
  */
-    
+
     const IndexType rows = 10;
     const IndexType cols = 10;
-    
+
     Matrix m;
     m.reset();
-    
+
     m.setDimensions( rows, cols );
-    
+
     typename Matrix::CompressedRowLengthsVector rowLengths;
     rowLengths.setSize( rows );
     rowLengths.setElement( 0, 4 );
@@ -292,29 +292,29 @@ void test_SetElement()
     rowLengths.setElement( 8, 10 );
     rowLengths.setElement( 9, 10 );
     m.setCompressedRowLengths( rowLengths );
-    
+
     RealType value = 1;
     for( IndexType i = 0; i < 4; i++ )
         m.setElement( 0, 2 * i, value++ );
-    
+
     for( IndexType i = 0; i < 3; i++ )
         m.setElement( 1, i, value++ );
-    
+
     for( IndexType i = 0; i < 8; i++ )
         m.setElement( 2, i, value++ );
-    
+
     for( IndexType i = 0; i < 2; i++ )
         m.setElement( 3, i, value++ );
-    
+
     for( IndexType i = 4; i < 8; i++ )
         m.setElement( i, 0, value++ );
-    
+
     for( IndexType j = 8; j < rows; j++)
     {
         for( IndexType i = 0; i < cols; i++ )
             m.setElement( j, i, value++ );
     }
-    
+
     EXPECT_EQ( m.getElement( 0, 0 ),  1 );
     EXPECT_EQ( m.getElement( 0, 1 ),  0 );
     EXPECT_EQ( m.getElement( 0, 2 ),  2 );
@@ -325,7 +325,7 @@ void test_SetElement()
     EXPECT_EQ( m.getElement( 0, 7 ),  0 );
     EXPECT_EQ( m.getElement( 0, 8 ),  0 );
     EXPECT_EQ( m.getElement( 0, 9 ),  0 );
-    
+
     EXPECT_EQ( m.getElement( 1, 0 ),  5 );
     EXPECT_EQ( m.getElement( 1, 1 ),  6 );
     EXPECT_EQ( m.getElement( 1, 2 ),  7 );
@@ -336,7 +336,7 @@ void test_SetElement()
     EXPECT_EQ( m.getElement( 1, 7 ),  0 );
     EXPECT_EQ( m.getElement( 1, 8 ),  0 );
     EXPECT_EQ( m.getElement( 1, 9 ),  0 );
-    
+
     EXPECT_EQ( m.getElement( 2, 0 ),  8 );
     EXPECT_EQ( m.getElement( 2, 1 ),  9 );
     EXPECT_EQ( m.getElement( 2, 2 ), 10 );
@@ -347,7 +347,7 @@ void test_SetElement()
     EXPECT_EQ( m.getElement( 2, 7 ), 15 );
     EXPECT_EQ( m.getElement( 2, 8 ),  0 );
     EXPECT_EQ( m.getElement( 2, 9 ),  0 );
-    
+
     EXPECT_EQ( m.getElement( 3, 0 ), 16 );
     EXPECT_EQ( m.getElement( 3, 1 ), 17 );
     EXPECT_EQ( m.getElement( 3, 2 ),  0 );
@@ -358,7 +358,7 @@ void test_SetElement()
     EXPECT_EQ( m.getElement( 3, 7 ),  0 );
     EXPECT_EQ( m.getElement( 3, 8 ),  0 );
     EXPECT_EQ( m.getElement( 3, 9 ),  0 );
-    
+
     EXPECT_EQ( m.getElement( 4, 0 ), 18 );
     EXPECT_EQ( m.getElement( 4, 1 ),  0 );
     EXPECT_EQ( m.getElement( 4, 2 ),  0 );
@@ -369,7 +369,7 @@ void test_SetElement()
     EXPECT_EQ( m.getElement( 4, 7 ),  0 );
     EXPECT_EQ( m.getElement( 4, 8 ),  0 );
     EXPECT_EQ( m.getElement( 4, 9 ),  0 );
-    
+
     EXPECT_EQ( m.getElement( 5, 0 ), 19 );
     EXPECT_EQ( m.getElement( 5, 1 ),  0 );
     EXPECT_EQ( m.getElement( 5, 2 ),  0 );
@@ -380,7 +380,7 @@ void test_SetElement()
     EXPECT_EQ( m.getElement( 5, 7 ),  0 );
     EXPECT_EQ( m.getElement( 5, 8 ),  0 );
     EXPECT_EQ( m.getElement( 5, 9 ),  0 );
-    
+
     EXPECT_EQ( m.getElement( 6, 0 ), 20 );
     EXPECT_EQ( m.getElement( 6, 1 ),  0 );
     EXPECT_EQ( m.getElement( 6, 2 ),  0 );
@@ -391,7 +391,7 @@ void test_SetElement()
     EXPECT_EQ( m.getElement( 6, 7 ),  0 );
     EXPECT_EQ( m.getElement( 6, 8 ),  0 );
     EXPECT_EQ( m.getElement( 6, 9 ),  0 );
-    
+
     EXPECT_EQ( m.getElement( 7, 0 ), 21 );
     EXPECT_EQ( m.getElement( 7, 1 ),  0 );
     EXPECT_EQ( m.getElement( 7, 2 ),  0 );
@@ -402,7 +402,7 @@ void test_SetElement()
     EXPECT_EQ( m.getElement( 7, 7 ),  0 );
     EXPECT_EQ( m.getElement( 7, 8 ),  0 );
     EXPECT_EQ( m.getElement( 7, 9 ),  0 );
-    
+
     EXPECT_EQ( m.getElement( 8, 0 ), 22 );
     EXPECT_EQ( m.getElement( 8, 1 ), 23 );
     EXPECT_EQ( m.getElement( 8, 2 ), 24 );
@@ -413,7 +413,7 @@ void test_SetElement()
     EXPECT_EQ( m.getElement( 8, 7 ), 29 );
     EXPECT_EQ( m.getElement( 8, 8 ), 30 );
     EXPECT_EQ( m.getElement( 8, 9 ), 31 );
-    
+
     EXPECT_EQ( m.getElement( 9, 0 ), 32 );
     EXPECT_EQ( m.getElement( 9, 1 ), 33 );
     EXPECT_EQ( m.getElement( 9, 2 ), 34 );
@@ -432,7 +432,7 @@ 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:
  *
@@ -443,10 +443,10 @@ void test_AddElement()
  *    |  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 );
@@ -454,61 +454,61 @@ void test_AddElement()
     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.
 
 /*
@@ -521,7 +521,7 @@ void test_AddElement()
  *    |  0 11  0  0  0 |
  *    \  0  0  0 12  0 /
  */
-    
+
 /*
  * The following setup results in the following 6x5 sparse matrix:
  *
@@ -532,57 +532,57 @@ void test_AddElement()
  *    |  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 );
@@ -596,7 +596,7 @@ 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:
  *
@@ -604,10 +604,10 @@ void test_SetRow()
  *    |  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 );
@@ -616,7 +616,7 @@ void test_SetRow()
     rowLengths.setValue( 6 );
     rowLengths.setElement( 1, 3 );
     m.setCompressedRowLengths( rowLengths );
-    
+
     RealType value = 1;
     for( IndexType i = 0; i < 3; i++ )
     {
@@ -624,19 +624,19 @@ void test_SetRow()
         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 );
@@ -644,7 +644,7 @@ void test_SetRow()
     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 );
@@ -652,7 +652,7 @@ void test_SetRow()
     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 );
@@ -669,7 +669,7 @@ void test_VectorProduct()
     using DeviceType = typename Matrix::DeviceType;
     using IndexType = typename Matrix::IndexType;
     using VectorType = TNL::Containers::Vector< RealType, DeviceType, IndexType >;
-    
+
 /*
  * Sets up the following 4x4 sparse matrix:
  *
@@ -678,10 +678,10 @@ void test_VectorProduct()
  *    |  0  4  0  0 |
  *    \  0  0  5  0 /
  */
-    
+
     const IndexType m_rows_1 = 4;
     const IndexType m_cols_1 = 4;
-    
+
     Matrix m_1;
     m_1.reset();
     m_1.setDimensions( m_rows_1, m_cols_1 );
@@ -692,37 +692,37 @@ void test_VectorProduct()
     rowLengths_1.setElement( 2, 1 );
     rowLengths_1.setElement( 3, 1 );
     m_1.setCompressedRowLengths( rowLengths_1 );
-    
+
     RealType value_1 = 1;
     m_1.setElement( 0, 0, value_1++ );      // 0th row
-    
+
     m_1.setElement( 1, 1, value_1++ );      // 1st row
     m_1.setElement( 1, 3, value_1++ );
-        
+
     m_1.setElement( 2, 1, value_1++ );      // 2nd row
-        
+
     m_1.setElement( 3, 2, value_1++ );      // 3rd row
-    
+
     VectorType inVector_1;
     inVector_1.setSize( m_cols_1 );
-    for( IndexType i = 0; i < inVector_1.getSize(); i++ )        
+    for( IndexType i = 0; i < inVector_1.getSize(); i++ )
         inVector_1.setElement( i, 2 );
 
-    VectorType outVector_1;  
+    VectorType outVector_1;
     outVector_1.setSize( m_rows_1 );
     for( IndexType j = 0; j < outVector_1.getSize(); j++ )
         outVector_1.setElement( j, 0 );
- 
-    
+
+
     m_1.vectorProduct( inVector_1, outVector_1 );
-    
-   
+
+
     EXPECT_EQ( outVector_1.getElement( 0 ),  2 );
     EXPECT_EQ( outVector_1.getElement( 1 ), 10 );
     EXPECT_EQ( outVector_1.getElement( 2 ),  8 );
     EXPECT_EQ( outVector_1.getElement( 3 ), 10 );
-    
-    
+
+
 /*
  * Sets up the following 4x4 sparse matrix:
  *
@@ -731,10 +731,10 @@ void test_VectorProduct()
  *    |  5  6  7  0 |
  *    \  0  8  0  0 /
  */
-    
+
     const IndexType m_rows_2 = 4;
     const IndexType m_cols_2 = 4;
-    
+
     Matrix m_2;
     m_2.reset();
     m_2.setDimensions( m_rows_2, m_cols_2 );
@@ -744,39 +744,39 @@ void test_VectorProduct()
     rowLengths_2.setElement( 1, 1 );
     rowLengths_2.setElement( 3, 1 );
     m_2.setCompressedRowLengths( rowLengths_2 );
-    
+
     RealType value_2 = 1;
     for( IndexType i = 0; i < 3; i++ )   // 0th row
         m_2.setElement( 0, i, value_2++ );
-    
+
     m_2.setElement( 1, 3, value_2++ );      // 1st row
-        
+
     for( IndexType i = 0; i < 3; i++ )   // 2nd row
         m_2.setElement( 2, i, value_2++ );
-        
+
     for( IndexType i = 1; i < 2; i++ )       // 3rd row
         m_2.setElement( 3, i, value_2++ );
-    
+
     VectorType inVector_2;
     inVector_2.setSize( m_cols_2 );
-    for( IndexType i = 0; i < inVector_2.getSize(); i++ )        
+    for( IndexType i = 0; i < inVector_2.getSize(); i++ )
         inVector_2.setElement( i, 2 );
 
-    VectorType outVector_2;  
+    VectorType outVector_2;
     outVector_2.setSize( m_rows_2 );
     for( IndexType j = 0; j < outVector_2.getSize(); j++ )
         outVector_2.setElement( j, 0 );
- 
-    
+
+
     m_2.vectorProduct( inVector_2, outVector_2 );
-    
-   
+
+
     EXPECT_EQ( outVector_2.getElement( 0 ), 12 );
     EXPECT_EQ( outVector_2.getElement( 1 ),  8 );
     EXPECT_EQ( outVector_2.getElement( 2 ), 36 );
     EXPECT_EQ( outVector_2.getElement( 3 ), 16 );
-    
-    
+
+
 /*
  * Sets up the following 4x4 sparse matrix:
  *
@@ -785,10 +785,10 @@ void test_VectorProduct()
  *    |  7  8  9  0 |
  *    \  0 10 11 12 /
  */
-    
+
     const IndexType m_rows_3 = 4;
     const IndexType m_cols_3 = 4;
-    
+
     Matrix m_3;
     m_3.reset();
     m_3.setDimensions( m_rows_3, m_cols_3 );
@@ -796,40 +796,40 @@ void test_VectorProduct()
     rowLengths_3.setSize( m_rows_3 );
     rowLengths_3.setValue( 3 );
     m_3.setCompressedRowLengths( rowLengths_3 );
-    
+
     RealType value_3 = 1;
     for( IndexType i = 0; i < 3; i++ )          // 0th row
         m_3.setElement( 0, i, value_3++ );
-    
+
     for( IndexType i = 1; i < 4; i++ )
         m_3.setElement( 1, i, value_3++ );      // 1st row
-        
+
     for( IndexType i = 0; i < 3; i++ )          // 2nd row
         m_3.setElement( 2, i, value_3++ );
-        
+
     for( IndexType i = 1; i < 4; i++ )          // 3rd row
         m_3.setElement( 3, i, value_3++ );
-    
+
     VectorType inVector_3;
     inVector_3.setSize( m_cols_3 );
-    for( IndexType i = 0; i < inVector_3.getSize(); i++ )        
+    for( IndexType i = 0; i < inVector_3.getSize(); i++ )
         inVector_3.setElement( i, 2 );
 
-    VectorType outVector_3;  
+    VectorType outVector_3;
     outVector_3.setSize( m_rows_3 );
     for( IndexType j = 0; j < outVector_3.getSize(); j++ )
         outVector_3.setElement( j, 0 );
- 
-    
+
+
     m_3.vectorProduct( inVector_3, outVector_3 );
-    
-   
+
+
     EXPECT_EQ( outVector_3.getElement( 0 ), 12 );
     EXPECT_EQ( outVector_3.getElement( 1 ), 30 );
     EXPECT_EQ( outVector_3.getElement( 2 ), 48 );
     EXPECT_EQ( outVector_3.getElement( 3 ), 66 );
-    
-    
+
+
 /*
  * Sets up the following 8x8 sparse matrix:
  *
@@ -842,10 +842,10 @@ void test_VectorProduct()
  *    | 26 27 28 29 30  0  0  0 |
  *    \ 31 32 33 34 35  0  0  0 /
  */
-    
+
     const IndexType m_rows_4 = 8;
     const IndexType m_cols_4 = 8;
-    
+
     Matrix m_4;
     m_4.reset();
     m_4.setDimensions( m_rows_4, m_cols_4 );
@@ -856,48 +856,48 @@ void test_VectorProduct()
     rowLengths_4.setElement( 6, 5 );
     rowLengths_4.setElement( 7, 5 );
     m_4.setCompressedRowLengths( rowLengths_4 );
-    
+
     RealType value_4 = 1;
     for( IndexType i = 0; i < 3; i++ )       // 0th row
         m_4.setElement( 0, i, value_4++ );
-    
+
     m_4.setElement( 0, 5, value_4++ );
-    
+
     for( IndexType i = 1; i < 5; i++ )       // 1st row
         m_4.setElement( 1, i, value_4++ );
-    
+
     for( IndexType i = 0; i < 5; i++ )       // 2nd row
         m_4.setElement( 2, i, value_4++ );
-    
+
     for( IndexType i = 1; i < 5; i++ )       // 3rd row
         m_4.setElement( 3, i, value_4++ );
-    
+
     for( IndexType i = 2; i < 6; i++ )       // 4th row
         m_4.setElement( 4, i, value_4++ );
-    
+
     for( IndexType i = 3; i < 7; i++ )       // 5th row
         m_4.setElement( 5, i, value_4++ );
-    
+
     for( IndexType i = 0; i < 5; i++ )       // 6th row
         m_4.setElement( 6, i, value_4++ );
-    
+
     for( IndexType i = 0; i < 5; i++ )       // 7th row
         m_4.setElement( 7, i, value_4++ );
-    
+
     VectorType inVector_4;
     inVector_4.setSize( m_cols_4 );
-    for( IndexType i = 0; i < inVector_4.getSize(); i++ )        
+    for( IndexType i = 0; i < inVector_4.getSize(); i++ )
         inVector_4.setElement( i, 2 );
 
-    VectorType outVector_4;  
+    VectorType outVector_4;
     outVector_4.setSize( m_rows_4 );
     for( IndexType j = 0; j < outVector_4.getSize(); j++ )
         outVector_4.setElement( j, 0 );
-    
-    
+
+
     m_4.vectorProduct( inVector_4, outVector_4 );
-    
-   
+
+
     EXPECT_EQ( outVector_4.getElement( 0 ),  20 );
     EXPECT_EQ( outVector_4.getElement( 1 ),  52 );
     EXPECT_EQ( outVector_4.getElement( 2 ), 110 );
@@ -906,8 +906,8 @@ void test_VectorProduct()
     EXPECT_EQ( outVector_4.getElement( 5 ), 188 );
     EXPECT_EQ( outVector_4.getElement( 6 ), 280 );
     EXPECT_EQ( outVector_4.getElement( 7 ), 330 );
-    
-  
+
+
    /*
     * Sets up the following 8x8 sparse matrix:
     *
@@ -976,7 +976,7 @@ void test_VectorProduct()
     for( IndexType i = 0; i < inVector_5.getSize(); i++ )
         inVector_5.setElement( i, 2 );
 
-    VectorType outVector_5;  
+    VectorType outVector_5;
     outVector_5.setSize( m_rows_5 );
     for( IndexType j = 0; j < outVector_5.getSize(); j++ )
         outVector_5.setElement( j, 0 );
@@ -1077,6 +1077,8 @@ void test_RowsReduction()
    };
    m.allRowsReduction( fetch, reduce, keep, 0 );
    EXPECT_EQ( rowsCapacities, rowLengths );
+   m.getCompressedRowLengths( rowLengths );
+   EXPECT_EQ( rowsCapacities, rowLengths );
 
    ////
    // Compute max norm
@@ -1102,7 +1104,7 @@ 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:
  *
@@ -1111,10 +1113,10 @@ void test_PerformSORIteration()
  *    |  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 );
@@ -1122,54 +1124,54 @@ void test_PerformSORIteration()
     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 );
@@ -1183,7 +1185,7 @@ void test_OperatorEquals()
    using RealType = typename Matrix::RealType;
    using DeviceType = typename Matrix::DeviceType;
    using IndexType = typename Matrix::IndexType;
-   
+
    if( std::is_same< DeviceType, TNL::Devices::Cuda >::value )
        return;
    else
@@ -1229,33 +1231,33 @@ void test_OperatorEquals()
 
         m_host.setElement( 0, 4, value++ );           // 0th row
         m_host.setElement( 0, 5, value++ );
-        
+
         m_host.setElement( 1, 1, value++ );           // 1st row
         m_host.setElement( 1, 3, value++ );
 
         for( IndexType i = 1; i < 3; i++ )            // 2nd row
             m_host.setElement( 2, i, value++ );
-        
+
         m_host.setElement( 2, 4, value++ );           // 2nd row
 
-        
+
         for( IndexType i = 1; i < 5; i++ )            // 3rd row
             m_host.setElement( 3, i, value++ );
 
         m_host.setElement( 4, 1, value++ );           // 4th row
-        
+
         for( IndexType i = 1; i < 7; i++ )            // 5th row
             m_host.setElement( 5, i, value++ );
-        
+
         for( IndexType i = 0; i < 7; i++ )            // 6th row
             m_host.setElement( 6, i, value++ );
-        
+
         for( IndexType i = 0; i < 8; i++ )            // 7th row
             m_host.setElement( 7, i, value++ );
-        
+
         for( IndexType i = 0; i < 7; i++ )            // 1s at the end or rows: 5, 6
             m_host.setElement( i, 7, 1);
-        
+
         EXPECT_EQ( m_host.getElement( 0, 0 ),  1 );
         EXPECT_EQ( m_host.getElement( 0, 1 ),  2 );
         EXPECT_EQ( m_host.getElement( 0, 2 ),  3 );
@@ -1264,7 +1266,7 @@ void test_OperatorEquals()
         EXPECT_EQ( m_host.getElement( 0, 5 ),  5 );
         EXPECT_EQ( m_host.getElement( 0, 6 ),  0 );
         EXPECT_EQ( m_host.getElement( 0, 7 ),  1 );
-        
+
         EXPECT_EQ( m_host.getElement( 1, 0 ),  0 );
         EXPECT_EQ( m_host.getElement( 1, 1 ),  6 );
         EXPECT_EQ( m_host.getElement( 1, 2 ),  0 );
@@ -1273,7 +1275,7 @@ void test_OperatorEquals()
         EXPECT_EQ( m_host.getElement( 1, 5 ),  0 );
         EXPECT_EQ( m_host.getElement( 1, 6 ),  0 );
         EXPECT_EQ( m_host.getElement( 1, 7 ),  1 );
-        
+
         EXPECT_EQ( m_host.getElement( 2, 0 ),  0 );
         EXPECT_EQ( m_host.getElement( 2, 1 ),  8 );
         EXPECT_EQ( m_host.getElement( 2, 2 ),  9 );
@@ -1282,7 +1284,7 @@ void test_OperatorEquals()
         EXPECT_EQ( m_host.getElement( 2, 5 ),  0 );
         EXPECT_EQ( m_host.getElement( 2, 6 ),  0 );
         EXPECT_EQ( m_host.getElement( 2, 7 ),  1 );
-        
+
         EXPECT_EQ( m_host.getElement( 3, 0 ),  0 );
         EXPECT_EQ( m_host.getElement( 3, 1 ), 11 );
         EXPECT_EQ( m_host.getElement( 3, 2 ), 12 );
@@ -1291,7 +1293,7 @@ void test_OperatorEquals()
         EXPECT_EQ( m_host.getElement( 3, 5 ),  0 );
         EXPECT_EQ( m_host.getElement( 3, 6 ),  0 );
         EXPECT_EQ( m_host.getElement( 3, 7 ),  1 );
-        
+
         EXPECT_EQ( m_host.getElement( 4, 0 ),  0 );
         EXPECT_EQ( m_host.getElement( 4, 1 ), 15 );
         EXPECT_EQ( m_host.getElement( 4, 2 ),  0 );
@@ -1300,7 +1302,7 @@ void test_OperatorEquals()
         EXPECT_EQ( m_host.getElement( 4, 5 ),  0 );
         EXPECT_EQ( m_host.getElement( 4, 6 ),  0 );
         EXPECT_EQ( m_host.getElement( 4, 7 ),  1 );
-        
+
         EXPECT_EQ( m_host.getElement( 5, 0 ),  0 );
         EXPECT_EQ( m_host.getElement( 5, 1 ), 16 );
         EXPECT_EQ( m_host.getElement( 5, 2 ), 17 );
@@ -1309,7 +1311,7 @@ void test_OperatorEquals()
         EXPECT_EQ( m_host.getElement( 5, 5 ), 20 );
         EXPECT_EQ( m_host.getElement( 5, 6 ), 21 );
         EXPECT_EQ( m_host.getElement( 5, 7 ),  1 );
-        
+
         EXPECT_EQ( m_host.getElement( 6, 0 ), 22 );
         EXPECT_EQ( m_host.getElement( 6, 1 ), 23 );
         EXPECT_EQ( m_host.getElement( 6, 2 ), 24 );
@@ -1318,7 +1320,7 @@ void test_OperatorEquals()
         EXPECT_EQ( m_host.getElement( 6, 5 ), 27 );
         EXPECT_EQ( m_host.getElement( 6, 6 ), 28 );
         EXPECT_EQ( m_host.getElement( 6, 7 ),  1 );
-        
+
         EXPECT_EQ( m_host.getElement( 7, 0 ), 29 );
         EXPECT_EQ( m_host.getElement( 7, 1 ), 30 );
         EXPECT_EQ( m_host.getElement( 7, 2 ), 31 );
@@ -1348,7 +1350,7 @@ void test_OperatorEquals()
         EXPECT_EQ( m_host.getElement( 0, 5 ),  5 );
         EXPECT_EQ( m_host.getElement( 0, 6 ),  0 );
         EXPECT_EQ( m_host.getElement( 0, 7 ),  1 );
-        
+
         EXPECT_EQ( m_host.getElement( 1, 0 ),  0 );
         EXPECT_EQ( m_host.getElement( 1, 1 ),  6 );
         EXPECT_EQ( m_host.getElement( 1, 2 ),  0 );
@@ -1357,7 +1359,7 @@ void test_OperatorEquals()
         EXPECT_EQ( m_host.getElement( 1, 5 ),  0 );
         EXPECT_EQ( m_host.getElement( 1, 6 ),  0 );
         EXPECT_EQ( m_host.getElement( 1, 7 ),  1 );
-        
+
         EXPECT_EQ( m_host.getElement( 2, 0 ),  0 );
         EXPECT_EQ( m_host.getElement( 2, 1 ),  8 );
         EXPECT_EQ( m_host.getElement( 2, 2 ),  9 );
@@ -1366,7 +1368,7 @@ void test_OperatorEquals()
         EXPECT_EQ( m_host.getElement( 2, 5 ),  0 );
         EXPECT_EQ( m_host.getElement( 2, 6 ),  0 );
         EXPECT_EQ( m_host.getElement( 2, 7 ),  1 );
-        
+
         EXPECT_EQ( m_host.getElement( 3, 0 ),  0 );
         EXPECT_EQ( m_host.getElement( 3, 1 ), 11 );
         EXPECT_EQ( m_host.getElement( 3, 2 ), 12 );
@@ -1375,7 +1377,7 @@ void test_OperatorEquals()
         EXPECT_EQ( m_host.getElement( 3, 5 ),  0 );
         EXPECT_EQ( m_host.getElement( 3, 6 ),  0 );
         EXPECT_EQ( m_host.getElement( 3, 7 ),  1 );
-        
+
         EXPECT_EQ( m_host.getElement( 4, 0 ),  0 );
         EXPECT_EQ( m_host.getElement( 4, 1 ), 15 );
         EXPECT_EQ( m_host.getElement( 4, 2 ),  0 );
@@ -1384,7 +1386,7 @@ void test_OperatorEquals()
         EXPECT_EQ( m_host.getElement( 4, 5 ),  0 );
         EXPECT_EQ( m_host.getElement( 4, 6 ),  0 );
         EXPECT_EQ( m_host.getElement( 4, 7 ),  1 );
-        
+
         EXPECT_EQ( m_host.getElement( 5, 0 ),  0 );
         EXPECT_EQ( m_host.getElement( 5, 1 ), 16 );
         EXPECT_EQ( m_host.getElement( 5, 2 ), 17 );
@@ -1393,7 +1395,7 @@ void test_OperatorEquals()
         EXPECT_EQ( m_host.getElement( 5, 5 ), 20 );
         EXPECT_EQ( m_host.getElement( 5, 6 ), 21 );
         EXPECT_EQ( m_host.getElement( 5, 7 ),  1 );
-        
+
         EXPECT_EQ( m_host.getElement( 6, 0 ), 22 );
         EXPECT_EQ( m_host.getElement( 6, 1 ), 23 );
         EXPECT_EQ( m_host.getElement( 6, 2 ), 24 );
@@ -1402,7 +1404,7 @@ void test_OperatorEquals()
         EXPECT_EQ( m_host.getElement( 6, 5 ), 27 );
         EXPECT_EQ( m_host.getElement( 6, 6 ), 28 );
         EXPECT_EQ( m_host.getElement( 6, 7 ),  1 );
-        
+
         EXPECT_EQ( m_host.getElement( 7, 0 ), 29 );
         EXPECT_EQ( m_host.getElement( 7, 1 ), 30 );
         EXPECT_EQ( m_host.getElement( 7, 2 ), 31 );
@@ -1411,22 +1413,22 @@ void test_OperatorEquals()
         EXPECT_EQ( m_host.getElement( 7, 5 ), 34 );
         EXPECT_EQ( m_host.getElement( 7, 6 ), 35 );
         EXPECT_EQ( m_host.getElement( 7, 7 ), 36 );
-        
+
         // Try vectorProduct with copied cuda matrix to see if it works correctly.
         using VectorType = TNL::Containers::Vector< RealType, TNL::Devices::Cuda, IndexType >;
-    
+
         VectorType inVector;
         inVector.setSize( m_cols );
-        for( IndexType i = 0; i < inVector.getSize(); i++ )        
+        for( IndexType i = 0; i < inVector.getSize(); i++ )
             inVector.setElement( i, 2 );
 
-        VectorType outVector;  
+        VectorType outVector;
         outVector.setSize( m_rows );
         for( IndexType j = 0; j < outVector.getSize(); j++ )
             outVector.setElement( j, 0 );
-        
+
         m_cuda.vectorProduct( inVector, outVector );
-        
+
         EXPECT_EQ( outVector.getElement( 0 ),  32 );
         EXPECT_EQ( outVector.getElement( 1 ),  28 );
         EXPECT_EQ( outVector.getElement( 2 ),  56 );
@@ -1444,7 +1446,7 @@ void test_SaveAndLoad( const char* filename )
    using RealType = typename Matrix::RealType;
    using DeviceType = typename Matrix::DeviceType;
    using IndexType = typename Matrix::IndexType;
-    
+
    /*
     * Sets up the following 4x4 sparse matrix:
     *
@@ -1453,10 +1455,10 @@ void test_SaveAndLoad( const char* filename )
     *    |  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 );
@@ -1464,22 +1466,22 @@ void test_SaveAndLoad( const char* filename )
     rowLengths.setSize( m_rows );
     rowLengths.setValue( 3 );
     savedMatrix.setCompressedRowLengths( rowLengths );
-    
+
     RealType value = 1;
     for( IndexType i = 0; i < m_cols - 1; i++ )   // 0th row
         savedMatrix.setElement( 0, i, value++ );
-        
+
     savedMatrix.setElement( 1, 1, value++ );
     savedMatrix.setElement( 1, 3, value++ );      // 1st row
-        
+
     for( IndexType i = 0; i < m_cols - 1; i++ )   // 2nd row
         savedMatrix.setElement( 2, i, value++ );
-        
+
     for( IndexType i = 1; i < m_cols; i++ )       // 3rd row
         savedMatrix.setElement( 3, i, value++ );
-        
+
     ASSERT_NO_THROW( savedMatrix.save( filename ) );
-    
+
     Matrix loadedMatrix;
     loadedMatrix.reset();
     loadedMatrix.setDimensions( m_rows, m_cols );
@@ -1487,51 +1489,51 @@ void test_SaveAndLoad( const char* filename )
     rowLengths2.setSize( m_rows );
     rowLengths2.setValue( 3 );
     loadedMatrix.setCompressedRowLengths( rowLengths2 );
-    
-    
+
+
     ASSERT_NO_THROW( loadedMatrix.load( filename ) );
-    
-    
+
+
     EXPECT_EQ( savedMatrix.getElement( 0, 0 ), loadedMatrix.getElement( 0, 0 ) );
     EXPECT_EQ( savedMatrix.getElement( 0, 1 ), loadedMatrix.getElement( 0, 1 ) );
     EXPECT_EQ( savedMatrix.getElement( 0, 2 ), loadedMatrix.getElement( 0, 2 ) );
     EXPECT_EQ( savedMatrix.getElement( 0, 3 ), loadedMatrix.getElement( 0, 3 ) );
-    
+
     EXPECT_EQ( savedMatrix.getElement( 1, 0 ), loadedMatrix.getElement( 1, 0 ) );
     EXPECT_EQ( savedMatrix.getElement( 1, 1 ), loadedMatrix.getElement( 1, 1 ) );
     EXPECT_EQ( savedMatrix.getElement( 1, 2 ), loadedMatrix.getElement( 1, 2 ) );
     EXPECT_EQ( savedMatrix.getElement( 1, 3 ), loadedMatrix.getElement( 1, 3 ) );
-    
+
     EXPECT_EQ( savedMatrix.getElement( 2, 0 ), loadedMatrix.getElement( 2, 0 ) );
     EXPECT_EQ( savedMatrix.getElement( 2, 1 ), loadedMatrix.getElement( 2, 1 ) );
     EXPECT_EQ( savedMatrix.getElement( 2, 2 ), loadedMatrix.getElement( 2, 2 ) );
     EXPECT_EQ( savedMatrix.getElement( 2, 3 ), loadedMatrix.getElement( 2, 3 ) );
-    
+
     EXPECT_EQ( savedMatrix.getElement( 3, 0 ), loadedMatrix.getElement( 3, 0 ) );
     EXPECT_EQ( savedMatrix.getElement( 3, 1 ), loadedMatrix.getElement( 3, 1 ) );
     EXPECT_EQ( savedMatrix.getElement( 3, 2 ), loadedMatrix.getElement( 3, 2 ) );
     EXPECT_EQ( savedMatrix.getElement( 3, 3 ), loadedMatrix.getElement( 3, 3 ) );
-    
+
     EXPECT_EQ( savedMatrix.getElement( 0, 0 ),  1 );
     EXPECT_EQ( savedMatrix.getElement( 0, 1 ),  2 );
     EXPECT_EQ( savedMatrix.getElement( 0, 2 ),  3 );
     EXPECT_EQ( savedMatrix.getElement( 0, 3 ),  0 );
-    
+
     EXPECT_EQ( savedMatrix.getElement( 1, 0 ),  0 );
     EXPECT_EQ( savedMatrix.getElement( 1, 1 ),  4 );
     EXPECT_EQ( savedMatrix.getElement( 1, 2 ),  0 );
     EXPECT_EQ( savedMatrix.getElement( 1, 3 ),  5 );
-    
+
     EXPECT_EQ( savedMatrix.getElement( 2, 0 ),  6 );
     EXPECT_EQ( savedMatrix.getElement( 2, 1 ),  7 );
     EXPECT_EQ( savedMatrix.getElement( 2, 2 ),  8 );
     EXPECT_EQ( savedMatrix.getElement( 2, 3 ),  0 );
-    
+
     EXPECT_EQ( savedMatrix.getElement( 3, 0 ),  0 );
     EXPECT_EQ( savedMatrix.getElement( 3, 1 ),  9 );
     EXPECT_EQ( savedMatrix.getElement( 3, 2 ), 10 );
     EXPECT_EQ( savedMatrix.getElement( 3, 3 ), 11 );
-    
+
     EXPECT_EQ( std::remove( filename ), 0 );
 }
 
@@ -1541,7 +1543,7 @@ 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:
  *
@@ -1551,10 +1553,10 @@ void test_Print()
  *    |  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 );
@@ -1562,40 +1564,40 @@ void test_Print()
     rowLengths.setSize( m_rows );
     rowLengths.setValue( 3 );
     m.setCompressedRowLengths( rowLengths );
-    
+
     RealType value = 1;
     for( IndexType i = 0; i < m_cols - 1; i++ )   // 0th row
         m.setElement( 0, i, value++ );
-    
+
     m.setElement( 1, 3, value++ );      // 1st row
-        
+
     for( IndexType i = 0; i < m_cols - 1; i++ )   // 2nd row
         m.setElement( 2, i, value++ );
-        
+
     for( IndexType i = 1; i < m_cols; i++ )       // 3rd row
         m.setElement( 3, i, value++ );
-        
+
     for( IndexType i = 2; i < m_cols; i++ )       // 4th row
         m.setElement( 4, i, value++ );
-    
+
     #include <sstream>
     std::stringstream printed;
     std::stringstream couted;
-    
+
     //change the underlying buffer and save the old buffer
-    auto old_buf = std::cout.rdbuf(printed.rdbuf()); 
+    auto old_buf = std::cout.rdbuf(printed.rdbuf());
 
     m.print( std::cout ); //all the std::cout goes to ss
 
     std::cout.rdbuf(old_buf); //reset
-    
+
     couted << "Row: 0 ->  Col:0->1	 Col:1->2	 Col:2->3\t\n"
                "Row: 1 ->  Col:3->4\t\n"
                "Row: 2 ->  Col:0->5	 Col:1->6	 Col:2->7\t\n"
                "Row: 3 ->  Col:1->8	 Col:2->9	 Col:3->10\t\n"
                "Row: 4 ->  Col:2->11	 Col:3->12\t\n";
-    
-    
+
+
     EXPECT_EQ( printed.str(), couted.str() );
 }
 
-- 
GitLab