diff --git a/src/Python/pytnl/tnl/SparseMatrix.h b/src/Python/pytnl/tnl/SparseMatrix.h
index 5d278ae8a17fc204f51b39652cc5dfe735c28dca..7f8aa040478eec2b3628d7fab6262c415d083ec0 100644
--- a/src/Python/pytnl/tnl/SparseMatrix.h
+++ b/src/Python/pytnl/tnl/SparseMatrix.h
@@ -99,7 +99,7 @@ void export_Matrix( py::module & m, const char* name )
         // TODO: these two don't work
         //.def("addMatrix",           &Matrix::addMatrix)
         //.def("getTransposition",    &Matrix::getTransposition)
-        .def("performSORIteration", &Matrix::template performSORIteration< VectorType >)
+        .def("performSORIteration", &Matrix::template performSORIteration< VectorType, VectorType >)
 //        .def("assign",              &Matrix::operator=)
         .def("assign", []( Matrix& matrix, const Matrix& other ) -> Matrix& {
                 return matrix = other;
diff --git a/src/TNL/Matrices/CSR.h b/src/TNL/Matrices/CSR.h
index c6d7dd5b7f48806635ddb9b49c9b3412e92f51cc..f6d4c6d31b895476e16bc2ce2616a9bc77fce505 100644
--- a/src/TNL/Matrices/CSR.h
+++ b/src/TNL/Matrices/CSR.h
@@ -164,10 +164,10 @@ public:
    void getTransposition( const CSR< Real2, Device, Index2 >& matrix,
                           const RealType& matrixMultiplicator = 1.0 );
 
-   template< typename Vector >
-   bool performSORIteration( const Vector& b,
+   template< typename Vector1, typename Vector2 >
+   bool performSORIteration( const Vector1& b,
                              const IndexType row,
-                             Vector& x,
+                             Vector2& x,
                              const RealType& omega = 1.0 ) const;
 
    // copy assignment
diff --git a/src/TNL/Matrices/CSR_impl.h b/src/TNL/Matrices/CSR_impl.h
index 9e22c5a9d6abd27be5523e1428aa318dcd802f42..1516e932231c900c6e56b5442b88133cd5267a1a 100644
--- a/src/TNL/Matrices/CSR_impl.h
+++ b/src/TNL/Matrices/CSR_impl.h
@@ -494,10 +494,10 @@ void CSR< Real, Device, Index >::getTransposition( const CSR< Real2, Device, Ind
 template< typename Real,
           typename Device,
           typename Index >
-   template< typename Vector >
-bool CSR< Real, Device, Index >::performSORIteration( const Vector& b,
+   template< typename Vector1, typename Vector2 >
+bool CSR< Real, Device, Index >::performSORIteration( const Vector1& b,
                                                       const IndexType row,
-                                                      Vector& x,
+                                                      Vector2& x,
                                                       const RealType& omega ) const
 {
    TNL_ASSERT( row >=0 && row < this->getRows(),
diff --git a/src/TNL/Matrices/ChunkedEllpack.h b/src/TNL/Matrices/ChunkedEllpack.h
index ba14092163815ca30e7f296705d00013ed84e124..8c4a47a320b6a38fff9ab7c9622f79c3caae6b4b 100644
--- a/src/TNL/Matrices/ChunkedEllpack.h
+++ b/src/TNL/Matrices/ChunkedEllpack.h
@@ -224,10 +224,10 @@ public:
    void getTransposition( const ChunkedEllpack< Real2, Device, Index2 >& matrix,
                           const RealType& matrixMultiplicator = 1.0 );
 
-   template< typename Vector >
-   bool performSORIteration( const Vector& b,
+   template< typename Vector1, typename Vector2 >
+   bool performSORIteration( const Vector1& b,
                              const IndexType row,
-                             Vector& x,
+                             Vector2& x,
                              const RealType& omega = 1.0 ) const;
 
    // copy assignment
diff --git a/src/TNL/Matrices/ChunkedEllpack_impl.h b/src/TNL/Matrices/ChunkedEllpack_impl.h
index 5c5c71543c82c10cc2da55728fd7352d83b943b0..1a47fe4e608fa2d5087eb404730356517316a913 100644
--- a/src/TNL/Matrices/ChunkedEllpack_impl.h
+++ b/src/TNL/Matrices/ChunkedEllpack_impl.h
@@ -1154,11 +1154,11 @@ void ChunkedEllpack< Real, Device, Index >::getTransposition( const ChunkedEllpa
 template< typename Real,
           typename Device,
           typename Index >
-   template< typename Vector >
-bool ChunkedEllpack< Real, Device, Index >::performSORIteration( const Vector& b,
-                                                                                    const IndexType row,
-                                                                                    Vector& x,
-                                                                                    const RealType& omega ) const
+   template< typename Vector1, typename Vector2 >
+bool ChunkedEllpack< Real, Device, Index >::performSORIteration( const Vector1& b,
+                                                                 const IndexType row,
+                                                                 Vector2& x,
+                                                                 const RealType& omega ) const
 {
    TNL_ASSERT( row >=0 && row < this->getRows(),
               std::cerr << "row = " << row
diff --git a/src/TNL/Matrices/Dense.h b/src/TNL/Matrices/Dense.h
index 3904f5c059b210ae9b61aa0f1455d4a2ca762964..bb7a156c200f90979f0d0a79ae886d3d555e3d90 100644
--- a/src/TNL/Matrices/Dense.h
+++ b/src/TNL/Matrices/Dense.h
@@ -177,10 +177,10 @@ public:
    void getTransposition( const Matrix& matrix,
                           const RealType& matrixMultiplicator = 1.0 );
 
-   template< typename Vector >
-   void performSORIteration( const Vector& b,
+   template< typename Vector1, typename Vector2 >
+   void performSORIteration( const Vector1& b,
                              const IndexType row,
-                             Vector& x,
+                             Vector2& x,
                              const RealType& omega = 1.0 ) const;
 
    // copy assignment
diff --git a/src/TNL/Matrices/Dense_impl.h b/src/TNL/Matrices/Dense_impl.h
index 64d154779d7eaa2cc2ab7da0298c126015472703..36610ac78d9c00536c340a085ce7549605b86485 100644
--- a/src/TNL/Matrices/Dense_impl.h
+++ b/src/TNL/Matrices/Dense_impl.h
@@ -847,11 +847,11 @@ void Dense< Real, Device, Index >::getTransposition( const Matrix& matrix,
 template< typename Real,
           typename Device,
           typename Index >
-   template< typename Vector >
-void Dense< Real, Device, Index >::performSORIteration( const Vector& b,
-                                                                 const IndexType row,
-                                                                 Vector& x,
-                                                                 const RealType& omega ) const
+   template< typename Vector1, typename Vector2 >
+void Dense< Real, Device, Index >::performSORIteration( const Vector1& b,
+                                                        const IndexType row,
+                                                        Vector2& x,
+                                                        const RealType& omega ) const
 {
    RealType sum( 0.0 ), diagonalValue;
    for( IndexType i = 0; i < this->getColumns(); i++ )
diff --git a/src/TNL/Matrices/Ellpack.h b/src/TNL/Matrices/Ellpack.h
index 3b4b00b3e0a30bb3fa3298149e7c37f02699b317..38333685bfbc59cd94dec2197463ca40557a57e5 100644
--- a/src/TNL/Matrices/Ellpack.h
+++ b/src/TNL/Matrices/Ellpack.h
@@ -161,10 +161,10 @@ public:
    void getTransposition( const Ellpack< Real2, Device, Index2 >& matrix,
                           const RealType& matrixMultiplicator = 1.0 );
 
-   template< typename Vector >
-   bool performSORIteration( const Vector& b,
+   template< typename Vector1, typename Vector2 >
+   bool performSORIteration( const Vector1& b,
                              const IndexType row,
-                             Vector& x,
+                             Vector2& x,
                              const RealType& omega = 1.0 ) const;
 
    template< typename Vector >
diff --git a/src/TNL/Matrices/Ellpack_impl.h b/src/TNL/Matrices/Ellpack_impl.h
index be72eb42e354f1b59676871a4426cd820d78eee7..9801b6bcac54bdff89337428d3b83968ebf3759a 100644
--- a/src/TNL/Matrices/Ellpack_impl.h
+++ b/src/TNL/Matrices/Ellpack_impl.h
@@ -540,11 +540,11 @@ void Ellpack< Real, Device, Index >::getTransposition( const Ellpack< Real2, Dev
 template< typename Real,
           typename Device,
           typename Index >
-   template< typename Vector >
-bool Ellpack< Real, Device, Index > :: performSORIteration( const Vector& b,
-                                                                           const IndexType row,
-                                                                           Vector& x,
-                                                                           const RealType& omega ) const
+   template< typename Vector1, typename Vector2 >
+bool Ellpack< Real, Device, Index > :: performSORIteration( const Vector1& b,
+                                                            const IndexType row,
+                                                            Vector2& x,
+                                                            const RealType& omega ) const
 {
    TNL_ASSERT( row >=0 && row < this->getRows(),
               std::cerr << "row = " << row
diff --git a/src/TNL/Matrices/Multidiagonal.h b/src/TNL/Matrices/Multidiagonal.h
index 28de74b8ac64b9fdf002aa71a4ca6a34a6ad4b26..9b8f18779ceb9afa15e3e27f3610a2b0fa23fde6 100644
--- a/src/TNL/Matrices/Multidiagonal.h
+++ b/src/TNL/Matrices/Multidiagonal.h
@@ -177,10 +177,10 @@ public:
    void getTransposition( const Multidiagonal< Real2, Device, Index2 >& matrix,
                           const RealType& matrixMultiplicator = 1.0 );
 
-   template< typename Vector >
-   bool performSORIteration( const Vector& b,
+   template< typename Vector1, typename Vector2 >
+   bool performSORIteration( const Vector1& b,
                              const IndexType row,
-                             Vector& x,
+                             Vector2& x,
                              const RealType& omega = 1.0 ) const;
 
    // copy assignment
diff --git a/src/TNL/Matrices/Multidiagonal_impl.h b/src/TNL/Matrices/Multidiagonal_impl.h
index 47d827b93df0904b568c234cc8bfbea602479e46..5f7228d698db1a47dbc62f2b540c08b1e3f9b86c 100644
--- a/src/TNL/Matrices/Multidiagonal_impl.h
+++ b/src/TNL/Matrices/Multidiagonal_impl.h
@@ -576,11 +576,11 @@ void Multidiagonal< Real, Device, Index >::getTransposition( const Multidiagonal
 template< typename Real,
           typename Device,
           typename Index >
-   template< typename Vector >
-bool Multidiagonal< Real, Device, Index > :: performSORIteration( const Vector& b,
-                                                                           const IndexType row,
-                                                                           Vector& x,
-                                                                           const RealType& omega ) const
+   template< typename Vector1, typename Vector2 >
+bool Multidiagonal< Real, Device, Index > :: performSORIteration( const Vector1& b,
+                                                                  const IndexType row,
+                                                                  Vector2& x,
+                                                                  const RealType& omega ) const
 {
    TNL_ASSERT( row >=0 && row < this->getRows(),
               std::cerr << "row = " << row
diff --git a/src/TNL/Matrices/SlicedEllpack.h b/src/TNL/Matrices/SlicedEllpack.h
index 0557d26ebbd94ee57ab2bcd51bb029f48cbd30c7..815728d7a58d588a1791c1aa80b84bcd81da8f4b 100644
--- a/src/TNL/Matrices/SlicedEllpack.h
+++ b/src/TNL/Matrices/SlicedEllpack.h
@@ -187,10 +187,10 @@ public:
    void getTransposition( const SlicedEllpack< Real2, Device, Index2 >& matrix,
                           const RealType& matrixMultiplicator = 1.0 );
 
-   template< typename Vector >
-   bool performSORIteration( const Vector& b,
+   template< typename Vector1, typename Vector2 >
+   bool performSORIteration( const Vector1& b,
                              const IndexType row,
-                             Vector& x,
+                             Vector2& x,
                              const RealType& omega = 1.0 ) const;
 
    // copy assignment
diff --git a/src/TNL/Matrices/SlicedEllpack_impl.h b/src/TNL/Matrices/SlicedEllpack_impl.h
index 52288c650fbdadb7c06e2aba8933616fcd89786a..2ff01b49c51943c7411626257a57f273c6880b05 100644
--- a/src/TNL/Matrices/SlicedEllpack_impl.h
+++ b/src/TNL/Matrices/SlicedEllpack_impl.h
@@ -562,11 +562,11 @@ template< typename Real,
           typename Device,
           typename Index,
           int SliceSize >
-   template< typename Vector >
-bool SlicedEllpack< Real, Device, Index, SliceSize >::performSORIteration( const Vector& b,
-                                                                                    const IndexType row,
-                                                                                    Vector& x,
-                                                                                    const RealType& omega ) const
+   template< typename Vector1, typename Vector2 >
+bool SlicedEllpack< Real, Device, Index, SliceSize >::performSORIteration( const Vector1& b,
+                                                                           const IndexType row,
+                                                                           Vector2& x,
+                                                                           const RealType& omega ) const
 {
    TNL_ASSERT( row >=0 && row < this->getRows(),
               std::cerr << "row = " << row
diff --git a/src/TNL/Matrices/Tridiagonal.h b/src/TNL/Matrices/Tridiagonal.h
index 0769ca83f2f0355ee3ad4c726c601bd16bae8c3f..472cadffcd4194270d5218e3c0ea1415b2c7ae5c 100644
--- a/src/TNL/Matrices/Tridiagonal.h
+++ b/src/TNL/Matrices/Tridiagonal.h
@@ -167,11 +167,11 @@ public:
    void getTransposition( const Tridiagonal< Real2, Device, Index2 >& matrix,
                           const RealType& matrixMultiplicator = 1.0 );
 
-   template< typename Vector >
+   template< typename Vector1, typename Vector2 >
    __cuda_callable__
-   void performSORIteration( const Vector& b,
+   void performSORIteration( const Vector1& b,
                              const IndexType row,
-                             Vector& x,
+                             Vector2& x,
                              const RealType& omega = 1.0 ) const;
 
    // copy assignment
diff --git a/src/TNL/Matrices/Tridiagonal_impl.h b/src/TNL/Matrices/Tridiagonal_impl.h
index f3c073cd0d65285643df34d8fc2885e802c77e36..66fe9d7e80a80f93a1d16f1e741008f4de6a0787 100644
--- a/src/TNL/Matrices/Tridiagonal_impl.h
+++ b/src/TNL/Matrices/Tridiagonal_impl.h
@@ -535,12 +535,12 @@ void Tridiagonal< Real, Device, Index >::getTransposition( const Tridiagonal< Re
 template< typename Real,
           typename Device,
           typename Index >
-   template< typename Vector >
+   template< typename Vector1, typename Vector2 >
 __cuda_callable__
-void Tridiagonal< Real, Device, Index >::performSORIteration( const Vector& b,
-                                                                       const IndexType row,
-                                                                       Vector& x,
-                                                                       const RealType& omega ) const
+void Tridiagonal< Real, Device, Index >::performSORIteration( const Vector1& b,
+                                                              const IndexType row,
+                                                              Vector2& x,
+                                                              const RealType& omega ) const
 {
    RealType sum( 0.0 );
    if( row > 0 )