From 73fd90eabf898b586a07b8832c574147f1cd420f Mon Sep 17 00:00:00 2001
From: Tomas Oberhuber <tomas.oberhuber@fjfi.cvut.cz>
Date: Fri, 20 May 2011 18:30:53 +0000
Subject: [PATCH] Debuging tnlAdaptiveRgCSRMatrix.

---
 src/matrix/tnlAdaptiveRgCSRMatrix.h           |  87 +++--
 tests/run-sparse-matrix-benchmark             |   1 +
 tests/sparse-matrix-benchmark.h               |   2 +
 tests/tnlSpmvBenchmark.h                      |   1 +
 tests/tnlSpmvBenchmarkAdaptiveRgCSRMatrix.h   |   4 +-
 tests/tnlSpmvBenchmarkHybridMatrix.h          |   2 +-
 tests/unit-tests/Makefile.am                  |   2 +-
 tests/unit-tests/Makefile.in                  |   2 +-
 .../matrix/tnlAdaptiveRgCSRMatrixTester.h     | 322 +++++++++++++++---
 tests/unit-tests/tnl-unit-tests.cpp           |   6 +-
 10 files changed, 350 insertions(+), 79 deletions(-)

diff --git a/src/matrix/tnlAdaptiveRgCSRMatrix.h b/src/matrix/tnlAdaptiveRgCSRMatrix.h
index 624c97151e..a2af863172 100644
--- a/src/matrix/tnlAdaptiveRgCSRMatrix.h
+++ b/src/matrix/tnlAdaptiveRgCSRMatrix.h
@@ -276,10 +276,10 @@ void tnlAdaptiveRgCSRMatrix< Real, Device, Index > :: tuneFormat( const Index cu
 }
 
 template< typename Real, tnlDevice Device, typename Index >
-bool tnlAdaptiveRgCSRMatrix< Real, Device, Index > :: copyFrom( const tnlCSRMatrix< Real, tnlHost, Index >& mat )
+bool tnlAdaptiveRgCSRMatrix< Real, Device, Index > :: copyFrom( const tnlCSRMatrix< Real, tnlHost, Index >& csrMatrix )
 {
 	dbgFunctionName( "tnlAdaptiveRgCSRMatrix< Real, tnlHost >", "copyFrom" );
-	if( ! this -> setSize( mat.getSize() ) )
+	if( ! this -> setSize( csrMatrix.getSize() ) )
 		return false;
 	
 	Index nonzerosInGroup( 0 );
@@ -304,7 +304,7 @@ bool tnlAdaptiveRgCSRMatrix< Real, Device, Index > :: copyFrom( const tnlCSRMatr
 		if( groupEnd > this -> getSize() )
 		   groupEnd = this -> getSize();
 
-		nonzerosInGroup = mat. row_offsets[ groupEnd ] - mat. row_offsets[ groupBegin ];
+		nonzerosInGroup = csrMatrix. row_offsets[ groupEnd ] - csrMatrix. row_offsets[ groupBegin ];
 		rowsInGroup = groupEnd - groupBegin;
 
 		if( nonzerosInGroup < targetNonzeroesPerGroup &&
@@ -326,7 +326,7 @@ bool tnlAdaptiveRgCSRMatrix< Real, Device, Index > :: copyFrom( const tnlCSRMatr
 		Index usedThreads = 0;
 		for( Index i = groupBegin; i < groupEnd; i++ )
 		{
-		   double nonzerosInRow = mat. getNonzeroElementsInRow( i );
+		   double nonzerosInRow = csrMatrix. getNonzeroElementsInRow( i );
 		   double nonzerosInRowRatio = nonzerosInRow / ( double ) nonzerosInGroup;
 		   usedThreads += threadsPerRow[ i - groupBegin ] = floor( freeThreads * nonzerosInRowRatio );
 		}
@@ -357,7 +357,7 @@ bool tnlAdaptiveRgCSRMatrix< Real, Device, Index > :: copyFrom( const tnlCSRMatr
 		Index rounds( 0 ), roundsFinal( 0 );
 		for( Index i = groupBegin; i < groupEnd; i++ )
 		{
-		   double nonzerosInRow = mat. getNonzeroElementsInRow( i );
+		   double nonzerosInRow = csrMatrix. getNonzeroElementsInRow( i );
 		   rounds = ceil( nonzerosInRow / ( double ) threadsPerRow[ i - groupBegin ] );
 		   roundsFinal = Max( rounds, roundsFinal );
 		}
@@ -392,7 +392,7 @@ bool tnlAdaptiveRgCSRMatrix< Real, Device, Index > :: copyFrom( const tnlCSRMatr
 	dbgCout( "Allocating " << numberOfStoredValues << " elements.");
 	if( ! setNonzeroElements( numberOfStoredValues ) )
 		return false;
-	artificial_zeros = numberOfStoredValues - mat. getNonzeroElements();
+	artificial_zeros = numberOfStoredValues - csrMatrix. getNonzeroElements();
 
 	last_nonzero_element = numberOfStoredValues;
 
@@ -413,7 +413,7 @@ bool tnlAdaptiveRgCSRMatrix< Real, Device, Index > :: copyFrom( const tnlCSRMatr
 			 */
 			for( Index j = 0; j < groupInfo[ i ]. numRows; j++ )
 			{
-			   NZperRow[ j ] = mat. getNonzeroElementsInRow( baseRow + j );
+			   NZperRow[ j ] = csrMatrix. getNonzeroElementsInRow( baseRow + j );
 			   if( j > 0 )
 			      threadsPerRow[ j ] = threads[ baseRow + j ] - threads[ baseRow + j - 1 ];
 			   else
@@ -429,10 +429,10 @@ bool tnlAdaptiveRgCSRMatrix< Real, Device, Index > :: copyFrom( const tnlCSRMatr
 					{
 						if( counters[ j ] < NZperRow[ j ] )
 						{
-						   Index pos = mat. row_offsets[ baseRow + j ] + counters[ j ];
+						   Index pos = csrMatrix. row_offsets[ baseRow + j ] + counters[ j ];
 						   dbgCout( "Inserting data from CSR format at position " << pos << " to AdaptiveRgCSR at " << index );
-							nonzeroElements[ index ] = mat. nonzero_elements[ pos ];
-							columns[ index ] = mat.columns[ pos ];
+							nonzeroElements[ index ] = csrMatrix. nonzero_elements[ pos ];
+							columns[ index ] = csrMatrix. columns[ pos ];
 						}
 						else
 						{
@@ -442,6 +442,7 @@ bool tnlAdaptiveRgCSRMatrix< Real, Device, Index > :: copyFrom( const tnlCSRMatr
 						}
 						counters[ j ] ++;
 						index ++;
+						tnlAssert( index < numberOfStoredValues, );
 					}
 		}
 	}
@@ -518,7 +519,7 @@ Real tnlAdaptiveRgCSRMatrix< Real, Device, Index > :: getElement( Index row,
                   if( columns[ pointer ] == column )
                      return nonzeroElements[ pointer ];
                   pointer ++;
-                  }
+               }
          }
       return 0.0;
    }
@@ -547,54 +548,76 @@ void tnlAdaptiveRgCSRMatrix< Real, Device, Index > :: vectorProduct( const tnlLo
 
 
    const Index TB_SIZE = 256;
-   const Index MAX_ROWS = 128;
+   const Index MAX_ROWS = 256;
    if( Device == tnlHost )
    {
-      Index idx[TB_SIZE];
-      Real psum[TB_SIZE];        //partial sums for each thread
-      Index limits[MAX_ROWS + 1];  //indices of first threads for each row + index of first unused thread
-      Real results[MAX_ROWS];
+      Index idx[ TB_SIZE ];
+      Real psum[ TB_SIZE ];        //partial sums for each thread
+      Index limits[ MAX_ROWS + 1 ];  //indices of first threads for each row + index of first unused thread
+      Real results[ MAX_ROWS ];
 
       /****
        * Go over all groups ...
        */
-      for( Index group = 0; group < this -> numberOfGroups; group ++ )
+      dbgExpr( this -> numberOfGroups );
+      for( Index groupId = 0; groupId < this -> numberOfGroups; groupId ++ )
       {
          /****
           * In each group compute partial sums of each thread
           */
+         dbgExpr( groupId );
+         dbgExpr( this -> usedThreadsInGroup[ groupId ] );
          for( Index thread = 0;
-              thread < this -> usedThreadsInGroup[ group ];
+              thread < this -> usedThreadsInGroup[ groupId ];
               thread ++ )
          {
-            idx[ thread ] = this -> groupInfo[ group ]. idxFirstValue + thread;
-            psum[thread] = 0;
-
+            idx[ thread ] = this -> groupInfo[ groupId ]. idxFirstValue + thread;
+            psum[ thread ] = 0;
             for( Index j = 0;
-                 j < this -> groupInfo[ group ]. numRounds;
+                 j < this -> groupInfo[ groupId ]. numRounds;
                  j ++ )
             {
-               psum[ thread ] += this -> nonzeroElements[ idx[ thread ] ] * vec[ this -> columns[ idx[ thread ] ] ];
-               idx[ j ] += TB_SIZE;
+               if( this -> columns[ idx[ thread ] ] != -1  )
+                  psum[ thread ] += this -> nonzeroElements[ idx[ thread ] ] * vec[ this -> columns[ idx[ thread ] ] ];
+               idx[ thread ] += TB_SIZE;
             }
+            dbgExpr( psum[ thread ] );
          }
 
+         for( Index row = groupInfo[ groupId ]. idxFirstRow;
+              row < groupInfo[ groupId ]. idxFirstRow + groupInfo[ groupId ]. numRows;
+              row ++ )
+         {
+            Index firstThreadInRow = 0;
+            if( row > 0 )
+               firstThreadInRow = threads[ row - 1 ];
+            Index lastThreadInRow = threads[ row ];
+
+            dbgCout( "Row: " << row << " firstThreadInRow: " << firstThreadInRow << " lastThreadInRow: " << lastThreadInRow );
+            result[ row ] = 0.0;
+            for( Index thread = firstThreadInRow; thread < lastThreadInRow; thread ++ )
+            {
+               result[ row ] += psum[ thread ];
+               dbgCout( "Thread: " << thread << " psum[ thread ]: " << psum[ thread ] << " result[ row ]: " << result[ row ] );
+            }
+         }
+#ifdef UNDEF
          /****
           * Compute local copy of thread indexes mapped to given row of the group.
           * (this is only to simulate copying data to the fast shared memory on GPU)
           */
          for( Index thread = 0;
-              thread < this -> groupInfo[ group ]. numRows;
+              thread < this -> groupInfo[ groupId ]. numRows;
               thread ++ )
-            limits[ thread ] = this -> threads[ this -> groupInfo[ group ]. idxFirstRow + thread ];
+            limits[ thread ] = this -> threads[ this -> groupInfo[ groupId ]. idxFirstRow + thread ];
          /****
           * For convenience, add the index of first unused row.
           */
-         limits[ this -> groupInfo[ group ]. numRows ] = this -> usedThreadsInGroup[ group ];
+         limits[ this -> groupInfo[ groupId ]. numRows ] = this -> usedThreadsInGroup[ groupId ];
 
          //reduction of partial sums and writing to the output
          for( Index thread = 0;
-              thread < this -> groupInfo[ group ]. numRows;
+              thread < this -> groupInfo[ groupId ]. numRows;
               thread ++)         //for threads corresponding to rows in group
          {
             results[ thread ] = 0;
@@ -602,8 +625,9 @@ void tnlAdaptiveRgCSRMatrix< Real, Device, Index > :: vectorProduct( const tnlLo
                  j < limits[ thread + 1 ];
                  j++ )              //sum up partial sums belonging to that row
                results[ thread ] += psum[ j ];
-            result[ this -> groupInfo[ group ]. idxFirstRow + thread ] = results[ thread ];
+            result[ this -> groupInfo[ groupId ]. idxFirstRow + thread ] = results[ thread ];
          }
+#endif
       }
    }
    if( Device == tnlCuda )
@@ -640,6 +664,11 @@ template< typename Real, tnlDevice Device, typename Index >
 void tnlAdaptiveRgCSRMatrix< Real, Device, Index > :: printOut( ostream& str,
 		                                                          const Index lines ) const
 {
+   /*****
+    * THIS IS NOT CORRECT
+    */
+   cerr << "tnlAdaptiveRgCSRMatrix< Real, Device, Index > :: printOut is not correct" << endl;
+   abort();
    str << "Structure of tnlAdaptiveRgCSRMatrix" << endl;
    str << "Matrix name:" << this -> getName() << endl;
    str << "Matrix size:" << this -> getSize() << endl;
diff --git a/tests/run-sparse-matrix-benchmark b/tests/run-sparse-matrix-benchmark
index c6feaa096a..be8869d2ca 100755
--- a/tests/run-sparse-matrix-benchmark
+++ b/tests/run-sparse-matrix-benchmark
@@ -7,6 +7,7 @@ IWD="$PWD"
 DEBUG="no"
 SPARSE_MATRIX_BENCHMARK="$IWD/sparse-matrix-benchmark"
 SPARSE_MATRIX_BENCHMARK_DBG="$IWD/.libs/sparse-matrix-benchmark-dbg"
+#SPARSE_MATRIX_BENCHMARK="$IWD/.libs/sparse-matrix-benchmark-dbg"
 STOP_TIME="1"
 
 source ../tnl-env-variables
diff --git a/tests/sparse-matrix-benchmark.h b/tests/sparse-matrix-benchmark.h
index b2a448c692..9bc61a6d56 100644
--- a/tests/sparse-matrix-benchmark.h
+++ b/tests/sparse-matrix-benchmark.h
@@ -221,6 +221,7 @@ bool benchmarkMatrix( const tnlString& input_file,
       }
    }
 
+#ifdef UNDEF
    /****
     * Row-Grouped CSR format
     */
@@ -371,6 +372,7 @@ bool benchmarkMatrix( const tnlString& input_file,
 
 
    }
+#endif
    csrMatrix. vectorProduct( refX, refB );
 
    /****
diff --git a/tests/tnlSpmvBenchmark.h b/tests/tnlSpmvBenchmark.h
index 86c891ec58..b412eff1f0 100644
--- a/tests/tnlSpmvBenchmark.h
+++ b/tests/tnlSpmvBenchmark.h
@@ -238,6 +238,7 @@ void tnlSpmvBenchmark< Real, Device, Index, Matrix > :: runBenchmark( const tnlL
          benchmarkWasSuccesful = false;
          firstErrorOccurence = j;
          //cerr << " b[ " << j << " ] = " << resB[ j ] << " while refB[ " << j << " ] = " << refB[ j ] << endl;
+         //abort();
       }
    }
 
diff --git a/tests/tnlSpmvBenchmarkAdaptiveRgCSRMatrix.h b/tests/tnlSpmvBenchmarkAdaptiveRgCSRMatrix.h
index a27b873188..d1b4773ebe 100644
--- a/tests/tnlSpmvBenchmarkAdaptiveRgCSRMatrix.h
+++ b/tests/tnlSpmvBenchmarkAdaptiveRgCSRMatrix.h
@@ -68,7 +68,7 @@ template< typename Real,
           typename Index>
 bool tnlSpmvBenchmarkAdaptiveRgCSRMatrix< Real, Device, Index > :: setup( const tnlCSRMatrix< Real, tnlHost, Index >& matrix )
 {
-   tnlAssert( this -> groupSize > 0, cerr << "groupSize = " << this -> groupSize );
+   //tnlAssert( this -> groupSize > 0, cerr << "groupSize = " << this -> groupSize );
    if( Device == tnlHost )
    {
       if( ! this -> matrix. copyFrom( matrix ) )
@@ -102,7 +102,7 @@ template< typename Real,
           typename Index >
 void tnlSpmvBenchmarkAdaptiveRgCSRMatrix< Real, Device, Index > :: writeProgress() const
 {
-   cout << left << setw( this -> formatColumnWidth - 15 ) << "Row-grouped CSR ";
+   cout << left << setw( this -> formatColumnWidth - 15 ) << "Adaptive Row-grouped CSR ";
    if( Device == tnlCuda )
       cout << setw( 5 ) << this -> groupSize
            << setw( 10 ) << this -> cudaBlockSize;
diff --git a/tests/tnlSpmvBenchmarkHybridMatrix.h b/tests/tnlSpmvBenchmarkHybridMatrix.h
index cdda4122d4..7505b04ec1 100644
--- a/tests/tnlSpmvBenchmarkHybridMatrix.h
+++ b/tests/tnlSpmvBenchmarkHybridMatrix.h
@@ -59,7 +59,7 @@ void tnlSpmvBenchmarkHybridMatrix< Real, Index > :: setFileName( const tnlString
 template< typename Real, typename Index>
 bool tnlSpmvBenchmarkHybridMatrix< Real, Index > :: setup( const tnlCSRMatrix< Real, tnlHost, Index >& matrix )
 {
-
+   return true;
 }
 
 template< typename Real,
diff --git a/tests/unit-tests/Makefile.am b/tests/unit-tests/Makefile.am
index b35b63fc99..b6db4a4361 100644
--- a/tests/unit-tests/Makefile.am
+++ b/tests/unit-tests/Makefile.am
@@ -67,7 +67,7 @@ endif
 #                                 ../../src/debug/libtnldebug-dbg-0.1.la                            
 #endif
 
-TESTS = tnl-unit-tests
+TESTS = tnl-unit-tests-dbg
 
 #if BUILD_CUDA
 #TESTS += tnl-cuda-unit-tests
diff --git a/tests/unit-tests/Makefile.in b/tests/unit-tests/Makefile.in
index d3a66eed46..c4de5fbc90 100644
--- a/tests/unit-tests/Makefile.in
+++ b/tests/unit-tests/Makefile.in
@@ -42,7 +42,7 @@ check_PROGRAMS = tnl-unit-tests$(EXEEXT) $(am__EXEEXT_1) \
 	$(am__EXEEXT_2)
 @BUILD_CUDA_TRUE@am__append_1 = tnl-cuda-unit-tests
 @BUILD_DBG_TRUE@am__append_2 = tnl-unit-tests-dbg
-TESTS = tnl-unit-tests$(EXEEXT)
+TESTS = tnl-unit-tests-dbg$(EXEEXT)
 subdir = tests/unit-tests
 DIST_COMMON = $(srcdir)/Makefile.am $(srcdir)/Makefile.in
 ACLOCAL_M4 = $(top_srcdir)/aclocal.m4
diff --git a/tests/unit-tests/matrix/tnlAdaptiveRgCSRMatrixTester.h b/tests/unit-tests/matrix/tnlAdaptiveRgCSRMatrixTester.h
index ae7b79e8c9..f92e28ea3e 100644
--- a/tests/unit-tests/matrix/tnlAdaptiveRgCSRMatrixTester.h
+++ b/tests/unit-tests/matrix/tnlAdaptiveRgCSRMatrixTester.h
@@ -38,59 +38,219 @@ template< class T > class tnlAdaptiveRgCSRMatrixTester : public CppUnit :: TestC
       CppUnit :: TestSuite* suiteOfTests = new CppUnit :: TestSuite( "tnlAdaptiveRgCSRMatrixTester" );
       CppUnit :: TestResult result;
       suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlAdaptiveRgCSRMatrixTester< T > >(
-                               "diagonalMatrixTest",
-                               & tnlAdaptiveRgCSRMatrixTester< T > :: diagonalMatrixTest )
+                               "ifDiagonalMatrixIsStoredProperly",
+                               & tnlAdaptiveRgCSRMatrixTester< T > :: ifDiagonalMatrixIsStoredProperly )
                              );
       suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlAdaptiveRgCSRMatrixTester< T > >(
-                               "triDiagonalMatrixTest",
-                               & tnlAdaptiveRgCSRMatrixTester< T > :: triDiagonalMatrixTest )
+                               "ifTriDiagonalMatrixIsStoredProperly",
+                               & tnlAdaptiveRgCSRMatrixTester< T > :: ifTriDiagonalMatrixIsStoredProperly )
                              );
       suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlAdaptiveRgCSRMatrixTester< T > >(
-                               "upperTriangularMatrixTest",
-                               & tnlAdaptiveRgCSRMatrixTester< T > :: upperTriangularMatrixTest )
+                               "ifUpperTriangularMatrixIsStoredProperly",
+                               & tnlAdaptiveRgCSRMatrixTester< T > :: ifUpperTriangularMatrixIsStoredProperly )
                              );
+      suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlAdaptiveRgCSRMatrixTester< T > >(
+                               "ifSpmvWithDiagonalMatrixWorks",
+                               & tnlAdaptiveRgCSRMatrixTester< T > :: ifSpmvWithDiagonalMatrixWorks )
+                             );
+      suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlAdaptiveRgCSRMatrixTester< T > >(
+                               "ifSpmvWithTriDiagonalMatrixWorks",
+                               & tnlAdaptiveRgCSRMatrixTester< T > :: ifSpmvWithTriDiagonalMatrixWorks )
+                             );
+
+      suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlAdaptiveRgCSRMatrixTester< T > >(
+                               "ifSpmvWithUpperTriangularMatrixWorks",
+                               & tnlAdaptiveRgCSRMatrixTester< T > :: ifSpmvWithUpperTriangularMatrixWorks  )
+                             );
+      suiteOfTests -> addTest( new CppUnit :: TestCaller< tnlAdaptiveRgCSRMatrixTester< T > >(
+                               "ifSpmvWithBcsstk20MatrixWorks",
+                               & tnlAdaptiveRgCSRMatrixTester< T > :: ifSpmvWithBcsstk20MatrixWorks )
+                             );
+
 
 
-            return suiteOfTests;
+      return suiteOfTests;
    }
 
-   void diagonalMatrixTest()
+   void setDiagonalMatrix( tnlCSRMatrix< T >& csrMatrix,
+                           tnlAdaptiveRgCSRMatrix< T >& argcsrMatrix,
+                           const int size )
    {
-      tnlCSRMatrix< T > csr_matrix( "test-matrix:Diagonal" );
-      tnlAdaptiveRgCSRMatrix< T > argcsrMatrix( "test-matrix:Diagonal" );
-      csr_matrix. setSize( 12 );
-      csr_matrix. setNonzeroElements( 12 );
-      for( int i = 0; i < 12; i ++ )
-         csr_matrix. setElement( i, i, T( i + 1 ) );
+      csrMatrix. setSize( size );
+      csrMatrix. setNonzeroElements( size );
+      for( int i = 0; i < size; i ++ )
+         csrMatrix. setElement( i, i, T( i + 1 ) );
       //cerr << "Copying data to coalesced CSR matrix." << endl;
-      argcsrMatrix. copyFrom( csr_matrix );
+      argcsrMatrix. copyFrom( csrMatrix );
       //argcsrMatrix. printOut( cout );
+   }
+
+   void setTridiagonalMatrix( tnlCSRMatrix< T >& csrMatrix,
+                              tnlAdaptiveRgCSRMatrix< T >& argcsrMatrix,
+                              const int size )
+   {
+      csrMatrix. setSize( size );
+      csrMatrix. setNonzeroElements( size * 3 - 2 );
+      T data[] = { -1.0, 2.0, -1.0 };
+      int offsets[] = { -1, 0, 1 };
+      for( int i = 0; i < size; i ++ )
+      {
+         csrMatrix. insertRow( i,      // row
+                               3,      // elements
+                               data,   // data
+                               i,      // first column
+                               offsets );
+      }
+      argcsrMatrix. copyFrom( csrMatrix );
+      //argcsr_matrix. printOut( cout );
+   }
+
+   void setUpperTriangularMatrix( tnlCSRMatrix< T >& csrMatrix,
+                                  tnlAdaptiveRgCSRMatrix< T >& argcsrMatrix,
+                                  const int size )
+   {
+      csrMatrix. setSize( size );
+      csrMatrix. setNonzeroElements( size * size );
+      for( int i = 0; i < size; i ++ )
+         for( int j = i; j < size; j ++ )
+            csrMatrix. setElement( i, j, 1.0 );
+
+      argcsrMatrix. copyFrom( csrMatrix );
+      //argcsr_matrix. printOut( cout );
+   }
+
+   void setBcsstk20Matrix( tnlCSRMatrix< T >& csrMatrix,
+                           tnlAdaptiveRgCSRMatrix< T >& argcsrMatrix )
+   {
+      /****
+       * Sets a real matrix from a matrix market
+       */
+      csrMatrix. setSize( 50 );
+      csrMatrix. setNonzeroElements( 100 );
+      
+      csrMatrix. setElement( 1, 1,  3.4651842369950e+08    );
+      csrMatrix. setElement( 3, 1,  3.0695529658430e+10    );
+      csrMatrix. setElement( 4, 1, -3.4651842369950e+08    );
+      csrMatrix. setElement( 6, 1,  3.0695529658430e+10    );
+      csrMatrix. setElement( 2, 2,  5.4150467497750e+08    );
+      csrMatrix. setElement( 5, 2, -5.4150467497750e+08    );
+      csrMatrix. setElement( 3, 3,  3.6254562588700e+12    );
+      csrMatrix. setElement( 4, 3, -3.0695529658430e+10    );
+      csrMatrix. setElement( 6, 3,  1.8127281294350e+12    );
+      csrMatrix. setElement( 4, 4,  5.8200610488580e+09    );
+      csrMatrix. setElement( 6, 4,  1.6324889406770e+11    );
+      csrMatrix. setElement( 7, 4, -5.4735426251580e+09    );
+      csrMatrix. setElement( 9, 4,  1.9394442372610e+11    );
+      csrMatrix. setElement( 5, 5,  1.7185252532270e+09    );
+      csrMatrix. setElement( 8, 5, -1.1770205782500e+09    );
+      csrMatrix. setElement( 6, 6,  1.2788184938840e+13    );
+      csrMatrix. setElement( 7, 6, -1.9394442372610e+11    );
+      csrMatrix. setElement( 9, 6,  4.5813643399870e+12    );
+      csrMatrix. setElement( 7, 7,  1.1034132683920e+10    );
+      csrMatrix. setElement( 9, 7,  3.0843578836850e+09    );
+      csrMatrix. setElement( 10, 7, -5.5605900587650e+09   );
+      csrMatrix. setElement( 12, 7,  1.9702878160980e+11   );
+      csrMatrix. setElement( 8, 8,  2.3681428407440e+09    );
+      csrMatrix. setElement( 11, 8, -1.1911222624940e+09   );
+      csrMatrix. setElement( 9, 9,  1.8471175055240e+13    );
+      csrMatrix. setElement( 10, 9, -1.9702878160980e+11   );
+      csrMatrix. setElement( 12, 9,  4.6542231876330e+12   );
+      csrMatrix. setElement( 10, 10,  1.1211709448480e+10  );
+      csrMatrix. setElement( 12, 10,  3.2077321990330e+09  );
+      csrMatrix. setElement( 13, 10, -5.6511193897150e+09  );
+      csrMatrix. setElement( 15, 10,  2.0023651380880e+11  );
+      csrMatrix. setElement( 11, 11,  2.3954060969490e+09  );
+      csrMatrix. setElement( 14, 11, -1.2042838344550e+09  );
+      csrMatrix. setElement( 12, 12,  1.8768439153640e+13  );
+      csrMatrix. setElement( 13, 12, -2.0023651380880e+11  );
+      csrMatrix. setElement( 15, 12,  4.7299963891850e+12  );
+      csrMatrix. setElement( 13, 13,  1.1392768110380e+10  );
+      csrMatrix. setElement( 15, 13,  3.2077321990300e+09  );
+      csrMatrix. setElement( 16, 13, -5.7416487206660e+09  );
+      csrMatrix. setElement( 18, 13,  2.0344424600780e+11  );
+      csrMatrix. setElement( 14, 14,  2.4226693531550e+09  );
+      csrMatrix. setElement( 17, 14, -1.2183855187000e+09  );
+      csrMatrix. setElement( 15, 15,  1.9071531959840e+13  );
+      csrMatrix. setElement( 16, 15, -2.0344424600780e+11  );
+      csrMatrix. setElement( 18, 15,  4.8057695907370e+12  );
+      csrMatrix. setElement( 16, 16,  1.1570344874940e+10  );
+      csrMatrix. setElement( 18, 16,  3.0843578836840e+09  );
+      csrMatrix. setElement( 19, 16, -5.8286961542720e+09  );
+      csrMatrix. setElement( 21, 16,  2.0652860389150e+11  );
+      csrMatrix. setElement( 17, 17,  2.4499326093600e+09  );
+      csrMatrix. setElement( 20, 17, -1.2315470906610e+09  );
+      csrMatrix. setElement( 18, 18,  1.9368796058240e+13  );
+      csrMatrix. setElement( 19, 18, -2.0652860389150e+11  );
+      csrMatrix. setElement( 21, 18,  4.8786284383830e+12  );
+      csrMatrix. setElement( 19, 19,  1.1747921639490e+10  );
+      csrMatrix. setElement( 21, 19,  3.2077321990320e+09  );
+      csrMatrix. setElement( 22, 19, -5.9192254852220e+09  );
+      csrMatrix. setElement( 24, 19,  2.0973633609060e+11  );
+      csrMatrix. setElement( 20, 20,  2.4771958655660e+09  );
+      csrMatrix. setElement( 23, 20, -1.2456487749050e+09  );
+      csrMatrix. setElement( 21, 21,  1.9666060156640e+13  );
+      csrMatrix. setElement( 22, 21, -2.0973633609060e+11  );
+      csrMatrix. setElement( 24, 21,  4.9544016399350e+12  );
+      csrMatrix. setElement( 22, 22,  1.1928980301390e+10  );
+      csrMatrix. setElement( 24, 22,  3.2077321990300e+09  );
+      csrMatrix. setElement( 25, 22, -6.0097548161730e+09  );
+      csrMatrix. setElement( 27, 22,  2.1294406828960e+11  );
+      csrMatrix. setElement( 23, 23,  2.5044591217710e+09  );
+      csrMatrix. setElement( 26, 23, -1.2588103468660e+09  );
+      csrMatrix. setElement( 24, 24,  1.9969152962840e+13  );
+      csrMatrix. setElement( 25, 24, -2.1294406828960e+11  );
+      csrMatrix. setElement( 27, 24,  5.0301748414870e+12  );
+      csrMatrix. setElement( 25, 25,  1.2106557065950e+10  );
+      csrMatrix. setElement( 27, 25,  3.0843578836860e+09  );
+      csrMatrix. setElement( 28, 25, -6.0968022497790e+09  );
+      csrMatrix. setElement( 30, 25,  2.1602842617330e+11  );
+      csrMatrix. setElement( 26, 26,  2.5317223779770e+09  );
+      csrMatrix. setElement( 29, 26, -1.2729120311100e+09  );
+      csrMatrix. setElement( 27, 27,  2.0266417061240e+13  );
+      csrMatrix. setElement( 28, 27, -2.1602842617330e+11  );
+      csrMatrix. setElement( 30, 27,  5.1030336891330e+12  );
+      csrMatrix. setElement( 28, 28,  1.2284133830500e+10  );
+      csrMatrix. setElement( 30, 28,  3.2077321988540e+09  );
+      csrMatrix. setElement( 31, 28, -6.1873315807220e+09  );
+      csrMatrix. setElement( 33, 28,  2.1923615837210e+11  );
+      csrMatrix. setElement( 29, 29,  2.5589856341820e+09  );
+      csrMatrix. setElement( 32, 29, -1.2860736030710e+09  );
+      csrMatrix. setElement( 30, 30,  2.0563681159630e+13  );
+      csrMatrix. setElement( 31, 30, -2.1923615837210e+11  );
+      csrMatrix. setElement( 33, 30,  5.1788068906830e+12  );
+      csrMatrix. setElement( 31, 31,  1.2465192492400e+10  );
+      csrMatrix. setElement( 33, 31,  3.2077321992110e+09  );
+      csrMatrix. setElement( 34, 31, -6.2778609116800e+09  );
+      csrMatrix. setElement( 36, 31,  2.2244389057130e+11  );
+      csrMatrix. setElement( 32, 32,  2.5862488903870e+09  );
+      csrMatrix. setElement( 35, 32, -1.3001752873160e+09  );
+      csrMatrix. setElement( 33, 33,  2.0866773965840e+13  );
+
+
+      argcsrMatrix. copyFrom( csrMatrix );
+
+   }
+
+   void ifDiagonalMatrixIsStoredProperly()
+   {
+      const int size = 12;
+      tnlCSRMatrix< T > csrMatrix( "test-matrix:Diagonal" );
+      tnlAdaptiveRgCSRMatrix< T > argcsrMatrix( "test-matrix:Diagonal" );
+      setDiagonalMatrix( csrMatrix, argcsrMatrix, size );
+
       bool error( false );
-      for( int i = 0; i < 12; i ++ )
+      for( int i = 0; i < size; i ++ )
          if( argcsrMatrix. getElement( i, i ) != T( i + 1 ) )
             error = true;
       CPPUNIT_ASSERT( ! error );
    };
 
-   void triDiagonalMatrixTest()
+   void ifTriDiagonalMatrixIsStoredProperly()
    {
       tnlCSRMatrix< T > csr_matrix( "test-matrix:Tridiagonal" );
       tnlAdaptiveRgCSRMatrix< T > argcsr_matrix( "test-matrix:Tridiagonal" );
       int size = 12;
-      csr_matrix. setSize( size );
-      csr_matrix. setNonzeroElements( size * 3 - 2 );
-      T data[] = { -1.0, 2.0, -1.0 };
-      int offsets[] = { -1, 0, 1 };
-      for( int i = 0; i < size; i ++ )
-      {
-         csr_matrix. insertRow( i,      // row
-                                3,      // elements
-                                data,   // data
-                                i,      // first column
-                                offsets );
-      }
-      argcsr_matrix. copyFrom( csr_matrix );
-      //argcsr_matrix. printOut( cout );
+      setTridiagonalMatrix( csr_matrix, argcsr_matrix, size );
       bool error( false );
       for( int i = 0; i < size; i ++ )
          for( int j = 0; j < size; j ++ )
@@ -99,20 +259,12 @@ template< class T > class tnlAdaptiveRgCSRMatrixTester : public CppUnit :: TestC
       CPPUNIT_ASSERT( ! error );
 
    }
-   void upperTriangularMatrixTest()
+   void ifUpperTriangularMatrixIsStoredProperly()
    {
       tnlCSRMatrix< T > csr_matrix( "test-matrix:upperTriangular" );
       tnlAdaptiveRgCSRMatrix< T > argcsr_matrix( "test-matrix:upperTriangular" );
-      int size = 12;
-      csr_matrix. setSize( size );
-      csr_matrix. setNonzeroElements( size * size );
-      for( int i = 0; i < size; i ++ )
-         for( int j = i; j < size; j ++ )
-            csr_matrix. setElement( i, j, 1.0 );
-
-      argcsr_matrix. copyFrom( csr_matrix );
-      //cerr << "----------------" << endl;
-      //argcsr_matrix. printOut( cout );
+      const int size = 12;
+      setUpperTriangularMatrix( csr_matrix, argcsr_matrix, size );
       bool error( false );
       for( int i = 0; i < size; i ++ )
          for( int j = 0; j < size; j ++ )
@@ -120,9 +272,95 @@ template< class T > class tnlAdaptiveRgCSRMatrixTester : public CppUnit :: TestC
                error = true;
 
       CPPUNIT_ASSERT( ! error );
+   }
+
+   void ifSpmvWithDiagonalMatrixWorks()
+   {
+      const int size = 12;
+      tnlCSRMatrix< T > csrMatrix( "test-matrix:Diagonal" );
+      tnlAdaptiveRgCSRMatrix< T > argcsrMatrix( "test-matrix:Diagonal" );
+      setDiagonalMatrix( csrMatrix, argcsrMatrix, size );
+
+      tnlLongVector< T > x( "x" ), b1( "b1" ), b2( "b2" );
+      x. setSize( size );
+      b1. setSize( size );
+      b2. setSize( size );
+      x. setValue( 1.0 );
+      csrMatrix. vectorProduct( x, b1 );
+      argcsrMatrix. vectorProduct( x, b2 );
+
+      cout << b1 << endl;
+      cout << b2 << endl;
 
+      CPPUNIT_ASSERT( b1 == b2 );
    }
 
+   void ifSpmvWithTriDiagonalMatrixWorks()
+   {
+      const int size = 12;
+      tnlCSRMatrix< T > csrMatrix( "test-matrix:TriDiagonal" );
+      tnlAdaptiveRgCSRMatrix< T > argcsrMatrix( "test-matrix:TriDiagonal" );
+      setTridiagonalMatrix( csrMatrix, argcsrMatrix, size );
+
+      tnlLongVector< T > x( "x" ), b1( "b1" ), b2( "b2" );
+      x. setSize( size );
+      b1. setSize( size );
+      b2. setSize( size );
+      x. setValue( 1.0 );
+      csrMatrix. vectorProduct( x, b1 );
+      argcsrMatrix. vectorProduct( x, b2 );
+
+      cout << b1 << endl;
+      cout << b2 << endl;
+
+      CPPUNIT_ASSERT( b1 == b2 );
+   }
+
+   void ifSpmvWithUpperTriangularMatrixWorks()
+   {
+      const int size = 12;
+      tnlCSRMatrix< T > csrMatrix( "test-matrix:TriDiagonal" );
+      tnlAdaptiveRgCSRMatrix< T > argcsrMatrix( "test-matrix:TriDiagonal" );
+      setUpperTriangularMatrix( csrMatrix, argcsrMatrix, size );
+
+      tnlLongVector< T > x( "x" ), b1( "b1" ), b2( "b2" );
+      x. setSize( size );
+      b1. setSize( size );
+      b2. setSize( size );
+      x. setValue( 1.0 );
+      csrMatrix. vectorProduct( x, b1 );
+      argcsrMatrix. vectorProduct( x, b2 );
+
+      cout << b1 << endl;
+      cout << b2 << endl;
+
+      CPPUNIT_ASSERT( b1 == b2 );
+   }
+
+   void ifSpmvWithBcsstk20MatrixWorks()
+   {
+      tnlCSRMatrix< T > csrMatrix( "test-matrix:TriDiagonal" );
+      tnlAdaptiveRgCSRMatrix< T > argcsrMatrix( "test-matrix:TriDiagonal" );
+      setBcsstk20Matrix( csrMatrix, argcsrMatrix );
+      const int size = csrMatrix. getSize();
+
+      tnlLongVector< T > x( "x" ), b1( "b1" ), b2( "b2" );
+      x. setSize( size );
+      b1. setSize( size );
+      b2. setSize( size );
+      x. setValue( 1.0 );
+      csrMatrix. vectorProduct( x, b1 );
+      argcsrMatrix. vectorProduct( x, b2 );
+
+      cout << b1 << endl;
+      cout << b2 << endl;
+
+      CPPUNIT_ASSERT( b1 == b2 );
+   }
+
+
+
+
 };
 
 #endif /* TNLAdaptiveRgCSRMATRIXTESTER_H_ */
diff --git a/tests/unit-tests/tnl-unit-tests.cpp b/tests/unit-tests/tnl-unit-tests.cpp
index e1e84498f2..77af728501 100644
--- a/tests/unit-tests/tnl-unit-tests.cpp
+++ b/tests/unit-tests/tnl-unit-tests.cpp
@@ -76,12 +76,12 @@ int main( int argc, char* argv[] )
    runner. addTest( tnlRgCSRMatrixTester< float > :: suite() );
    runner. addTest( tnlRgCSRMatrixTester< double > :: suite() );
 
-   runner. addTest( tnlAdaptiveRgCSRMatrixTester< float > :: suite() );
-   runner. addTest( tnlAdaptiveRgCSRMatrixTester< double > :: suite() );
-
    runner. addTest( tnlEllpackMatrixTester< float > :: suite() );
    runner. addTest( tnlEllpackMatrixTester< double > :: suite() );
 
+   runner. addTest( tnlAdaptiveRgCSRMatrixTester< float > :: suite() );
+   runner. addTest( tnlAdaptiveRgCSRMatrixTester< double > :: suite() );
+
    runner. addTest( tnlMPIMeshTester< float > :: suite() );
 
    //runner. addTest( tnlSharedMemoryTester< tnlHost > :: suite() );
-- 
GitLab