diff --git a/src/TNL/Matrices/CSR.h b/src/TNL/Matrices/CSR.h
index ef7ba5d6f925d2e56e0df64b2951fe1752a7f84f..1ce7d330bbf17aec36373169389889f4b61773ee 100644
--- a/src/TNL/Matrices/CSR.h
+++ b/src/TNL/Matrices/CSR.h
@@ -75,6 +75,8 @@ public:
 
    __cuda_callable__
    IndexType getRowLengthFast( const IndexType row ) const;
+   
+   IndexType getNonZeroRowLength( const IndexType row ) const;
 
    template< typename Real2, typename Device2, typename Index2 >
    void setLike( const CSR< Real2, Device2, Index2 >& matrix );
diff --git a/src/TNL/Matrices/CSR_impl.h b/src/TNL/Matrices/CSR_impl.h
index b4dff85470bf86021c69478138eb3be86f74d593..5849bbaa76aa1d1af8f87c3b2e18eb1c05c33530 100644
--- a/src/TNL/Matrices/CSR_impl.h
+++ b/src/TNL/Matrices/CSR_impl.h
@@ -131,6 +131,15 @@ Index CSR< Real, Device, Index >::getRowLengthFast( const IndexType row ) const
    return this->rowPointers[ row + 1 ] - this->rowPointers[ row ];
 }
 
+template< typename Real,
+          typename Device,
+          typename Index >
+Index CSR< Real, Device, Index >::getNonZeroRowLength( const IndexType row ) const
+{
+    ConstMatrixRow matrixRow = getRow( row );
+    return matrixRow.getNonZeroElementsCount();
+}
+
 template< typename Real,
           typename Device,
           typename Index >
diff --git a/src/TNL/Matrices/ChunkedEllpack.h b/src/TNL/Matrices/ChunkedEllpack.h
index 35bbfa89799eff2b248283cda4ef141bcf7eb039..ff889a49fe0d7b71c4e1fb1e6e1165e7d8970b86 100644
--- a/src/TNL/Matrices/ChunkedEllpack.h
+++ b/src/TNL/Matrices/ChunkedEllpack.h
@@ -104,6 +104,8 @@ public:
 
    __cuda_callable__
    IndexType getRowLengthFast( const IndexType row ) const;
+   
+   IndexType getNonZeroRowLength( const IndexType row ) const;
 
    template< typename Real2, typename Device2, typename Index2 >
    void setLike( const ChunkedEllpack< Real2, Device2, Index2 >& matrix );
diff --git a/src/TNL/Matrices/ChunkedEllpack_impl.h b/src/TNL/Matrices/ChunkedEllpack_impl.h
index 00cee63bb26298526a3d8255775d24980f6b7d25..56a51149125de61fc345efff836d96c694f39c4e 100644
--- a/src/TNL/Matrices/ChunkedEllpack_impl.h
+++ b/src/TNL/Matrices/ChunkedEllpack_impl.h
@@ -308,6 +308,15 @@ Index ChunkedEllpack< Real, Device, Index >::getRowLengthFast( const IndexType r
    return rowPointers[ row + 1 ] - rowPointers[ row ];
 }
 
+template< typename Real,
+          typename Device,
+          typename Index >
+Index ChunkedEllpack< Real, Device, Index >::getNonZeroRowLength( const IndexType row ) const
+{
+    ConstMatrixRow matrixRow = getRow( row );
+    return matrixRow.getNonZeroElementsCount();
+}
+
 template< typename Real,
           typename Device,
           typename Index >
@@ -979,10 +988,10 @@ getRow( const IndexType rowIndex ) const
 {
    const IndexType rowOffset = this->rowPointers[ rowIndex ];
    const IndexType rowLength = this->rowPointers[ rowIndex + 1 ] - rowOffset;
-   return MatrixRow( &this->columnIndexes[ rowOffset ],
-                     &this->values[ rowOffset ],
-                     rowLength,
-                     1 );
+   return ConstMatrixRow( &this->columnIndexes[ rowOffset ],
+                          &this->values[ rowOffset ],
+                          rowLength,
+                          1 );
 }
 
 
diff --git a/src/TNL/Matrices/Ellpack.h b/src/TNL/Matrices/Ellpack.h
index 1646db1c5c8b37bb635af0ee3501afd3fce6e431..7d17ff07e8649b7c3db872f08672cf0f36fb92ce 100644
--- a/src/TNL/Matrices/Ellpack.h
+++ b/src/TNL/Matrices/Ellpack.h
@@ -67,6 +67,8 @@ public:
 
    __cuda_callable__
    IndexType getRowLengthFast( const IndexType row ) const;
+   
+   IndexType getNonZeroRowLength( const IndexType row ) const;
 
    template< typename Real2, typename Device2, typename Index2 >
    void setLike( const Ellpack< Real2, Device2, Index2 >& matrix );
diff --git a/src/TNL/Matrices/Ellpack_impl.h b/src/TNL/Matrices/Ellpack_impl.h
index 6186206439474d97d8edced12b4671b257b6f0ed..bf3063cd60a37e0352570e01258e793c082c1a5b 100644
--- a/src/TNL/Matrices/Ellpack_impl.h
+++ b/src/TNL/Matrices/Ellpack_impl.h
@@ -123,6 +123,15 @@ Index Ellpack< Real, Device, Index >::getRowLengthFast( const IndexType row ) co
    return this->rowLengths;
 }
 
+template< typename Real,
+          typename Device,
+          typename Index >
+Index Ellpack< Real, Device, Index >::getNonZeroRowLength( const IndexType row ) const
+{
+    ConstMatrixRow matrixRow = getRow( row );
+    return matrixRow.getNonZeroElementsCount();
+}
+
 template< typename Real,
           typename Device,
           typename Index >
diff --git a/src/TNL/Matrices/SlicedEllpack.h b/src/TNL/Matrices/SlicedEllpack.h
index 6f68f2fa8aea4979b8f4685d2ee25d3039653ea7..0fc9ccb0b2ca42bd2ec6441f9fe76cd661679f8e 100644
--- a/src/TNL/Matrices/SlicedEllpack.h
+++ b/src/TNL/Matrices/SlicedEllpack.h
@@ -95,6 +95,8 @@ public:
 
    __cuda_callable__
    IndexType getRowLengthFast( const IndexType row ) const;
+   
+   IndexType getNonZeroRowLength( const IndexType row ) const;
 
    template< typename Real2, typename Device2, typename Index2 >
    void setLike( const SlicedEllpack< Real2, Device2, Index2, SliceSize >& matrix );
diff --git a/src/TNL/Matrices/SlicedEllpack_impl.h b/src/TNL/Matrices/SlicedEllpack_impl.h
index d186bc0471ee3c4f59043f9b6466a5270045b0d5..59b548ade4a572d6b87ea0a0bff3d2cc9783d27b 100644
--- a/src/TNL/Matrices/SlicedEllpack_impl.h
+++ b/src/TNL/Matrices/SlicedEllpack_impl.h
@@ -121,6 +121,16 @@ Index SlicedEllpack< Real, Device, Index, SliceSize >::getRowLengthFast( const I
    return this->sliceCompressedRowLengths[ slice ];
 }
 
+template< typename Real,
+          typename Device,
+          typename Index ,
+          int SliceSize >
+Index SlicedEllpack< Real, Device, Index, SliceSize >::getNonZeroRowLength( const IndexType row ) const
+{
+    ConstMatrixRow matrixRow = getRow( row );
+    return matrixRow.getNonZeroElementsCount();
+}
+
 template< typename Real,
           typename Device,
           typename Index,
diff --git a/src/TNL/Matrices/SparseRow.h b/src/TNL/Matrices/SparseRow.h
index e7547ee679e89cef5e2c607720d13d538ba3c31a..4f8efbdb50c97dac58b51ac2012336f7913c738a 100644
--- a/src/TNL/Matrices/SparseRow.h
+++ b/src/TNL/Matrices/SparseRow.h
@@ -51,6 +51,9 @@ class SparseRow
 
       __cuda_callable__
       Index getLength() const;
+      
+      __cuda_callable__
+      Index getNonZeroElementsCount() const;
 
       void print( std::ostream& str ) const;
 
diff --git a/src/TNL/Matrices/SparseRow_impl.h b/src/TNL/Matrices/SparseRow_impl.h
index f6921b15bbf8f282f0dee1a2f6c842a230cf0dd3..2f0d87d5eb2c354e0fce554f7836fa65cc551a77 100644
--- a/src/TNL/Matrices/SparseRow_impl.h
+++ b/src/TNL/Matrices/SparseRow_impl.h
@@ -11,6 +11,7 @@
 #pragma once
 
 #include <TNL/Matrices/SparseRow.h>
+#include <TNL/ParallelFor.h>
 
 namespace TNL {
 namespace Matrices {   
@@ -107,6 +108,32 @@ getLength() const
    return length;
 }
 
+template< typename Real, typename Index >
+__cuda_callable__
+Index
+SparseRow< Real, Index >::
+getNonZeroElementsCount() const
+{
+    using NonConstIndex = typename std::remove_const< Index >::type;
+    
+    NonConstIndex elementCount ( 0 );
+    
+//    auto computeNonzeros = [this, &elementCount] /*__cuda_callable__*/ ( NonConstIndex i ) mutable
+//    {
+//        if( getElementValue( i ) != ( Real ) 0 )
+//            elementCount++;
+//    };
+   
+//    ParallelFor< Device >::exec( ( NonConstIndex ) 0, length, computeNonzeros );
+//    The ParallelFor::exec() function needs a < DeviceType >, how to get this into SparseRow?
+   
+    for( NonConstIndex i = 0; i < length; i++ )
+        if( getElementValue( i ) != 0 ) // This returns the same amount of elements in a row as does getRowLength(). WHY?
+            elementCount++;
+    
+    return elementCount;
+}
+
 template< typename Real, typename Index >
 void
 SparseRow< Real, Index >::
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_impl.h b/src/UnitTests/Matrices/SparseMatrixTest_impl.h
index 0c1ccbbdd993caa56766eca389a5735ad6cedb6b..d2dae912e609f2e1870aff0b0e893348d5e706b4 100644
--- a/src/UnitTests/Matrices/SparseMatrixTest_impl.h
+++ b/src/UnitTests/Matrices/SparseMatrixTest_impl.h
@@ -136,66 +136,108 @@ void test_SetCompressedRowLengths()
     
     m.setCompressedRowLengths( rowLengths );
     
-    
-    if( m.getType() == TNL::String( TNL::String( "Matrices::CSR< ") +
-                       TNL::String( TNL::getType< RealType >() ) +
-                       TNL::String( ", " ) +
-                       TNL::String( Matrix::DeviceType::getDeviceType() ) +
-                       TNL::String( ", " ) +
-                       TNL::String( TNL::getType< IndexType >() ) +
-                       TNL::String( " >" ) )
-      )
-    {
-        EXPECT_EQ( m.getRowLength( 0 ), 3 );
-        EXPECT_EQ( m.getRowLength( 1 ), 3 );
-        EXPECT_EQ( m.getRowLength( 2 ), 1 );
-        EXPECT_EQ( m.getRowLength( 3 ), 2 );
-        EXPECT_EQ( m.getRowLength( 4 ), 3 );
-        EXPECT_EQ( m.getRowLength( 5 ), 4 );
-        EXPECT_EQ( m.getRowLength( 6 ), 5 );
-        EXPECT_EQ( m.getRowLength( 7 ), 6 );
-        EXPECT_EQ( m.getRowLength( 8 ), 7 );
-        EXPECT_EQ( m.getRowLength( 9 ), 8 );
-    }
-    else if( m.getType() == TNL::String( TNL::String( "Matrices::Ellpack< ") +
-                            TNL::String( TNL::getType< RealType >() ) +
-                            TNL::String( ", " ) +
-                            TNL::String( Matrix::DeviceType::getDeviceType() ) +
-                            TNL::String( ", " ) +
-                            TNL::String( TNL::getType< IndexType >() ) +
-                            TNL::String( " >" ) ) 
-                            || 
-             m.getType() == TNL::String( TNL::String( "Matrices::SlicedEllpack< ") +
-                            TNL::String( TNL::getType< RealType >() ) +
-                            TNL::String( ", " ) +
-                            TNL::String( Matrix::DeviceType::getDeviceType() ) +
-                            TNL::String( " >" ) )
-           )
-    {
-        EXPECT_EQ( m.getRowLength( 0 ), 8 );
-        EXPECT_EQ( m.getRowLength( 1 ), 8 );
-        EXPECT_EQ( m.getRowLength( 2 ), 8 );
-        EXPECT_EQ( m.getRowLength( 3 ), 8 );
-        EXPECT_EQ( m.getRowLength( 4 ), 8 );
-        EXPECT_EQ( m.getRowLength( 5 ), 8 );
-        EXPECT_EQ( m.getRowLength( 6 ), 8 );
-        EXPECT_EQ( m.getRowLength( 7 ), 8 );
-        EXPECT_EQ( m.getRowLength( 8 ), 8 );
-        EXPECT_EQ( m.getRowLength( 9 ), 8 );
-    }
-    else
-    {
-        EXPECT_EQ( m.getRowLength( 0 ), 3 );
-        EXPECT_EQ( m.getRowLength( 1 ), 3 );
-        EXPECT_EQ( m.getRowLength( 2 ), 1 );
-        EXPECT_EQ( m.getRowLength( 3 ), 2 );
-        EXPECT_EQ( m.getRowLength( 4 ), 3 );
-        EXPECT_EQ( m.getRowLength( 5 ), 4 );
-        EXPECT_EQ( m.getRowLength( 6 ), 5 );
-        EXPECT_EQ( m.getRowLength( 7 ), 6 );
-        EXPECT_EQ( m.getRowLength( 8 ), 7 );
-        EXPECT_EQ( m.getRowLength( 9 ), 8 );
-    }
+    RealType realValue = 1; // Do this for every individual row, to assure that non-zero values are not assigned where they're not supposed to be, aka, outside of compressed Row Length
+    for( IndexType i = 0; i < rows; i++ )
+        for( IndexType j = 0; j < cols; j++ )
+            m.setElement( i, j, realValue++ );
+    
+    
+    EXPECT_EQ( m.getNonZeroRowLength( 0 ), 3 );
+    EXPECT_EQ( m.getNonZeroRowLength( 1 ), 3 );
+    EXPECT_EQ( m.getNonZeroRowLength( 2 ), 1 );
+    EXPECT_EQ( m.getNonZeroRowLength( 3 ), 2 );
+    EXPECT_EQ( m.getNonZeroRowLength( 4 ), 3 );
+    EXPECT_EQ( m.getNonZeroRowLength( 5 ), 4 );
+    EXPECT_EQ( m.getNonZeroRowLength( 6 ), 5 );
+    EXPECT_EQ( m.getNonZeroRowLength( 7 ), 6 );
+    EXPECT_EQ( m.getNonZeroRowLength( 8 ), 7 );
+    EXPECT_EQ( m.getNonZeroRowLength( 9 ), 8 );
+    
+//    if( m.getType() == TNL::String( TNL::String( "Matrices::CSR< ") +
+//                       TNL::String( TNL::getType< RealType >() ) +
+//                       TNL::String( ", " ) +
+//                       TNL::String( Matrix::DeviceType::getDeviceType() ) +
+//                       //TNL::String( ", " ) +
+//                       //TNL::String( TNL::getType< IndexType >() ) +
+//                       TNL::String( " >" ) )
+//      )
+//    {
+//        EXPECT_EQ( m.getRowLength( 0 ), 3 );
+//        EXPECT_EQ( m.getRowLength( 1 ), 3 );
+//        EXPECT_EQ( m.getRowLength( 2 ), 1 );
+//        EXPECT_EQ( m.getRowLength( 3 ), 2 );
+//        EXPECT_EQ( m.getRowLength( 4 ), 3 );
+//        EXPECT_EQ( m.getRowLength( 5 ), 4 );
+//        EXPECT_EQ( m.getRowLength( 6 ), 5 );
+//        EXPECT_EQ( m.getRowLength( 7 ), 6 );
+//        EXPECT_EQ( m.getRowLength( 8 ), 7 );
+//        EXPECT_EQ( m.getRowLength( 9 ), 8 );
+//    }
+//    else if( m.getType() == TNL::String( TNL::String( "Matrices::AdEllpack< ") +
+//                            TNL::String( TNL::getType< RealType >() ) +
+//                            TNL::String( ", " ) +
+//                            TNL::String( Matrix::DeviceType::getDeviceType() ) +
+//                            TNL::String( ", " ) +
+//                            TNL::String( TNL::getType< IndexType >() ) +
+//                            TNL::String( " >" ) ) 
+//                            || 
+//             m.getType() == TNL::String( TNL::String( "Matrices::SlicedEllpack< ") +
+//                            TNL::String( TNL::getType< RealType >() ) +
+//                            TNL::String( ", " ) +
+//                            TNL::String( Matrix::DeviceType::getDeviceType() ) +
+//                            TNL::String( " >" ) )
+//           )
+//    {
+//        EXPECT_EQ( m.getRowLength( 0 ), 8 );
+//        EXPECT_EQ( m.getRowLength( 1 ), 8 );
+//        EXPECT_EQ( m.getRowLength( 2 ), 8 );
+//        EXPECT_EQ( m.getRowLength( 3 ), 8 );
+//        EXPECT_EQ( m.getRowLength( 4 ), 8 );
+//        EXPECT_EQ( m.getRowLength( 5 ), 8 );
+//        EXPECT_EQ( m.getRowLength( 6 ), 8 );
+//        EXPECT_EQ( m.getRowLength( 7 ), 8 );
+//        EXPECT_EQ( m.getRowLength( 8 ), 8 );
+//        EXPECT_EQ( m.getRowLength( 9 ), 8 );
+//    }
+//    else if( m.getType() == TNL::String( TNL::String( "Matrices::Ellpack< ") +
+//                            TNL::String( TNL::getType< RealType >() ) +
+//                            TNL::String( ", " ) +
+//                            TNL::String( Matrix::DeviceType::getDeviceType() ) +
+//                            TNL::String( ", " ) +
+//                            TNL::String( TNL::getType< IndexType >() ) +
+//                            TNL::String( " >" ) ) 
+//                            ||
+//             m.getType() == TNL::String( TNL::String( "Matrices::ChunkedEllpack< ") +
+//                            TNL::String( TNL::getType< RealType >() ) +
+//                            TNL::String( ", " ) +
+//                            TNL::String( Matrix::DeviceType::getDeviceType() ) +
+//                            TNL::String( " >" ) )
+//           )
+//    {
+//        EXPECT_EQ( m.getNonZeroRowLength( 0 ), 3 );
+//        EXPECT_EQ( m.getNonZeroRowLength( 1 ), 3 );
+//        EXPECT_EQ( m.getNonZeroRowLength( 2 ), 1 );
+//        EXPECT_EQ( m.getNonZeroRowLength( 3 ), 2 );
+//        EXPECT_EQ( m.getNonZeroRowLength( 4 ), 3 );
+//        EXPECT_EQ( m.getNonZeroRowLength( 5 ), 4 );
+//        EXPECT_EQ( m.getNonZeroRowLength( 6 ), 5 );
+//        EXPECT_EQ( m.getNonZeroRowLength( 7 ), 6 );
+//        EXPECT_EQ( m.getNonZeroRowLength( 8 ), 7 );
+//        EXPECT_EQ( m.getNonZeroRowLength( 9 ), 8 );
+//    }
+//    else
+//    {
+//        EXPECT_EQ( m.getRowLength( 0 ), 3 );
+//        EXPECT_EQ( m.getRowLength( 1 ), 3 );
+//        EXPECT_EQ( m.getRowLength( 2 ), 1 );
+//        EXPECT_EQ( m.getRowLength( 3 ), 2 );
+//        EXPECT_EQ( m.getRowLength( 4 ), 3 );
+//        EXPECT_EQ( m.getRowLength( 5 ), 4 );
+//        EXPECT_EQ( m.getRowLength( 6 ), 5 );
+//        EXPECT_EQ( m.getRowLength( 7 ), 6 );
+//        EXPECT_EQ( m.getRowLength( 8 ), 7 );
+//        EXPECT_EQ( m.getRowLength( 9 ), 8 );
+//    }
 }
 
 template< typename Matrix1, typename Matrix2 >