From 999d1ba3889ad36752524364791712ba07c31233 Mon Sep 17 00:00:00 2001
From: Tomas Oberhuber <tomas.oberhuber@fjfi.cvut.cz>
Date: Wed, 25 Dec 2013 19:50:48 +0100
Subject: [PATCH] Debuging the matrix formats. Implementing the matrix formats
 tests.

---
 src/core/arrays/tnlArray.h                    |   2 +-
 src/core/tnlCuda.h                            |   3 +
 src/core/tnlHost.h                            |   3 +
 src/implementation/core/tnlCuda.cpp           |   5 +
 src/implementation/core/tnlHost.cpp           |  11 +-
 src/implementation/matrices/CMakeLists.txt    |   1 +
 .../matrices/tnlCSRMatrix_impl.h              |  66 ++---------
 .../matrices/tnlChunkedEllpackMatrix_impl.h   | 111 +++++++-----------
 .../matrices/tnlDenseMatrix_impl.h            |  80 +++++++++++--
 .../matrices/tnlEllpackMatrix_impl.h          |  80 +++----------
 .../matrices/tnlMatrixReader_impl.h           |  49 ++++++--
 src/implementation/matrices/tnlMatrix_impl.h  |  81 +++++++++++++
 .../matrices/tnlMultidiagonalMatrix_impl.h    |  47 +++-----
 .../matrices/tnlSlicedEllpackMatrix_impl.h    |  74 ++----------
 .../matrices/tnlSparseMatrix_impl.h           | 102 ++++++++++++++++
 .../matrices/tnlTridiagonalMatrix_impl.h      |  60 ++++++----
 src/legacy/matrices/tnlCusparseCSRMatrix.h    |   1 -
 src/matrices/CMakeLists.txt                   |   1 +
 src/matrices/tnlCSRMatrix.h                   |  19 +--
 src/matrices/tnlChunkedEllpackMatrix.h        |  19 +--
 src/matrices/tnlDenseMatrix.h                 |  30 ++++-
 src/matrices/tnlEllpackMatrix.h               |  20 +---
 src/matrices/tnlMatrix.h                      |  27 +++--
 src/matrices/tnlMatrixReader.h                |   6 +
 src/matrices/tnlMultidiagonalMatrix.h         |  14 +--
 src/matrices/tnlSlicedEllpackMatrix.h         |  24 ++--
 src/matrices/tnlSparseMatrix.h                |  59 ++++++++++
 src/matrices/tnlTridiagonalMatrix.h           |  21 ++--
 tests/benchmarks/CMakeLists.txt               |  14 +--
 tests/benchmarks/matrix-solvers-benchmark.h   |   4 +-
 tests/benchmarks/sparse-matrix-benchmark.h    | 105 ++++-------------
 .../tnlSpmvBenchmarkAdaptiveRgCSRMatrix.h     |   1 -
 .../benchmarks/tnlSpmvBenchmarkRgCSRMatrix.h  |   1 -
 tests/data/tnl-test-matrix-formats.cfg.desc   |   8 +-
 .../matrix-formats-test.h                     |  65 +++++++++-
 .../tnl-run-matrix-formats-test               |   2 +-
 .../matrices/tnlTridiagonalMatrixTester.h     |  18 +--
 37 files changed, 715 insertions(+), 519 deletions(-)
 create mode 100644 src/implementation/matrices/tnlSparseMatrix_impl.h
 create mode 100644 src/matrices/tnlSparseMatrix.h

diff --git a/src/core/arrays/tnlArray.h b/src/core/arrays/tnlArray.h
index d9e0c6f33e..b511849598 100644
--- a/src/core/arrays/tnlArray.h
+++ b/src/core/arrays/tnlArray.h
@@ -30,7 +30,7 @@ class tnlSharedArray;
 template< typename Element,
           typename Device = tnlHost,
           typename Index = int >
-class tnlArray : public tnlObject
+class tnlArray : public virtual tnlObject
 {
    public:
 
diff --git a/src/core/tnlCuda.h b/src/core/tnlCuda.h
index 48069908c0..89c67df3d3 100644
--- a/src/core/tnlCuda.h
+++ b/src/core/tnlCuda.h
@@ -19,6 +19,7 @@
 #define TNLCUDA_H_
 
 #include <iostream>
+#include <unistd.h>
 #include <core/tnlDevice.h>
 #include <core/tnlString.h>
 #include <core/tnlAssert.h>
@@ -43,6 +44,8 @@ class tnlCuda
 
    static bool checkDevice( const char* file_name, int line );
 
+   static size_t getFreeMemory();
+
    protected:
 
    static int maxGridSize, maxBlockSize;
diff --git a/src/core/tnlHost.h b/src/core/tnlHost.h
index 8adcd3ed6a..16ad2fa932 100644
--- a/src/core/tnlHost.h
+++ b/src/core/tnlHost.h
@@ -18,6 +18,7 @@
 #ifndef TNLHOST_H_
 #define TNLHOST_H_
 
+#include <unistd.h>
 #include <core/tnlDevice.h>
 #include <core/tnlString.h>
 
@@ -28,6 +29,8 @@ class tnlHost
    static tnlString getDeviceType();
 
    static tnlDeviceEnum getDevice();
+
+   static size_t getFreeMemory();
 };
 
 #endif /* TNLHOST_H_ */
diff --git a/src/implementation/core/tnlCuda.cpp b/src/implementation/core/tnlCuda.cpp
index 36856fe5c2..a06f45ed37 100644
--- a/src/implementation/core/tnlCuda.cpp
+++ b/src/implementation/core/tnlCuda.cpp
@@ -421,3 +421,8 @@ bool tnlCuda::checkDevice( const char* file_name, int line )
 #endif
 }
 
+size_t tnlCuda::getFreeMemory()
+{
+
+}
+
diff --git a/src/implementation/core/tnlHost.cpp b/src/implementation/core/tnlHost.cpp
index ea839b9621..837323be29 100644
--- a/src/implementation/core/tnlHost.cpp
+++ b/src/implementation/core/tnlHost.cpp
@@ -20,14 +20,21 @@
 
 #include <core/tnlHost.h>
 
-tnlString tnlHost :: getDeviceType()
+tnlString tnlHost::getDeviceType()
 {
    return tnlString( "tnlHost" );
 };
 
-tnlDeviceEnum tnlHost :: getDevice()
+tnlDeviceEnum tnlHost::getDevice()
 {
    return tnlHostDevice;
 };
 
+size_t tnlHost::getFreeMemory()
+{
+   long pages = sysconf(_SC_PHYS_PAGES);
+   long page_size = sysconf(_SC_PAGE_SIZE);
+   return pages * page_size;
+};
+
 #endif /* TNLHOST_H_ */
diff --git a/src/implementation/matrices/CMakeLists.txt b/src/implementation/matrices/CMakeLists.txt
index e31fb4bb42..b4a222734d 100755
--- a/src/implementation/matrices/CMakeLists.txt
+++ b/src/implementation/matrices/CMakeLists.txt
@@ -2,6 +2,7 @@ SET( headers tnlMatrix_impl.h
              tnlDenseMatrix_impl.h
              tnlTridiagonalMatrix_impl.h
              tnlMultidiagonalMatrix_impl.h 
+             tnlSparseMatrix_impl.h
              tnlEllpackMatrix_impl.h
              tnlSlicedEllpackMatrix_impl.h
              tnlChunkedEllpackMatrix_impl.h
diff --git a/src/implementation/matrices/tnlCSRMatrix_impl.h b/src/implementation/matrices/tnlCSRMatrix_impl.h
index 76d5e963a1..580eb7a8e2 100644
--- a/src/implementation/matrices/tnlCSRMatrix_impl.h
+++ b/src/implementation/matrices/tnlCSRMatrix_impl.h
@@ -26,8 +26,6 @@ template< typename Real,
           typename Device,
           typename Index >
 tnlCSRMatrix< Real, Device, Index >::tnlCSRMatrix()
-: rows( 0 ),
-  columns( 0 )
 {
 };
 
@@ -55,14 +53,10 @@ template< typename Real,
           typename Device,
           typename Index >
 bool tnlCSRMatrix< Real, Device, Index >::setDimensions( const IndexType rows,
-                                                                              const IndexType columns )
+                                                         const IndexType columns )
 {
-   tnlAssert( rows > 0 && columns > 0,
-              cerr << "rows = " << rows
-                   << " columns = " << columns << endl );
-   this->rows = rows;
-   this->columns = columns;
-   if( ! this->rowPointers.setSize( this->rows + 1 ) )
+   if( ! tnlSparseMatrix< Real, Device, Index >::setDimensions( rows, columns ) ||
+       ! this->rowPointers.setSize( this->rows + 1 ) )
       return false;
    this->rowPointers.setValue( 0 );
    return true;
@@ -71,8 +65,7 @@ bool tnlCSRMatrix< Real, Device, Index >::setDimensions( const IndexType rows,
 template< typename Real,
           typename Device,
           typename Index >
-   template< typename Vector >
-bool tnlCSRMatrix< Real, Device, Index >::setRowLengths( const Vector& rowLengths )
+bool tnlCSRMatrix< Real, Device, Index >::setRowLengths( const RowLengthsVector& rowLengths )
 {
    /****
     * Compute the rows pointers. The last one is
@@ -103,50 +96,21 @@ template< typename Real,
              typename Index2 >
 bool tnlCSRMatrix< Real, Device, Index >::setLike( const tnlCSRMatrix< Real2, Device2, Index2 >& matrix )
 {
-   if( ! this->setDimensions( matrix.getRows(), matrix.getColumns() ) ||
-       ! this->values.setLike( matrix.values ) ||
-       ! this->columnIndexes.setLike( matrix.columnIndexes ) ||
+   if( ! tnlSparseMatrix< Real, Device, Index >::setLike( matrix ) ||
        ! this->rowPointers.setLike( matrix.rowPointers ) )
       return false;
    return true;
 }
 
-template< typename Real,
-          typename Device,
-          typename Index >
-Index tnlCSRMatrix< Real, Device, Index >::getNumberOfAllocatedElements() const
-{
-   return this->values.getSize();
-}
-
 template< typename Real,
           typename Device,
           typename Index >
 void tnlCSRMatrix< Real, Device, Index >::reset()
 {
-   this->columns = 0;
-   this->rows = 0;
-   this->values.reset();
-   this->columnIndexes.reset();
+   tnlSparseMatrix< Real, Device, Index >::reset();
    this->rowPointers.reset();
 }
 
-template< typename Real,
-          typename Device,
-          typename Index >
-Index tnlCSRMatrix< Real, Device, Index >::getRows() const
-{
-   return this->rows;
-}
-
-template< typename Real,
-          typename Device,
-          typename Index >
-Index tnlCSRMatrix< Real, Device, Index >::getColumns() const
-{
-   return this->columns;
-}
-
 template< typename Real,
           typename Device,
           typename Index >
@@ -197,7 +161,7 @@ Real tnlCSRMatrix< Real, Device, Index >::getElement( const IndexType row,
    const IndexType rowEnd = this->rowPointers[ row + 1 ];
    while( elementPtr < rowEnd && this->columnIndexes[ elementPtr ] < column )
       elementPtr++;
-   if( this->columnIndexes[ elementPtr ] == column )
+   if( elementPtr < rowEnd && this->columnIndexes[ elementPtr ] == column )
       return this->values[ elementPtr ];
    return 0.0;
 }
@@ -347,11 +311,9 @@ template< typename Real,
           typename Index >
 bool tnlCSRMatrix< Real, Device, Index >::save( tnlFile& file ) const
 {
-   if( ! file.write( &this->rows ) ) return false;
-   if( ! file.write( &this->columns ) ) return false;
-   if( ! this->values.save( file ) ) return false;
-   if( ! this->columnIndexes.save( file ) ) return false;
-   if( ! this->rowPointers.save( file ) ) return false;
+   if( ! tnlSparseMatrix< Real, Device, Index >::save( file ) ||
+       ! this->rowPointers.save( file ) )
+      return false;
    return true;
 }
 
@@ -360,11 +322,9 @@ template< typename Real,
           typename Index >
 bool tnlCSRMatrix< Real, Device, Index >::load( tnlFile& file )
 {
-   if( ! file.read( &this->rows ) ) return false;
-   if( ! file.read( &this->columns ) ) return false;
-   if( ! this->values.load( file ) ) return false;
-   if( ! this->columnIndexes.load( file ) ) return false;
-   if( ! this->rowPointers.load( file ) ) return false;
+   if( ! tnlSparseMatrix< Real, Device, Index >::load( file ) ||
+       ! this->rowPointers.load( file ) )
+      return false;
    return true;
 }
 
diff --git a/src/implementation/matrices/tnlChunkedEllpackMatrix_impl.h b/src/implementation/matrices/tnlChunkedEllpackMatrix_impl.h
index c44679a058..46ab5408e3 100644
--- a/src/implementation/matrices/tnlChunkedEllpackMatrix_impl.h
+++ b/src/implementation/matrices/tnlChunkedEllpackMatrix_impl.h
@@ -27,13 +27,11 @@ template< typename Real,
           typename Device,
           typename Index >
 tnlChunkedEllpackMatrix< Real, Device, Index >::tnlChunkedEllpackMatrix()
-: rows( 0 ),
-  columns( 0 ),
-  chunksInSlice( 256 ),
+: chunksInSlice( 256 ),
   desiredChunkSize( 16 )
 {
-   values.setName( "tnlChunkedEllpackMatrix::values" );
-   columnIndexes.setName( "tnlChunkedEllpackMatrix::columnIndexes" );
+   this->values.setName( "tnlChunkedEllpackMatrix::values" );
+   this->columnIndexes.setName( "tnlChunkedEllpackMatrix::columnIndexes" );
    chunksToRowsMapping.setName( "tnlChunkedEllpackMatrix::chunksToRowsMapping" );
    slicesToRowsMapping.setName( "tnlChunkedEllpackMatrix::slicesToRowsMapping" );
    rowPointers.setName( "tnlChunkedEllpackMatrix::rowPointers" );
@@ -69,8 +67,8 @@ bool tnlChunkedEllpackMatrix< Real, Device, Index >::setDimensions( const IndexT
    tnlAssert( rows > 0 && columns > 0,
               cerr << "rows = " << rows
                    << " columns = " << columns << endl );
-   this->rows = rows;
-   this->columns = columns;
+   if( ! tnlSparseMatrix< Real, Device, Index >::setDimensions( rows, columns ) )
+      return false;
 
    /****
     * Allocate slice info array. Note that there cannot be
@@ -87,8 +85,7 @@ bool tnlChunkedEllpackMatrix< Real, Device, Index >::setDimensions( const IndexT
 template< typename Real,
           typename Device,
           typename Index >
-   template< typename Vector >
-bool tnlChunkedEllpackMatrix< Real, Device, Index >::setRowLengths( const Vector& rowLengths )
+bool tnlChunkedEllpackMatrix< Real, Device, Index >::setRowLengths( const RowLengthsVector& rowLengths )
 {
    /****
     * Iterate over rows and allocate slices so that each slice has
@@ -137,9 +134,14 @@ bool tnlChunkedEllpackMatrix< Real, Device, Index >::setRowLengths( const Vector
             RealType rowRatio( 0.0 );
             if( allocatedElementsInSlice != 0 )
                rowRatio = ( RealType ) rowLengths[ i ] / ( RealType ) allocatedElementsInSlice;
-            freeChunks -= this->chunksToRowsMapping[ i ] = ceil( freeChunks * rowRatio );
+            const IndexType addedChunks = ceil( freeChunks * rowRatio );
+            freeChunks -= addedChunks;
+            this->chunksToRowsMapping[ i ] += addedChunks;
+            tnlAssert( chunksToRowsMapping[ i ] > 0,
+                       cerr << " chunksToRowsMapping[ i ] = " << chunksToRowsMapping[ i ] << endl );
          }
          tnlAssert( freeChunks >= 0, );
+         //cout << "freeChunks = " << freeChunks << endl;
       }
 
       /****
@@ -150,6 +152,9 @@ bool tnlChunkedEllpackMatrix< Real, Device, Index >::setRowLengths( const Vector
          maxChunkInSlice = Max( maxChunkInSlice,
                              ceil( ( RealType ) rowLengths[ i ] /
                                    ( RealType ) this->chunksToRowsMapping[ i ] ) );
+      tnlAssert( maxChunkInSlice > 0,
+                 cerr << " maxChunkInSlice = " << maxChunkInSlice << endl );
+
 
       /****
        * Set-up the slice info.
@@ -165,8 +170,14 @@ bool tnlChunkedEllpackMatrix< Real, Device, Index >::setRowLengths( const Vector
       sliceIndex++;
 
       for( IndexType i = sliceBegin; i < sliceEnd; i++ )
+      {
          this->rowPointers[ i + 1 ] =
             this->rowPointers[ i ] + maxChunkInSlice*chunksToRowsMapping[ i ];
+         tnlAssert( this->rowPointers[ i ] >= 0,
+                    cerr << "this->rowPointers[ i ] = " << this->rowPointers[ i ] );
+         tnlAssert( this->rowPointers[ i + 1 ] >= 0,
+                    cerr << "this->rowPointers[ i + 1 ] = " << this->rowPointers[ i + 1 ] );
+      }
 
       /****
        * Finish the chunks to rows mapping by computing the prefix sum.
@@ -177,18 +188,12 @@ bool tnlChunkedEllpackMatrix< Real, Device, Index >::setRowLengths( const Vector
        * Proceed to the next row
        */
       sliceSize = 0;
+      allocatedElementsInSlice = 0;
       if( row < this->rows ) continue;
       else break;
    }
 
-   /****
-    * Allocate values and column indexes
-    */
-   if( ! this->values.setSize( elementsToAllocation ) ||
-       ! this->columnIndexes.setSize( elementsToAllocation ) )
-      return false;
-   this->columnIndexes.setValue( this->columns );
-   return true;
+   return tnlSparseMatrix< Real, Device, Index >::allocateMatrixElements( elementsToAllocation );
 }
 
 template< typename Real,
@@ -201,52 +206,23 @@ bool tnlChunkedEllpackMatrix< Real, Device, Index >::setLike( const tnlChunkedEl
 {
    this->chunksInSlice = matrix.chunksInSlice;
    this->desiredChunkSize = matrix.desiredChunkSize;
-   if( ! this->setDimensions( matrix.getRows(), matrix.getColumns() ) ||
-       ! this->values.setLike( matrix.values ) ||
-       ! this->columnIndexes.setLike( matrix.columnIndexes ) ||
+   if( ! tnlSparseMatrix< Real, Device, Index >::setLike( matrix ) ||
        ! this->chunksToRowsMapping.setLike( matrix.chunksToRowsMapping ) ||
        ! this->slicesToRowsMapping.setLike( matrix.slicesToRowsMapping ) )
       return false;
    return true;
 }
 
-template< typename Real,
-          typename Device,
-          typename Index >
-Index tnlChunkedEllpackMatrix< Real, Device, Index >::getNumberOfAllocatedElements() const
-{
-   return this->values.getSize();
-}
-
 template< typename Real,
           typename Device,
           typename Index >
 void tnlChunkedEllpackMatrix< Real, Device, Index >::reset()
 {
-   this->columns = 0;
-   this->rows = 0;
-   this->values.reset();
-   this->columnIndexes.reset();
+   tnlSparseMatrix< Real, Device, Index >::reset();
    this->chunksToRowsMapping.reset();
    this->slicesToRowsMapping.reset();
 }
 
-template< typename Real,
-          typename Device,
-          typename Index >
-Index tnlChunkedEllpackMatrix< Real, Device, Index >::getRows() const
-{
-   return this->rows;
-}
-
-template< typename Real,
-          typename Device,
-          typename Index >
-Index tnlChunkedEllpackMatrix< Real, Device, Index >::getColumns() const
-{
-   return this->columns;
-}
-
 template< typename Real,
           typename Device,
           typename Index >
@@ -330,9 +306,11 @@ Real tnlChunkedEllpackMatrix< Real, Device, Index >::getElement( const IndexType
    const IndexType& chunkSize = slices.getElement( sliceIndex ).chunkSize;
    IndexType elementPtr = rowPointers[ row ];
    const IndexType rowEnd = rowPointers[ row + 1 ];
+   tnlAssert( rowEnd <= this->columnIndexes.getSize(),
+            cerr << "rowEnd = " << rowEnd << " this->columnIndexes.getSize() = " << this->columnIndexes.getSize() );
    while( elementPtr < rowEnd && this->columnIndexes[ elementPtr ] < column )
       elementPtr++;
-   if( this->columnIndexes[ elementPtr ] == column )
+   if( elementPtr < rowEnd && this->columnIndexes[ elementPtr ] == column )
       return this->values[ elementPtr ];
    return 0.0;
 }
@@ -358,6 +336,11 @@ bool tnlChunkedEllpackMatrix< Real, Device, Index >::addElement( const IndexType
    IndexType elementPtr = rowPointers[ row ];
    const IndexType rowEnd = rowPointers[ row + 1 ];
 
+   tnlAssert( elementPtr >= 0,
+            cerr << "elementPtr = " << elementPtr );
+   tnlAssert( rowEnd <= this->columnIndexes.getSize(),
+            cerr << "rowEnd = " << rowEnd << " this->columnIndexes.getSize() = " << this->columnIndexes.getSize() );
+
    while( elementPtr < rowEnd && this->columnIndexes[ elementPtr ] < column ) elementPtr++;
    if( elementPtr == rowEnd )
       return false;
@@ -496,14 +479,12 @@ template< typename Real,
           typename Index >
 bool tnlChunkedEllpackMatrix< Real, Device, Index >::save( tnlFile& file ) const
 {
-   if( ! file.write( &this->rows ) ) return false;
-   if( ! file.write( &this->columns ) ) return false;
-   if( ! this->values.save( file ) ) return false;
-   if( ! this->columnIndexes.save( file ) ) return false;
-   if( ! this->chunksToRowsMapping.save( file ) ) return false;
-   if( ! this->slicesToRowsMapping.save( file ) ) return false;
-   if( ! this->rowPointers.save( file ) ) return false;
-   if( ! this->slices.save( file ) ) return false;
+   if( ! tnlSparseMatrix< Real, Device, Index >::save( file ) ||
+       ! this->chunksToRowsMapping.save( file ) ||
+       ! this->slicesToRowsMapping.save( file ) ||
+       ! this->rowPointers.save( file ) ||
+       ! this->slices.save( file ) )
+      return false;
    return true;
 }
 
@@ -512,14 +493,12 @@ template< typename Real,
           typename Index >
 bool tnlChunkedEllpackMatrix< Real, Device, Index >::load( tnlFile& file )
 {
-   if( ! file.read( &this->rows ) ) return false;
-   if( ! file.read( &this->columns ) ) return false;
-   if( ! this->values.load( file ) ) return false;
-   if( ! this->columnIndexes.load( file ) ) return false;
-   if( ! this->chunksToRowsMapping.load( file ) ) return false;
-   if( ! this->slicesToRowsMapping.load( file ) ) return false;
-   if( ! this->rowPointers.load( file ) ) return false;
-   if( ! this->slices.load( file ) ) return false;
+   if( ! tnlSparseMatrix< Real, Device, Index >::load( file ) ||
+       ! this->chunksToRowsMapping.load( file ) ||
+       ! this->slicesToRowsMapping.load( file ) ||
+       ! this->rowPointers.load( file ) ||
+       ! this->slices.load( file ) )
+      return false;
    return true;
 }
 
diff --git a/src/implementation/matrices/tnlDenseMatrix_impl.h b/src/implementation/matrices/tnlDenseMatrix_impl.h
index 5983c6f00b..c610061bab 100644
--- a/src/implementation/matrices/tnlDenseMatrix_impl.h
+++ b/src/implementation/matrices/tnlDenseMatrix_impl.h
@@ -53,15 +53,28 @@ template< typename Real,
 bool tnlDenseMatrix< Real, Device, Index >::setDimensions( const IndexType rows,
                                                            const IndexType columns )
 {
-   return tnlMultiArray< 2, Real, Device, Index >::setDimensions( columns, rows );
+   if( ! tnlMatrix< Real, Device, Index >::setDimensions( rows, columns ) ||
+       ! tnlMultiArray< 2, Real, Device, Index >::setDimensions( rows, columns ) )
+     return false;
    tnlMultiArray< 2, Real, Device, Index >::setValue( 0.0 );
+   return true;
 }
 
 template< typename Real,
           typename Device,
           typename Index >
-   template< typename Vector >
-bool tnlDenseMatrix< Real, Device, Index >::setRowLengths( const Vector& rowLengths )
+   template< typename Real2,
+             typename Device2,
+             typename Index2 >
+bool tnlDenseMatrix< Real, Device, Index >::setLike( const tnlDenseMatrix< Real2, Device2, Index2 >& matrix )
+{
+   return this->setDimensions( matrix.getRows(), matrix.getColumns() );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+bool tnlDenseMatrix< Real, Device, Index >::setRowLengths( const RowLengthsVector& rowLengths )
 {
    return true;
 }
@@ -69,7 +82,7 @@ bool tnlDenseMatrix< Real, Device, Index >::setRowLengths( const Vector& rowLeng
 template< typename Real,
           typename Device,
           typename Index >
-Index tnlDenseMatrix< Real, Device, Index >::getNumberOfAllocatedElements() const
+Index tnlDenseMatrix< Real, Device, Index >::getNumberOfMatrixElements() const
 {
    return this->getRows() * this->getColumns();
 }
@@ -77,17 +90,30 @@ Index tnlDenseMatrix< Real, Device, Index >::getNumberOfAllocatedElements() cons
 template< typename Real,
           typename Device,
           typename Index >
-Index tnlDenseMatrix< Real, Device, Index >::getRows() const
+void tnlDenseMatrix< Real, Device, Index >::reset()
 {
-   return this->getDimensions().x();
+   tnlMatrix< Real, Device, Index >::reset();
+   tnlMultiArray< 2, Real, Device, Index >::reset();
 }
 
 template< typename Real,
           typename Device,
           typename Index >
-Index tnlDenseMatrix< Real, Device, Index >::getColumns() const
+bool tnlDenseMatrix< Real, Device, Index >::setElement( const IndexType row,
+                                                        const IndexType column,
+                                                        const RealType& value )
 {
-   return this->getDimensions().y();
+   tnlMultiArray< 2, Real, Device, Index >::setElement( row, column, value );
+   return true;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+Real tnlDenseMatrix< Real, Device, Index >::getElement( const IndexType row,
+                                                        const IndexType column ) const
+{
+   return tnlMultiArray< 2, Real, Device, Index >::getElement( row, column );
 }
 
 template< typename Real,
@@ -293,6 +319,44 @@ void tnlDenseMatrix< Real, Device, Index >::performSORIteration( const Vector& b
    x[ row ] += omega / this->operator()( row, row )( b[ row ] - sum );
 }
 
+template< typename Real,
+          typename Device,
+          typename Index >
+bool tnlDenseMatrix< Real, Device, Index >::save( const tnlString& fileName ) const
+{
+   return tnlObject::save( fileName );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+bool tnlDenseMatrix< Real, Device, Index >::load( const tnlString& fileName )
+{
+   return tnlObject::load( fileName );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+bool tnlDenseMatrix< Real, Device, Index >::save( tnlFile& file ) const
+{
+   if( ! tnlMatrix< Real, Device, Index >::save( file ) ||
+       ! tnlMultiArray< 2, Real, Device, Index >::save( file ) )
+      return false;
+   return true;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+bool tnlDenseMatrix< Real, Device, Index >::load( tnlFile& file )
+{
+   if( ! tnlMatrix< Real, Device, Index >::load( file ) ||
+       ! tnlMultiArray< 2, Real, Device, Index >::load( file ) )
+      return false;
+   return true;
+}
+
 template< typename Real,
           typename Device,
           typename Index >
diff --git a/src/implementation/matrices/tnlEllpackMatrix_impl.h b/src/implementation/matrices/tnlEllpackMatrix_impl.h
index 61a3fa6f74..611d0d83ec 100644
--- a/src/implementation/matrices/tnlEllpackMatrix_impl.h
+++ b/src/implementation/matrices/tnlEllpackMatrix_impl.h
@@ -26,9 +26,7 @@ template< typename Real,
           typename Device,
           typename Index >
 tnlEllpackMatrix< Real, Device, Index > :: tnlEllpackMatrix()
-: rows( 0 ),
-  columns( 0 ),
-  rowLengths( 0 )
+: rowLengths( 0 )
 {
 };
 
@@ -55,8 +53,8 @@ tnlString tnlEllpackMatrix< Real, Device, Index >::getTypeVirtual() const
 template< typename Real,
           typename Device,
           typename Index >
-bool tnlEllpackMatrix< Real, Device, Index > :: setDimensions( const IndexType rows,
-                                                               const IndexType columns )
+bool tnlEllpackMatrix< Real, Device, Index >::setDimensions( const IndexType rows,
+                                                             const IndexType columns )
 {
    tnlAssert( rows > 0 && columns > 0,
               cerr << "rows = " << rows
@@ -71,8 +69,7 @@ bool tnlEllpackMatrix< Real, Device, Index > :: setDimensions( const IndexType r
 template< typename Real,
           typename Device,
           typename Index >
-   template< typename Vector >          
-bool tnlEllpackMatrix< Real, Device, Index > :: setRowLengths( const Vector& rowLengths )
+bool tnlEllpackMatrix< Real, Device, Index >::setRowLengths( const RowLengthsVector& rowLengths )
 {
    this->rowLengths = 0;
    for( IndexType i = 0; i < rowLengths.getSize(); i++ )
@@ -87,7 +84,7 @@ bool tnlEllpackMatrix< Real, Device, Index > :: setRowLengths( const Vector& row
 template< typename Real,
           typename Device,
           typename Index >
-bool tnlEllpackMatrix< Real, Device, Index > :: setConstantRowLengths( const IndexType& rowLengths )
+bool tnlEllpackMatrix< Real, Device, Index >::setConstantRowLengths( const IndexType& rowLengths )
 {
    tnlAssert( rowLengths > 0,
               cerr << " rowLengths = " << rowLengths );
@@ -103,49 +100,21 @@ template< typename Real,
    template< typename Real2,
              typename Device2,
              typename Index2 >
-bool tnlEllpackMatrix< Real, Device, Index > :: setLike( const tnlEllpackMatrix< Real2, Device2, Index2 >& matrix )
+bool tnlEllpackMatrix< Real, Device, Index >::setLike( const tnlEllpackMatrix< Real2, Device2, Index2 >& matrix )
 {
-   this->rowLengths = 0;
-   this->setDimensions( matrix.getRows(), matrix.getColumns() );
-   if( ! this->setConstantRowLengths( matrix.rowLengths ) )
+   if( ! tnlSparseMatrix< Real, Device, Index >::setLike( matrix ) )
       return false;
+   this->rowLengths = matrix.rowLengths;
    return true;
 }
 
-template< typename Real,
-          typename Device,
-          typename Index >
-Index tnlEllpackMatrix< Real, Device, Index > :: getNumberOfAllocatedElements() const
-{
-   return this->values.getSize();
-}
-
 template< typename Real,
           typename Device,
           typename Index >
 void tnlEllpackMatrix< Real, Device, Index > :: reset()
 {
-   this->rows = 0;
-   this->columns = 0;
+   tnlSparseMatrix< Real, Device, Index >::reset();
    this->rowLengths = 0;
-   this->values.reset();
-   this->columnIndexes.reset();
-}
-
-template< typename Real,
-          typename Device,
-          typename Index >
-Index tnlEllpackMatrix< Real, Device, Index >::getRows() const
-{
-   return this->rows;
-}
-
-template< typename Real,
-          typename Device,
-          typename Index >
-Index tnlEllpackMatrix< Real, Device, Index >::getColumns() const
-{
-   return this->columns;
 }
 
 template< typename Real,
@@ -194,12 +163,12 @@ template< typename Real,
 Real tnlEllpackMatrix< Real, Device, Index >::getElement( const IndexType row,
                                                           const IndexType column ) const
 {
-   IndexType i( row * this->rowLengths );
-   const IndexType rowEnd( i + this->rowLengths );
-   while( i < rowEnd && this->columnIndexes[ i ] < column ) i++;
-   if( i == rowEnd || this->columnIndexes[ i ] != column )
-      return 0.0;
-   return this->values.getElement( i );
+   IndexType elementPtr( row * this->rowLengths );
+   const IndexType rowEnd( elementPtr + this->rowLengths );
+   while( elementPtr < rowEnd && this->columnIndexes[ elementPtr ] < column ) elementPtr++;
+   if( elementPtr < rowEnd && this->columnIndexes[ elementPtr ] == column )
+      return this->values.getElement( elementPtr );
+   return 0.0;
 }
 
 template< typename Real,
@@ -339,11 +308,8 @@ template< typename Real,
           typename Index >
 bool tnlEllpackMatrix< Real, Device, Index >::save( tnlFile& file ) const
 {
-   if( ! file.write( &this->rows ) ) return false;
-   if( ! file.write( &this->columns ) ) return false;
+   if( ! tnlSparseMatrix< Real, Device, Index >::save( file) ) return false;
    if( ! file.write( &this->rowLengths ) ) return false;
-   if( ! this->values.save( file ) ) return false;
-   if( ! this->columnIndexes.save( file ) ) return false;
    return true;
 }
 
@@ -352,11 +318,8 @@ template< typename Real,
           typename Index >
 bool tnlEllpackMatrix< Real, Device, Index >::load( tnlFile& file )
 {
-   if( ! file.read( &this->rows ) ) return false;
-   if( ! file.read( &this->columns ) ) return false;
+   if( ! tnlSparseMatrix< Real, Device, Index >::load( file) ) return false;
    if( ! file.read( &this->rowLengths ) ) return false;
-   if( ! this->values.load( file ) ) return false;
-   if( ! this->columnIndexes.load( file ) ) return false;
    return true;
 }
 
@@ -401,15 +364,8 @@ template< typename Real,
           typename Index >
 bool tnlEllpackMatrix< Real, Device, Index >::allocateElements()
 {
-   if( ! this->values.setSize( this->rows * this->rowLengths ) ||
-       ! this->columnIndexes.setSize( this->rows * this->rowLengths ) )
+   if( ! tnlSparseMatrix< Real, Device, Index >::allocateMatrixElements( this->rows * this->rowLengths ) )
       return false;
-
-   /****
-    * Setting a column index to this->columns means that the
-    * index is undefined.
-    */
-   this->columnIndexes.setValue( this->columns );
    return true;
 }
 
diff --git a/src/implementation/matrices/tnlMatrixReader_impl.h b/src/implementation/matrices/tnlMatrixReader_impl.h
index ecfc088673..a03a386764 100644
--- a/src/implementation/matrices/tnlMatrixReader_impl.h
+++ b/src/implementation/matrices/tnlMatrixReader_impl.h
@@ -78,7 +78,7 @@ bool tnlMatrixReader< Matrix >::verifyMtxFile( std::istream& file,
          dimensionsLine = true;
          continue;
       }
-      IndexType row, column;
+      IndexType row( 1 ), column( 1 );
       RealType value;
       if( ! parseMtxLineWithElement( line, row, column, value ) )
          return false;
@@ -95,17 +95,52 @@ bool tnlMatrixReader< Matrix >::verifyMtxFile( std::istream& file,
       if( symmetricMatrix && row != column )
          processedElements++;
       if( verbose )
-         cout << " Verifying the matrix elements ... " << processedElements << " / " << matrix.getNumberOfAllocatedElements() << "                       \r" << flush;
+         cout << " Verifying the matrix elements ... " << processedElements << " / " << matrix.getNumberOfMatrixElements() << "                       \r" << flush;
    }
    file.clear();
    long int fileSize = file.tellg();
    if( verbose )
-      cout << " Verifying the matrix elements ... " << processedElements << " / " << matrix.getNumberOfAllocatedElements()
+      cout << " Verifying the matrix elements ... " << processedElements << " / " << matrix.getNumberOfMatrixElements()
            << " -> " << timer.GetTime()
            << " sec. i.e. " << fileSize / ( timer.GetTime() * ( 1 << 20 ))  << "MB/s." << endl;
    return true;
 }
 
+template< typename Matrix >
+bool tnlMatrixReader< Matrix >::findLineByElement( std::istream& file,
+                                                   const IndexType& row,
+                                                   const IndexType& column,
+                                                   tnlString& line,
+                                                   IndexType& lineNumber )
+{
+   file.clear();
+   file.seekg( 0, ios::beg );
+   bool symmetricMatrix( false );
+   bool dimensionsLine( false );
+   lineNumber = 0;
+   tnlTimerRT timer;
+   IndexType currentRow, currentColumn;
+   RealType value;
+   while( line.getLine( file ) )
+   {
+      lineNumber++;
+      if( line[ 0 ] == '%' ) continue;
+      if( ! dimensionsLine )
+      {
+         dimensionsLine = true;
+         continue;
+      }
+      IndexType currentRow( 1 ), currentColumn( 1 );
+      RealType value;
+      if( ! parseMtxLineWithElement( line, currentRow, currentColumn, value ) )
+         return false;
+      if( ( currentRow == row + 1 && currentColumn == column + 1 ) ||
+          ( symmetricMatrix && currentRow == column + 1 && currentColumn == row + 1 ) )
+         return true;
+   }
+   return false;
+}
+
 template< typename Matrix >
 bool tnlMatrixReader< Matrix >::checkMtxHeader( const tnlString& header,
                                                 bool& symmetric )
@@ -221,7 +256,7 @@ bool tnlMatrixReader< Matrix >::computeRowLengthsFromMtxFile( std::istream& file
          dimensionsLine = true;
          continue;
       }
-      IndexType row, column;
+      IndexType row( 1 ), column( 1 );
       RealType value;
       if( ! parseMtxLineWithElement( line, row, column, value ) )
          return false;
@@ -261,7 +296,7 @@ bool tnlMatrixReader< Matrix >::readMatrixElementsFromMtxFile( std::istream& fil
          dimensionsLine = true;
          continue;
       }
-      IndexType row, column;
+      IndexType row( 1 ), column( 1 );
       RealType value;
       if( ! parseMtxLineWithElement( line, row, column, value ) )
          return false;
@@ -273,12 +308,12 @@ bool tnlMatrixReader< Matrix >::readMatrixElementsFromMtxFile( std::istream& fil
          processedElements++;
       }
       if( verbose )
-         cout << " Reading the matrix elements ... " << processedElements << " / " << matrix.getNumberOfAllocatedElements() << "                       \r" << flush;
+         cout << " Reading the matrix elements ... " << processedElements << " / " << matrix.getNumberOfMatrixElements() << "                       \r" << flush;
    }
    file.clear();
    long int fileSize = file.tellg();
    if( verbose )
-      cout << " Reading the matrix elements ... " << processedElements << " / " << matrix.getNumberOfAllocatedElements()
+      cout << " Reading the matrix elements ... " << processedElements << " / " << matrix.getNumberOfMatrixElements()
               << " -> " << timer.GetTime()
               << " sec. i.e. " << fileSize / ( timer.GetTime() * ( 1 << 20 ))  << "MB/s." << endl;
    return true;
diff --git a/src/implementation/matrices/tnlMatrix_impl.h b/src/implementation/matrices/tnlMatrix_impl.h
index 346eb0fce8..0b8fa5f6ad 100644
--- a/src/implementation/matrices/tnlMatrix_impl.h
+++ b/src/implementation/matrices/tnlMatrix_impl.h
@@ -18,7 +18,88 @@
 #ifndef TNLMATRIX_IMPL_H_
 #define TNLMATRIX_IMPL_H_
 
+#include <matrices/tnlMatrix.h>
 
+template< typename Real,
+          typename Device,
+          typename Index >
+tnlMatrix< Real, Device, Index >::tnlMatrix()
+: rows( 0 ),
+  columns( 0 )
+{
+}
 
+template< typename Real,
+          typename Device,
+          typename Index >
+ bool tnlMatrix< Real, Device, Index >::setDimensions( const IndexType rows,
+                                                       const IndexType columns )
+{
+   tnlAssert( rows > 0 && columns > 0,
+            cerr << " rows = " << rows << " columns = " << columns );
+   this->rows = rows;
+   this->columns = columns;
+   return true;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+   template< typename Real2,
+             typename Device2,
+             typename Index2 >
+bool tnlMatrix< Real, Device, Index >::setLike( const tnlMatrix< Real2, Device2, Index2 >& matrix )
+{
+   return setDimensions( matrix.getRows(), matrix.getColumns() );
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+Index tnlMatrix< Real, Device, Index >::getRows() const
+{
+   return this->rows;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+Index tnlMatrix< Real, Device, Index >::getColumns() const
+{
+   return this->columns;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+void tnlMatrix< Real, Device, Index >::reset()
+{
+   this->rows = 0;
+   this->columns = 0;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+bool tnlMatrix< Real, Device, Index >::save( tnlFile& file ) const
+{
+   if( ! tnlObject::save( file ) ||
+       ! file.write( &this->rows ) ||
+       ! file.write( &this->columns ) )
+      return false;
+   return true;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+bool tnlMatrix< Real, Device, Index >::load( tnlFile& file )
+{
+   if( ! tnlObject::load( file ) ||
+       ! file.read( &this->rows ) ||
+       ! file.read( &this->columns ) )
+      return false;
+   return true;
+}
 
 #endif /* TNLMATRIX_IMPL_H_ */
diff --git a/src/implementation/matrices/tnlMultidiagonalMatrix_impl.h b/src/implementation/matrices/tnlMultidiagonalMatrix_impl.h
index 78735c1542..512f12ff5b 100644
--- a/src/implementation/matrices/tnlMultidiagonalMatrix_impl.h
+++ b/src/implementation/matrices/tnlMultidiagonalMatrix_impl.h
@@ -26,8 +26,6 @@ template< typename Real,
           typename Device,
           typename Index >
 tnlMultidiagonalMatrix< Real, Device, Index > :: tnlMultidiagonalMatrix()
-: rows( 0 ),
-  columns( 0 )
 {
 };
 
@@ -54,14 +52,14 @@ tnlString tnlMultidiagonalMatrix< Real, Device, Index >::getTypeVirtual() const
 template< typename Real,
           typename Device,
           typename Index >
-bool tnlMultidiagonalMatrix< Real, Device, Index > :: setDimensions( const IndexType rows,
-                                                                     const IndexType columns )
+bool tnlMultidiagonalMatrix< Real, Device, Index >::setDimensions( const IndexType rows,
+                                                                   const IndexType columns )
 {
    tnlAssert( rows > 0 && columns > 0,
               cerr << "rows = " << rows
                    << " columns = " << columns << endl );
-   this->rows = rows;
-   this->columns = columns;
+   if( ! tnlMatrix< Real, Device, Index >::setDimensions( rows, columns ) )
+      return false;
    if( this->diagonalsShift.getSize() != 0 )
    {
       if( ! this->values.setSize( Min( this->rows, this->columns ) * this->diagonalsShift.getSize() ) )
@@ -71,6 +69,17 @@ bool tnlMultidiagonalMatrix< Real, Device, Index > :: setDimensions( const Index
    return true;
 }
 
+template< typename Real,
+          typename Device,
+          typename Index >
+bool tnlMultidiagonalMatrix< Real, Device, Index >::setRowLengths( const RowLengthsVector& rowLengths )
+{
+   /****
+    * TODO: implement some check here similar to the one in the tridiagonal matrix
+    */
+   return true;
+}
+
 template< typename Real,
           typename Device,
           typename Index >
@@ -116,7 +125,7 @@ bool tnlMultidiagonalMatrix< Real, Device, Index > :: setLike( const tnlMultidia
 template< typename Real,
           typename Device,
           typename Index >
-Index tnlMultidiagonalMatrix< Real, Device, Index > :: getNumberOfAllocatedElements() const
+Index tnlMultidiagonalMatrix< Real, Device, Index > :: getNumberOfMatrixElements() const
 {
    return this->values.getSize();
 }
@@ -131,22 +140,6 @@ void tnlMultidiagonalMatrix< Real, Device, Index > :: reset()
    this->values.reset();
 }
 
-template< typename Real,
-          typename Device,
-          typename Index >
-Index tnlMultidiagonalMatrix< Real, Device, Index >::getRows() const
-{
-   return this->rows;
-}
-
-template< typename Real,
-          typename Device,
-          typename Index >
-Index tnlMultidiagonalMatrix< Real, Device, Index >::getColumns() const
-{
-   return this->columns;
-}
-
 template< typename Real,
           typename Device,
           typename Index >
@@ -318,8 +311,7 @@ template< typename Real,
           typename Index >
 bool tnlMultidiagonalMatrix< Real, Device, Index >::save( tnlFile& file ) const
 {
-   if( ! file.write( &this->rows ) ) return false;
-   if( ! file.write( &this->columns ) ) return false;
+   if( ! tnlMatrix< Real, Device, Index >::save( file ) ) return false;
    if( ! this->values.save( file ) ) return false;
    if( ! this->diagonalsShift.save( file ) ) return false;
    return true;
@@ -330,8 +322,7 @@ template< typename Real,
           typename Index >
 bool tnlMultidiagonalMatrix< Real, Device, Index >::load( tnlFile& file )
 {
-   if( ! file.read( &this->rows ) ) return false;
-   if( ! file.read( &this->columns ) ) return false;
+   if( ! tnlMatrix< Real, Device, Index >::load( file ) ) return false;
    if( ! this->values.load( file ) ) return false;
    if( ! this->diagonalsShift.load( file ) ) return false;
    return true;
@@ -364,7 +355,7 @@ void tnlMultidiagonalMatrix< Real, Device, Index >::print( ostream& str ) const
       for( IndexType i = 0; i < this->diagonalsShift.getSize(); i++ )
       {
          const IndexType column = row + diagonalsShift[ i ];
-         if( column >=0 && columns < this-columns )
+         if( column >=0 && column < this->columns )
             str << " Col:" << column << "->" << this->operator()( row, column ) << "\t";
       }
       str << endl;
diff --git a/src/implementation/matrices/tnlSlicedEllpackMatrix_impl.h b/src/implementation/matrices/tnlSlicedEllpackMatrix_impl.h
index 55e2fc7e49..66553809b4 100644
--- a/src/implementation/matrices/tnlSlicedEllpackMatrix_impl.h
+++ b/src/implementation/matrices/tnlSlicedEllpackMatrix_impl.h
@@ -27,8 +27,6 @@ template< typename Real,
           typename Index,
           int SliceSize >
 tnlSlicedEllpackMatrix< Real, Device, Index, SliceSize >::tnlSlicedEllpackMatrix()
-: rows( 0 ),
-  columns( 0 )
 {
 };
 
@@ -64,17 +62,14 @@ bool tnlSlicedEllpackMatrix< Real, Device, Index, SliceSize >::setDimensions( co
    tnlAssert( rows > 0 && columns > 0,
               cerr << "rows = " << rows
                    << " columns = " << columns << endl );
-   this->rows = rows;
-   this->columns = columns;
-   return true;
+   return tnlSparseMatrix< Real, Device, Index >::setDimensions( rows, columns );
 }
 
 template< typename Real,
           typename Device,
           typename Index,
           int SliceSize >
-   template< typename Vector >
-bool tnlSlicedEllpackMatrix< Real, Device, Index, SliceSize >::setRowLengths( const Vector& rowLengths )
+bool tnlSlicedEllpackMatrix< Real, Device, Index, SliceSize >::setRowLengths( const RowLengthsVector& rowLengths )
 {
    const IndexType slices = roundUpDivision( this->rows, SliceSize );
    if( ! this->sliceRowLengths.setSize( slices ) ||
@@ -107,14 +102,7 @@ bool tnlSlicedEllpackMatrix< Real, Device, Index, SliceSize >::setRowLengths( co
    this->slicePointers.setElement( slices, 0 );
    this->slicePointers.computeExclusivePrefixSum();
 
-   /****
-    * Allocate values and column indexes
-    */
-   if( ! this->values.setSize( this->slicePointers[ slices ] ) ||
-       ! this->columnIndexes.setSize( this->slicePointers[ slices ] ) )
-      return false;
-   this->columnIndexes.setValue( this->columns );
-   return true;
+   return this->allocateMatrixElements( this->slicePointers[ slices ] );
 }
 
 template< typename Real,
@@ -126,56 +114,24 @@ template< typename Real,
              typename Index2 >
 bool tnlSlicedEllpackMatrix< Real, Device, Index, SliceSize >::setLike( const tnlSlicedEllpackMatrix< Real2, Device2, Index2, SliceSize >& matrix )
 {
-   if( ! this->setDimensions( matrix.getRows(), matrix.getColumns() ) ||
-       ! this->values.setLike( matrix.values ) ||
-       ! this->columnIndexes.setLike( matrix.columnIndexes ) ||
+   if( !tnlSparseMatrix< Real, Device, Index >::setLike( matrix ) ||
        ! this->slicePointers.setLike( matrix.slicePointers ) ||
        ! this->sliceRowLengths.setLike( matrix.sliceRowLengths ) )
       return false;
    return true;
 }
 
-template< typename Real,
-          typename Device,
-          typename Index,
-          int SliceSize >
-Index tnlSlicedEllpackMatrix< Real, Device, Index, SliceSize >::getNumberOfAllocatedElements() const
-{
-   return this->values.getSize();
-}
-
 template< typename Real,
           typename Device,
           typename Index,
           int SliceSize >
 void tnlSlicedEllpackMatrix< Real, Device, Index, SliceSize >::reset()
 {
-   this->columns = 0;
-   this->rows = 0;
-   this->values.reset();
-   this->columnIndexes.reset();
+   tnlSparseMatrix< Real, Device, Index >::reset();
    this->slicePointers.reset();
    this->sliceRowLengths.reset();
 }
 
-template< typename Real,
-          typename Device,
-          typename Index,
-          int SliceSize >
-Index tnlSlicedEllpackMatrix< Real, Device, Index, SliceSize >::getRows() const
-{
-   return this->rows;
-}
-
-template< typename Real,
-          typename Device,
-          typename Index,
-          int SliceSize >
-Index tnlSlicedEllpackMatrix< Real, Device, Index, SliceSize >::getColumns() const
-{
-   return this->columns;
-}
-
 template< typename Real,
           typename Device,
           typename Index,
@@ -233,7 +189,7 @@ Real tnlSlicedEllpackMatrix< Real, Device, Index, SliceSize >::getElement( const
    const IndexType rowEnd = elementPtr + rowLength;
    while( elementPtr < rowEnd && this->columnIndexes[ elementPtr ] < column )
       elementPtr++;
-   if( this->columnIndexes[ elementPtr ] == column )
+   if( elementPtr < rowEnd && this->columnIndexes[ elementPtr ] == column )
       return this->values[ elementPtr ];
    return 0.0;
 }
@@ -391,12 +347,9 @@ template< typename Real,
           int SliceSize >
 bool tnlSlicedEllpackMatrix< Real, Device, Index, SliceSize >::save( tnlFile& file ) const
 {
-   if( ! file.write( &this->rows ) ) return false;
-   if( ! file.write( &this->columns ) ) return false;
-   if( ! this->values.save( file ) ) return false;
-   if( ! this->columnIndexes.save( file ) ) return false;
-   if( ! this->slicePointers.save( file ) ) return false;
-   if( ! this->sliceRowLengths.save( file ) ) return false;
+   if( ! tnlSparseMatrix< Real, Device, Index >::save( file ) ||
+       ! this->slicePointers.save( file ) ||
+       ! this->sliceRowLengths.save( file ) )
    return true;
 }
 
@@ -406,12 +359,9 @@ template< typename Real,
           int SliceSize >
 bool tnlSlicedEllpackMatrix< Real, Device, Index, SliceSize >::load( tnlFile& file )
 {
-   if( ! file.read( &this->rows ) ) return false;
-   if( ! file.read( &this->columns ) ) return false;
-   if( ! this->values.load( file ) ) return false;
-   if( ! this->columnIndexes.load( file ) ) return false;
-   if( ! this->slicePointers.load( file ) ) return false;
-   if( ! this->sliceRowLengths.load( file ) ) return false;
+   if( ! tnlSparseMatrix< Real, Device, Index >::load( file ) ||
+       ! this->slicePointers.load( file ) ||
+       ! this->sliceRowLengths.load( file ) )
    return true;
 }
 
diff --git a/src/implementation/matrices/tnlSparseMatrix_impl.h b/src/implementation/matrices/tnlSparseMatrix_impl.h
new file mode 100644
index 0000000000..eb71680da0
--- /dev/null
+++ b/src/implementation/matrices/tnlSparseMatrix_impl.h
@@ -0,0 +1,102 @@
+/***************************************************************************
+                          tnlSparseMatrix_impl.h  -  description
+                             -------------------
+    begin                : Dec 21, 2013
+    copyright            : (C) 2013 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef TNLSPARSEMATRIX_IMPL_H_
+#define TNLSPARSEMATRIX_IMPL_H_
+
+template< typename Real,
+          typename Device,
+          typename Index >
+tnlSparseMatrix< Real, Device, Index >::tnlSparseMatrix()
+{
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+   template< typename Real2,
+             typename Device2,
+             typename Index2 >
+bool tnlSparseMatrix< Real, Device, Index >::setLike( const tnlSparseMatrix< Real2, Device2, Index2 >& matrix )
+{
+   if( ! tnlMatrix< Real, Device, Index >::setLike( matrix ) ||
+       ! this->allocateMatrixElements( matrix.getNumberOfMatrixElements() ) )
+      return false;
+   return true;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+Index tnlSparseMatrix< Real, Device, Index >::getNumberOfMatrixElements() const
+{
+   return this->values.getSize();
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+void tnlSparseMatrix< Real, Device, Index >::reset()
+{
+   tnlMatrix< Real, Device, Index >::reset();
+   this->values.reset();
+   this->columnIndexes.reset();
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+bool tnlSparseMatrix< Real, Device, Index >::save( tnlFile& file ) const
+{
+   if( ! tnlMatrix< Real, Device, Index >::save( file ) ||
+       ! this->values.save( file ) ||
+       ! this->columnIndexes.save( file ) )
+      return false;
+   return true;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+bool tnlSparseMatrix< Real, Device, Index >::load( tnlFile& file )
+{
+   if( ! tnlMatrix< Real, Device, Index >::load( file ) ||
+       ! this->values.load( file ) ||
+       ! this->columnIndexes.load( file ) )
+      return false;
+   return true;
+}
+
+template< typename Real,
+          typename Device,
+          typename Index >
+bool tnlSparseMatrix< Real, Device, Index >::allocateMatrixElements( const IndexType& numberOfMatrixElements )
+{
+   if( ! this->values.setSize( numberOfMatrixElements ) ||
+       ! this->columnIndexes.setSize( numberOfMatrixElements ) )
+      return false;
+
+   /****
+    * Setting a column index to this->columns means that the
+    * index is undefined.
+    */
+   this->columnIndexes.setValue( this->columns );
+   return true;
+}
+
+
+#endif /* TNLSPARSEMATRIX_IMPL_H_ */
diff --git a/src/implementation/matrices/tnlTridiagonalMatrix_impl.h b/src/implementation/matrices/tnlTridiagonalMatrix_impl.h
index 704e35988b..bd051a3a83 100644
--- a/src/implementation/matrices/tnlTridiagonalMatrix_impl.h
+++ b/src/implementation/matrices/tnlTridiagonalMatrix_impl.h
@@ -25,7 +25,6 @@ template< typename Real,
           typename Device,
           typename Index >
 tnlTridiagonalMatrix< Real, Device, Index >::tnlTridiagonalMatrix()
-: rows( 0 )
 {
 }
 
@@ -51,55 +50,63 @@ tnlString tnlTridiagonalMatrix< Real, Device, Index >::getTypeVirtual() const
 template< typename Real,
           typename Device,
           typename Index >
-bool tnlTridiagonalMatrix< Real, Device, Index >::setDimensions( const IndexType rows )
+bool tnlTridiagonalMatrix< Real, Device, Index >::setDimensions( const IndexType rows,
+                                                                 const IndexType columns )
 {
-   if( ! values.setSize( 3*rows - 2 ) )
+   if( ! tnlMatrix< Real, Device, Index >::setDimensions( rows, columns ) )
+      return false;
+   if( ! values.setSize( 3*Min( rows, columns ) ) )
       return false;
-   this->rows = rows;
-   this->columns = rows;
    this->values.setValue( 0.0 );
 }
 
 template< typename Real,
           typename Device,
           typename Index >
-   template< typename Real2, typename Device2, typename Index2 >
-bool tnlTridiagonalMatrix< Real, Device, Index >::setLike( const tnlTridiagonalMatrix< Real2, Device2, Index2 >& m )
+bool tnlTridiagonalMatrix< Real, Device, Index >::setRowLengths( const RowLengthsVector& rowLengths )
 {
-   return this->setDimensions( m.getRows() );
-}
-
-template< typename Real,
-          typename Device,
-          typename Index >
-Index tnlTridiagonalMatrix< Real, Device, Index >::getNumberOfAllocatedElements() const
-{
-   return 3 * rows - 2;
+   if( rowLengths[ 0 ] > 2 )
+      return false;
+   const IndexType diagonalLength = Min( this->getRows(), this->getColumns() );
+   for( Index i = 1; i < diagonalLength-1; i++ )
+      if( rowLengths[ i ] > 3 )
+         return false;
+   if( this->getRows() > this->getColumns() )
+      if( rowLengths[ this->getRows()-1 ] > 1 )
+         return false;
+   if( this->getRows() == this->getColumns() )
+      if( rowLengths[ this->getRows()-1 ] > 2 )
+         return false;
+   if( this->getRows() < this->getColumns() )
+      if( rowLengths[ this->getRows()-1 ] > 3 )
+         return false;
+   return true;
 }
 
 template< typename Real,
           typename Device,
           typename Index >
-void tnlTridiagonalMatrix< Real, Device, Index >::reset()
+   template< typename Real2, typename Device2, typename Index2 >
+bool tnlTridiagonalMatrix< Real, Device, Index >::setLike( const tnlTridiagonalMatrix< Real2, Device2, Index2 >& m )
 {
-   this->rows = 0;
-   this->values.reset();
+   return this->setDimensions( m.getRows(), m.getColumns() );
 }
 
 template< typename Real,
           typename Device,
           typename Index >
-Index tnlTridiagonalMatrix< Real, Device, Index >::getRows() const
+Index tnlTridiagonalMatrix< Real, Device, Index >::getNumberOfMatrixElements() const
 {
-   return this->rows;
+   return 3 * Min( this->getRows(), this->getColumns() );
 }
 
 template< typename Real,
           typename Device,
           typename Index >
-Index tnlTridiagonalMatrix< Real, Device, Index >::getColumns() const
+void tnlTridiagonalMatrix< Real, Device, Index >::reset()
 {
-   return this->rows;
+   tnlMatrix< Real, Device, Index >::reset();
+   this->values.reset();
 }
 
 template< typename Real,
@@ -131,11 +138,12 @@ void tnlTridiagonalMatrix< Real, Device, Index >::setValue( const RealType& v )
 template< typename Real,
           typename Device,
           typename Index >
-void tnlTridiagonalMatrix< Real, Device, Index >::setElement( const IndexType row,
+bool tnlTridiagonalMatrix< Real, Device, Index >::setElement( const IndexType row,
                                                               const IndexType column,
                                                               const RealType& value )
 {
    this->values.setElement( this->getElementIndex( row, column ), value );
+   return true;
 }
 
 template< typename Real,
@@ -332,7 +340,7 @@ template< typename Real,
           typename Index >
 bool tnlTridiagonalMatrix< Real, Device, Index >::save( tnlFile& file ) const
 {
-   if( ! file.write( &this->rows ) ||
+   if( ! tnlMatrix< Real, Device, Index >::save( file ) ||
        ! this->values.save( file ) )
    {
       cerr << "Unable to save the tridiagonal matrix " << this->getName() << "." << endl;
@@ -346,7 +354,7 @@ template< typename Real,
           typename Index >
 bool tnlTridiagonalMatrix< Real, Device, Index >::load( tnlFile& file )
 {
-   if( ! file.read( &this->rows ) ||
+   if( ! tnlMatrix< Real, Device, Index >::load( file ) ||
        ! this->values.load( file ) )
    {
       cerr << "Unable to save the tridiagonal matrix " << this->getName() << "." << endl;
diff --git a/src/legacy/matrices/tnlCusparseCSRMatrix.h b/src/legacy/matrices/tnlCusparseCSRMatrix.h
index f40d59ef29..80b1b4d82f 100644
--- a/src/legacy/matrices/tnlCusparseCSRMatrix.h
+++ b/src/legacy/matrices/tnlCusparseCSRMatrix.h
@@ -26,7 +26,6 @@
 #include <core/tnlAssert.h>
 #include <core/mfuncs.h>
 #include <matrices/tnlCSRMatrix.h>
-#include <matrices/tnlRgCSRMatrix.h>
 #include <debug/tnlDebug.h>
 
 #ifdef HAVE_CUSPARSE
diff --git a/src/matrices/CMakeLists.txt b/src/matrices/CMakeLists.txt
index 0cdd303654..00caba1fb8 100755
--- a/src/matrices/CMakeLists.txt
+++ b/src/matrices/CMakeLists.txt
@@ -2,6 +2,7 @@ SET( headers tnlMatrix.h
              tnlDenseMatrix.h
              tnlTridiagonalMatrix.h
              tnlMultidiagonalMatrix.h
+             tnlSparseMatrix.h
              tnlEllpackMatrix.h
              tnlSlicedEllpackMatrix.h
              tnlChunkedEllpackMatrix.h
diff --git a/src/matrices/tnlCSRMatrix.h b/src/matrices/tnlCSRMatrix.h
index 3382173a0d..11fdffd9d1 100644
--- a/src/matrices/tnlCSRMatrix.h
+++ b/src/matrices/tnlCSRMatrix.h
@@ -18,16 +18,18 @@
 #ifndef TNLCSRMATRIX_H_
 #define TNLCSRMATRIX_H_
 
+#include <matrices/tnlSparseMatrix.h>
 #include <core/vectors/tnlVector.h>
 
 template< typename Real, typename Device = tnlHost, typename Index = int >
-class tnlCSRMatrix : public tnlObject
+class tnlCSRMatrix : public tnlSparseMatrix< Real, Device, Index >
 {
    public:
 
    typedef Real RealType;
    typedef Device DeviceType;
    typedef Index IndexType;
+   typedef typename tnlSparseMatrix< RealType, DeviceType, IndexType >:: RowLengthsVector RowLengthsVector;
 
    tnlCSRMatrix();
 
@@ -38,20 +40,13 @@ class tnlCSRMatrix : public tnlObject
    bool setDimensions( const IndexType rows,
                        const IndexType columns );
 
-   template< typename Vector >
-   bool setRowLengths( const Vector& rowLengths );
+   bool setRowLengths( const RowLengthsVector& rowLengths );
 
    template< typename Real2, typename Device2, typename Index2 >
    bool setLike( const tnlCSRMatrix< Real2, Device2, Index2 >& matrix );
 
-   IndexType getNumberOfAllocatedElements() const;
-
    void reset();
 
-   IndexType getRows() const;
-
-   IndexType getColumns() const;
-
    template< typename Real2, typename Device2, typename Index2 >
    bool operator == ( const tnlCSRMatrix< Real2, Device2, Index2 >& matrix ) const;
 
@@ -110,11 +105,7 @@ class tnlCSRMatrix : public tnlObject
 
    protected:
 
-   IndexType rows, columns;
-
-   tnlVector< Real, Device, Index > values;
-
-   tnlVector< Index, Device, Index > columnIndexes, rowPointers;
+   tnlVector< Index, Device, Index > rowPointers;
 
 };
 
diff --git a/src/matrices/tnlChunkedEllpackMatrix.h b/src/matrices/tnlChunkedEllpackMatrix.h
index c906b3560b..e78b958166 100644
--- a/src/matrices/tnlChunkedEllpackMatrix.h
+++ b/src/matrices/tnlChunkedEllpackMatrix.h
@@ -18,16 +18,18 @@
 #ifndef TNLCHUNKEDELLPACKMATRIX_H_
 #define TNLCHUNKEDELLPACKMATRIX_H_
 
+#include <matrices/tnlSparseMatrix.h>
 #include <core/vectors/tnlVector.h>
 
 template< typename Real, typename Device = tnlHost, typename Index = int >
-class tnlChunkedEllpackMatrix : public tnlObject
+class tnlChunkedEllpackMatrix : public tnlSparseMatrix< Real, Device, Index >
 {
    public:
 
    typedef Real RealType;
    typedef Device DeviceType;
    typedef Index IndexType;
+   typedef typename tnlSparseMatrix< RealType, DeviceType, IndexType >:: RowLengthsVector RowLengthsVector;
 
    tnlChunkedEllpackMatrix();
 
@@ -38,20 +40,13 @@ class tnlChunkedEllpackMatrix : public tnlObject
    bool setDimensions( const IndexType rows,
                        const IndexType columns );
 
-   template< typename Vector >
-   bool setRowLengths( const Vector& rowLengths );
+   bool setRowLengths( const RowLengthsVector& rowLengths );
 
    template< typename Real2, typename Device2, typename Index2 >
    bool setLike( const tnlChunkedEllpackMatrix< Real2, Device2, Index2 >& matrix );
 
-   IndexType getNumberOfAllocatedElements() const;
-
    void reset();
 
-   IndexType getRows() const;
-
-   IndexType getColumns() const;
-
    template< typename Real2, typename Device2, typename Index2 >
    bool operator == ( const tnlChunkedEllpackMatrix< Real2, Device2, Index2 >& matrix ) const;
 
@@ -129,13 +124,9 @@ class tnlChunkedEllpackMatrix : public tnlObject
       { return tnlString( "tnlChunkedEllpackSliceInfo" ); };
    };
 
-   IndexType rows, columns;
-
    IndexType chunksInSlice, desiredChunkSize;
 
-   tnlVector< Real, Device, Index > values;
-
-   tnlVector< Index, Device, Index > columnIndexes, chunksToRowsMapping, slicesToRowsMapping, rowPointers;
+   tnlVector< Index, Device, Index > chunksToRowsMapping, slicesToRowsMapping, rowPointers;
 
    tnlArray< tnlChunkedEllpackSliceInfo, Device, Index > slices;
 
diff --git a/src/matrices/tnlDenseMatrix.h b/src/matrices/tnlDenseMatrix.h
index fbdc7cc755..0c22c13011 100644
--- a/src/matrices/tnlDenseMatrix.h
+++ b/src/matrices/tnlDenseMatrix.h
@@ -19,18 +19,21 @@
 #define TNLDENSEMATRIX_H_
 
 #include <core/tnlHost.h>
+#include <matrices/tnlMatrix.h>
 #include <core/arrays/tnlMultiArray.h>
 
 template< typename Real = double,
           typename Device = tnlHost,
           typename Index = int >
-class tnlDenseMatrix : public tnlMultiArray< 2, Real, Device, Index >
+class tnlDenseMatrix : public tnlMatrix< Real, Device, Index >,
+                       public tnlMultiArray< 2, Real, Device, Index >
 {
    public:
 
    typedef Real RealType;
    typedef Device DeviceType;
    typedef Index IndexType;
+   typedef typename tnlMatrix< Real, Device, Index >::RowLengthsVector RowLengthsVector;
 
    tnlDenseMatrix();
 
@@ -41,17 +44,24 @@ class tnlDenseMatrix : public tnlMultiArray< 2, Real, Device, Index >
    bool setDimensions( const IndexType rows,
                        const IndexType columns );
 
+   template< typename Real2, typename Device2, typename Index2 >
+   bool setLike( const tnlDenseMatrix< Real2, Device2, Index2 >& matrix );
+
    /****
     * This method is only for the compatibility with the sparse matrices.
     */
-   template< typename Vector >
-   bool setRowLengths( const Vector& rowLengths );
+   bool setRowLengths( const RowLengthsVector& rowLengths );
 
-   IndexType getNumberOfAllocatedElements() const;
+   IndexType getNumberOfMatrixElements() const;
 
-   IndexType getRows() const;
+   void reset();
+
+   bool setElement( const IndexType row,
+                    const IndexType column,
+                    const RealType& value );
 
-   IndexType getColumns() const;
+   Real getElement( const IndexType row,
+                    const IndexType column ) const;
 
    bool addElement( const IndexType row,
                     const IndexType column,
@@ -98,6 +108,14 @@ class tnlDenseMatrix : public tnlMultiArray< 2, Real, Device, Index >
                              Vector& x,
                              const RealType& omega = 1.0 ) const;
 
+   bool save( const tnlString& fileName ) const;
+
+   bool load( const tnlString& fileName );
+
+   bool save( tnlFile& file ) const;
+
+   bool load( tnlFile& file );
+
    void print( ostream& str ) const;
 };
 
diff --git a/src/matrices/tnlEllpackMatrix.h b/src/matrices/tnlEllpackMatrix.h
index ae848ad19e..194b3971e6 100644
--- a/src/matrices/tnlEllpackMatrix.h
+++ b/src/matrices/tnlEllpackMatrix.h
@@ -18,16 +18,18 @@
 #ifndef TNLELLPACKMATRIX_H_
 #define TNLELLPACKMATRIX_H_
 
+#include <matrices/tnlSparseMatrix.h>
 #include <core/vectors/tnlVector.h>
 
 template< typename Real, typename Device = tnlHost, typename Index = int >
-class tnlEllpackMatrix : public tnlObject
+class tnlEllpackMatrix : public tnlSparseMatrix< Real, Device, Index >
 {
    public:
 
    typedef Real RealType;
    typedef Device DeviceType;
    typedef Index IndexType;
+   typedef typename tnlSparseMatrix< RealType, DeviceType, IndexType >:: RowLengthsVector RowLengthsVector;
 
    tnlEllpackMatrix();
 
@@ -38,22 +40,15 @@ class tnlEllpackMatrix : public tnlObject
    bool setDimensions( const IndexType rows,
                        const IndexType columns );
 
-   template< typename Vector >
-   bool setRowLengths( const Vector& rowLengths );
+   bool setRowLengths( const RowLengthsVector& rowLengths );
 
    bool setConstantRowLengths( const IndexType& rowLengths );
 
    template< typename Real2, typename Device2, typename Index2 >
    bool setLike( const tnlEllpackMatrix< Real2, Device2, Index2 >& matrix );
 
-   IndexType getNumberOfAllocatedElements() const;
-
    void reset();
 
-   IndexType getRows() const;
-
-   IndexType getColumns() const;
-
    template< typename Real2, typename Device2, typename Index2 >
    bool operator == ( const tnlEllpackMatrix< Real2, Device2, Index2 >& matrix ) const;
 
@@ -110,12 +105,7 @@ class tnlEllpackMatrix : public tnlObject
 
    bool allocateElements();
 
-   IndexType rows, columns, rowLengths;
-
-   tnlVector< Real, Device, Index > values;
-
-   tnlVector< Index, Device, Index > columnIndexes;
-
+   IndexType rowLengths;
 };
 
 #include <implementation/matrices/tnlEllpackMatrix_impl.h>
diff --git a/src/matrices/tnlMatrix.h b/src/matrices/tnlMatrix.h
index 32a0369eb8..498c18289c 100644
--- a/src/matrices/tnlMatrix.h
+++ b/src/matrices/tnlMatrix.h
@@ -20,49 +20,60 @@
 
 #include <core/tnlObject.h>
 #include <core/tnlHost.h>
+#include <core/vectors/tnlVector.h>
 
 template< typename Real = double,
           typename Device = tnlHost,
           typename Index = int >
-class tnlMatrix : public tnlObject
+class tnlMatrix : public virtual tnlObject
 {
    public:
 
    typedef Real RealType;
    typedef Device DeviceType;
    typedef Index IndexType;
+   typedef tnlVector< IndexType, DeviceType, IndexType > RowLengthsVector;
+
+   tnlMatrix();
 
    virtual bool setDimensions( const IndexType rows,
                                const IndexType columns );
 
+   virtual bool setRowLengths( const RowLengthsVector& rowLengths ) = 0;
+
+   template< typename Real2, typename Device2, typename Index2 >
+   bool setLike( const tnlMatrix< Real2, Device2, Index2 >& matrix );
+
+   virtual IndexType getNumberOfMatrixElements() const = 0;
+
+   void reset();
+
    IndexType getRows() const;
 
    IndexType getColumns() const;
 
    virtual bool setElement( const IndexType row,
                             const IndexType column,
-                            const RealType& value );
+                            const RealType& value ) = 0;
 
    virtual bool addElement( const IndexType row,
                             const IndexType column,
                             const RealType& value,
-                            const RealType& thisElementMultiplicator = 1.0 );
+                            const RealType& thisElementMultiplicator = 1.0 ) = 0;
 
+   virtual Real getElement( const IndexType row,
+                            const IndexType column ) const = 0;
 
    virtual bool save( tnlFile& file ) const;
 
    virtual bool load( tnlFile& file );
 
-   virtual bool save( const tnlString& fileName ) const;
-
-   virtual bool load( const tnlString& fileName );
-
    protected:
 
    IndexType rows, columns;
 
 };
 
-#include <implementation/matrices/tnlMatrix.h>
+#include <implementation/matrices/tnlMatrix_impl.h>
 
 #endif /* TNLMATRIX_H_ */
diff --git a/src/matrices/tnlMatrixReader.h b/src/matrices/tnlMatrixReader.h
index 3d3b8d401f..5f2ee48cc5 100644
--- a/src/matrices/tnlMatrixReader.h
+++ b/src/matrices/tnlMatrixReader.h
@@ -37,6 +37,12 @@ class tnlMatrixReader
    static bool verifyMtxFile( std::istream& file,
                               const Matrix& matrix,
                               bool verbose = false );
+
+   static bool findLineByElement( std::istream& file,
+                                  const IndexType& row,
+                                  const IndexType& column,
+                                  tnlString& line,
+                                  IndexType& lineNumber );
    protected:
 
    static bool checkMtxHeader( const tnlString& header,
diff --git a/src/matrices/tnlMultidiagonalMatrix.h b/src/matrices/tnlMultidiagonalMatrix.h
index 09d91bc96b..f853096173 100644
--- a/src/matrices/tnlMultidiagonalMatrix.h
+++ b/src/matrices/tnlMultidiagonalMatrix.h
@@ -18,16 +18,18 @@
 #ifndef TNLMULTIDIAGONALMATRIX_H_
 #define TNLMULTIDIAGONALMATRIX_H_
 
+#include <matrices/tnlMatrix.h>
 #include <core/vectors/tnlVector.h>
 
 template< typename Real, typename Device = tnlHost, typename Index = int >
-class tnlMultidiagonalMatrix : public tnlObject
+class tnlMultidiagonalMatrix : public tnlMatrix< Real, Device, Index >
 {
    public:
 
    typedef Real RealType;
    typedef Device DeviceType;
    typedef Index IndexType;
+   typedef typename tnlMatrix< Real, Device, Index >::RowLengthsVector RowLengthsVector;
 
    tnlMultidiagonalMatrix();
 
@@ -38,6 +40,8 @@ class tnlMultidiagonalMatrix : public tnlObject
    bool setDimensions( const IndexType rows,
                        const IndexType columns );
 
+   bool setRowLengths( const RowLengthsVector& rowLengths );
+
    bool setDiagonals( const IndexType diagonalsNumber,
                       const IndexType* diagonalsShift );
 
@@ -46,14 +50,10 @@ class tnlMultidiagonalMatrix : public tnlObject
    template< typename Real2, typename Device2, typename Index2 >
    bool setLike( const tnlMultidiagonalMatrix< Real2, Device2, Index2 >& matrix );
 
-   IndexType getNumberOfAllocatedElements() const;
+   IndexType getNumberOfMatrixElements() const;
 
    void reset();
 
-   IndexType getRows() const;
-
-   IndexType getColumns() const;
-
    template< typename Real2, typename Device2, typename Index2 >
    bool operator == ( const tnlMultidiagonalMatrix< Real2, Device2, Index2 >& matrix ) const;
 
@@ -109,8 +109,6 @@ class tnlMultidiagonalMatrix : public tnlObject
                          const IndexType column,
                          IndexType& index ) const;
 
-   IndexType rows, columns;
-
    tnlVector< Real, Device, Index > values;
 
    tnlVector< Index, Device, Index > diagonalsShift;
diff --git a/src/matrices/tnlSlicedEllpackMatrix.h b/src/matrices/tnlSlicedEllpackMatrix.h
index 49de206492..17ad1b9508 100644
--- a/src/matrices/tnlSlicedEllpackMatrix.h
+++ b/src/matrices/tnlSlicedEllpackMatrix.h
@@ -18,16 +18,21 @@
 #ifndef TNLSLICEDELLPACKMATRIX_H_
 #define TNLSLICEDELLPACKMATRIX_H_
 
+#include <matrices/tnlSparseMatrix.h>
 #include <core/vectors/tnlVector.h>
 
-template< typename Real, typename Device = tnlHost, typename Index = int, int SliceSize = 32 >
-class tnlSlicedEllpackMatrix : public tnlObject
+template< typename Real = double,
+          typename Device = tnlHost,
+          typename Index = int,
+          int SliceSize = 32 >
+class tnlSlicedEllpackMatrix : public tnlSparseMatrix< Real, Device, Index >
 {
    public:
 
    typedef Real RealType;
    typedef Device DeviceType;
    typedef Index IndexType;
+   typedef typename tnlSparseMatrix< RealType, DeviceType, IndexType >:: RowLengthsVector RowLengthsVector;
 
    tnlSlicedEllpackMatrix();
 
@@ -38,20 +43,13 @@ class tnlSlicedEllpackMatrix : public tnlObject
    bool setDimensions( const IndexType rows,
                        const IndexType columns );
 
-   template< typename Vector >
-   bool setRowLengths( const Vector& rowLengths );
+   bool setRowLengths( const RowLengthsVector& rowLengths );
 
    template< typename Real2, typename Device2, typename Index2 >
    bool setLike( const tnlSlicedEllpackMatrix< Real2, Device2, Index2, SliceSize >& matrix );
 
-   IndexType getNumberOfAllocatedElements() const;
-
    void reset();
 
-   IndexType getRows() const;
-
-   IndexType getColumns() const;
-
    template< typename Real2, typename Device2, typename Index2 >
    bool operator == ( const tnlSlicedEllpackMatrix< Real2, Device2, Index2 >& matrix ) const;
 
@@ -106,11 +104,7 @@ class tnlSlicedEllpackMatrix : public tnlObject
 
    protected:
 
-   IndexType rows, columns;
-
-   tnlVector< Real, Device, Index > values;
-
-   tnlVector< Index, Device, Index > columnIndexes, slicePointers, sliceRowLengths;
+   tnlVector< Index, Device, Index > slicePointers, sliceRowLengths;
 
 };
 
diff --git a/src/matrices/tnlSparseMatrix.h b/src/matrices/tnlSparseMatrix.h
new file mode 100644
index 0000000000..b9ccc26e0c
--- /dev/null
+++ b/src/matrices/tnlSparseMatrix.h
@@ -0,0 +1,59 @@
+/***************************************************************************
+                          tnlSparseMatrix.h  -  description
+                             -------------------
+    begin                : Dec 21, 2013
+    copyright            : (C) 2013 by Tomas Oberhuber
+    email                : tomas.oberhuber@fjfi.cvut.cz
+ ***************************************************************************/
+
+/***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef TNLSPARSEMATRIX_H_
+#define TNLSPARSEMATRIX_H_
+
+#include <matrices/tnlMatrix.h>
+
+template< typename Real,
+          typename Device,
+          typename Index >
+class tnlSparseMatrix : public tnlMatrix< Real, Device, Index >
+{
+   public:
+
+   typedef Real RealType;
+   typedef Device DeviceType;
+   typedef Index IndexType;
+   typedef typename tnlMatrix< RealType, DeviceType, IndexType >:: RowLengthsVector RowLengthsVector;
+
+   tnlSparseMatrix();
+
+   template< typename Real2, typename Device2, typename Index2 >
+   bool setLike( const tnlSparseMatrix< Real2, Device2, Index2 >& matrix );
+
+   IndexType getNumberOfMatrixElements() const;
+
+   void reset();
+
+   bool save( tnlFile& file ) const;
+
+   bool load( tnlFile& file );
+
+   protected:
+
+   bool allocateMatrixElements( const IndexType& numberOfMatrixElements );
+
+   tnlVector< Index, Device, Index > columnIndexes;
+
+   tnlVector< Real, Device, Index > values;
+};
+
+#include <implementation/matrices/tnlSparseMatrix_impl.h>
+
+#endif /* TNLSPARSEMATRIX_H_ */
diff --git a/src/matrices/tnlTridiagonalMatrix.h b/src/matrices/tnlTridiagonalMatrix.h
index 1b6bb14728..5d7328d229 100644
--- a/src/matrices/tnlTridiagonalMatrix.h
+++ b/src/matrices/tnlTridiagonalMatrix.h
@@ -18,20 +18,20 @@
 #ifndef TNLTRIDIAGONALMATRIX_H_
 #define TNLTRIDIAGONALMATRIX_H_
 
-#include <core/tnlObject.h>
-#include <core/tnlHost.h>
+#include <matrices/tnlMatrix.h>
 #include <core/vectors/tnlVector.h>
 
 template< typename Real = double,
           typename Device = tnlHost,
           typename Index = int >
-class tnlTridiagonalMatrix : public tnlObject
+class tnlTridiagonalMatrix : public tnlMatrix< Real, Device, Index >
 {
    public:
 
    typedef Real RealType;
    typedef Device DeviceType;
    typedef Index IndexType;
+   typedef typename tnlMatrix< Real, Device, Index >::RowLengthsVector RowLengthsVector;
 
    tnlTridiagonalMatrix();
 
@@ -39,19 +39,18 @@ class tnlTridiagonalMatrix : public tnlObject
 
    tnlString getTypeVirtual() const;
 
-   bool setDimensions( const IndexType rows );
+   bool setDimensions( const IndexType rows,
+                       const IndexType columns );
+
+   bool setRowLengths( const RowLengthsVector& rowLengths );
 
    template< typename Real2, typename Device2, typename Index2 >
    bool setLike( const tnlTridiagonalMatrix< Real2, Device2, Index2 >& m );
 
-   IndexType getNumberOfAllocatedElements() const;
+   IndexType getNumberOfMatrixElements() const;
 
    void reset();
 
-   IndexType getRows() const;
-
-   IndexType getColumns() const;
-
    template< typename Real2, typename Device2, typename Index2 >
    bool operator == ( const tnlTridiagonalMatrix< Real2, Device2, Index2 >& matrix ) const;
 
@@ -60,7 +59,7 @@ class tnlTridiagonalMatrix : public tnlObject
 
    void setValue( const RealType& v );
 
-   void setElement( const IndexType row,
+   bool setElement( const IndexType row,
                     const IndexType column,
                     const RealType& value );
 
@@ -127,8 +126,6 @@ class tnlTridiagonalMatrix : public tnlObject
    IndexType getElementIndex( const IndexType row,
                               const IndexType column ) const;
 
-   IndexType rows, columns;
-
    tnlVector< RealType, DeviceType, IndexType > values;
 };
 
diff --git a/tests/benchmarks/CMakeLists.txt b/tests/benchmarks/CMakeLists.txt
index 462ad789e9..19e8a854bb 100755
--- a/tests/benchmarks/CMakeLists.txt
+++ b/tests/benchmarks/CMakeLists.txt
@@ -9,13 +9,13 @@ SET( tnlSpmvBenchmark_headers sparse-matrix-benchmark.h
                               tnlSpmvBenchmarkRgCSRMatrix.h )
      
 
-IF( BUILD_CUDA )
-    #CUDA_ADD_EXECUTABLE( tnl-sparse-matrix-benchmark${debugExt} sparse-matrix-benchmark.cu )
-    #SET_TARGET_PROPERTIES( tnl-sparse-matrix-benchmark${debugExt} PROPERTIES CUDA_COMPILE_FLAGS "${CXX_OPTIMIZE_FLAGS}" )        
-ELSE()
-    #ADD_EXECUTABLE( tnl-sparse-matrix-benchmark${debugExt} sparse-matrix-benchmark.cpp )
-    #SET_TARGET_PROPERTIES( tnl-sparse-matrix-benchmark${debugExt} PROPERTIES COMPILE_FLAGS "${CXX_OPTIMIZE_FLAGS}" )
-ENDIF()
+#IF( BUILD_CUDA )
+#    CUDA_ADD_EXECUTABLE( tnl-sparse-matrix-benchmark${debugExt} sparse-matrix-benchmark.cu )
+#    SET_TARGET_PROPERTIES( tnl-sparse-matrix-benchmark${debugExt} PROPERTIES CUDA_COMPILE_FLAGS "${CXX_OPTIMIZE_FLAGS}" )        
+#ELSE()
+#    ADD_EXECUTABLE( tnl-sparse-matrix-benchmark${debugExt} sparse-matrix-benchmark.cpp )
+#    SET_TARGET_PROPERTIES( tnl-sparse-matrix-benchmark${debugExt} PROPERTIES COMPILE_FLAGS "${CXX_OPTIMIZE_FLAGS}" )
+#ENDIF()
 #TARGET_LINK_LIBRARIES( tnl-sparse-matrix-benchmark${debugExt} tnl${debugExt}-${tnlVersion}
 #                                                              ${CUSPARSE_LIBRARY} )
 
diff --git a/tests/benchmarks/matrix-solvers-benchmark.h b/tests/benchmarks/matrix-solvers-benchmark.h
index 213e7ef01f..d8f0d9ee7f 100644
--- a/tests/benchmarks/matrix-solvers-benchmark.h
+++ b/tests/benchmarks/matrix-solvers-benchmark.h
@@ -72,7 +72,7 @@ bool benchmarkSolver( const tnlParameterContainer&  parameters,
 
    const RealType& maxResidue = parameters. GetParameter< double >( "max-residue" );
    const IndexType& size = matrix. getRows();
-   const IndexType nonZeros = matrix. getNumberOfAllocatedElements();
+   const IndexType nonZeros = matrix. getNumberOfMatrixElements();
    //const IndexType maxIterations = size * ( ( double ) size * size / ( double ) nonZeros );
    const IndexType maxIterations = size;
    cout << "Setting max. number of iterations to " << maxIterations << endl;
@@ -277,7 +277,7 @@ bool benchmarkMatrix( const tnlParameterContainer&  parameters )
          return false;
       }
       matrixStatsFile << "             <td> " << csrMatrix. getRows() << " </td> " << endl
-                      << "             <td> " << csrMatrix. getNumberOfAllocatedElements() << " </td> " << endl;
+                      << "             <td> " << csrMatrix. getNumberOfMatrixElements() << " </td> " << endl;
       matrixStatsFile. close();
    }
 
diff --git a/tests/benchmarks/sparse-matrix-benchmark.h b/tests/benchmarks/sparse-matrix-benchmark.h
index 6b5743233f..104ff9a294 100644
--- a/tests/benchmarks/sparse-matrix-benchmark.h
+++ b/tests/benchmarks/sparse-matrix-benchmark.h
@@ -23,12 +23,11 @@
 #include <iomanip>
 #include <config/tnlConfigDescription.h>
 #include <config/tnlParameterContainer.h>
-#include <legacy/matrices/tnlFullMatrix.h>
-#include <legacy/matrices/tnlFastCSRMatrix.h>
-#include <legacy/matrices/tnlFastRgCSRMatrix.h>
-#include <legacy/matrices/tnlFastRgCSRMatrixCUDA.h>
-#include <legacy/matrices/tnlEllpackMatrix.h>
-#include <legacy/matrices/tnlEllpackMatrixCUDA.h>
+#include <matrices/tnlDenseMatrix.h>
+#include <matrices/tnlEllpackMatrix.h>
+#include <matrices/tnlSlicedEllpackMatrix.h>
+#include <matrices/tnlChunkedEllpackMatrix.h>
+#include <matrices/tnlCSRMatrix.h>
 #include <core/mfuncs.h>
 #include "tnlSpmvBenchmarkCSRMatrix.h"
 #include "tnlSpmvBenchmarkCusparseCSRMatrix.h"
@@ -102,38 +101,25 @@ void benchmarkRgCSRFormat( const tnlCSRMatrix< Real, tnlHost, int >& csrMatrix,
    }
 }
 
-template< class Real >
-bool benchmarkMatrix( const tnlString& inputFile,
-                      const tnlString& inputMtxFile,
-                      const tnlString& pdfFile,
-                      const tnlString& logFileName,
-                      bool formatTest,
-                      int maxIterations,
-                      int verbose )
+template< typename RealType >
+bool benchmarkMatrix( const tnlParameterContainer& parameters )
 {
    /****
     * Read the CSR matrix ...
     */
-   tnlCSRMatrix< Real > csrMatrix;
-   csrMatrix.setName( "csr-matrix" );
-   tnlString inputMtxSortedFile( inputMtxFile );
-   inputMtxSortedFile += tnlString( ".sort" );
-   tnlFile binaryFile;
-   if( ! binaryFile. open( inputFile, tnlReadMode ) )
-   {
-      cerr << "I am not able to open the file " << inputFile << "." << endl;
-      return 1;
-   }
-   if( verbose )
-      cout << "Reading the CSR matrix ... " << flush;
-   if( ! csrMatrix. load( binaryFile ) )
+   typedef tnlCSRMatrix< RealType, tnlHost, int > CsrMatrix;
+   CsrMatrix csrMatrix;
+
+   const tnlString& inputFileName = parameters.GetParameter< tnlString >( "input-file" );
+   fstream inputFile;
+   inputFile.open( inputFileName.getString(), ios::in );
+   if( ! inputFile )
    {
-      cerr << "Unable to restore the CSR matrix." << endl;
+      cerr << "I am not able to open the file " << inputFileName << "." << endl;
       return false;
    }
-   if( verbose )
-      cout << " OK." << endl;
-   binaryFile. close();
+   if( ! tnlMatrixReader< CsrMatrix >::readMtxFile( inputFile, csrMatrix ) )
+      return false;
 
    /****
     * Check the number of the non-zero elements
@@ -436,57 +422,14 @@ int main( int argc, char* argv[] )
       conf_desc. PrintUsage( argv[ 0 ] );
       return 1;
    }
-   tnlString inputFile = parameters. GetParameter< tnlString >( "input-file" );
-   tnlString inputMtxFile = parameters. GetParameter< tnlString >( "input-mtx-file" );
-   tnlString pdfFile = parameters. GetParameter< tnlString >( "pdf-file" );
-   tnlString logFileName = parameters. GetParameter< tnlString >( "log-file" );
-   double stop_time = parameters. GetParameter< double >( "stop-time" );
-   bool formatTest = parameters. GetParameter< bool >( "format-test" );
-   int maxIterations = parameters. GetParameter< int >( "max-iterations" );
-   int verbose = parameters. GetParameter< int >( "verbose");
-
-
-   tnlFile binaryFile;
-   if( ! binaryFile. open( inputFile, tnlReadMode ) )
-   {
-      cerr << "I am not able to open the file " << inputFile << "." << endl;
-      return 1;
-   }
-   tnlString object_type;
-   if( ! getObjectType( binaryFile, object_type ) )
-   {
-      cerr << "Unknown object ... SKIPPING!" << endl;
-      return EXIT_FAILURE;
-   }
-   if( verbose )
-      cout << object_type << " detected ... " << endl;
-   binaryFile. close();
-
-   if( object_type == "tnlCSRMatrix< float, tnlHost >")
-      benchmarkMatrix< float >( inputFile,
-                                inputMtxFile,
-                                pdfFile,
-                                logFileName,
-                                formatTest,
-                                maxIterations,
-                                verbose );
-
-   if( object_type == "tnlCSRMatrix< double, tnlHost >" )
-   {
-      benchmarkMatrix< double >( inputFile,
-                                 inputMtxFile,
-                                 pdfFile,
-                                 logFileName,
-                                 formatTest,
-                                 maxIterations,
-                                 verbose );
-   }
-   //binaryFile. close();
-
-
-
+   const tnlString& precision = parameters.GetParameter< tnlString >( "precision" );
+   if( precision == "float" )
+      if( ! benchmarkMatrix< float >( parameters ) )
+         return EXIT_FAILURE;
+   if( precision == "double" )
+      if( ! benchmarkMatrix< double >( parameters ) )
+         return EXIT_FAILURE;
    return EXIT_SUCCESS;
 }
 
-
 #endif /* SPARSEMATRIXBENCHMARK_H_ */
diff --git a/tests/benchmarks/tnlSpmvBenchmarkAdaptiveRgCSRMatrix.h b/tests/benchmarks/tnlSpmvBenchmarkAdaptiveRgCSRMatrix.h
index 81e17151b8..5d51a27231 100644
--- a/tests/benchmarks/tnlSpmvBenchmarkAdaptiveRgCSRMatrix.h
+++ b/tests/benchmarks/tnlSpmvBenchmarkAdaptiveRgCSRMatrix.h
@@ -19,7 +19,6 @@
 #define TNLSPMVBENCHMARKADAPTIVERGCSRMATRIX_H_
 
 #include "tnlSpmvBenchmark.h"
-#include <matrices/tnlAdaptiveRgCSRMatrix.h>
 #include <core/tnlAssert.h>
 
 template< typename Real, typename Device, typename Index>
diff --git a/tests/benchmarks/tnlSpmvBenchmarkRgCSRMatrix.h b/tests/benchmarks/tnlSpmvBenchmarkRgCSRMatrix.h
index 67be6da562..5209fededd 100644
--- a/tests/benchmarks/tnlSpmvBenchmarkRgCSRMatrix.h
+++ b/tests/benchmarks/tnlSpmvBenchmarkRgCSRMatrix.h
@@ -19,7 +19,6 @@
 #define TNLSPMVBENCHMARKRGCSRMATRIX_H_
 
 #include "tnlSpmvBenchmark.h"
-#include <matrices/tnlRgCSRMatrix.h>
 
 template< typename Real, typename Device, typename Index>
 class tnlSpmvBenchmarkRgCSRMatrix : public tnlSpmvBenchmark< Real, Device, Index, tnlRgCSRMatrix >
diff --git a/tests/data/tnl-test-matrix-formats.cfg.desc b/tests/data/tnl-test-matrix-formats.cfg.desc
index e166bc61db..bb1ceab808 100644
--- a/tests/data/tnl-test-matrix-formats.cfg.desc
+++ b/tests/data/tnl-test-matrix-formats.cfg.desc
@@ -1,7 +1,9 @@
 group Arguments
 {
-   string input-file            [Input file name.];
-   string matrix-format         [Matrix format. It can be: dense, ellpack, sliced-ellpack, chunked-ellpack and csr.];
-   bool verbose                 [Verbose mode.];
+   string input-file              [Input file name.];
+   string matrix-format           [Matrix format. It can be: dense, ellpack, sliced-ellpack, chunked-ellpack and csr.];
+   bool hard-test( no )           [Comparison against the dense matrix.];
+   bool multiplication-test( no ) [Matrix-vector multiplication test.];
+   bool verbose                   [Verbose mode.];
 },[Arguments of the matrix format test.];
 
diff --git a/tests/long-time-unit-tests/matrix-formats-test.h b/tests/long-time-unit-tests/matrix-formats-test.h
index 0d7591f130..5bdf13ac67 100644
--- a/tests/long-time-unit-tests/matrix-formats-test.h
+++ b/tests/long-time-unit-tests/matrix-formats-test.h
@@ -40,6 +40,10 @@ template< typename Matrix >
 bool testMatrix( const tnlParameterContainer& parameters )
 {
    Matrix matrix;
+   typedef typename Matrix::RealType RealType;
+   typedef typename Matrix::DeviceType DeviceType;
+   typedef typename Matrix::IndexType IndexType;
+
    const tnlString& fileName = parameters.GetParameter< tnlString >( "input-file" );
    bool verbose = parameters.GetParameter< bool >( "verbose" );
    fstream file;
@@ -50,16 +54,65 @@ bool testMatrix( const tnlParameterContainer& parameters )
       return false;
    }
    if( ! tnlMatrixReader< Matrix >::readMtxFile( file, matrix, verbose ) )
-   {
-      file.close();
       return false;
-   }
    if( ! tnlMatrixReader< Matrix >::verifyMtxFile( file, matrix, verbose ) )
-   {
-      file.close();
       return false;
+   if( parameters.GetParameter< bool >( "hard-test" ) )
+   {
+      typedef tnlDenseMatrix< RealType, DeviceType, IndexType > DenseMatrix;
+      DenseMatrix denseMatrix;
+      if( ! tnlMatrixReader< DenseMatrix >::readMtxFile( file, denseMatrix, verbose ) )
+         return false;
+      if( ! tnlMatrixReader< DenseMatrix >::verifyMtxFile( file, denseMatrix, verbose ) )
+         return false;
+      //matrix.print( cout );
+      //denseMatrix.print( cout );
+      for( IndexType i = 0; i < matrix.getRows(); i++ )
+      {
+         for( IndexType j = 0; j < matrix.getColumns(); j++ )
+            if( matrix.getElement( i, j ) != denseMatrix.getElement( i, j ) )
+            {
+               cerr << "The matrices differ at position " << i << ", " << j << "." << endl
+                    << " The values are " << matrix.getElement( i, j ) << " (sparse) and "
+                    << denseMatrix.getElement( i, j ) << " (dense)." << endl;
+               tnlString line;
+               IndexType lineNumber;
+               if( tnlMatrixReader< Matrix >::findLineByElement( file, i, j, line, lineNumber ) )
+                  cerr << "The mtx file says ( line " << lineNumber << " ): " << line << endl;
+               else
+                  cerr << "The element is missing in the file. Should be zero therefore." << endl;
+               return false;
+            }
+         if( verbose )
+            cout << " Comparing the sparse matrix with the dense matrix ... " << i << " / " << matrix.getRows() << "             \r" << flush;
+      }
+      if( verbose )
+         cout << " Comparing the sparse matrix with the dense matrix ... OK.           " << endl;
+   }
+   if( parameters.GetParameter< bool >( "multiplication-test" ) )
+   {
+      tnlVector< RealType, DeviceType, IndexType > x, b;
+      x.setSize( matrix.getColumns() );
+      b.setSize( matrix.getRows() );
+      for( IndexType i = 0; i < x.getSize(); i++ )
+      {
+         x.setValue( 0 );
+         x.setElement( i, 1.0 );
+         matrix.vectorProduct( x, b );
+         for( IndexType j = 0; j < b.getSize(); j++ )
+            if( b.getElement( j ) != matrix.getElement( j, i ) )
+            {
+               cerr << "The matrix-vector multiplication gives wrong result at positions "
+                    << j << ", " << i << ". The result is " << b.getElement( j ) << " and it should be "
+                    << matrix.getElement( j, i ) << "." << endl;
+               return false;
+            }
+         if( verbose )
+            cerr << " Testing the matrix-vector multiplication ... " << i << " / " << matrix.getRows() << "            \r" << flush;
+      }
+      if( verbose )
+         cerr << " Testing the matrix-vector multiplication ...  OK.                                       " << endl;
    }
-   file.close();
    return true;
 }
 
diff --git a/tests/long-time-unit-tests/tnl-run-matrix-formats-test b/tests/long-time-unit-tests/tnl-run-matrix-formats-test
index 563ee83be7..74c74700d3 100644
--- a/tests/long-time-unit-tests/tnl-run-matrix-formats-test
+++ b/tests/long-time-unit-tests/tnl-run-matrix-formats-test
@@ -18,7 +18,7 @@ test_matrix()
    # $1 - matrix name
    FILESIZE=$(stat -c%s "$1")
    echo "Testing $1 ( `echo \"$FILESIZE / 1024 / 1024 \" | bc`  MB)."
-   $TNL_MATRIX_TEST --input-file $1 --matrix-format csr --verbose $VERBOSE
+   $TNL_MATRIX_TEST --input-file $1 --matrix-format chunked-ellpack --hard-test yes --verbose $VERBOSE
    echo "-------------------------------------------------------------------------------------"   
 }
 
diff --git a/tests/unit-tests/matrices/tnlTridiagonalMatrixTester.h b/tests/unit-tests/matrices/tnlTridiagonalMatrixTester.h
index 13d3854fd9..ae232530df 100644
--- a/tests/unit-tests/matrices/tnlTridiagonalMatrixTester.h
+++ b/tests/unit-tests/matrices/tnlTridiagonalMatrixTester.h
@@ -62,7 +62,7 @@ class tnlTridiagonalMatrixTester : public CppUnit :: TestCase
    void setDimensionsTest()
    {
       MatrixType m;
-      m.setDimensions( 10 );
+      m.setDimensions( 10, 10 );
       CPPUNIT_ASSERT( m.getRows() == 10 );
       CPPUNIT_ASSERT( m.getColumns() == 10 );
    }
@@ -70,7 +70,7 @@ class tnlTridiagonalMatrixTester : public CppUnit :: TestCase
    void setLikeTest()
    {
       MatrixType m1, m2;
-      m1.setDimensions( 10 );
+      m1.setDimensions( 10, 10 );
       m2.setLike( m1 );
       CPPUNIT_ASSERT( m1.getRows() == m2.getRows() );
    }
@@ -79,7 +79,7 @@ class tnlTridiagonalMatrixTester : public CppUnit :: TestCase
    {
       const int size( 10 );
       MatrixType m;
-      m.setDimensions( size );
+      m.setDimensions( size, size );
       m.setValue( 1.0 );
       for( int i = 0; i < size; i++ )
          for( int j = 0; j < size; j++ )
@@ -92,7 +92,7 @@ class tnlTridiagonalMatrixTester : public CppUnit :: TestCase
    void setElementTest()
    {
       MatrixType m;
-      m.setDimensions( 10 );
+      m.setDimensions( 10, 10 );
       for( int i = 0; i < 10; i++ )
          m.setElement( i, i, i );
       for( int i = 0; i < 10; i++ )
@@ -107,7 +107,7 @@ class tnlTridiagonalMatrixTester : public CppUnit :: TestCase
    void addElementTest()
    {
       MatrixType m;
-      m.setDimensions( 10 );
+      m.setDimensions( 10, 10 );
       for( int i = 0; i < 10; i++ )
          m.setElement( i, i, i );
       for( int i = 0; i < 10; i++ )
@@ -129,7 +129,7 @@ class tnlTridiagonalMatrixTester : public CppUnit :: TestCase
    void setRowTest()
    {
       MatrixType m;
-      m.setDimensions( 10 );
+      m.setDimensions( 10, 10 );
       m.setValue( 0.0 );
       for( int i = 0; i < 10; i++ )
          m.setElement( i, i, i );
@@ -165,7 +165,7 @@ class tnlTridiagonalMatrixTester : public CppUnit :: TestCase
       v.setSize( size );
       w.setSize( size );
       MatrixType m;
-      m.setDimensions( size );
+      m.setDimensions( size, size );
       for( int i = 0; i < size; i++ )
       {
          v.setElement( i, i );
@@ -181,7 +181,7 @@ class tnlTridiagonalMatrixTester : public CppUnit :: TestCase
    {
       const int size = 10;
       MatrixType m;
-      m.setDimensions( 10);
+      m.setDimensions( 10, 10 );
       for( int i = 0; i < size; i++ )
          for( int j = 0; j < size; j++ )
             if( abs( i - j ) <= 1 )
@@ -213,7 +213,7 @@ class tnlTridiagonalMatrixTester : public CppUnit :: TestCase
    {
       const int size = 10;
       MatrixType m;
-      m.setDimensions( 10 );
+      m.setDimensions( 10, 10 );
       for( int i = 0; i < size; i++ )
          for( int j = 0; j < size; j++ )
             if( abs( i - j ) <= 1 )
-- 
GitLab