diff --git a/src/Benchmarks/SpMV/spmv.h b/src/Benchmarks/SpMV/spmv.h
index 26ef145c9b80e9d678cea910ca36622bcc4211a5..8a1b0614e1dd0bce16f2bd1fdb27a833bec37429 100644
--- a/src/Benchmarks/SpMV/spmv.h
+++ b/src/Benchmarks/SpMV/spmv.h
@@ -27,6 +27,7 @@
 #include <TNL/Matrices/MatrixReader.h>
 
 #include <TNL/Matrices/SparseMatrix.h>
+#include <TNL/Matrices/MatrixType.h>
 #include <TNL/Containers/Segments/CSR.h>
 #include <TNL/Containers/Segments/Ellpack.h>
 #include <TNL/Containers/Segments/SlicedEllpack.h>
@@ -43,19 +44,19 @@ using SlicedEllpackAlias = Matrices::SlicedEllpack< Real, Device, Index >;
 
 // Segments based sparse matrix aliases
 template< typename Real, typename Device, typename Index >
-using SparseMatrix_CSR = Matrices::SparseMatrix< Real, Containers::Segments::CSR, Device, Index >;
+using SparseMatrix_CSR = Matrices::SparseMatrix< Real, Device, Index, Matrices::GeneralMatrix, Containers::Segments::CSR >;
 
 template< typename Device, typename Index, typename IndexAllocator >
 using EllpackSegments = Containers::Segments::Ellpack< Device, Index, IndexAllocator >;
 
 template< typename Real, typename Device, typename Index >
-using SparseMatrix_Ellpack = Matrices::SparseMatrix< Real, EllpackSegments, Device, Index >;
+using SparseMatrix_Ellpack = Matrices::SparseMatrix< Real, Device, Index, Matrices::GeneralMatrix, EllpackSegments >;
 
 template< typename Device, typename Index, typename IndexAllocator >
 using SlicedEllpackSegments = Containers::Segments::SlicedEllpack< Device, Index, IndexAllocator >;
 
 template< typename Real, typename Device, typename Index >
-using SparseMatrix_SlicedEllpack = Matrices::SparseMatrix< Real, SlicedEllpackSegments, Device, Index >;
+using SparseMatrix_SlicedEllpack = Matrices::SparseMatrix< Real, Device, Index, Matrices::GeneralMatrix, SlicedEllpackSegments >;
 
 // Get the name (with extension) of input matrix file
 std::string getMatrixFileName( const String& InputFileName )
diff --git a/src/TNL/Matrices/SparseMatrix.h b/src/TNL/Matrices/SparseMatrix.h
index a5effce936d4294ed70089dd10919278ef92d7b6..0d8527daf34d98ddaf70222ffa2719b06defdfd2 100644
--- a/src/TNL/Matrices/SparseMatrix.h
+++ b/src/TNL/Matrices/SparseMatrix.h
@@ -11,15 +11,18 @@
 #pragma once
 
 #include <TNL/Matrices/Matrix.h>
+#include <TNL/Matrices/MatrixType.h>
 #include <TNL/Allocators/Default.h>
+#include <TNL/Containers/Segments/CSR.h>
 
 namespace TNL {
 namespace Matrices {
 
 template< typename Real,
-          template< typename Device_, typename Index_, typename IndexAllocator_ > class Segments,
           typename Device = Devices::Host,
           typename Index = int,
+          typename MatrixType = GeneralMatrix,
+          template< typename Device_, typename Index_, typename IndexAllocator_ > class Segments = Containers::Segments::CSR,
           typename RealAllocator = typename Allocators::Default< Device >::template Allocator< Real >,
           typename IndexAllocator = typename Allocators::Default< Device >::template Allocator< Index > >
 class SparseMatrix : public Matrix< Real, Device, Index, RealAllocator >
@@ -45,6 +48,7 @@ class SparseMatrix : public Matrix< Real, Device, Index, RealAllocator >
       typedef Containers::VectorView< IndexType, DeviceType, IndexType > CompressedRowLengthsVectorView;
       typedef typename CompressedRowLengthsVectorView::ConstViewType ConstCompressedRowLengthsVectorView;
 
+      static constexpr bool isSymmetric() { return MatrixType::isSymmetric(); };
 
       SparseMatrix( const RealAllocatorType& realAllocator = RealAllocatorType(),
                     const IndexAllocatorType& indexAllocator = IndexAllocatorType() );
@@ -83,8 +87,8 @@ class SparseMatrix : public Matrix< Real, Device, Index, RealAllocator >
       __cuda_callable__
       IndexType getNonZeroRowLengthFast( const IndexType row ) const;
 
-      template< typename Real2, template< typename, typename, typename > class Segments2, typename Device2, typename Index2, typename RealAllocator2, typename IndexAllocator2 >
-      void setLike( const SparseMatrix< Real2, Segments2, Device2, Index2, RealAllocator2, IndexAllocator2 >& matrix );
+      template< typename Real2, typename Device2, typename Index2, typename MatrixType2, template< typename, typename, typename > class Segments2, typename RealAllocator2, typename IndexAllocator2 >
+      void setLike( const SparseMatrix< Real2, Device2, Index2, MatrixType2, Segments2, RealAllocator2, IndexAllocator2 >& matrix );
 
       IndexType getNumberOfNonzeroMatrixElements() const;
 
@@ -197,12 +201,13 @@ class SparseMatrix : public Matrix< Real, Device, Index, RealAllocator >
 
       // cross-device copy assignment
       template< typename Real2,
-                template< typename, typename, typename > class Segments2,
                 typename Device2,
                 typename Index2,
+                typename MatrixType2,
+                template< typename, typename, typename > class Segments2,
                 typename RealAllocator2,
                 typename IndexAllocator2 >
-      SparseMatrix& operator=( const SparseMatrix< Real2, Segments2, Device2, Index2, RealAllocator2, IndexAllocator2 >& matrix );
+      SparseMatrix& operator=( const SparseMatrix< Real2, Device2, Index2, MatrixType2, Segments2, RealAllocator2, IndexAllocator2 >& matrix );
 
       void save( File& file ) const;
 
diff --git a/src/TNL/Matrices/SparseMatrix.hpp b/src/TNL/Matrices/SparseMatrix.hpp
index 68f33b93e3c3134d20a3c42015e6b13418e7f2c7..b8091d307899e33c853ea0467a3ce0111f6ec852 100644
--- a/src/TNL/Matrices/SparseMatrix.hpp
+++ b/src/TNL/Matrices/SparseMatrix.hpp
@@ -17,13 +17,14 @@
 namespace TNL {
 namespace Matrices {
 
-   template< typename Real,
-          template< typename, typename, typename > class Segments,
+template< typename Real,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 SparseMatrix( const RealAllocatorType& realAllocator,
               const IndexAllocatorType& indexAllocator )
    : Matrix< Real, Device, Index, RealAllocator >( realAllocator ), columnIndexes( indexAllocator )
@@ -31,36 +32,39 @@ SparseMatrix( const RealAllocatorType& realAllocator,
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 SparseMatrix( const SparseMatrix& m )
    : Matrix< Real, Device, Index, RealAllocator >( m ), columnIndexes( m.columnIndexes )
 {
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 SparseMatrix( const SparseMatrix&& m )
    : Matrix< Real, Device, Index, RealAllocator >( std::move( m ) ), columnIndexes( std::move( m.columnIndexes ) )
 {
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 SparseMatrix( const IndexType rows,
               const IndexType columns,
               const RealAllocatorType& realAllocator,
@@ -70,13 +74,14 @@ SparseMatrix( const IndexType rows,
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 String
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 getSerializationType()
 {
    return String( "Matrices::SparseMatrix< " ) +
@@ -86,27 +91,29 @@ getSerializationType()
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 String
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 getSerializationTypeVirtual() const
 {
    return this->getSerializationType();
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
    template< typename RowsCapacitiesVector >
 void
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 setCompressedRowLengths( const RowsCapacitiesVector& rowsCapacities )
 {
    TNL_ASSERT_EQ( rowsCapacities.getSize(), this->getRows(), "Number of matrix rows does not fit with rowLengths vector size." );
@@ -126,14 +133,15 @@ setCompressedRowLengths( const RowsCapacitiesVector& rowsCapacities )
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
    template< typename Vector >
 void
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 getCompressedRowLengths( Vector& rowLengths ) const
 {
    rowLengths.setSize( this->getRows() );
@@ -152,81 +160,87 @@ getCompressedRowLengths( Vector& rowLengths ) const
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 Index
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 getRowLength( const IndexType row ) const
 {
 
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 __cuda_callable__
 Index
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 getRowLengthFast( const IndexType row ) const
 {
 
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 Index
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 getNonZeroRowLength( const IndexType row ) const
 {
 
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 __cuda_callable__
 Index
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 getNonZeroRowLengthFast( const IndexType row ) const
 {
 
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
-   template< typename Real2, template< typename, typename, typename > class Segments2,  typename Device2, typename Index2, typename RealAllocator2, typename IndexAllocator2 >
+   template< typename Real2, typename Device2, typename Index2, typename MatrixType2, template< typename, typename, typename > class Segments2, typename RealAllocator2, typename IndexAllocator2 >
 void
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
-setLike( const SparseMatrix< Real2, Segments2, Device2, Index2, RealAllocator2, IndexAllocator2 >& matrix )
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
+setLike( const SparseMatrix< Real2, Device2, Index2, MatrixType2, Segments2, RealAllocator2, IndexAllocator2 >& matrix )
 {
    Matrix< Real, Device, Index, RealAllocator >::setLike( matrix );
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 Index
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 getNumberOfNonzeroMatrixElements() const
 {
    const auto columns_view = this->columnIndexes.getConstView();
@@ -238,13 +252,14 @@ getNumberOfNonzeroMatrixElements() const
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 void
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 reset()
 {
    Matrix< Real, Device, Index >::reset();
@@ -253,14 +268,15 @@ reset()
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 __cuda_callable__
 bool
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 setElementFast( const IndexType row,
                 const IndexType column,
                 const RealType& value )
@@ -269,13 +285,14 @@ setElementFast( const IndexType row,
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 bool
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 setElement( const IndexType row,
             const IndexType column,
             const RealType& value )
@@ -284,14 +301,15 @@ setElement( const IndexType row,
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 __cuda_callable__
 bool
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 addElementFast( const IndexType row,
                 const IndexType column,
                 const RealType& value,
@@ -301,13 +319,14 @@ addElementFast( const IndexType row,
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 bool
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 addElement( const IndexType row,
             const IndexType column,
             const RealType& value,
@@ -367,14 +386,15 @@ addElement( const IndexType row,
 
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 __cuda_callable__
 bool
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 setRowFast( const IndexType row,
             const IndexType* columnIndexes,
             const RealType* values,
@@ -383,13 +403,14 @@ setRowFast( const IndexType row,
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 bool
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 setRow( const IndexType row,
         const IndexType* columnIndexes,
         const RealType* values,
@@ -412,14 +433,15 @@ setRow( const IndexType row,
 
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 __cuda_callable__
 bool
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 addRowFast( const IndexType row,
             const IndexType* columns,
             const RealType* values,
@@ -430,13 +452,14 @@ addRowFast( const IndexType row,
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 bool
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 addRow( const IndexType row,
         const IndexType* columns,
         const RealType* values,
@@ -446,16 +469,16 @@ addRow( const IndexType row,
 
 }
 
-
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 __cuda_callable__
 Real
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 getElementFast( const IndexType row,
                 const IndexType column ) const
 {
@@ -463,13 +486,14 @@ getElementFast( const IndexType row,
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 Real
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 getElement( const IndexType row,
             const IndexType column ) const
 {
@@ -486,14 +510,15 @@ getElement( const IndexType row,
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 __cuda_callable__
 void
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 getRowFast( const IndexType row,
             IndexType* columns,
             RealType* values ) const
@@ -502,15 +527,16 @@ getRowFast( const IndexType row,
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
    template< typename Vector >
 __cuda_callable__
 typename Vector::RealType
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 rowVectorProduct( const IndexType row,
                   const Vector& vector ) const
 {
@@ -518,15 +544,16 @@ rowVectorProduct( const IndexType row,
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 template< typename InVector,
        typename OutVector >
 void
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 vectorProduct( const InVector& inVector,
                OutVector& outVector,
                const RealType& matrixMultiplicator,
@@ -553,14 +580,15 @@ vectorProduct( const InVector& inVector,
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
    template< typename Fetch, typename Reduce, typename Keep, typename FetchValue >
 void
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 rowsReduction( IndexType first, IndexType last, Fetch& fetch, Reduce& reduce, Keep& keep, const FetchValue& zero ) const
 {
    const auto columns_view = this->columnIndexes.getConstView();
@@ -576,28 +604,30 @@ rowsReduction( IndexType first, IndexType last, Fetch& fetch, Reduce& reduce, Ke
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
    template< typename Fetch, typename Reduce, typename Keep, typename FetchReal >
 void
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 allRowsReduction( Fetch& fetch, Reduce& reduce, Keep& keep, const FetchReal& zero ) const
 {
    this->rowsReduction( 0, this->getRows(), fetch, reduce, keep, zero );
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
    template< typename Function >
 void
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 forRows( IndexType first, IndexType last, Function& function ) const
 {
    const auto columns_view = this->columnIndexes.getConstView();
@@ -612,14 +642,15 @@ forRows( IndexType first, IndexType last, Function& function ) const
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
    template< typename Function >
 void
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 forAllRows( Function& function ) const
 {
    this->forRows( 0, this->getRows(), function );
@@ -633,7 +664,7 @@ forAllRows( Function& function ) const
           typename IndexAllocator >
 template< typename Real2, template< typename, typename > class Segments2, typename Index2, typename RealAllocator2, typename IndexAllocator2 >
 void
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 addMatrix( const SparseMatrix< Real2, Segments2, Device, Index2, RealAllocator2, IndexAllocator2 >& matrix,
            const RealType& matrixMultiplicator,
            const RealType& thisMatrixMultiplicator )
@@ -649,7 +680,7 @@ template< typename Real,
           typename IndexAllocator >
 template< typename Real2, typename Index2 >
 void
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 getTransposition( const SparseMatrix< Real2, Device, Index2 >& matrix,
                   const RealType& matrixMultiplicator )
 {
@@ -657,14 +688,15 @@ getTransposition( const SparseMatrix< Real2, Device, Index2 >& matrix,
 }*/
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 template< typename Vector1, typename Vector2 >
 bool
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 performSORIteration( const Vector1& b,
                      const IndexType row,
                      Vector2& x,
@@ -675,13 +707,14 @@ performSORIteration( const Vector1& b,
 
 // copy assignment
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >&
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >&
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 operator=( const SparseMatrix& matrix )
 {
    Matrix< Real, Device, Index >::operator=( matrix );
@@ -692,22 +725,24 @@ operator=( const SparseMatrix& matrix )
 
 // cross-device copy assignment
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
    template< typename Real2,
-             template< typename, typename, typename > class Segments2,
              typename Device2,
              typename Index2,
+             typename MatrixType2,
+             template< typename, typename, typename > class Segments2,
              typename RealAllocator2,
              typename IndexAllocator2 >
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >&
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
-operator=( const SparseMatrix< Real2, Segments2, Device2, Index2, RealAllocator2, IndexAllocator2 >& matrix )
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >&
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
+operator=( const SparseMatrix< Real2, Device2, Index2, MatrixType2, Segments2, RealAllocator2, IndexAllocator2 >& matrix )
 {
-   using RHSMatrixType = SparseMatrix< Real2, Segments2, Device2, Index2, RealAllocator2, IndexAllocator2 >;
+   using RHSMatrixType = SparseMatrix< Real2, Device2, Index2, MatrixType2, Segments2, RealAllocator2, IndexAllocator2 >;
    typename RHSMatrixType::RowsCapacitiesType rowLengths;
    matrix.getCompressedRowLengths( rowLengths );
    this->setDimensions( matrix.getRows(), matrix.getColumns() );
@@ -800,13 +835,14 @@ operator=( const SparseMatrix< Real2, Segments2, Device2, Index2, RealAllocator2
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 void
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 save( File& file ) const
 {
    Matrix< RealType, DeviceType, IndexType >::save( file );
@@ -815,13 +851,14 @@ save( File& file ) const
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 void
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 load( File& file )
 {
    Matrix< RealType, DeviceType, IndexType >::load( file );
@@ -830,39 +867,42 @@ load( File& file )
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 void
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 save( const String& fileName ) const
 {
    Object::save( fileName );
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 void
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 load( const String& fileName )
 {
    Object::load( fileName );
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 void
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 print( std::ostream& str ) const
 {
    for( IndexType row = 0; row < this->getRows(); row++ )
@@ -882,14 +922,15 @@ print( std::ostream& str ) const
 }
 
 template< typename Real,
-          template< typename, typename, typename > class Segments,
           typename Device,
           typename Index,
+          typename MatrixType,
+          template< typename, typename, typename > class Segments,
           typename RealAllocator,
           typename IndexAllocator >
 __cuda_callable__
 Index
-SparseMatrix< Real, Segments, Device, Index, RealAllocator, IndexAllocator >::
+SparseMatrix< Real, Device, Index, MatrixType, Segments, RealAllocator, IndexAllocator >::
 getPaddingIndex() const
 {
    return -1;
diff --git a/src/UnitTests/Matrices/SparseMatrixCopyTest.h b/src/UnitTests/Matrices/SparseMatrixCopyTest.h
index 34ffd600ded107c5c12d510e6354cb3adf36812b..d100bb939826ea243041ea1fa165d3b2b1e4c4f0 100644
--- a/src/UnitTests/Matrices/SparseMatrixCopyTest.h
+++ b/src/UnitTests/Matrices/SparseMatrixCopyTest.h
@@ -13,6 +13,7 @@
 #include <TNL/Matrices/SlicedEllpack.h>
 
 #include <TNL/Matrices/SparseMatrix.h>
+#include <TNL/Matrices/MatrixType.h>
 #include <TNL/Containers/Segments/CSR.h>
 #include <TNL/Containers/Segments/Ellpack.h>
 #include <TNL/Containers/Segments/SlicedEllpack.h>
@@ -30,12 +31,12 @@ using EllpackSegments = TNL::Containers::Segments::Ellpack< Device, Index, Index
 template< typename Device, typename Index, typename IndexAllocator >
 using SlicedEllpackSegments = TNL::Containers::Segments::SlicedEllpack< Device, Index, IndexAllocator >;
 
-using CSR_host = TNL::Matrices::SparseMatrix< int, TNL::Containers::Segments::CSR, TNL::Devices::Host, int >;
-using CSR_cuda = TNL::Matrices::SparseMatrix< int, TNL::Containers::Segments::CSR, TNL::Devices::Cuda, int >;
-using E_host   = TNL::Matrices::SparseMatrix< int, EllpackSegments, TNL::Devices::Host, int >;
-using E_cuda   = TNL::Matrices::SparseMatrix< int, EllpackSegments, TNL::Devices::Cuda, int >;
-using SE_host  = TNL::Matrices::SparseMatrix< int, SlicedEllpackSegments, TNL::Devices::Host, int >;
-using SE_cuda  = TNL::Matrices::SparseMatrix< int, SlicedEllpackSegments, TNL::Devices::Cuda, int >;
+using CSR_host = TNL::Matrices::SparseMatrix< int, TNL::Devices::Host, int, TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >;
+using CSR_cuda = TNL::Matrices::SparseMatrix< int, TNL::Devices::Cuda, int, TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >;
+using E_host   = TNL::Matrices::SparseMatrix< int, TNL::Devices::Host, int, TNL::Matrices::GeneralMatrix, EllpackSegments >;
+using E_cuda   = TNL::Matrices::SparseMatrix< int, TNL::Devices::Cuda, int, TNL::Matrices::GeneralMatrix, EllpackSegments >;
+using SE_host  = TNL::Matrices::SparseMatrix< int, TNL::Devices::Host, int, TNL::Matrices::GeneralMatrix, SlicedEllpackSegments >;
+using SE_cuda  = TNL::Matrices::SparseMatrix< int, TNL::Devices::Cuda, int, TNL::Matrices::GeneralMatrix, SlicedEllpackSegments >;
 
 
 #ifdef HAVE_GTEST 
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_CSR_segments.h b/src/UnitTests/Matrices/SparseMatrixTest_CSR_segments.h
index 0718e3a69a8978d322288d9ee65ec1a471b6ee22..353dcdbb0e1721f54803a20b04c417ca04bcca3f 100644
--- a/src/UnitTests/Matrices/SparseMatrixTest_CSR_segments.h
+++ b/src/UnitTests/Matrices/SparseMatrixTest_CSR_segments.h
@@ -29,31 +29,31 @@ protected:
 // types for which MatrixTest is instantiated
 using CSRMatrixTypes = ::testing::Types
 <
-    TNL::Matrices::SparseMatrix< int,     TNL::Containers::Segments::CSR, TNL::Devices::Host, short >,
-    TNL::Matrices::SparseMatrix< long,    TNL::Containers::Segments::CSR, TNL::Devices::Host, short >,
-    TNL::Matrices::SparseMatrix< float,   TNL::Containers::Segments::CSR, TNL::Devices::Host, short >,
-    TNL::Matrices::SparseMatrix< double,  TNL::Containers::Segments::CSR, TNL::Devices::Host, short >,
-    TNL::Matrices::SparseMatrix< int,     TNL::Containers::Segments::CSR, TNL::Devices::Host, int   >,
-    TNL::Matrices::SparseMatrix< long,    TNL::Containers::Segments::CSR, TNL::Devices::Host, int   >,
-    TNL::Matrices::SparseMatrix< float,   TNL::Containers::Segments::CSR, TNL::Devices::Host, int   >,
-    TNL::Matrices::SparseMatrix< double,  TNL::Containers::Segments::CSR, TNL::Devices::Host, int   >,
-    TNL::Matrices::SparseMatrix< int,     TNL::Containers::Segments::CSR, TNL::Devices::Host, long  >,
-    TNL::Matrices::SparseMatrix< long,    TNL::Containers::Segments::CSR, TNL::Devices::Host, long  >,
-    TNL::Matrices::SparseMatrix< float,   TNL::Containers::Segments::CSR, TNL::Devices::Host, long  >,
-    TNL::Matrices::SparseMatrix< double,  TNL::Containers::Segments::CSR, TNL::Devices::Host, long  >
+    TNL::Matrices::SparseMatrix< int,     TNL::Devices::Host, short, TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >,
+    TNL::Matrices::SparseMatrix< long,    TNL::Devices::Host, short, TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >,
+    TNL::Matrices::SparseMatrix< float,   TNL::Devices::Host, short, TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >,
+    TNL::Matrices::SparseMatrix< double,  TNL::Devices::Host, short, TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >,
+    TNL::Matrices::SparseMatrix< int,     TNL::Devices::Host, int,   TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >,
+    TNL::Matrices::SparseMatrix< long,    TNL::Devices::Host, int,   TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >,
+    TNL::Matrices::SparseMatrix< float,   TNL::Devices::Host, int,   TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >,
+    TNL::Matrices::SparseMatrix< double,  TNL::Devices::Host, int,   TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >,
+    TNL::Matrices::SparseMatrix< int,     TNL::Devices::Host, long,  TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >,
+    TNL::Matrices::SparseMatrix< long,    TNL::Devices::Host, long,  TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >,
+    TNL::Matrices::SparseMatrix< float,   TNL::Devices::Host, long,  TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >,
+    TNL::Matrices::SparseMatrix< double,  TNL::Devices::Host, long,  TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >
 #ifdef HAVE_CUDA
-   ,TNL::Matrices::SparseMatrix< int,     TNL::Containers::Segments::CSR, TNL::Devices::Cuda, short >,
-    TNL::Matrices::SparseMatrix< long,    TNL::Containers::Segments::CSR, TNL::Devices::Cuda, short >,
-    TNL::Matrices::SparseMatrix< float,   TNL::Containers::Segments::CSR, TNL::Devices::Cuda, short >,
-    TNL::Matrices::SparseMatrix< double,  TNL::Containers::Segments::CSR, TNL::Devices::Cuda, short >,
-    TNL::Matrices::SparseMatrix< int,     TNL::Containers::Segments::CSR, TNL::Devices::Cuda, int   >,
-    TNL::Matrices::SparseMatrix< long,    TNL::Containers::Segments::CSR, TNL::Devices::Cuda, int   >,
-    TNL::Matrices::SparseMatrix< float,   TNL::Containers::Segments::CSR, TNL::Devices::Cuda, int   >,
-    TNL::Matrices::SparseMatrix< double,  TNL::Containers::Segments::CSR, TNL::Devices::Cuda, int   >,
-    TNL::Matrices::SparseMatrix< int,     TNL::Containers::Segments::CSR, TNL::Devices::Cuda, long  >,
-    TNL::Matrices::SparseMatrix< long,    TNL::Containers::Segments::CSR, TNL::Devices::Cuda, long  >,
-    TNL::Matrices::SparseMatrix< float,   TNL::Containers::Segments::CSR, TNL::Devices::Cuda, long  >,
-    TNL::Matrices::SparseMatrix< double,  TNL::Containers::Segments::CSR, TNL::Devices::Cuda, long  >
+   ,TNL::Matrices::SparseMatrix< int,     TNL::Devices::Cuda, short, TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >,
+    TNL::Matrices::SparseMatrix< long,    TNL::Devices::Cuda, short, TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >,
+    TNL::Matrices::SparseMatrix< float,   TNL::Devices::Cuda, short, TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >,
+    TNL::Matrices::SparseMatrix< double,  TNL::Devices::Cuda, short, TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >,
+    TNL::Matrices::SparseMatrix< int,     TNL::Devices::Cuda, int,   TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >,
+    TNL::Matrices::SparseMatrix< long,    TNL::Devices::Cuda, int,   TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >,
+    TNL::Matrices::SparseMatrix< float,   TNL::Devices::Cuda, int,   TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >,
+    TNL::Matrices::SparseMatrix< double,  TNL::Devices::Cuda, int,   TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >,
+    TNL::Matrices::SparseMatrix< int,     TNL::Devices::Cuda, long,  TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >,
+    TNL::Matrices::SparseMatrix< long,    TNL::Devices::Cuda, long,  TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >,
+    TNL::Matrices::SparseMatrix< float,   TNL::Devices::Cuda, long,  TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >,
+    TNL::Matrices::SparseMatrix< double,  TNL::Devices::Cuda, long,  TNL::Matrices::GeneralMatrix, TNL::Containers::Segments::CSR >
 #endif
 >;
 
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_Ellpack_segments.h b/src/UnitTests/Matrices/SparseMatrixTest_Ellpack_segments.h
index 16c22d9cad957dd62cfaec4fe35303f41a5c002b..b7dc338345e4300786d6844f31ab70bf4637eeaa 100644
--- a/src/UnitTests/Matrices/SparseMatrixTest_Ellpack_segments.h
+++ b/src/UnitTests/Matrices/SparseMatrixTest_Ellpack_segments.h
@@ -40,31 +40,31 @@ using ColumnMajorEllpack = TNL::Containers::Segments::Ellpack< Device, Index, In
 // types for which MatrixTest is instantiated
 using EllpackMatrixTypes = ::testing::Types
 <
-    TNL::Matrices::SparseMatrix< int,     RowMajorEllpack, TNL::Devices::Host, short >,
-    TNL::Matrices::SparseMatrix< long,    RowMajorEllpack, TNL::Devices::Host, short >,
-    TNL::Matrices::SparseMatrix< float,   RowMajorEllpack, TNL::Devices::Host, short >,
-    TNL::Matrices::SparseMatrix< double,  RowMajorEllpack, TNL::Devices::Host, short >,
-    TNL::Matrices::SparseMatrix< int,     RowMajorEllpack, TNL::Devices::Host, int   >,
-    TNL::Matrices::SparseMatrix< long,    RowMajorEllpack, TNL::Devices::Host, int   >,
-    TNL::Matrices::SparseMatrix< float,   RowMajorEllpack, TNL::Devices::Host, int   >,
-    TNL::Matrices::SparseMatrix< double,  RowMajorEllpack, TNL::Devices::Host, int   >,
-    TNL::Matrices::SparseMatrix< int,     RowMajorEllpack, TNL::Devices::Host, long  >,
-    TNL::Matrices::SparseMatrix< long,    RowMajorEllpack, TNL::Devices::Host, long  >,
-    TNL::Matrices::SparseMatrix< float,   RowMajorEllpack, TNL::Devices::Host, long  >,
-    TNL::Matrices::SparseMatrix< double,  RowMajorEllpack, TNL::Devices::Host, long  >
+    TNL::Matrices::SparseMatrix< int,     TNL::Devices::Host, short, TNL::Matrices::GeneralMatrix, RowMajorEllpack >,
+    TNL::Matrices::SparseMatrix< long,    TNL::Devices::Host, short, TNL::Matrices::GeneralMatrix, RowMajorEllpack >,
+    TNL::Matrices::SparseMatrix< float,   TNL::Devices::Host, short, TNL::Matrices::GeneralMatrix, RowMajorEllpack >,
+    TNL::Matrices::SparseMatrix< double,  TNL::Devices::Host, short, TNL::Matrices::GeneralMatrix, RowMajorEllpack >,
+    TNL::Matrices::SparseMatrix< int,     TNL::Devices::Host, int,   TNL::Matrices::GeneralMatrix, RowMajorEllpack >,
+    TNL::Matrices::SparseMatrix< long,    TNL::Devices::Host, int,   TNL::Matrices::GeneralMatrix, RowMajorEllpack >,
+    TNL::Matrices::SparseMatrix< float,   TNL::Devices::Host, int,   TNL::Matrices::GeneralMatrix, RowMajorEllpack >,
+    TNL::Matrices::SparseMatrix< double,  TNL::Devices::Host, int,   TNL::Matrices::GeneralMatrix, RowMajorEllpack >,
+    TNL::Matrices::SparseMatrix< int,     TNL::Devices::Host, long,  TNL::Matrices::GeneralMatrix, RowMajorEllpack >,
+    TNL::Matrices::SparseMatrix< long,    TNL::Devices::Host, long,  TNL::Matrices::GeneralMatrix, RowMajorEllpack >,
+    TNL::Matrices::SparseMatrix< float,   TNL::Devices::Host, long,  TNL::Matrices::GeneralMatrix, RowMajorEllpack >,
+    TNL::Matrices::SparseMatrix< double,  TNL::Devices::Host, long,  TNL::Matrices::GeneralMatrix, RowMajorEllpack >
 #ifdef HAVE_CUDA
-   ,TNL::Matrices::SparseMatrix< int,     ColumnMajorEllpack, TNL::Devices::Cuda, short >,
-    TNL::Matrices::SparseMatrix< long,    ColumnMajorEllpack, TNL::Devices::Cuda, short >,
-    TNL::Matrices::SparseMatrix< float,   ColumnMajorEllpack, TNL::Devices::Cuda, short >,
-    TNL::Matrices::SparseMatrix< double,  ColumnMajorEllpack, TNL::Devices::Cuda, short >,
-    TNL::Matrices::SparseMatrix< int,     ColumnMajorEllpack, TNL::Devices::Cuda, int   >,
-    TNL::Matrices::SparseMatrix< long,    ColumnMajorEllpack, TNL::Devices::Cuda, int   >,
-    TNL::Matrices::SparseMatrix< float,   ColumnMajorEllpack, TNL::Devices::Cuda, int   >,
-    TNL::Matrices::SparseMatrix< double,  ColumnMajorEllpack, TNL::Devices::Cuda, int   >,
-    TNL::Matrices::SparseMatrix< int,     ColumnMajorEllpack, TNL::Devices::Cuda, long  >,
-    TNL::Matrices::SparseMatrix< long,    ColumnMajorEllpack, TNL::Devices::Cuda, long  >,
-    TNL::Matrices::SparseMatrix< float,   ColumnMajorEllpack, TNL::Devices::Cuda, long  >,
-    TNL::Matrices::SparseMatrix< double,  ColumnMajorEllpack, TNL::Devices::Cuda, long  >
+   ,TNL::Matrices::SparseMatrix< int,     TNL::Devices::Cuda, short, TNL::Matrices::GeneralMatrix, ColumnMajorEllpack >,
+    TNL::Matrices::SparseMatrix< long,    TNL::Devices::Cuda, short, TNL::Matrices::GeneralMatrix, ColumnMajorEllpack >,
+    TNL::Matrices::SparseMatrix< float,   TNL::Devices::Cuda, short, TNL::Matrices::GeneralMatrix, ColumnMajorEllpack >,
+    TNL::Matrices::SparseMatrix< double,  TNL::Devices::Cuda, short, TNL::Matrices::GeneralMatrix, ColumnMajorEllpack >,
+    TNL::Matrices::SparseMatrix< int,     TNL::Devices::Cuda, int,   TNL::Matrices::GeneralMatrix, ColumnMajorEllpack >,
+    TNL::Matrices::SparseMatrix< long,    TNL::Devices::Cuda, int,   TNL::Matrices::GeneralMatrix, ColumnMajorEllpack >,
+    TNL::Matrices::SparseMatrix< float,   TNL::Devices::Cuda, int,   TNL::Matrices::GeneralMatrix, ColumnMajorEllpack >,
+    TNL::Matrices::SparseMatrix< double,  TNL::Devices::Cuda, int,   TNL::Matrices::GeneralMatrix, ColumnMajorEllpack >,
+    TNL::Matrices::SparseMatrix< int,     TNL::Devices::Cuda, long,  TNL::Matrices::GeneralMatrix, ColumnMajorEllpack >,
+    TNL::Matrices::SparseMatrix< long,    TNL::Devices::Cuda, long,  TNL::Matrices::GeneralMatrix, ColumnMajorEllpack >,
+    TNL::Matrices::SparseMatrix< float,   TNL::Devices::Cuda, long,  TNL::Matrices::GeneralMatrix, ColumnMajorEllpack >,
+    TNL::Matrices::SparseMatrix< double,  TNL::Devices::Cuda, long,  TNL::Matrices::GeneralMatrix, ColumnMajorEllpack >
 #endif
 >;
 
diff --git a/src/UnitTests/Matrices/SparseMatrixTest_SlicedEllpack_segments.h b/src/UnitTests/Matrices/SparseMatrixTest_SlicedEllpack_segments.h
index 8597121e463a1d8349c1d8f93a9dd71a88f576fd..b2404fe68cc01c163757b4568bc2706508c24c50 100644
--- a/src/UnitTests/Matrices/SparseMatrixTest_SlicedEllpack_segments.h
+++ b/src/UnitTests/Matrices/SparseMatrixTest_SlicedEllpack_segments.h
@@ -10,6 +10,7 @@
 
 #include <TNL/Containers/Segments/SlicedEllpack.h>
 #include <TNL/Matrices/SparseMatrix.h>
+#include <TNL/Matrices/MatrixType.h>
 
 
 #include "SparseMatrixTest.hpp"
@@ -40,31 +41,31 @@ using ColumnMajorSlicedEllpack = TNL::Containers::Segments::SlicedEllpack< Devic
 // types for which MatrixTest is instantiated
 using SlicedEllpackMatrixTypes = ::testing::Types
 <
-    TNL::Matrices::SparseMatrix< int,     RowMajorSlicedEllpack, TNL::Devices::Host, short >,
-    TNL::Matrices::SparseMatrix< long,    RowMajorSlicedEllpack, TNL::Devices::Host, short >,
-    TNL::Matrices::SparseMatrix< float,   RowMajorSlicedEllpack, TNL::Devices::Host, short >,
-    TNL::Matrices::SparseMatrix< double,  RowMajorSlicedEllpack, TNL::Devices::Host, short >,
-    TNL::Matrices::SparseMatrix< int,     RowMajorSlicedEllpack, TNL::Devices::Host, int   >,
-    TNL::Matrices::SparseMatrix< long,    RowMajorSlicedEllpack, TNL::Devices::Host, int   >,
-    TNL::Matrices::SparseMatrix< float,   RowMajorSlicedEllpack, TNL::Devices::Host, int   >,
-    TNL::Matrices::SparseMatrix< double,  RowMajorSlicedEllpack, TNL::Devices::Host, int   >,
-    TNL::Matrices::SparseMatrix< int,     RowMajorSlicedEllpack, TNL::Devices::Host, long  >,
-    TNL::Matrices::SparseMatrix< long,    RowMajorSlicedEllpack, TNL::Devices::Host, long  >,
-    TNL::Matrices::SparseMatrix< float,   RowMajorSlicedEllpack, TNL::Devices::Host, long  >,
-    TNL::Matrices::SparseMatrix< double,  RowMajorSlicedEllpack, TNL::Devices::Host, long  >
+    TNL::Matrices::SparseMatrix< int,     TNL::Devices::Host, short, TNL::Matrices::GeneralMatrix, RowMajorSlicedEllpack >,
+    TNL::Matrices::SparseMatrix< long,    TNL::Devices::Host, short, TNL::Matrices::GeneralMatrix, RowMajorSlicedEllpack >,
+    TNL::Matrices::SparseMatrix< float,   TNL::Devices::Host, short, TNL::Matrices::GeneralMatrix, RowMajorSlicedEllpack >,
+    TNL::Matrices::SparseMatrix< double,  TNL::Devices::Host, short, TNL::Matrices::GeneralMatrix, RowMajorSlicedEllpack >,
+    TNL::Matrices::SparseMatrix< int,     TNL::Devices::Host, int,   TNL::Matrices::GeneralMatrix, RowMajorSlicedEllpack >,
+    TNL::Matrices::SparseMatrix< long,    TNL::Devices::Host, int,   TNL::Matrices::GeneralMatrix, RowMajorSlicedEllpack >,
+    TNL::Matrices::SparseMatrix< float,   TNL::Devices::Host, int,   TNL::Matrices::GeneralMatrix, RowMajorSlicedEllpack >,
+    TNL::Matrices::SparseMatrix< double,  TNL::Devices::Host, int,   TNL::Matrices::GeneralMatrix, RowMajorSlicedEllpack >,
+    TNL::Matrices::SparseMatrix< int,     TNL::Devices::Host, long,  TNL::Matrices::GeneralMatrix, RowMajorSlicedEllpack >,
+    TNL::Matrices::SparseMatrix< long,    TNL::Devices::Host, long,  TNL::Matrices::GeneralMatrix, RowMajorSlicedEllpack >,
+    TNL::Matrices::SparseMatrix< float,   TNL::Devices::Host, long,  TNL::Matrices::GeneralMatrix, RowMajorSlicedEllpack >,
+    TNL::Matrices::SparseMatrix< double,  TNL::Devices::Host, long,  TNL::Matrices::GeneralMatrix, RowMajorSlicedEllpack >
 #ifdef HAVE_CUDA
-   ,TNL::Matrices::SparseMatrix< int,     ColumnMajorSlicedEllpack, TNL::Devices::Cuda, short >,
-    TNL::Matrices::SparseMatrix< long,    ColumnMajorSlicedEllpack, TNL::Devices::Cuda, short >,
-    TNL::Matrices::SparseMatrix< float,   ColumnMajorSlicedEllpack, TNL::Devices::Cuda, short >,
-    TNL::Matrices::SparseMatrix< double,  ColumnMajorSlicedEllpack, TNL::Devices::Cuda, short >,
-    TNL::Matrices::SparseMatrix< int,     ColumnMajorSlicedEllpack, TNL::Devices::Cuda, int   >,
-    TNL::Matrices::SparseMatrix< long,    ColumnMajorSlicedEllpack, TNL::Devices::Cuda, int   >,
-    TNL::Matrices::SparseMatrix< float,   ColumnMajorSlicedEllpack, TNL::Devices::Cuda, int   >,
-    TNL::Matrices::SparseMatrix< double,  ColumnMajorSlicedEllpack, TNL::Devices::Cuda, int   >,
-    TNL::Matrices::SparseMatrix< int,     ColumnMajorSlicedEllpack, TNL::Devices::Cuda, long  >,
-    TNL::Matrices::SparseMatrix< long,    ColumnMajorSlicedEllpack, TNL::Devices::Cuda, long  >,
-    TNL::Matrices::SparseMatrix< float,   ColumnMajorSlicedEllpack, TNL::Devices::Cuda, long  >,
-    TNL::Matrices::SparseMatrix< double,  ColumnMajorSlicedEllpack, TNL::Devices::Cuda, long  >
+   ,TNL::Matrices::SparseMatrix< int,     TNL::Devices::Cuda, short, TNL::Matrices::GeneralMatrix, ColumnMajorSlicedEllpack >,
+    TNL::Matrices::SparseMatrix< long,    TNL::Devices::Cuda, short, TNL::Matrices::GeneralMatrix, ColumnMajorSlicedEllpack >,
+    TNL::Matrices::SparseMatrix< float,   TNL::Devices::Cuda, short, TNL::Matrices::GeneralMatrix, ColumnMajorSlicedEllpack >,
+    TNL::Matrices::SparseMatrix< double,  TNL::Devices::Cuda, short, TNL::Matrices::GeneralMatrix, ColumnMajorSlicedEllpack >,
+    TNL::Matrices::SparseMatrix< int,     TNL::Devices::Cuda, int,   TNL::Matrices::GeneralMatrix, ColumnMajorSlicedEllpack >,
+    TNL::Matrices::SparseMatrix< long,    TNL::Devices::Cuda, int,   TNL::Matrices::GeneralMatrix, ColumnMajorSlicedEllpack >,
+    TNL::Matrices::SparseMatrix< float,   TNL::Devices::Cuda, int,   TNL::Matrices::GeneralMatrix, ColumnMajorSlicedEllpack >,
+    TNL::Matrices::SparseMatrix< double,  TNL::Devices::Cuda, int,   TNL::Matrices::GeneralMatrix, ColumnMajorSlicedEllpack >,
+    TNL::Matrices::SparseMatrix< int,     TNL::Devices::Cuda, long,  TNL::Matrices::GeneralMatrix, ColumnMajorSlicedEllpack >,
+    TNL::Matrices::SparseMatrix< long,    TNL::Devices::Cuda, long,  TNL::Matrices::GeneralMatrix, ColumnMajorSlicedEllpack >,
+    TNL::Matrices::SparseMatrix< float,   TNL::Devices::Cuda, long,  TNL::Matrices::GeneralMatrix, ColumnMajorSlicedEllpack >,
+    TNL::Matrices::SparseMatrix< double,  TNL::Devices::Cuda, long,  TNL::Matrices::GeneralMatrix, ColumnMajorSlicedEllpack >
 #endif
 >;