diff --git a/src/TNL/Matrices/Matrix.h b/src/TNL/Matrices/Matrix.h
index a9b458d7b735cdf8b5fda2217bfe32cfd8cbbf28..7813fa9625ce25b7d4d0890d04df13ffb7a4eeba 100644
--- a/src/TNL/Matrices/Matrix.h
+++ b/src/TNL/Matrices/Matrix.h
@@ -64,7 +64,7 @@ public:
    template< typename Matrix_ >
    void setLike( const Matrix_& matrix );
 
-   IndexType getNumberOfMatrixElements() const;
+   IndexType getAllocatedElementsCount() const;
 
    virtual IndexType getNumberOfNonzeroMatrixElements() const = 0;
 
diff --git a/src/TNL/Matrices/Matrix.hpp b/src/TNL/Matrices/Matrix.hpp
index 29226cb00fa72d968481bb565bdc6bd8302e4083..efd26e1fa6368d807a86b58803800a4f239dca14 100644
--- a/src/TNL/Matrices/Matrix.hpp
+++ b/src/TNL/Matrices/Matrix.hpp
@@ -91,7 +91,7 @@ template< typename Real,
           typename Device,
           typename Index,
           typename RealAllocator >
-Index Matrix< Real, Device, Index, RealAllocator >::getNumberOfMatrixElements() const
+Index Matrix< Real, Device, Index, RealAllocator >::getAllocatedElementsCount() const
 {
    return this->values.getSize();
 }
diff --git a/src/TNL/Matrices/MatrixView.h b/src/TNL/Matrices/MatrixView.h
index 76965e511291bfae5b178970f869b1d50814fec6..b8adfd791cb9b8f069568aed8d31a6b94e0b45ee 100644
--- a/src/TNL/Matrices/MatrixView.h
+++ b/src/TNL/Matrices/MatrixView.h
@@ -57,12 +57,10 @@ public:
 
    virtual void getCompressedRowLengths( CompressedRowLengthsVectorView rowLengths ) const;
 
-   IndexType getNumberOfMatrixElements() const;
+   IndexType getAllocatedElementsCount() const;
 
    virtual IndexType getNumberOfNonzeroMatrixElements() const;
 
-   void reset();
-
    __cuda_callable__
    IndexType getRows() const;
 
@@ -91,6 +89,15 @@ public:
 
    ValuesView& getValues();
 
+   /**
+    * \brief Shallow copy of the matrix view.
+    *
+    * @param view
+    * @return 
+    */
+   __cuda_callable__
+   MatrixView& operator=( const MatrixView& view );
+   
    // TODO: parallelize and optimize for sparse matrices
    template< typename Matrix >
    bool operator == ( const Matrix& matrix ) const;
diff --git a/src/TNL/Matrices/MatrixView.hpp b/src/TNL/Matrices/MatrixView.hpp
index bd3d9beaee73a8d55bc73c1e9de870a4c14b330c..b2739ae1dc9c458ddfe093700159ee4328054e8f 100644
--- a/src/TNL/Matrices/MatrixView.hpp
+++ b/src/TNL/Matrices/MatrixView.hpp
@@ -64,7 +64,7 @@ void MatrixView< Real, Device, Index >::getCompressedRowLengths( CompressedRowLe
 template< typename Real,
           typename Device,
           typename Index >
-Index MatrixView< Real, Device, Index >::getNumberOfMatrixElements() const
+Index MatrixView< Real, Device, Index >::getAllocatedElementsCount() const
 {
    return this->values.getSize();
 }
@@ -118,15 +118,17 @@ getValues()
 {
    return this->values;
 }
-
 template< typename Real,
           typename Device,
           typename Index >
-void MatrixView< Real, Device, Index >::reset()
+__cuda_callable__
+MatrixView< Real, Device, Index >& 
+MatrixView< Real, Device, Index >::
+operator=( const MatrixView& view )
 {
-   this->rows = 0;
-   this->columns = 0;
-   this->values.reset();
+   rows = view.rows;
+   columns = view.columns;
+   values.copy( view.values );
 }
 
 template< typename Real,
diff --git a/src/TNL/Matrices/Tridiagonal.h b/src/TNL/Matrices/Tridiagonal.h
index 51e05c8992d2e7657bf488c10e096eebc91ada23..d2827015695da2183c16c152e10d1d711149b51e 100644
--- a/src/TNL/Matrices/Tridiagonal.h
+++ b/src/TNL/Matrices/Tridiagonal.h
@@ -79,12 +79,8 @@ class Tridiagonal : public Matrix< Real, Device, Index, RealAllocator >
       template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_, typename RealAllocator_ >
       void setLike( const Tridiagonal< Real_, Device_, Index_, RowMajorOrder_, RealAllocator_ >& m );
 
-      IndexType getNumberOfMatrixElements() const;
-
       IndexType getNumberOfNonzeroMatrixElements() const;
 
-      IndexType getMaxRowlength() const;
-
       void reset();
 
       template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_, typename RealAllocator_ >
@@ -179,6 +175,8 @@ class Tridiagonal : public Matrix< Real, Device, Index, RealAllocator >
                                  const IndexType localIdx ) const;
 
       IndexerType indexer;
+
+      ViewType view;
 };
 
 } // namespace Matrices
diff --git a/src/TNL/Matrices/Tridiagonal.hpp b/src/TNL/Matrices/Tridiagonal.hpp
index a7178f86ed7311f6a623b97a5f9c7a3b4d643515..c6d359d3be42c618eaf28c33a0b320f5d0dee1f7 100644
--- a/src/TNL/Matrices/Tridiagonal.hpp
+++ b/src/TNL/Matrices/Tridiagonal.hpp
@@ -10,6 +10,7 @@
 
 #pragma once
 
+#include <sstream>
 #include <TNL/Assert.h>
 #include <TNL/Matrices/Tridiagonal.h>
 #include <TNL/Exceptions/NotImplementedError.h>
@@ -105,6 +106,7 @@ setDimensions( const IndexType rows, const IndexType columns )
    this->indexer.setDimensions( rows, columns );
    this->values.setSize( this->indexer.getStorageSize() );
    this->values = 0.0;
+   this->view = this->getView();
 }
 
 template< typename Real,
@@ -138,20 +140,25 @@ template< typename Real,
           typename Index,
           bool RowMajorOrder,
           typename RealAllocator >
-Index
+   template< typename Vector >
+void
 Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
-getRowLength( const IndexType row ) const
+getCompressedRowLengths( Vector& rowLengths ) const
 {
-   const IndexType diagonalLength = min( this->getRows(), this->getColumns() );
-   if( row == 0 )
-      return 2;
-   if( row > 0 && row < diagonalLength - 1 )
-      return 3;
-   if( this->getRows() > this->getColumns() )
-      return 1;
-   if( this->getRows() == this->getColumns() )
-      return 2;
-   return 3;
+   return this->view.getCompressedRowLengths( rowLengths );
+   /*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,
@@ -161,9 +168,10 @@ template< typename Real,
           typename RealAllocator >
 Index
 Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
-getMaxRowLength() const
+getRowLength( const IndexType row ) const
 {
-   return 3;
+   return this->view.getRowLength( row );
+   //return this->indexer.getRowSize( row );
 }
 
 template< typename Real,
@@ -171,12 +179,11 @@ template< typename Real,
           typename Index,
           bool RowMajorOrder,
           typename RealAllocator >
-   template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_, typename RealAllocator_ >
-void
+Index
 Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
-setLike( const Tridiagonal< Real_, Device_, Index_, RowMajorOrder_, RealAllocator_ >& m )
+getMaxRowLength() const
 {
-   this->setDimensions( m.getRows(), m.getColumns() );
+   return this->view.getMaxRowLength();
 }
 
 template< typename Real,
@@ -184,11 +191,12 @@ template< typename Real,
           typename Index,
           bool RowMajorOrder,
           typename RealAllocator >
-Index
+   template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_, typename RealAllocator_ >
+void
 Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
-getNumberOfMatrixElements() const
+setLike( const Tridiagonal< Real_, Device_, Index_, RowMajorOrder_, RealAllocator_ >& m )
 {
-   return 3 * min( this->getRows(), this->getColumns() );
+   this->setDimensions( m.getRows(), m.getColumns() );
 }
 
 template< typename Real,
@@ -200,23 +208,12 @@ Index
 Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
 getNumberOfNonzeroMatrixElements() const
 {
-   const auto values_view = this->values.getConstView();
+   return this->view.getNumberOfNonzeroMatrixElements();
+   /*const auto values_view = this->values.getConstView();
    auto fetch = [=] __cuda_callable__ ( const IndexType i ) -> IndexType {
       return ( values_view[ i ] != 0.0 );
    };
-   return Algorithms::Reduction< DeviceType >::reduce( this->values.getSize(), std::plus<>{}, fetch, 0 );
-}
-
-template< typename Real,
-          typename Device,
-          typename Index,
-          bool RowMajorOrder,
-          typename RealAllocator >
-Index
-Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
-getMaxRowlength() const
-{
-   return 3;
+   return Algorithms::Reduction< DeviceType >::reduce( this->values.getSize(), std::plus<>{}, fetch, 0 );*/
 }
 
 template< typename Real,
@@ -272,7 +269,7 @@ void
 Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
 setValue( const RealType& v )
 {
-   this->values = v;
+   this->view.setValue( v );
 }
 
 template< typename Real,
@@ -285,7 +282,8 @@ auto
 Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
 getRow( const IndexType& rowIdx ) const -> const RowView
 {
-   return RowView( this->values.getView(), this->indexer );
+   return this->view.getRow( rowIdx );
+   //return RowView( this->values.getView(), this->indexer );
 }
 
 template< typename Real,
@@ -298,7 +296,8 @@ auto
 Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
 getRow( const IndexType& rowIdx ) -> RowView
 {
-   return RowView( this->values.getView(), this->indexer );
+   return this->view.getRow( rowIdx );
+   //return RowView( this->values.getView(), this->indexer );
 }
 
 template< typename Real,
@@ -310,14 +309,19 @@ bool
 Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
 setElement( const IndexType row, const IndexType column, const RealType& value )
 {
-   TNL_ASSERT_GE( row, 0, "" );
+   return this->view.setElement( row, column, value );
+   /*TNL_ASSERT_GE( row, 0, "" );
    TNL_ASSERT_LT( row, this->getRows(), "" );
    TNL_ASSERT_GE( column, 0, "" );
    TNL_ASSERT_LT( column, this->getColumns(), "" );
    if( abs( row - column ) > 1 )
-      throw std::logic_error( "Wrong matrix element coordinates in tridiagonal matrix." );
+   {
+      std::stringstream msg;
+      msg << "Wrong matrix element coordinates ( "  << row << ", " << column << " ) in tridiagonal matrix.";
+      throw std::logic_error( msg.str() );
+   }
    this->values.setElement( this->getElementIndex( row, column ), value );
-   return true;
+   return true;*/
 }
 
 template< typename Real,
@@ -332,15 +336,20 @@ addElement( const IndexType row,
             const RealType& value,
             const RealType& thisElementMultiplicator )
 {
-   TNL_ASSERT_GE( row, 0, "" );
+   return this->view.addElement( row, column, value, thisElementMultiplicator );
+   /*TNL_ASSERT_GE( row, 0, "" );
    TNL_ASSERT_LT( row, this->getRows(), "" );
    TNL_ASSERT_GE( column, 0, "" );
    TNL_ASSERT_LT( column, this->getColumns(), "" );
    if( abs( row - column ) > 1 )
-      throw std::logic_error( "Wrong matrix element coordinates in tridiagonal matrix." );
+   {
+      std::stringstream msg;
+      msg << "Wrong matrix element coordinates ( "  << row << ", " << column << " ) in tridiagonal matrix.";
+      throw std::logic_error( msg.str() );
+   }
    const Index i = this->getElementIndex( row, column );
    this->values.setElement( i, thisElementMultiplicator * this->values.getElement( i ) + value );
-   return true;
+   return true;*/
 }
 
 template< typename Real,
@@ -352,14 +361,15 @@ Real
 Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
 getElement( const IndexType row, const IndexType column ) const
 {
-   TNL_ASSERT_GE( row, 0, "" );
+   return this->view.getElement( row, column );
+   /*TNL_ASSERT_GE( row, 0, "" );
    TNL_ASSERT_LT( row, this->getRows(), "" );
    TNL_ASSERT_GE( column, 0, "" );
    TNL_ASSERT_LT( column, this->getColumns(), "" );
 
    if( abs( column - row ) > 1 )
       return 0.0;
-   return this->values.getElement( this->getElementIndex( row, column ) );
+   return this->values.getElement( this->getElementIndex( row, column ) );*/
 }
 
 template< typename Real,
@@ -372,46 +382,40 @@ void
 Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
 rowsReduction( IndexType first, IndexType last, Fetch& fetch, Reduce& reduce, Keep& keep, const FetchReal& zero ) const
 {
+   this->view.rowsReduction( first, last, fetch, reduce, keep, zero );
+   /*using Real_ = decltype( fetch( IndexType(), IndexType(), RealType() ) );
    const auto values_view = this->values.getConstView();
-   const auto indexer_ = this->indexer;
-   const auto rows = this->getRows();
-   const auto columns = this->getColumns();
-   const auto size = this->size;
+   const auto indexer = this->indexer;
+   const auto zero = zero_;
    auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
-      //bool compute;
+      Real_ sum( zero );
       if( rowIdx == 0 )
       {
-         IndexType i_0 = indexer.getGlobalIndex( 0, 0 );
-         IndexType i_1 = indexer.getGlobalIndex( 0, 1 );
-         keep( 0, reduce( fetch( 0, 0, i_0, values_view[ i_0 ] ),
-                          fetch( 0, 1, i_1, values_view[ i_1 ] ) ) );
+         reduce( sum, fetch( 0, 0, values_view[ indexer.getGlobalIndex( 0, 0 ) ] ) );
+         reduce( sum, fetch( 0, 1, values_view[ indexer.getGlobalIndex( 0, 1 ) ] ) );
+         keep( 0, sum );
          return;
       }
-      if( rowIdx < size || columns > rows )
+      if( rowIdx < indexer.getSize() || indexer.getColumns() > indexer.getRows() )
       {
-         IndexType i_0 = indexer.getGlobalIndex( rowIdx, 0 );
-         IndexType i_1 = indexer.getGlobalIndex( rowIdx, 1 );
-         IndexType i_2 = indexer.getGlobalIndex( rowIdx, 2 );
-
-         keep( rowIdx, reduce( reduce( fetch( rowIdx, rowIdx - 1, i_0, values_view[ i_0 ] ),
-                                       fetch( rowIdx, rowIdx, i_1, values_view[ i_1 ] ) ),
-                               fetch( rowIdx, rowIdx + 1, i_2, values_view[ i_2] ) ) );
+         reduce( sum, fetch( rowIdx, rowIdx - 1, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ] ) );
+         reduce( sum, fetch( rowIdx, rowIdx,     values_view[ indexer.getGlobalIndex( rowIdx, 1 ) ] ) );
+         reduce( sum, fetch( rowIdx, rowIdx + 1, values_view[ indexer.getGlobalIndex( rowIdx, 2 ) ] ) );
+         keep( rowIdx, sum );
          return;
       }
-      if( rows == columns )
+      if( indexer.getRows() == indexer.getColumns() )
       {
-         IndexType i_0 = indexer.getGlobalIndex( rowIdx, 0 );
-         IndexType i_1 = indexer.getGlobalIndex( rowIdx, 1 );
-         keep( rowIdx, reduce( fetch( rowIdx, rowIdx - 1, i_0, values_view[ i_0 ] ),
-                               fetch( rowIdx, rowIdx, i_1, values_view[ i_1 ] ) ) );
+         reduce( sum, fetch( rowIdx, rowIdx - 1, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ] ) );
+         reduce( sum, fetch( rowIdx, rowIdx,     values_view[ indexer.getGlobalIndex( rowIdx, 1 ) ] ) );
+         keep( rowIdx, sum );
       }
       else
       {
-         IndexType i_0 = indexer.getGlobalIndex( rowIdx, 0 );
-         keep( rowIdx, fetch( rowIdx, rowIdx, i_0, values_view[ i_0 ] ) );
+         keep( rowIdx, fetch( rowIdx, rowIdx, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ] ) );
       }
    };
-   Algorithms::ParallelFor< DeviceType >::exec( first, last, f );
+   Algorithms::ParallelFor< DeviceType >::exec( first, last, f );*/
 }
 
 template< typename Real,
@@ -424,7 +428,7 @@ void
 Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
 allRowsReduction( Fetch& fetch, Reduce& reduce, Keep& keep, const FetchReal& zero ) const
 {
-   this->rowsReduction( 0, this->getRows(), fetch, reduce, keep, zero );
+   this->view.rowsReduction( 0, this->getRows(), fetch, reduce, keep, zero );
 }
 
 template< typename Real,
@@ -437,7 +441,8 @@ void
 Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
 forRows( IndexType first, IndexType last, Function& function ) const
 {
-   const auto values_view = this->values.getConstView();
+   this->view.forRows( first, last, function );
+   /*const auto values_view = this->values.getConstView();
    const auto indexer_ = this->indexer;
    const auto rows = this->getRows();
    const auto columns = this->getColumns();
@@ -475,7 +480,7 @@ forRows( IndexType first, IndexType last, Function& function ) const
          function( rowIdx, 0, rowIdx, values_view[ i_0 ] );
       }
    };
-   Algorithms::ParallelFor< DeviceType >::exec( first, last, f );
+   Algorithms::ParallelFor< DeviceType >::exec( first, last, f );*/
 }
 
 template< typename Real,
@@ -488,7 +493,8 @@ void
 Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
 forRows( IndexType first, IndexType last, Function& function )
 {
-   const auto values_view = this->values.getConstView();
+   this->view.forRows( first, last, function );
+   /*const auto values_view = this->values.getConstView();
    const auto indexer_ = this->indexer;
    const auto rows = this->getRows();
    const auto columns = this->getColumns();
@@ -526,7 +532,7 @@ forRows( IndexType first, IndexType last, Function& function )
          function( rowIdx, 0, rowIdx, values_view[ i_0 ] );
       }
    };
-   Algorithms::ParallelFor< DeviceType >::exec( first, last, f );
+   Algorithms::ParallelFor< DeviceType >::exec( first, last, f );*/
 }
 
 template< typename Real,
@@ -539,7 +545,7 @@ void
 Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
 forAllRows( Function& function ) const
 {
-   this->forRows( 0, this->getRows(), function );
+   this->view.forRows( 0, this->getRows(), function );
 }
 
 template< typename Real,
@@ -552,7 +558,7 @@ void
 Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
 forAllRows( Function& function )
 {
-   this->forRows( 0, this->getRows(), function );
+   this->view.forRows( 0, this->getRows(), function );
 }
 
 template< typename Real,
@@ -566,11 +572,12 @@ typename Vector::RealType
 Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
 rowVectorProduct( const IndexType row, const Vector& vector ) const
 {
-   return TridiagonalDeviceDependentCode< Device >::
+   return this->view.rowVectorProduct();
+   /*return TridiagonalDeviceDependentCode< Device >::
              rowVectorProduct( this->rows,
                                this->values,
                                row,
-                               vector );
+                               vector );*/
 }
 
 template< typename Real,
@@ -584,12 +591,13 @@ void
 Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
 vectorProduct( const InVector& inVector, OutVector& outVector ) const
 {
-   TNL_ASSERT( this->getColumns() == inVector.getSize(),
+   this->view.vectorProduct( inVector, outVector );
+   /*TNL_ASSERT( this->getColumns() == inVector.getSize(),
             std::cerr << "Matrix columns: " << this->getColumns() << std::endl
                  << "Vector size: " << inVector.getSize() << std::endl );
    TNL_ASSERT( this->getRows() == outVector.getSize(),
                std::cerr << "Matrix rows: " << this->getRows() << std::endl
-                    << "Vector size: " << outVector.getSize() << std::endl );
+                    << "Vector size: " << outVector.getSize() << std::endl );*/
 
    //DeviceDependentCode::vectorProduct( *this, inVector, outVector );
 }
@@ -815,12 +823,15 @@ template< typename Real,
           typename RealAllocator >
 __cuda_callable__
 Index Tridiagonal< Real, Device, Index, RowMajorOrder, RealAllocator >::
-getElementIndex( const IndexType row, const IndexType localIdx ) const
+getElementIndex( const IndexType row, const IndexType column ) const
 {
-   TNL_ASSERT_GE( row, 0, "" );
-   TNL_ASSERT_LT( row, this->getRows(), "" );
+   IndexType localIdx = column - row;
+   if( row > 0 )
+      localIdx++;
+
    TNL_ASSERT_GE( localIdx, 0, "" );
    TNL_ASSERT_LT( localIdx, 3, "" );
+
    return this->indexer.getGlobalIndex( row, localIdx );
 }
 
diff --git a/src/TNL/Matrices/TridiagonalMatrixView.h b/src/TNL/Matrices/TridiagonalMatrixView.h
index 05f7663c9f134c7ae468562b861bf8f56edc8a7f..78593acf511153ab4051be71f5895210cdafb93a 100644
--- a/src/TNL/Matrices/TridiagonalMatrixView.h
+++ b/src/TNL/Matrices/TridiagonalMatrixView.h
@@ -70,12 +70,8 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
 
       IndexType getMaxRowLength() const;
 
-      IndexType getNumberOfMatrixElements() const;
-
       IndexType getNumberOfNonzeroMatrixElements() const;
 
-      IndexType getMaxRowlength() const;
-
       template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_ >
       bool operator == ( const TridiagonalMatrixView< Real_, Device_, Index_, RowMajorOrder_ >& matrix ) const;
 
@@ -144,13 +140,6 @@ class TridiagonalMatrixView : public MatrixView< Real, Device, Index >
                                 Vector2& x,
                                 const RealType& omega = 1.0 ) const;
 
-      // copy assignment
-      TridiagonalMatrixView& operator=( const TridiagonalMatrixView& matrix );
-
-      // cross-device copy assignment
-      template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_ >
-      TridiagonalMatrixView& operator=( const TridiagonalMatrixView< Real_, Device_, Index_, RowMajorOrder_ >& matrix );
-
       void save( File& file ) const;
 
       void save( const String& fileName ) const;
diff --git a/src/TNL/Matrices/TridiagonalMatrixView.hpp b/src/TNL/Matrices/TridiagonalMatrixView.hpp
index ef893295e8396158798883e3de91006c07642548..83ff6035dd36b628f02e0ac4d324b890dfb05b82 100644
--- a/src/TNL/Matrices/TridiagonalMatrixView.hpp
+++ b/src/TNL/Matrices/TridiagonalMatrixView.hpp
@@ -87,31 +87,36 @@ template< typename Real,
           typename Device,
           typename Index,
           bool RowMajorOrder >
-Index
+   template< typename Vector >
+void
 TridiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
-getRowLength( const IndexType row ) const
+getCompressedRowLengths( Vector& rowLengths ) const
 {
-   const IndexType diagonalLength = min( this->getRows(), this->getColumns() );
-   if( row == 0 )
-      return 2;
-   if( row > 0 && row < diagonalLength - 1 )
-      return 3;
-   if( this->getRows() > this->getColumns() )
-      return 1;
-   if( this->getRows() == this->getColumns() )
-      return 2;
-   return 3;
+   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,
           typename Device,
           typename Index,
           bool RowMajorOrder >
 Index
 TridiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
-getMaxRowLength() const
+getRowLength( const IndexType row ) const
 {
-   return 3;
+   return this->indexer.getRowSize( row );
 }
 
 template< typename Real,
@@ -120,9 +125,9 @@ template< typename Real,
           bool RowMajorOrder >
 Index
 TridiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
-getNumberOfMatrixElements() const
+getMaxRowLength() const
 {
-   return 3 * min( this->getRows(), this->getColumns() );
+   return 3;
 }
 
 template< typename Real,
@@ -140,17 +145,6 @@ getNumberOfNonzeroMatrixElements() const
    return Algorithms::Reduction< DeviceType >::reduce( this->values.getSize(), std::plus<>{}, fetch, 0 );
 }
 
-template< typename Real,
-          typename Device,
-          typename Index,
-          bool RowMajorOrder >
-Index
-TridiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
-getMaxRowlength() const
-{
-   return 3;
-}
-
 template< typename Real,
           typename Device,
           typename Index,
@@ -228,7 +222,11 @@ setElement( const IndexType row, const IndexType column, const RealType& value )
    TNL_ASSERT_GE( column, 0, "" );
    TNL_ASSERT_LT( column, this->getColumns(), "" );
    if( abs( row - column ) > 1 )
-      throw std::logic_error( "Wrong matrix element coordinates in tridiagonal matrix." );
+   {
+      std::stringstream msg;
+      msg << "Wrong matrix element coordinates ( "  << row << ", " << column << " ) in tridiagonal matrix.";
+      throw std::logic_error( msg.str() );
+   }
    this->values.setElement( this->getElementIndex( row, column ), value );
    return true;
 }
@@ -249,7 +247,11 @@ addElement( const IndexType row,
    TNL_ASSERT_GE( column, 0, "" );
    TNL_ASSERT_LT( column, this->getColumns(), "" );
    if( abs( row - column ) > 1 )
-      throw std::logic_error( "Wrong matrix element coordinates in tridiagonal matrix." );
+   {
+      std::stringstream msg;
+      msg << "Wrong matrix element coordinates ( "  << row << ", " << column << " ) in tridiagonal matrix.";
+      throw std::logic_error( msg.str() );
+   }
    const Index i = this->getElementIndex( row, column );
    this->values.setElement( i, thisElementMultiplicator * this->values.getElement( i ) + value );
    return true;
@@ -280,45 +282,38 @@ template< typename Real,
    template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >
 void
 TridiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
-rowsReduction( IndexType first, IndexType last, Fetch& fetch, Reduce& reduce, Keep& keep, const FetchReal& zero ) const
+rowsReduction( IndexType first, IndexType last, Fetch& fetch, Reduce& reduce, Keep& keep, const FetchReal& zero_ ) const
 {
+   using Real_ = decltype( fetch( IndexType(), IndexType(), RealType() ) );
    const auto values_view = this->values.getConstView();
-   const auto indexer_ = this->indexer;
-   const auto rows = this->getRows();
-   const auto columns = this->getColumns();
-   const auto size = this->size;
+   const auto indexer = this->indexer;
+   const auto zero = zero_;
    auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
-      //bool compute;
+      Real_ sum( zero );
       if( rowIdx == 0 )
       {
-         IndexType i_0 = indexer.getGlobalIndex( 0, 0 );
-         IndexType i_1 = indexer.getGlobalIndex( 0, 1 );
-         keep( 0, reduce( fetch( 0, 0, i_0, values_view[ i_0 ] ),
-                          fetch( 0, 1, i_1, values_view[ i_1 ] ) ) );
+         reduce( sum, fetch( 0, 0, values_view[ indexer.getGlobalIndex( 0, 0 ) ] ) );
+         reduce( sum, fetch( 0, 1, values_view[ indexer.getGlobalIndex( 0, 1 ) ] ) );
+         keep( 0, sum );
          return;
       }
-      if( rowIdx < size || columns > rows )
+      if( rowIdx < indexer.getSize() || indexer.getColumns() > indexer.getRows() )
       {
-         IndexType i_0 = indexer.getGlobalIndex( rowIdx, 0 );
-         IndexType i_1 = indexer.getGlobalIndex( rowIdx, 1 );
-         IndexType i_2 = indexer.getGlobalIndex( rowIdx, 2 );
-
-         keep( rowIdx, reduce( reduce( fetch( rowIdx, rowIdx - 1, i_0, values_view[ i_0 ] ),
-                                       fetch( rowIdx, rowIdx, i_1, values_view[ i_1 ] ) ),
-                               fetch( rowIdx, rowIdx + 1, i_2, values_view[ i_2] ) ) );
+         reduce( sum, fetch( rowIdx, rowIdx - 1, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ] ) );
+         reduce( sum, fetch( rowIdx, rowIdx,     values_view[ indexer.getGlobalIndex( rowIdx, 1 ) ] ) );
+         reduce( sum, fetch( rowIdx, rowIdx + 1, values_view[ indexer.getGlobalIndex( rowIdx, 2 ) ] ) );
+         keep( rowIdx, sum );
          return;
       }
-      if( rows == columns )
+      if( indexer.getRows() == indexer.getColumns() )
       {
-         IndexType i_0 = indexer.getGlobalIndex( rowIdx, 0 );
-         IndexType i_1 = indexer.getGlobalIndex( rowIdx, 1 );
-         keep( rowIdx, reduce( fetch( rowIdx, rowIdx - 1, i_0, values_view[ i_0 ] ),
-                               fetch( rowIdx, rowIdx, i_1, values_view[ i_1 ] ) ) );
+         reduce( sum, fetch( rowIdx, rowIdx - 1, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ] ) );
+         reduce( sum, fetch( rowIdx, rowIdx,     values_view[ indexer.getGlobalIndex( rowIdx, 1 ) ] ) );
+         keep( rowIdx, sum );
       }
       else
       {
-         IndexType i_0 = indexer.getGlobalIndex( rowIdx, 0 );
-         keep( rowIdx, fetch( rowIdx, rowIdx, i_0, values_view[ i_0 ] ) );
+         keep( rowIdx, fetch( rowIdx, rowIdx, values_view[ indexer.getGlobalIndex( rowIdx, 0 ) ] ) );
       }
    };
    Algorithms::ParallelFor< DeviceType >::exec( first, last, f );
@@ -613,41 +608,6 @@ performSORIteration( const Vector1& b,
 }
 
 
-// copy assignment
-template< typename Real,
-          typename Device,
-          typename Index,
-          bool RowMajorOrder >
-TridiagonalMatrixView< Real, Device, Index, RowMajorOrder >&
-TridiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
-operator=( const TridiagonalMatrixView& matrix )
-{
-   this->setLike( matrix );
-   this->values = matrix.values;
-   return *this;
-}
-
-// cross-device copy assignment
-template< typename Real,
-          typename Device,
-          typename Index,
-          bool RowMajorOrder >
-   template< typename Real_, typename Device_, typename Index_, bool RowMajorOrder_ >
-TridiagonalMatrixView< Real, Device, Index, RowMajorOrder >&
-TridiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
-operator=( const TridiagonalMatrixView< Real_, Device_, Index_, RowMajorOrder_ >& matrix )
-{
-   static_assert( std::is_same< Device, Devices::Host >::value || std::is_same< Device, Devices::Cuda >::value,
-                  "unknown device" );
-   static_assert( std::is_same< Device_, Devices::Host >::value || std::is_same< Device_, Devices::Cuda >::value,
-                  "unknown device" );
-
-   this->setLike( matrix );
-
-   throw Exceptions::NotImplementedError("Cross-device assignment for the Tridiagonal format is not implemented yet.");
-}
-
-
 template< typename Real,
           typename Device,
           typename Index,
@@ -690,12 +650,15 @@ template< typename Real,
           bool RowMajorOrder >
 __cuda_callable__
 Index TridiagonalMatrixView< Real, Device, Index, RowMajorOrder >::
-getElementIndex( const IndexType row, const IndexType localIdx ) const
+getElementIndex( const IndexType row, const IndexType column ) const
 {
-   TNL_ASSERT_GE( row, 0, "" );
-   TNL_ASSERT_LT( row, this->getRows(), "" );
+   IndexType localIdx = column - row;
+   if( row > 0 )
+      localIdx++;
+
    TNL_ASSERT_GE( localIdx, 0, "" );
    TNL_ASSERT_LT( localIdx, 3, "" );
+
    return this->indexer.getGlobalIndex( row, localIdx );
 }
 
diff --git a/src/TNL/Matrices/details/TridiagonalMatrixIndexer.h b/src/TNL/Matrices/details/TridiagonalMatrixIndexer.h
index 2f245c38f6593f7542677f44964311d6cc547685..d9fdd0c23cc5bc93ea58b62f31130b9a3674adff 100644
--- a/src/TNL/Matrices/details/TridiagonalMatrixIndexer.h
+++ b/src/TNL/Matrices/details/TridiagonalMatrixIndexer.h
@@ -26,21 +26,21 @@ class TridiagonalMatrixIndexer
 
       __cuda_callable__
       TridiagonalMatrixIndexer()
-      : rows( 0 ), columns( 0 ), size( 0 ){};
+      : rows( 0 ), columns( 0 ), nonEmptyRows( 0 ){};
 
       __cuda_callable__
       TridiagonalMatrixIndexer( const IndexType& rows, const IndexType& columns )
-      : rows( rows ), columns( columns ), size( TNL::min( rows, columns ) ) {};
+      : rows( rows ), columns( columns ), nonEmptyRows( TNL::min( rows, columns ) + ( rows > columns ) ) {};
 
       __cuda_callable__
       TridiagonalMatrixIndexer( const TridiagonalMatrixIndexer& indexer )
-      : rows( indexer.rows ), columns( indexer.columns ), size( indexer.size ) {};
+      : rows( indexer.rows ), columns( indexer.columns ), nonEmptyRows( indexer.nonEmptyRows ) {};
 
       void setDimensions( const IndexType& rows, const IndexType& columns )
       {
          this->rows = rows;
          this->columns = columns;
-         this->size = min( rows, columns );
+         this->nonEmptyRows = min( rows, columns ) + ( rows > columns );
       };
 
       __cuda_callable__
@@ -59,13 +59,15 @@ class TridiagonalMatrixIndexer
       };
 
       __cuda_callable__
-      IndexType getRows() const { return this->rows; };
+      const IndexType& getRows() const { return this->rows; };
 
       __cuda_callable__
-      IndexType getColumns() const { return this->rows; };
+      const IndexType& getColumns() const { return this->columns; };
 
       __cuda_callable__
-      IndexType getStorageSize() const { return 3 * this->size; };
+      const IndexType& getSize() const { return this->nonEmptyRows; };
+      __cuda_callable__
+      IndexType getStorageSize() const { return 3 * this->nonEmptyRows; };
 
       __cuda_callable__
       IndexType getGlobalIndex( const Index rowIdx, const Index localIdx ) const
@@ -78,12 +80,12 @@ class TridiagonalMatrixIndexer
          if( RowMajorOrder )
             return 3 * rowIdx + localIdx;
          else
-            return localIdx * size + rowIdx;
+            return localIdx * nonEmptyRows + rowIdx;
       };
 
       protected:
 
-         IndexType rows, columns, size;
+         IndexType rows, columns, nonEmptyRows;
 };
       } //namespace details
    } // namespace Materices
diff --git a/src/UnitTests/Matrices/DenseMatrixTest.h b/src/UnitTests/Matrices/DenseMatrixTest.h
index 183783ea3ef20cd778d66fd6e1674bea5e34dfb0..0f715801007af25b110584aa17514eb7b06e268f 100644
--- a/src/UnitTests/Matrices/DenseMatrixTest.h
+++ b/src/UnitTests/Matrices/DenseMatrixTest.h
@@ -559,9 +559,9 @@ void test_SetRow()
    /*
     * 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 /
+    *    / 11 11 11 11 11  6  7 \
+    *    | 22 22 22 22 22 13 14 |
+    *    \ 15 16 33 33 33 33 33 /
     */
    const IndexType rows = 3;
    const IndexType cols = 7;
diff --git a/src/UnitTests/Matrices/MultidiagonalMatrixTest.h b/src/UnitTests/Matrices/MultidiagonalMatrixTest.h
index 01ae4a518e022300801c1a9eeccafbbc9e42b7dc..abe6b64c555fcb54068a5e53df968cfe8e11d10a 100644
--- a/src/UnitTests/Matrices/MultidiagonalMatrixTest.h
+++ b/src/UnitTests/Matrices/MultidiagonalMatrixTest.h
@@ -1,7 +1,7 @@
 /***************************************************************************
                           MultidiagonalMatrixTest.h -  description
                              -------------------
-    begin                : Jan 9, 2020
+    begin                : Jan 8, 2020
     copyright            : (C) 2020 by Tomas Oberhuber et al.
     email                : tomas.oberhuber@fjfi.cvut.cz
  ***************************************************************************/
diff --git a/src/UnitTests/Matrices/TridiagonalMatrixTest.h b/src/UnitTests/Matrices/TridiagonalMatrixTest.h
index 962f8c82df9a3cb9a22a888cf63b92c2c682dac1..dcd14302a446b32171574bf805921ba4b63d322a 100644
--- a/src/UnitTests/Matrices/TridiagonalMatrixTest.h
+++ b/src/UnitTests/Matrices/TridiagonalMatrixTest.h
@@ -46,42 +46,42 @@ void test_GetSerializationType()
 template< typename Matrix >
 void test_SetDimensions()
 {
-    using RealType = typename Matrix::RealType;
-    using DeviceType = typename Matrix::DeviceType;
-    using IndexType = typename Matrix::IndexType;
+   using RealType = typename Matrix::RealType;
+   using DeviceType = typename Matrix::DeviceType;
+   using IndexType = typename Matrix::IndexType;
 
-    const IndexType rows = 9;
-    const IndexType cols = 8;
+   const IndexType rows = 9;
+   const IndexType cols = 8;
 
-    Matrix m;
-    m.setDimensions( rows, cols );
+   Matrix m;
+   m.setDimensions( rows, cols );
 
-    EXPECT_EQ( m.getRows(), 9 );
-    EXPECT_EQ( m.getColumns(), 8 );
+   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;
+   using RealType = typename Matrix1::RealType;
+   using DeviceType = typename Matrix1::DeviceType;
+   using IndexType = typename Matrix1::IndexType;
 
-    const IndexType rows = 8;
-    const IndexType cols = 7;
+   const IndexType rows = 8;
+   const IndexType cols = 7;
 
-    Matrix1 m1;
-    m1.reset();
-    m1.setDimensions( rows + 1, cols + 2 );
+   Matrix1 m1;
+   m1.reset();
+   m1.setDimensions( rows + 1, cols + 2 );
 
-    Matrix2 m2;
-    m2.reset();
-    m2.setDimensions( rows, cols );
+   Matrix2 m2;
+   m2.reset();
+   m2.setDimensions( rows, cols );
 
-    m1.setLike( m2 );
+   m1.setLike( m2 );
 
-    EXPECT_EQ( m1.getRows(), m2.getRows() );
-    EXPECT_EQ( m1.getColumns(), m2.getColumns() );
+   EXPECT_EQ( m1.getRows(), m2.getRows() );
+   EXPECT_EQ( m1.getColumns(), m2.getColumns() );
 }
 
 template< typename Matrix >
@@ -94,459 +94,464 @@ void test_GetCompressedRowLengths()
    const IndexType rows = 10;
    const IndexType cols = 11;
 
-    Matrix m( rows, cols );
+   Matrix m( rows, cols );
 
-    // Insert values into the rows.
-    RealType value = 1;
+   // 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 < 2; i++ )  // 0th row -> 2 elements
+      m.setElement( 0, i, value++ );
 
-    for( IndexType i = 0; i < 3; i++ )      // 1st row
-        m.setElement( 1, i, value++ );
+   for( IndexType i = 0; i < 3; i++ )  // 1st row -> 3 elements
+      m.setElement( 1, i, value++ );
 
-    for( IndexType i = 0; i < 1; i++ )      // 2nd row
-        m.setElement( 2, i, value++ );
+   for( IndexType i = 1; i < 3; i++ )  // 2nd row -> 2 elements
+      m.setElement( 2, i, value++ );
 
-    for( IndexType i = 0; i < 2; i++ )      // 3rd row
-        m.setElement( 3, i, value++ );
+   for( IndexType i = 2; i < 5; i++ )  // 3rd row -> 3 elements
+      m.setElement( 3, i, value++ );
 
-    for( IndexType i = 0; i < 3; i++ )      // 4th row
-        m.setElement( 4, i, value++ );
+   for( IndexType i = 3; i < 6; i++ )  // 4th row -> 3 elements
+      m.setElement( 4, i, value++ );
 
-    for( IndexType i = 0; i < 4; i++ )      // 5th row
-        m.setElement( 5, i, value++ );
+   for( IndexType i = 4; i < 6; i++ )  // 5th row -> 2 elements
+      m.setElement( 5, i, value++ );
 
-    for( IndexType i = 0; i < 5; i++ )      // 6th row
-        m.setElement( 6, i, value++ );
+   for( IndexType i = 5; i < 8; i++ )  // 6th row -> 3 elements
+      m.setElement( 6, i, value++ );
 
-    for( IndexType i = 0; i < 6; i++ )      // 7th row
-        m.setElement( 7, i, value++ );
+   for( IndexType i = 6; i < 8; i++ )  // 7th row -> 2 elements
+      m.setElement( 7, i, value++ );
 
-    for( IndexType i = 0; i < 7; i++ )      // 8th row
-        m.setElement( 8, i, value++ );
+   for( IndexType i = 7; i < 10; i++ ) // 8th row -> 3 elements
+      m.setElement( 8, i, value++ );
 
-    for( IndexType i = 0; i < 8; i++ )      // 9th row
-        m.setElement( 9, i, value++ );
+   for( IndexType i = 8; i < 11; i++ ) // 9th row -> 3 elements
+      m.setElement( 9, i, value++ );
 
-   typename Matrix::CompressedRowLengthsVector rowLengths;
+   typename Matrix::CompressedRowLengthsVector rowLengths( rows );
    rowLengths = 0;
    m.getCompressedRowLengths( rowLengths );
-   typename Matrix::CompressedRowLengthsVector correctRowLengths{ 3, 3, 1, 2, 3, 4, 5, 6, 7, 8 };
+   typename Matrix::CompressedRowLengthsVector correctRowLengths{ 2, 3, 2, 3, 3, 2, 3, 2, 3, 3 };
    EXPECT_EQ( rowLengths, correctRowLengths );
 }
 
 template< typename Matrix >
 void test_GetRowLength()
 {
-    using RealType = typename Matrix::RealType;
-    using DeviceType = typename Matrix::DeviceType;
-    using IndexType = typename Matrix::IndexType;
+   using RealType = typename Matrix::RealType;
+   using DeviceType = typename Matrix::DeviceType;
+   using IndexType = typename Matrix::IndexType;
 
-    const IndexType rows = 8;
-    const IndexType cols = 7;
+   const IndexType rows = 8;
+   const IndexType cols = 7;
 
-    Matrix m;
-    m.reset();
-    m.setDimensions( rows, cols );
+   Matrix m( 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 );
+   EXPECT_EQ( m.getRowLength( 0 ), 2 );
+   EXPECT_EQ( m.getRowLength( 1 ), 3 );
+   EXPECT_EQ( m.getRowLength( 2 ), 3 );
+   EXPECT_EQ( m.getRowLength( 3 ), 3 );
+   EXPECT_EQ( m.getRowLength( 4 ), 3 );
+   EXPECT_EQ( m.getRowLength( 5 ), 3 );
+   EXPECT_EQ( m.getRowLength( 6 ), 2 );
+   EXPECT_EQ( m.getRowLength( 7 ), 1 );
 }
 
 template< typename Matrix >
-void test_GetNumberOfMatrixElements()
+void test_GetAllocatedElementsCount()
 {
-    using RealType = typename Matrix::RealType;
-    using DeviceType = typename Matrix::DeviceType;
-    using IndexType = typename Matrix::IndexType;
+   using RealType = typename Matrix::RealType;
+   using DeviceType = typename Matrix::DeviceType;
+   using IndexType = typename Matrix::IndexType;
 
-    const IndexType rows = 7;
-    const IndexType cols = 6;
+   const IndexType rows = 7;
+   const IndexType cols = 6;
 
-    Matrix m;
-    m.reset();
-    m.setDimensions( rows, cols );
+   Matrix m( rows, cols );
 
-    EXPECT_EQ( m.getNumberOfMatrixElements(), 42 );
+   EXPECT_EQ( m.getAllocatedElementsCount(), 21 );
 }
 
 template< typename Matrix >
 void test_GetNumberOfNonzeroMatrixElements()
 {
-    using RealType = typename Matrix::RealType;
-    using DeviceType = typename Matrix::DeviceType;
-    using IndexType = typename Matrix::IndexType;
+   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;
+   /*
+    * Sets up the following 7x6 dense matrix:
+    *
+    *    /  0  1  0  0  0  0 \
+    *    |  2  3  4  0  0  0 |
+    *    |  0  5  6  7  0  0 |
+    *    |  0  0  8  9 10  0 |
+    *    |  0  0  0 11 12 13 |
+    *    |  0  0  0  0 14  0 |
+    *    \  0  0  0  0  0 16 /
+    */
+   const IndexType rows = 7;
+   const IndexType cols = 6;
 
-    Matrix m;
-    m.reset();
-    m.setDimensions( rows, cols );
+   Matrix m( 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 value = 0;
+   for( IndexType i = 0; i < rows; i++ )
+      for( IndexType j = TNL::max( 0, i - 1 ); j < TNL::min( cols, i + 2 ); 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.
+   m.setElement( 5, 5, 0);
 
-    EXPECT_EQ( m.getNumberOfNonzeroMatrixElements(), 40 );
+   EXPECT_EQ( m.getNumberOfNonzeroMatrixElements(), 15 );
 }
 
 template< typename Matrix >
 void test_Reset()
 {
-    using RealType = typename Matrix::RealType;
-    using DeviceType = typename Matrix::DeviceType;
-    using IndexType = typename Matrix::IndexType;
+   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;
+   /*
+    * 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 );
+   Matrix m( rows, cols );
 
-    m.reset();
+   m.reset();
 
-    EXPECT_EQ( m.getRows(), 0 );
-    EXPECT_EQ( m.getColumns(), 0 );
+   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;
+   using RealType = typename Matrix::RealType;
+   using DeviceType = typename Matrix::DeviceType;
+   using IndexType = typename Matrix::IndexType;
 
-    Matrix m;
-    m.reset();
-    m.setDimensions( rows, cols );
+   /*
+    * Sets up the following 7x6 dense matrix:
+    *
+    *    /  0  1  0  0  0  0 \
+    *    |  2  3  4  0  0  0 |
+    *    |  0  5  6  7  0  0 |
+    *    |  0  0  8  9 10  0 |
+    *    |  0  0  0 11 12 13 |
+    *    |  0  0  0  0 14  0 |
+    *    \  0  0  0  0  0 16 /
+    */
+   const IndexType rows = 7;
+   const IndexType cols = 6;
 
-    RealType value = 1;
-    for( IndexType i = 0; i < rows; i++ )
-        for( IndexType j = 0; j < cols; j++ )
-            m.setElement( i, j, value++ );
+   Matrix m( rows, cols );
+
+   RealType value = 0;
+   for( IndexType i = 0; i < rows; i++ )
+      for( IndexType j = TNL::max( 0, i - 1 ); j < TNL::min( cols, i + 2 ); 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 );
+   m.setElement( 5, 5, 0);
+
+   EXPECT_EQ( m.getElement( 0, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 0, 1 ),  1 );
+   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 ),  2 );
+   EXPECT_EQ( m.getElement( 1, 1 ),  3 );
+   EXPECT_EQ( m.getElement( 1, 2 ),  4 );
+   EXPECT_EQ( m.getElement( 1, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 1, 4 ),  0 );
+   EXPECT_EQ( m.getElement( 1, 5 ),  0 );
+
+   EXPECT_EQ( m.getElement( 2, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 2, 1 ),  5 );
+   EXPECT_EQ( m.getElement( 2, 2 ),  6 );
+   EXPECT_EQ( m.getElement( 2, 3 ),  7 );
+   EXPECT_EQ( m.getElement( 2, 4 ),  0 );
+   EXPECT_EQ( m.getElement( 2, 5 ),  0 );
+
+   EXPECT_EQ( m.getElement( 3, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 3, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 3, 2 ),  8 );
+   EXPECT_EQ( m.getElement( 3, 3 ),  9 );
+   EXPECT_EQ( m.getElement( 3, 4 ), 10 );
+   EXPECT_EQ( m.getElement( 3, 5 ),  0 );
+
+   EXPECT_EQ( m.getElement( 4, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 4, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 4, 2 ),  0 );
+   EXPECT_EQ( m.getElement( 4, 3 ), 11 );
+   EXPECT_EQ( m.getElement( 4, 4 ), 12 );
+   EXPECT_EQ( m.getElement( 4, 5 ), 13 );
+
+   EXPECT_EQ( m.getElement( 5, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 2 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 4 ), 14 );
+   EXPECT_EQ( m.getElement( 5, 5 ),  0 );
+
+   EXPECT_EQ( m.getElement( 6, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 6, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 6, 2 ),  0 );
+   EXPECT_EQ( m.getElement( 6, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 6, 4 ),  0 );
+   EXPECT_EQ( m.getElement( 6, 5 ), 16 );
+
+   // 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 ),  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 ), 42 );
+   EXPECT_EQ( m.getElement( 1, 1 ), 42 );
+   EXPECT_EQ( m.getElement( 1, 2 ), 42 );
+   EXPECT_EQ( m.getElement( 1, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 1, 4 ),  0 );
+   EXPECT_EQ( m.getElement( 1, 5 ),  0 );
+
+   EXPECT_EQ( m.getElement( 2, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 2, 1 ), 42 );
+   EXPECT_EQ( m.getElement( 2, 2 ), 42 );
+   EXPECT_EQ( m.getElement( 2, 3 ), 42 );
+   EXPECT_EQ( m.getElement( 2, 4 ),  0 );
+   EXPECT_EQ( m.getElement( 2, 5 ),  0 );
+
+   EXPECT_EQ( m.getElement( 3, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 3, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 3, 2 ), 42 );
+   EXPECT_EQ( m.getElement( 3, 3 ), 42 );
+   EXPECT_EQ( m.getElement( 3, 4 ), 42 );
+   EXPECT_EQ( m.getElement( 3, 5 ),  0 );
+
+   EXPECT_EQ( m.getElement( 4, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 4, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 4, 2 ),  0 );
+   EXPECT_EQ( m.getElement( 4, 3 ), 42 );
+   EXPECT_EQ( m.getElement( 4, 4 ), 42 );
+   EXPECT_EQ( m.getElement( 4, 5 ), 42 );
+
+   EXPECT_EQ( m.getElement( 5, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 2 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 4 ), 42 );
+   EXPECT_EQ( m.getElement( 5, 5 ), 42 );
+
+   EXPECT_EQ( m.getElement( 6, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 6, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 6, 2 ),  0 );
+   EXPECT_EQ( m.getElement( 6, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 6, 4 ),  0 );
+   EXPECT_EQ( m.getElement( 6, 5 ), 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;
+   using RealType = typename Matrix::RealType;
+   using DeviceType = typename Matrix::DeviceType;
+   using IndexType = typename Matrix::IndexType;
 
-    Matrix m;
-    m.reset();
-    m.setDimensions( rows, cols );
+   /*
+    * Sets up the following 5x5 dense matrix:
+    *
+    *    /  1  2  0  0  0 \
+    *    |  6  7  8  0  0 |
+    *    |  0 12 13 14  0 |
+    *    |  0  0 18 19 20 |
+    *    \  0  0  0 24 25 /
+    */
+   const IndexType rows = 5;
+   const IndexType cols = 5;
 
-    RealType value = 1;
-    for( IndexType i = 0; i < rows; i++ )
-        for( IndexType j = 0; j < cols; j++ )
+   Matrix m( rows, cols );
+
+   RealType value = 1;
+   for( IndexType i = 0; i < rows; i++ )
+      for( IndexType j = 0; j < cols; j++ )
+      {
+         if( abs( i - j ) > 1 )
+         {
+            EXPECT_THROW( m.setElement( i, j, value++ ), std::logic_error );
+         }
+         else
             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 ),  0 );
+   EXPECT_EQ( m.getElement( 0, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 0, 4 ),  0 );
+
+   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 ),  0 );
+   EXPECT_EQ( m.getElement( 1, 4 ),  0 );
+
+   EXPECT_EQ( m.getElement( 2, 0 ),  0 );
+   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 ),  0 );
+
+   EXPECT_EQ( m.getElement( 3, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 3, 1 ),  0 );
+   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( 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( 4, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 4, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 4, 2 ),  0 );
+   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;
+   using RealType = typename Matrix::RealType;
+   using DeviceType = typename Matrix::DeviceType;
+   using IndexType = typename Matrix::IndexType;
 
-    Matrix m;
-    m.reset();
-    m.setDimensions( rows, cols );
+   /*
+    * Sets up the following 6x5 dense matrix:
+    *
+    *    /  1  2  0  0  0 \
+    *    |  6  7  8  0  0 |
+    *    |  0 12 13 14  0 |
+    *    |  0  0 18 19 20 |
+    *    |  0  0  0 24 25 |
+    *    \  0  0  0  0 30 /
+    */
+
+   const IndexType rows = 6;
+   const IndexType cols = 5;
+
+   Matrix m( rows, cols );
 
     RealType value = 1;
     for( IndexType i = 0; i < rows; i++ )
         for( IndexType j = 0; j < cols; j++ )
-            m.setElement( i, j, value++ );
+        {
+           if( abs( i - j ) <= 1 )
+               m.setElement( i, j, value );
+           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++ )
+   // Check the added elements
+   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( 1, 0 ),  6 );
+   EXPECT_EQ( m.getElement( 1, 1 ),  7 );
+   EXPECT_EQ( m.getElement( 1, 2 ),  8 );
+   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 ), 12 );
+   EXPECT_EQ( m.getElement( 2, 2 ), 13 );
+   EXPECT_EQ( m.getElement( 2, 3 ), 14 );
+   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 ), 18 );
+   EXPECT_EQ( m.getElement( 3, 3 ), 19 );
+   EXPECT_EQ( m.getElement( 3, 4 ), 20 );
+
+   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 ), 24 );
+   EXPECT_EQ( m.getElement( 4, 4 ), 25 );
+
+   EXPECT_EQ( m.getElement( 5, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 2 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 4 ), 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:
+    *
+    *     /  1  2  0  0  0 \    /  1  2  0  0  0 \   /  3  6  0  0  0 \
+    *     |  6  7  8  0  0 |    |  3  4  5  0  0 |   | 15 18 21  0  0 |
+    * 2 * |  0 12 13 14  0 |  + |  0  6  7  8  0 | = |  0 30 33 36  0 |
+    *     |  0  0 18 19 20 |    |  0  0  9 10 11 |   |  0  0 45 48 51 |
+    *     |  0  0  0 24 25 |    |  0  0  0 12 13 |   |  0  0  0 60 63 |
+    *     \  0  0  0  0 30 /    \  0  0  0  0 14 /   \  0  0  0  0 74 /
+    */
+
+   RealType newValue = 1;
+   RealType multiplicator = 2;
+   for( IndexType i = 0; i < rows; i++ )
+      for( IndexType j = 0; j < cols; j++ )
+         if( abs( i - j ) <= 1 )
             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 );
+   EXPECT_EQ( m.getElement( 0, 0 ),  3 );
+   EXPECT_EQ( m.getElement( 0, 1 ),  6 );
+   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 ), 15 );
+   EXPECT_EQ( m.getElement( 1, 1 ), 18 );
+   EXPECT_EQ( m.getElement( 1, 2 ), 21 );
+   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 ), 30 );
+   EXPECT_EQ( m.getElement( 2, 2 ), 33 );
+   EXPECT_EQ( m.getElement( 2, 3 ), 36 );
+   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 ), 45  );
+   EXPECT_EQ( m.getElement( 3, 3 ), 48 );
+   EXPECT_EQ( m.getElement( 3, 4 ), 51 );
+
+   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 ), 60 );
+   EXPECT_EQ( m.getElement( 4, 4 ), 63 );
+
+   EXPECT_EQ( m.getElement( 5, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 2 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 4 ), 74 );
 }
 
 template< typename Matrix >
@@ -559,61 +564,54 @@ void test_SetRow()
    /*
     * 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 /
+    *    /  1  2  0  0  0  0  0 \
+    *    |  8  9 10  0  0  0  0 |
+    *    \  0 16 17 18  0  0  0 /
     */
    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++ );
+   Matrix m( rows, cols );
 
    auto matrix_view = m.getView();
    auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
-      RealType values[ 3 ][ 5 ] {
-         { 11, 11, 11, 11, 11 },
-         { 22, 22, 22, 22, 22 },
-         { 33, 33, 33, 33, 33 } };
-      IndexType columnIndexes[ 3 ][ 5 ] {
-         { 0, 1, 2, 3, 4 },
-         { 0, 1, 2, 3, 4 },
-         { 2, 3, 4, 5, 6 } };
+      RealType values[ 3 ][ 3 ] {
+         {  1,  2,  0 },
+         {  8,  9, 10 },
+         { 16, 17, 18 } };
       auto row = matrix_view.getRow( rowIdx );
-      for( IndexType i = 0; i < 5; i++ )
-        row.setElement( i, values[ rowIdx ][ i ] );
+      for( IndexType i = 0; i < 3; i++ )
+      {
+         if( rowIdx == 0 && i > 1 )
+            break;
+         row.setElement( i, values[ rowIdx ][ i ] );
+      }
    };
    TNL::Algorithms::ParallelFor< DeviceType >::exec( 0, 3, f );
 
-   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( 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( 0, 6 ),  0 );
+
+   EXPECT_EQ( m.getElement( 1, 0 ),  8 );
+   EXPECT_EQ( m.getElement( 1, 1 ),  9 );
+   EXPECT_EQ( m.getElement( 1, 2 ), 10 );
+   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 ), 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 );
+   EXPECT_EQ( m.getElement( 2, 2 ), 17 );
+   EXPECT_EQ( m.getElement( 2, 3 ), 18 );
+   EXPECT_EQ( m.getElement( 2, 4 ),  0 );
+   EXPECT_EQ( m.getElement( 2, 5 ),  0 );
+   EXPECT_EQ( m.getElement( 2, 6 ),  0 );
 }
 
 template< typename Matrix >
@@ -625,12 +623,12 @@ void test_AddRow()
    /*
     * 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 /
+    *    /  1  2  0  0  0 \
+    *    |  6  7  8  0  0 |
+    *    |  0 12 13 14  0 |
+    *    |  0  0 18 19 20 |
+    *    |  0  0  0 24 25 |
+    *    \  0  0  0  0 30 /
     */
 
    const IndexType rows = 6;
@@ -641,68 +639,72 @@ void test_AddRow()
    RealType value = 1;
    for( IndexType i = 0; i < rows; i++ )
       for( IndexType j = 0; j < cols; j++ )
-         m.setElement( i, j, value++ );
+      {
+         if( abs( i - j ) <= 1 )
+            m.setElement( i, j, value );
+         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( 0, 2 ),  0 );
+   EXPECT_EQ( m.getElement( 0, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 0, 4 ),  0 );
 
    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( 1, 3 ),  0 );
+   EXPECT_EQ( m.getElement( 1, 4 ),  0 );
 
-   EXPECT_EQ( m.getElement( 2, 0 ), 11 );
+   EXPECT_EQ( m.getElement( 2, 0 ),  0 );
    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( 2, 4 ),  0 );
 
-   EXPECT_EQ( m.getElement( 3, 0 ), 16 );
-   EXPECT_EQ( m.getElement( 3, 1 ), 17 );
+   EXPECT_EQ( m.getElement( 3, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 3, 1 ),  0 );
    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, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 4, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 4, 2 ),  0 );
    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, 0 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 1 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 2 ),  0 );
+   EXPECT_EQ( m.getElement( 5, 3 ),  0 );
    EXPECT_EQ( m.getElement( 5, 4 ), 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 /
+    *  / 0  0  0  0  0  0 \   /  1  2  0  0  0 \   / 11 11  0  0  0 \   / 11  11  0   0   0 \
+    *  | 0  1  0  0  0  0 |   |  6  7  8  0  0 |   | 22 22 22  0  0 |   | 28  29 30   0   0 |
+    *  | 0  0  2  0  0  0 | * |  0 12 13 14  0 | + |  0 33 33 33  0 | = |  0  57 59  61   0 |
+    *  | 0  0  0  3  0  0 |   |  0  0 18 19 20 |   |  0  0 44 44 44 |   |  0   0 98 101 104 |
+    *  | 0  0  0  0  4  0 |   |  0  0  0 24 25 |   |  0  0  0 55 55 |   |  0   0  0 151 155 |
+    *  \ 0  0  0  0  0  5 /   \  0  0  0  0 30 /   \  0  0  0  0 66 /   \  0   0  0   0 216 /
     */
 
    auto matrix_view = m.getView();
    auto f = [=] __cuda_callable__ ( IndexType rowIdx ) mutable {
-      RealType values[ 6 ][ 5 ] {
-         { 11, 11, 11, 11, 0 },
-         { 22, 22, 22, 22, 0 },
-         { 33, 33, 33, 33, 0 },
-         { 44, 44, 44, 44, 0 },
-         { 55, 55, 55, 55, 0 },
-         { 66, 66, 66, 66, 0 } };
+      RealType values[ 6 ][ 3 ] {
+         { 11, 11,  0 },
+         { 22, 22, 22 },
+         { 33, 33, 33 },
+         { 44, 44, 44 },
+         { 55, 55, 55 },
+         { 66, 66, 66 } };
       auto row = matrix_view.getRow( rowIdx );
-      for( IndexType i = 0; i < 5; i++ )
+      for( IndexType i = 0; i < 3; i++ )
       {
          RealType& val = row.getValue( i );
          val = rowIdx * val + values[ rowIdx ][ i ];
@@ -711,90 +713,86 @@ void test_AddRow()
    TNL::Algorithms::ParallelFor< DeviceType >::exec( 0, 6, f );
 
 
-    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 );
+   EXPECT_EQ( m.getElement( 0, 0 ),  11 );
+   EXPECT_EQ( m.getElement( 0, 1 ),  11 );
+   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 ),  28 );
+   EXPECT_EQ( m.getElement( 1, 1 ),  29 );
+   EXPECT_EQ( m.getElement( 1, 2 ),  30 );
+   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 ),  57 );
+   EXPECT_EQ( m.getElement( 2, 2 ),  59 );
+   EXPECT_EQ( m.getElement( 2, 3 ),  61 );
+   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 ),  98 );
+   EXPECT_EQ( m.getElement( 3, 3 ), 101 );
+   EXPECT_EQ( m.getElement( 3, 4 ), 104 );
+
+   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 ), 151 );
+   EXPECT_EQ( m.getElement( 4, 4 ), 155 );
+
+   EXPECT_EQ( m.getElement( 5, 0 ),   0 );
+   EXPECT_EQ( m.getElement( 5, 1 ),   0 );
+   EXPECT_EQ( m.getElement( 5, 2 ),   0 );
+   EXPECT_EQ( m.getElement( 5, 3 ),   0 );
+   EXPECT_EQ( m.getElement( 5, 4 ), 216 );
 }
 
 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;
+   using RealType = typename Matrix::RealType;
+   using DeviceType = typename Matrix::DeviceType;
+   using IndexType = typename Matrix::IndexType;
 
-    Matrix m;
-    m.reset();
-    m.setDimensions( rows, cols );
+   /*
+    * Sets up the following 5x4 dense matrix:
+    *
+    *    /  1  2  0  0 \
+    *    |  5  6  7  0 |
+    *    |  0 10 11 12 |
+    *    |  0  0 15 16 |
+    *    \  0  0  0 20 /
+    */
+   const IndexType rows = 5;
+   const IndexType cols = 4;
 
-    RealType value = 1;
-    for( IndexType i = 0; i < rows; i++ )
-        for( IndexType j = 0; j < cols; j++)
-            m.setElement( i, j, value++ );
+   Matrix m( rows, cols );
 
-    using VectorType = TNL::Containers::Vector< RealType, DeviceType, IndexType >;
+   RealType value = 1;
+   for( IndexType i = 0; i < rows; i++ )
+      for( IndexType j = 0; j < cols; j++)
+         if( abs( i - j ) <= 1 )
+            m.setElement( i, j, value++ );
 
-    VectorType inVector;
-    inVector.setSize( 4 );
-    for( IndexType i = 0; i < inVector.getSize(); i++ )
-        inVector.setElement( i, 2 );
+   using VectorType = TNL::Containers::Vector< RealType, DeviceType, IndexType >;
 
-    VectorType outVector;
-    outVector.setSize( 5 );
-    for( IndexType j = 0; j < outVector.getSize(); j++ )
-        outVector.setElement( j, 0 );
+   VectorType inVector( 4 );
+   inVector = 2;
 
+   VectorType outVector( 5 );
+   outVector = 0;
 
-    m.vectorProduct( inVector, outVector);
+   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 );
+   std::cerr << outVector << std::endl;
+   EXPECT_EQ( outVector.getElement( 0 ),  6 );
+   EXPECT_EQ( outVector.getElement( 1 ), 36 );
+   EXPECT_EQ( outVector.getElement( 2 ), 66 );
+   EXPECT_EQ( outVector.getElement( 3 ), 62 );
+   EXPECT_EQ( outVector.getElement( 4 ), 40 );
 }
 
 template< typename Matrix >
@@ -1388,6 +1386,13 @@ TYPED_TEST( MatrixTest, setLikeTest )
     test_SetLike< MatrixType, MatrixType >();
 }
 
+TYPED_TEST( MatrixTest, getCompressedRowLengthTest )
+{
+    using MatrixType = typename TestFixture::MatrixType;
+
+    test_GetCompressedRowLengths< MatrixType >();
+}
+
 TYPED_TEST( MatrixTest, getRowLengthTest )
 {
     using MatrixType = typename TestFixture::MatrixType;
@@ -1395,11 +1400,11 @@ TYPED_TEST( MatrixTest, getRowLengthTest )
     test_GetRowLength< MatrixType >();
 }
 
-TYPED_TEST( MatrixTest, getNumberOfMatrixElementsTest )
+TYPED_TEST( MatrixTest, getAllocatedElementsCountTest )
 {
     using MatrixType = typename TestFixture::MatrixType;
 
-    test_GetNumberOfMatrixElements< MatrixType >();
+    test_GetAllocatedElementsCount< MatrixType >();
 }
 
 TYPED_TEST( MatrixTest, getNumberOfNonzeroMatrixElementsTest )