diff --git a/src/Benchmarks/BLAS/vector-operations.h b/src/Benchmarks/BLAS/vector-operations.h
index bb672a7ec464c7cfd1150f6fd0c986c0ede9ddce..f6006104f369e84ef31311ddd67b443b6dc071ab 100644
--- a/src/Benchmarks/BLAS/vector-operations.h
+++ b/src/Benchmarks/BLAS/vector-operations.h
@@ -117,10 +117,18 @@ benchmarkVectorOperations( Benchmark & benchmark,
    auto minCuda = [&]() {
       resultDevice = deviceVector.min();
    };
+   auto minHostET = [&]() {
+      resultHost = min( hostView );
+   };
+   auto minCudaET = [&]() {
+      resultDevice = min( deviceView );
+   };
    benchmark.setOperation( "min", datasetSize );
    benchmark.time< Devices::Host >( reset1, "CPU", minHost );
+   benchmark.time< Devices::Host >( reset1, "CPU", minHostET );
 #ifdef HAVE_CUDA
    benchmark.time< Devices::Cuda >( reset1, "GPU", minCuda );
+   benchmark.time< Devices::Cuda >( reset1, "GPU", minCudaET );
 #endif
 
 
@@ -136,7 +144,6 @@ benchmarkVectorOperations( Benchmark & benchmark,
    auto absMaxCudaET = [&]() {
       resultDevice = max( abs( deviceView ) );
    };
-   
 #ifdef HAVE_CUDA
    auto absMaxCublas = [&]() {
       int index = 0;
@@ -162,6 +169,12 @@ benchmarkVectorOperations( Benchmark & benchmark,
    auto absMinCuda = [&]() {
       resultDevice = deviceVector.absMin();
    };
+   auto absMinHostET = [&]() {
+      resultHost = min( abs( hostView ) );
+   };
+   auto absMinCudaET = [&]() {
+      resultDevice = min( abs( deviceView ) );
+   };
 #ifdef HAVE_CUDA
    auto absMinCublas = [&]() {
       int index = 0;
@@ -173,8 +186,10 @@ benchmarkVectorOperations( Benchmark & benchmark,
 #endif
    benchmark.setOperation( "absMin", datasetSize );
    benchmark.time< Devices::Host >( reset1, "CPU", absMinHost );
+   benchmark.time< Devices::Host >( reset1, "CPU", absMinHostET );
 #ifdef HAVE_CUDA
    benchmark.time< Devices::Cuda >( reset1, "GPU", absMinCuda );
+   benchmark.time< Devices::Cuda >( reset1, "GPU", absMinCudaET );
    benchmark.time< Devices::Cuda >( reset1, "cuBLAS", absMinCublas );
 #endif
 
@@ -185,10 +200,18 @@ benchmarkVectorOperations( Benchmark & benchmark,
    auto sumCuda = [&]() {
       resultDevice = deviceVector.sum();
    };
+   auto sumHostET = [&]() {
+      resultHost = sum( hostView );
+   };
+   auto sumCudaET = [&]() {
+      resultDevice = sum( deviceView );
+   };
    benchmark.setOperation( "sum", datasetSize );
    benchmark.time< Devices::Host >( reset1, "CPU", sumHost );
+   benchmark.time< Devices::Host >( reset1, "CPU", sumHostET );
 #ifdef HAVE_CUDA
    benchmark.time< Devices::Cuda >( reset1, "GPU", sumCuda );
+   benchmark.time< Devices::Cuda >( reset1, "GPU", sumCudaET );
 #endif
 
 
@@ -198,6 +221,12 @@ benchmarkVectorOperations( Benchmark & benchmark,
    auto l1normCuda = [&]() {
       resultDevice = deviceVector.lpNorm( 1.0 );
    };
+   auto l1normHostET = [&]() {
+      resultHost = lpNorm( hostView, 1.0 );
+   };
+   auto l1normCudaET = [&]() {
+      resultDevice = lpNorm( deviceView, 1.0 );
+   };
 #ifdef HAVE_CUDA
    auto l1normCublas = [&]() {
       cublasGasum( cublasHandle, size,
@@ -207,8 +236,10 @@ benchmarkVectorOperations( Benchmark & benchmark,
 #endif
    benchmark.setOperation( "l1 norm", datasetSize );
    benchmark.time< Devices::Host >( reset1, "CPU", l1normHost );
+   benchmark.time< Devices::Host >( reset1, "CPU", l1normHostET );
 #ifdef HAVE_CUDA
    benchmark.time< Devices::Cuda >( reset1, "GPU", l1normCuda );
+   benchmark.time< Devices::Cuda >( reset1, "GPU", l1normCudaET );
    benchmark.time< Devices::Cuda >( reset1, "cuBLAS", l1normCublas );
 #endif
 
@@ -219,6 +250,12 @@ benchmarkVectorOperations( Benchmark & benchmark,
    auto l2normCuda = [&]() {
       resultDevice = deviceVector.lpNorm( 2.0 );
    };
+   auto l2normHostET = [&]() {
+      resultHost = lpNorm( hostView, 2.0 );
+   };
+   auto l2normCudaET = [&]() {
+      resultDevice = lpNorm( deviceView, 2.0 );
+   };
 #ifdef HAVE_CUDA
    auto l2normCublas = [&]() {
       cublasGnrm2( cublasHandle, size,
@@ -228,8 +265,10 @@ benchmarkVectorOperations( Benchmark & benchmark,
 #endif
    benchmark.setOperation( "l2 norm", datasetSize );
    benchmark.time< Devices::Host >( reset1, "CPU", l2normHost );
+   benchmark.time< Devices::Host >( reset1, "CPU", l2normHostET );
 #ifdef HAVE_CUDA
    benchmark.time< Devices::Cuda >( reset1, "GPU", l2normCuda );
+   benchmark.time< Devices::Cuda >( reset1, "GPU", l2normCudaET );
    benchmark.time< Devices::Cuda >( reset1, "cuBLAS", l2normCublas );
 #endif
 
@@ -240,10 +279,19 @@ benchmarkVectorOperations( Benchmark & benchmark,
    auto l3normCuda = [&]() {
       resultDevice = deviceVector.lpNorm( 3.0 );
    };
+   auto l3normHostET = [&]() {
+      resultHost = lpNorm( hostView, 3.0 );
+   };
+   auto l3normCudaET = [&]() {
+      resultDevice = lpNorm( deviceView, 3.0 );
+   };
+
    benchmark.setOperation( "l3 norm", datasetSize );
    benchmark.time< Devices::Host >( reset1, "CPU", l3normHost );
+   benchmark.time< Devices::Host >( reset1, "CPU", l3normHostET );
 #ifdef HAVE_CUDA
    benchmark.time< Devices::Cuda >( reset1, "GPU", l3normCuda );
+   benchmark.time< Devices::Cuda >( reset1, "GPU", l3normCudaET );
 #endif
 
 
diff --git a/src/TNL/Containers/Expressions/ExpressionTemplates.h b/src/TNL/Containers/Expressions/ExpressionTemplates.h
index 280b2817e99e5c0278ec03ba50cc570a14564389..adfab7b95442be32cceb382a62b848ee3c1f2a14 100644
--- a/src/TNL/Containers/Expressions/ExpressionTemplates.h
+++ b/src/TNL/Containers/Expressions/ExpressionTemplates.h
@@ -278,7 +278,6 @@ template< typename L1,
           typename R1,
           typename R2,
           template< typename, typename > class ROperation >
-__cuda_callable__
 const Containers::Expressions::BinaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >,
@@ -295,7 +294,6 @@ operator + ( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LO
 template< typename T1,
           typename T2,
           template< typename, typename > class Operation >
-__cuda_callable__
 const Containers::Expressions::BinaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >,
    typename Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >::RealType,
@@ -314,7 +312,6 @@ template< typename L1,
           typename R1,
           typename R2,
           template< typename, typename > class ROperation >
-__cuda_callable__
 const Containers::Expressions::BinaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    typename Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >,
@@ -333,7 +330,6 @@ template< typename L1,
           template< typename, typename > class LOperation,
           typename R1,
           template< typename > class ROperation >
-__cuda_callable__
 const Containers::Expressions::BinaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    typename Containers::Expressions::UnaryExpressionTemplate< R1, ROperation >,
@@ -355,7 +351,6 @@ template< typename L1,
           typename R1,
           typename R2,
           template< typename, typename > class ROperation >
-__cuda_callable__
 const Containers::Expressions::BinaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >,
@@ -372,7 +367,6 @@ operator - ( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LO
 template< typename T1,
           typename T2,
           template< typename, typename > class Operation >
-__cuda_callable__
 const Containers::Expressions::BinaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >,
    typename Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >::RealType,
@@ -391,7 +385,6 @@ template< typename L1,
           typename R1,
           typename R2,
           template< typename, typename > class ROperation >
-__cuda_callable__
 const Containers::Expressions::BinaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    typename Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >,
@@ -410,7 +403,6 @@ template< typename L1,
           template< typename, typename > class LOperation,
           typename R1,
           template< typename > class ROperation >
-__cuda_callable__
 const Containers::Expressions::BinaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    typename Containers::Expressions::UnaryExpressionTemplate< R1, ROperation >,
@@ -432,7 +424,6 @@ template< typename L1,
           typename R1,
           typename R2,
           template< typename, typename > class ROperation >
-__cuda_callable__
 const Containers::Expressions::BinaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >,
@@ -449,7 +440,6 @@ operator * ( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LO
 template< typename T1,
           typename T2,
           template< typename, typename > class Operation >
-__cuda_callable__
 const Containers::Expressions::BinaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >,
    typename Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >::RealType,
@@ -468,7 +458,6 @@ template< typename L1,
           typename R1,
           typename R2,
           template< typename, typename > class ROperation >
-__cuda_callable__
 const Containers::Expressions::BinaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    typename Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >,
@@ -487,7 +476,6 @@ template< typename L1,
           template< typename, typename > class LOperation,
           typename R1,
           template< typename > class ROperation >
-__cuda_callable__
 const Containers::Expressions::BinaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    typename Containers::Expressions::UnaryExpressionTemplate< R1, ROperation >,
@@ -509,7 +497,6 @@ template< typename L1,
           typename R1,
           typename R2,
           template< typename, typename > class ROperation >
-__cuda_callable__
 const Containers::Expressions::BinaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >,
@@ -526,7 +513,6 @@ operator / ( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LO
 template< typename T1,
           typename T2,
           template< typename, typename > class Operation >
-__cuda_callable__
 const Containers::Expressions::BinaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >,
    typename Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >::RealType,
@@ -545,7 +531,6 @@ template< typename L1,
           typename R1,
           typename R2,
           template< typename, typename > class ROperation >
-__cuda_callable__
 const Containers::Expressions::BinaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    typename Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >,
@@ -564,7 +549,6 @@ template< typename L1,
           template< typename, typename > class LOperation,
           typename R1,
           template< typename > class ROperation >
-__cuda_callable__
 const Containers::Expressions::BinaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    typename Containers::Expressions::UnaryExpressionTemplate< R1, ROperation >,
@@ -586,7 +570,6 @@ template< typename L1,
           typename R1,
           typename R2,
           template< typename, typename > class ROperation >
-__cuda_callable__
 const Containers::Expressions::BinaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >,
@@ -603,7 +586,6 @@ min ( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperatio
 template< typename T1,
           typename T2,
           template< typename, typename > class Operation >
-__cuda_callable__
 const Containers::Expressions::BinaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >,
    typename Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >::RealType,
@@ -622,7 +604,6 @@ template< typename L1,
           typename R1,
           typename R2,
           template< typename, typename > class ROperation >
-__cuda_callable__
 const Containers::Expressions::BinaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    typename Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >,
@@ -641,7 +622,6 @@ template< typename L1,
           template< typename, typename > class LOperation,
           typename R1,
           template< typename > class ROperation >
-__cuda_callable__
 const Containers::Expressions::BinaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    typename Containers::Expressions::UnaryExpressionTemplate< R1, ROperation >,
@@ -663,7 +643,6 @@ template< typename L1,
           typename R1,
           typename R2,
           template< typename, typename > class ROperation >
-__cuda_callable__
 const Containers::Expressions::BinaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >,
@@ -680,7 +659,6 @@ max( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation
 template< typename T1,
           typename T2,
           template< typename, typename > class Operation >
-__cuda_callable__
 const Containers::Expressions::BinaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >,
    typename Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >::RealType,
@@ -699,7 +677,6 @@ template< typename L1,
           typename R1,
           typename R2,
           template< typename, typename > class ROperation >
-__cuda_callable__
 const Containers::Expressions::BinaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    typename Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >,
@@ -718,7 +695,6 @@ template< typename L1,
           template< typename, typename > class LOperation,
           typename R1,
           template< typename > class ROperation >
-__cuda_callable__
 const Containers::Expressions::BinaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    typename Containers::Expressions::UnaryExpressionTemplate< R1, ROperation >,
@@ -740,7 +716,6 @@ template< typename L1,
           typename R1,
           typename R2,
           template< typename, typename > class ROperation >
-__cuda_callable__
 bool
 operator == ( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >& a,
              const Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >& b )
@@ -751,7 +726,6 @@ operator == ( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, L
 template< typename T1,
           typename T2,
           template< typename, typename > class Operation >
-__cuda_callable__
 bool
 operator == ( const Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >& a,
               const typename Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >::RealType& b )
@@ -764,7 +738,6 @@ template< typename L1,
           typename R1,
           typename R2,
           template< typename, typename > class ROperation >
-__cuda_callable__
 bool
 operator == ( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >& a,
               const typename Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >& b )
@@ -777,7 +750,6 @@ template< typename L1,
           template< typename, typename > class LOperation,
           typename R1,
           template< typename > class ROperation >
-__cuda_callable__
 bool
 operator == ( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >& a,
              const typename Containers::Expressions::UnaryExpressionTemplate< R1,ROperation >& b )
@@ -793,7 +765,6 @@ template< typename L1,
           typename R1,
           typename R2,
           template< typename, typename > class ROperation >
-__cuda_callable__
 bool
 operator != ( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >& a,
              const Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >& b )
@@ -804,7 +775,6 @@ operator != ( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, L
 template< typename T1,
           typename T2,
           template< typename, typename > class Operation >
-__cuda_callable__
 bool
 operator != ( const Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >& a,
               const typename Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >::RealType& b )
@@ -817,7 +787,6 @@ template< typename L1,
           typename R1,
           typename R2,
           template< typename, typename > class ROperation >
-__cuda_callable__
 bool
 operator != ( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >& a,
               const typename Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >& b )
@@ -830,7 +799,6 @@ template< typename L1,
           template< typename, typename > class LOperation,
           typename R1,
           template< typename > class ROperation >
-__cuda_callable__
 bool
 operator != ( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >& a,
               const typename Containers::Expressions::UnaryExpressionTemplate< R1, ROperation >& b )
@@ -846,7 +814,6 @@ template< typename L1,
           typename R1,
           typename R2,
           template< typename, typename > class ROperation >
-__cuda_callable__
 bool
 operator < ( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >& a,
              const Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >& b )
@@ -857,7 +824,6 @@ operator < ( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LO
 template< typename T1,
           typename T2,
           template< typename, typename > class Operation >
-__cuda_callable__
 bool
 operator < ( const Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >& a,
              const typename Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >::RealType& b )
@@ -870,7 +836,6 @@ template< typename L1,
           typename R1,
           typename R2,
           template< typename, typename > class ROperation >
-__cuda_callable__
 bool
 operator < ( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >& a,
               const typename Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >& b )
@@ -883,7 +848,6 @@ template< typename L1,
           template< typename, typename > class LOperation,
           typename R1,
           template< typename > class ROperation >
-__cuda_callable__
 bool
 operator < ( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >& a,
               const typename Containers::Expressions::UnaryExpressionTemplate< R1, ROperation >& b )
@@ -899,7 +863,6 @@ template< typename L1,
           typename R1,
           typename R2,
           template< typename, typename > class ROperation >
-__cuda_callable__
 bool
 operator <= ( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >& a,
              const Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >& b )
@@ -910,7 +873,6 @@ operator <= ( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, L
 template< typename T1,
           typename T2,
           template< typename, typename > class Operation >
-__cuda_callable__
 bool
 operator <= ( const Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >& a,
              const typename Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >::RealType& b )
@@ -923,7 +885,6 @@ template< typename L1,
           typename R1,
           typename R2,
           template< typename, typename > class ROperation >
-__cuda_callable__
 bool
 operator <= ( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >& a,
               const typename Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >& b )
@@ -936,7 +897,6 @@ template< typename L1,
           template< typename, typename > class LOperation,
           typename R1,
           template< typename > class ROperation >
-__cuda_callable__
 bool
 operator <= ( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >& a,
               const typename Containers::Expressions::UnaryExpressionTemplate< R1, ROperation >& b )
@@ -952,7 +912,6 @@ template< typename L1,
           typename R1,
           typename R2,
           template< typename, typename > class ROperation >
-__cuda_callable__
 bool
 operator > ( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >& a,
              const Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >& b )
@@ -963,7 +922,6 @@ operator > ( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LO
 template< typename T1,
           typename T2,
           template< typename, typename > class Operation >
-__cuda_callable__
 bool
 operator > ( const Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >& a,
              const typename Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >::RealType& b )
@@ -976,7 +934,6 @@ template< typename L1,
           typename R1,
           typename R2,
           template< typename, typename > class ROperation >
-__cuda_callable__
 bool
 operator > ( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >& a,
              const typename Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >& b )
@@ -989,7 +946,6 @@ template< typename L1,
           template< typename, typename > class LOperation,
           typename R1,
           template< typename > class ROperation >
-__cuda_callable__
 bool
 operator > ( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >& a,
              const typename Containers::Expressions::UnaryExpressionTemplate< R1, ROperation >& b )
@@ -1005,7 +961,6 @@ template< typename L1,
           typename R1,
           typename R2,
           template< typename, typename > class ROperation >
-__cuda_callable__
 bool
 operator >= ( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >& a,
               const Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >& b )
@@ -1016,7 +971,6 @@ operator >= ( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, L
 template< typename T1,
           typename T2,
           template< typename, typename > class Operation >
-__cuda_callable__
 bool
 operator >= ( const Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >& a,
               const typename Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >::RealType& b )
@@ -1029,7 +983,6 @@ template< typename L1,
           typename R1,
           typename R2,
           template< typename, typename > class ROperation >
-__cuda_callable__
 bool
 operator >= ( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >& a,
               const typename Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >& b )
@@ -1042,7 +995,6 @@ template< typename L1,
           template< typename, typename > class LOperation,
           typename R1,
           template< typename > class ROperation >
-__cuda_callable__
 bool
 operator >= ( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >& a,
               const typename Containers::Expressions::UnaryExpressionTemplate< R1, ROperation >& b )
@@ -1059,7 +1011,6 @@ operator >= ( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, L
 template< typename L1,
           typename L2,
           template< typename, typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::Minus >
@@ -1072,7 +1023,6 @@ operator -( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOp
 
 template< typename L1,
           template< typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    Containers::Expressions::Abs >
@@ -1088,7 +1038,6 @@ operator -( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperati
 template< typename L1,
           typename L2,
           template< typename, typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::Abs >
@@ -1101,7 +1050,6 @@ abs( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation
 
 template< typename L1,
           template< typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    Containers::Expressions::Abs >
@@ -1117,7 +1065,6 @@ abs( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >& a
 template< typename L1,
           typename L2,
           template< typename, typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::Sin >
@@ -1130,7 +1077,6 @@ sin( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation
 
 template< typename L1,
           template< typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    Containers::Expressions::Sin >
@@ -1146,7 +1092,6 @@ sin( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >& a
 template< typename L1,
           typename L2,
           template< typename, typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::Cos >
@@ -1159,7 +1104,6 @@ cos( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation
 
 template< typename L1,
           template< typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    Containers::Expressions::Cos >
@@ -1175,7 +1119,6 @@ cos( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >& a
 template< typename L1,
           typename L2,
           template< typename, typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::Tan >
@@ -1188,7 +1131,6 @@ tan( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation
 
 template< typename L1,
           template< typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    Containers::Expressions::Tan >
@@ -1204,7 +1146,6 @@ tan( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >& a
 template< typename L1,
           typename L2,
           template< typename, typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::Sqrt >
@@ -1217,7 +1158,6 @@ sqrt( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperatio
 
 template< typename L1,
           template< typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    Containers::Expressions::Sqrt >
@@ -1233,7 +1173,6 @@ sqrt( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >&
 template< typename L1,
           typename L2,
           template< typename, typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::Cbrt >
@@ -1246,7 +1185,6 @@ cbrt( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperatio
 
 template< typename L1,
           template< typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    Containers::Expressions::Cbrt >
@@ -1263,7 +1201,6 @@ template< typename L1,
           typename L2,
           template< typename, typename > class LOperation,
           typename Real >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::Pow >
@@ -1279,7 +1216,6 @@ pow( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation
 template< typename L1,
           template< typename > class LOperation,
           typename Real >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    Containers::Expressions::Pow >
@@ -1297,7 +1233,6 @@ pow( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >& a
 template< typename L1,
           typename L2,
           template< typename, typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::Sin >
@@ -1310,7 +1245,6 @@ floor( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperati
 
 template< typename L1,
           template< typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    Containers::Expressions::Floor >
@@ -1326,7 +1260,6 @@ floor( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >&
 template< typename L1,
           typename L2,
           template< typename, typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::Ceil >
@@ -1339,7 +1272,6 @@ ceil( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperatio
 
 template< typename L1,
           template< typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    Containers::Expressions::Ceil >
@@ -1355,7 +1287,6 @@ sin( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >& a
 template< typename L1,
           typename L2,
           template< typename, typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::Asin >
@@ -1368,7 +1299,6 @@ asin( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperatio
 
 template< typename L1,
           template< typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    Containers::Expressions::Asin >
@@ -1384,7 +1314,6 @@ asin( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >&
 template< typename L1,
           typename L2,
           template< typename, typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::Acos >
@@ -1397,7 +1326,6 @@ cos( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation
 
 template< typename L1,
           template< typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    Containers::Expressions::Acos >
@@ -1413,7 +1341,6 @@ acos( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >&
 template< typename L1,
           typename L2,
           template< typename, typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::Atan >
@@ -1426,7 +1353,6 @@ tan( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation
 
 template< typename L1,
           template< typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    Containers::Expressions::Atan >
@@ -1442,7 +1368,6 @@ atan( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >&
 template< typename L1,
           typename L2,
           template< typename, typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::Sinh >
@@ -1455,7 +1380,6 @@ sinh( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperatio
 
 template< typename L1,
           template< typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    Containers::Expressions::Sinh >
@@ -1471,7 +1395,6 @@ sinh( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >&
 template< typename L1,
           typename L2,
           template< typename, typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::Cosh >
@@ -1484,7 +1407,6 @@ cosh( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperatio
 
 template< typename L1,
           template< typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    Containers::Expressions::Cosh >
@@ -1500,7 +1422,6 @@ cosh( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >&
 template< typename L1,
           typename L2,
           template< typename, typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::Tanh >
@@ -1513,7 +1434,6 @@ cosh( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperatio
 
 template< typename L1,
           template< typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    Containers::Expressions::Tanh >
@@ -1529,7 +1449,6 @@ tanh( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >&
 template< typename L1,
           typename L2,
           template< typename, typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::Log >
@@ -1542,7 +1461,6 @@ log( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation
 
 template< typename L1,
           template< typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    Containers::Expressions::Log >
@@ -1558,7 +1476,6 @@ log( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >& a
 template< typename L1,
           typename L2,
           template< typename, typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::Log10 >
@@ -1571,7 +1488,6 @@ log10( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperati
 
 template< typename L1,
           template< typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    Containers::Expressions::Log10 >
@@ -1587,7 +1503,6 @@ log10( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >&
 template< typename L1,
           typename L2,
           template< typename, typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::Log2 >
@@ -1600,7 +1515,6 @@ log2( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperatio
 
 template< typename L1,
           template< typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    Containers::Expressions::Log2 >
@@ -1616,7 +1530,6 @@ log2( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >&
 template< typename L1,
           typename L2,
           template< typename, typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >,
    Containers::Expressions::Exp >
@@ -1629,7 +1542,6 @@ exp( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation
 
 template< typename L1,
           template< typename > class LOperation >
-__cuda_callable__
 const Containers::Expressions::UnaryExpressionTemplate<
    Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >,
    Containers::Expressions::Exp >
@@ -1770,6 +1682,60 @@ binaryAnd( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperatio
    return ExpressionBinaryAnd( a );
 }
 
+
+////
+// Scalar product
+template< typename L1,
+          typename L2,
+          template< typename, typename > class LOperation,
+          typename R1,
+          typename R2,
+          template< typename, typename > class ROperation >
+auto
+operator,( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >& a,
+           const Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >& b )
+-> decltype( TNL::sum( a * b ) )
+{
+   return TNL::sum( a * b );
+}
+
+template< typename T1,
+          typename T2,
+          template< typename, typename > class Operation >
+auto
+operator,( const Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >& a,
+           const typename Containers::Expressions::BinaryExpressionTemplate< T1, T2, Operation >::RealType& b )
+-> decltype( TNL::sum( a * b ) )
+{
+   return TNL::sum( a * b );
+}
+
+template< typename L1,
+          template< typename > class LOperation,
+          typename R1,
+          typename R2,
+          template< typename, typename > class ROperation >
+auto
+operator,( const Containers::Expressions::UnaryExpressionTemplate< L1, LOperation >& a,
+           const typename Containers::Expressions::BinaryExpressionTemplate< R1, R2, ROperation >& b )
+-> decltype( TNL::sum( a * b ) )
+{
+   return TNL::sum( a * b );
+}
+
+template< typename L1,
+          typename L2,
+          template< typename, typename > class LOperation,
+          typename R1,
+          template< typename > class ROperation >
+auto
+operator,( const Containers::Expressions::BinaryExpressionTemplate< L1, L2, LOperation >& a,
+           const typename Containers::Expressions::UnaryExpressionTemplate< R1,ROperation >& b )
+-> decltype( TNL::sum( a * b ) )
+{
+   return TNL::sum( a * b );
+}
+
 ////
 // Output stream
 template< typename T1,
diff --git a/src/TNL/Containers/Expressions/StaticExpressionTemplates.h b/src/TNL/Containers/Expressions/StaticExpressionTemplates.h
index a1d6ae63974d68ee4830e6b03424acf5a093958e..6a97948bccc8e2147f89a87fdedd8ba740d821b4 100644
--- a/src/TNL/Containers/Expressions/StaticExpressionTemplates.h
+++ b/src/TNL/Containers/Expressions/StaticExpressionTemplates.h
@@ -1824,4 +1824,61 @@ binaryOr( const Containers::Expressions::StaticUnaryExpressionTemplate< L1, LOpe
    return StaticExpressionBinaryOr( a );
 }
 
+////
+// Scalar product
+template< typename L1,
+          typename L2,
+          template< typename, typename > class LOperation,
+          typename R1,
+          typename R2,
+          template< typename, typename > class ROperation >
+__cuda_callable__
+auto
+operator,( const Containers::Expressions::StaticBinaryExpressionTemplate< L1, L2, LOperation >& a,
+           const Containers::Expressions::StaticBinaryExpressionTemplate< R1, R2, ROperation >& b )
+-> decltype( TNL::sum( a * b ) )
+{
+   return TNL::sum( a * b );
+}
+
+template< typename T1,
+          typename T2,
+          template< typename, typename > class Operation >
+__cuda_callable__
+auto
+operator,( const Containers::Expressions::StaticBinaryExpressionTemplate< T1, T2, Operation >& a,
+           const typename Containers::Expressions::StaticBinaryExpressionTemplate< T1, T2, Operation >::RealType& b )
+-> decltype( TNL::sum( a * b ) )
+{
+   return TNL::sum( a * b );
+}
+
+template< typename L1,
+          template< typename > class LOperation,
+          typename R1,
+          typename R2,
+          template< typename, typename > class ROperation >
+__cuda_callable__
+auto
+operator,( const Containers::Expressions::StaticUnaryExpressionTemplate< L1, LOperation >& a,
+           const typename Containers::Expressions::StaticBinaryExpressionTemplate< R1, R2, ROperation >& b )
+-> decltype( TNL::sum( a * b ) )
+{
+   return TNL::sum( a * b );
+}
+
+template< typename L1,
+          typename L2,
+          template< typename, typename > class LOperation,
+          typename R1,
+          template< typename > class ROperation >
+__cuda_callable__
+auto
+operator,( const Containers::Expressions::StaticBinaryExpressionTemplate< L1, L2, LOperation >& a,
+           const typename Containers::Expressions::StaticUnaryExpressionTemplate< R1,ROperation >& b )
+-> decltype( TNL::sum( a * b ) )
+{
+   return TNL::sum( a * b );
+}
+
 } // namespace TNL
diff --git a/src/TNL/Containers/StaticVectorExpressions.h b/src/TNL/Containers/StaticVectorExpressions.h
index 311cc07d22447bb62de49eb5e33e14d0487cd355..9b6b78d64b29961579d77cb7248be1e3b16f90fc 100644
--- a/src/TNL/Containers/StaticVectorExpressions.h
+++ b/src/TNL/Containers/StaticVectorExpressions.h
@@ -588,6 +588,32 @@ binaryAnd( const Containers::StaticVector< Size, Real >& a )
    return Containers::Expressions::StaticExpressionBinaryAnd( a );
 }
 
+////
+// Scalar product
+template< int Size, typename Real, typename ET >
+__cuda_callable__
+auto operator,( const Containers::StaticVector< Size, Real >& a, const ET& b )
+->decltype( TNL::sum( a * b ) )
+{
+   return TNL::sum( a * b );
+}
+
+template< typename ET, int Size, typename Real >
+__cuda_callable__
+auto operator,( const ET& a, const Containers::StaticVector< Size, Real >& b )
+->decltype( TNL::sum( a * b ) )
+{
+   return TNL::sum( a * b );
+}
+
+template< typename Real1, int Size, typename Real2 >
+__cuda_callable__
+auto operator,( const Containers::StaticVector< Size, Real1 >& a, const Containers::StaticVector< Size, Real2 >& b )
+->decltype( TNL::sum( a * b ) )
+{
+   return TNL::sum( a * b );
+}
+
 ////
 // TODO: Replace this with multiplication when its safe
 template< int Size, typename Real, typename ET >
diff --git a/src/TNL/Containers/Vector.h b/src/TNL/Containers/Vector.h
index 457a556e2508a9855f85448bdf1dbaefae141ac6..08405560780dc484388eaaa7cec14fdb7aecad12 100644
--- a/src/TNL/Containers/Vector.h
+++ b/src/TNL/Containers/Vector.h
@@ -18,7 +18,7 @@ namespace Containers {
 
 template< typename Real, typename Device, typename Index >
 class VectorView;
-   
+
 /**
  * \brief This class extends TNL::Array with algebraic operations.
  *
@@ -180,13 +180,13 @@ public:
    __cuda_callable__ const Real& operator[]( const Index& i ) const;
 
    Vector& operator = ( const Vector& v );
-   
+
    template< typename Real_, typename Device_, typename Index_ >
    Vector& operator = ( const Vector< Real_, Device_, Index_ >& v );
 
    template< typename Real_, typename Device_, typename Index_ >
    Vector& operator = ( const VectorView< Real_, Device_, Index_ >& v );
-   
+
    template< typename VectorExpression >
    Vector& operator = ( const VectorExpression& expression );
 
@@ -226,6 +226,14 @@ public:
    template< typename Scalar >
    Vector& operator /= ( const Scalar c );
 
+   /**
+    * \brief Scalar product
+    * @param v
+    * @return
+    */
+   template< typename Vector_ >
+   Real operator, ( const Vector_& v ) const;
+
    /**
     * \brief Returns the maximum value out of all vector elements.
     */
diff --git a/src/TNL/Containers/Vector.hpp b/src/TNL/Containers/Vector.hpp
index dc6fce8b26ce4be41cc17e960b5e0e42c845f0b8..26b34d8039fe25c89ccfbc4c4395cb81a1cad112 100644
--- a/src/TNL/Containers/Vector.hpp
+++ b/src/TNL/Containers/Vector.hpp
@@ -329,6 +329,17 @@ Real Vector< Real, Device, Index >::min() const
 }
 
 
+template< typename Real,
+          typename Device,
+          typename Index >
+   template< typename Vector_ >
+Real Vector< Real, Device, Index >::
+operator,( const Vector_& v ) const
+{
+   static_assert( std::is_same< DeviceType, typename Vector_::DeviceType >::value, "Cannot compute product of vectors allocated on different devices." );
+   return Algorithms::VectorOperations< Device >::getScalarProduct( *this, v );
+}
+
 template< typename Real,
           typename Device,
           typename Index >
diff --git a/src/TNL/Containers/VectorView.h b/src/TNL/Containers/VectorView.h
index 6b86263976ee1cdb19283bee2931037901e89895..59dc29f83f5f5a74cf0a2e872ebf1f8852b27df3 100644
--- a/src/TNL/Containers/VectorView.h
+++ b/src/TNL/Containers/VectorView.h
@@ -127,6 +127,14 @@ public:
    template< typename Real_, typename Device_, typename Index_ >
    bool operator!=( const VectorView< Real_, Device_, Index_ >& v );
 
+   /**
+    * \brief Scalar product
+    * @param v
+    * @return
+    */
+   template< typename Vector_ >
+   NonConstReal operator, ( const Vector_& v ) const;
+
    NonConstReal max() const;
 
    NonConstReal min() const;
diff --git a/src/TNL/Containers/VectorViewExpressions.h b/src/TNL/Containers/VectorViewExpressions.h
index 077a6651a8990ccf9347ef493e3483f0f1cc005c..f83641f2f294ad504bff6c12622b865e0c0727f3 100644
--- a/src/TNL/Containers/VectorViewExpressions.h
+++ b/src/TNL/Containers/VectorViewExpressions.h
@@ -515,36 +515,6 @@ sign( const Containers::VectorView< Real, Device, Index >& a )
    return Containers::Expressions::UnaryExpressionTemplate< Containers::VectorView< Real, Device, Index >, Containers::Expressions::Sign >( a );
 }
 
-
-////
-// TODO: Replace this with multiplication when its safe
-template< typename Real, typename Device, typename Index, typename ET >
-__cuda_callable__
-Containers::VectorView< Real, Device, Index >
-Scale( const Containers::VectorView< Real, Device, Index >& a, const ET& b )
-{
-   Containers::VectorView< Real, Device, Index > result = Containers::Expressions::BinaryExpressionTemplate< Containers::VectorView< Real, Device, Index >, ET, Containers::Expressions::Multiplication >( a, b );
-   return result;
-}
-
-template< typename ET, typename Real, typename Device, typename Index >
-__cuda_callable__
-Containers::Expressions::BinaryExpressionTemplate< ET, Containers::VectorView< Real, Device, Index >, Containers::Expressions::Multiplication >
-Scale( const ET& a, const Containers::VectorView< Real, Device, Index >& b )
-{
-   Containers::VectorView< Real, Device, Index > result =  Containers::Expressions::BinaryExpressionTemplate< ET, Containers::VectorView< Real, Device, Index >, Containers::Expressions::Multiplication >( a, b );
-   return result;
-}
-
-template< typename Real1, typename Real2, typename Device, typename Index >
-__cuda_callable__
-Containers::Expressions::BinaryExpressionTemplate< Containers::VectorView< Real1, Device, Index >, Containers::VectorView< Real2, Device, Index >, Containers::Expressions::Multiplication >
-Scale( const Containers::VectorView< Real1, Device, Index >& a, const Containers::VectorView< Real2, Device, Index >& b )
-{
-   Containers::VectorView< Real1, Device, Index > result =  Containers::Expressions::BinaryExpressionTemplate< Containers::VectorView< Real1, Device, Index >, Containers::VectorView< Real2, Device, Index >, Containers::Expressions::Multiplication >( a, b );
-   return result;
-}
-
 ////
 // Vertical operations - min
 template< typename Real,
@@ -629,4 +599,57 @@ binaryAnd( const Containers::VectorView< Real, Device, Index >& a )
    return Containers::Expressions::ExpressionBinaryAnd( a );
 }
 
+////
+// Scalar product
+template< typename Real, typename Device, typename Index, typename ET >
+__cuda_callable__
+Real operator,( const Containers::VectorView< Real, Device, Index >& a, const ET& b )
+{
+   return TNL::sum( a * b );
+}
+
+template< typename ET, typename Real, typename Device, typename Index >
+__cuda_callable__
+Real operator,( const ET& a, const Containers::VectorView< Real, Device, Index >& b )
+{
+   return TNL::sum( a * b );
+}
+
+template< typename Real1, typename Real2, typename Device, typename Index >
+__cuda_callable__
+auto operator,( const Containers::VectorView< Real1, Device, Index >& a, const Containers::VectorView< Real2, Device, Index >& b )
+->decltype( TNL::sum( a * b ) )
+{
+   return TNL::sum( a * b );
+}
+
+////
+// TODO: Replace this with multiplication when its safe
+template< typename Real, typename Device, typename Index, typename ET >
+__cuda_callable__
+Containers::VectorView< Real, Device, Index >
+Scale( const Containers::VectorView< Real, Device, Index >& a, const ET& b )
+{
+   Containers::VectorView< Real, Device, Index > result = Containers::Expressions::BinaryExpressionTemplate< Containers::VectorView< Real, Device, Index >, ET, Containers::Expressions::Multiplication >( a, b );
+   return result;
+}
+
+template< typename ET, typename Real, typename Device, typename Index >
+__cuda_callable__
+Containers::Expressions::BinaryExpressionTemplate< ET, Containers::VectorView< Real, Device, Index >, Containers::Expressions::Multiplication >
+Scale( const ET& a, const Containers::VectorView< Real, Device, Index >& b )
+{
+   Containers::VectorView< Real, Device, Index > result =  Containers::Expressions::BinaryExpressionTemplate< ET, Containers::VectorView< Real, Device, Index >, Containers::Expressions::Multiplication >( a, b );
+   return result;
+}
+
+template< typename Real1, typename Real2, typename Device, typename Index >
+__cuda_callable__
+Containers::Expressions::BinaryExpressionTemplate< Containers::VectorView< Real1, Device, Index >, Containers::VectorView< Real2, Device, Index >, Containers::Expressions::Multiplication >
+Scale( const Containers::VectorView< Real1, Device, Index >& a, const Containers::VectorView< Real2, Device, Index >& b )
+{
+   Containers::VectorView< Real1, Device, Index > result =  Containers::Expressions::BinaryExpressionTemplate< Containers::VectorView< Real1, Device, Index >, Containers::VectorView< Real2, Device, Index >, Containers::Expressions::Multiplication >( a, b );
+   return result;
+}
+
 } // namespace TNL
diff --git a/src/TNL/Containers/VectorView_impl.h b/src/TNL/Containers/VectorView_impl.h
index d8b927d519bcc5056f48a4db5528a2683753b6cc..2433721b7fc33f9f9be34683d1701698632e2154 100644
--- a/src/TNL/Containers/VectorView_impl.h
+++ b/src/TNL/Containers/VectorView_impl.h
@@ -197,6 +197,18 @@ operator!=( const VectorView< Real_, Device_, Index_ >& v )
    return !ArrayView< Real, Device, Index >::operator ==( v );
 }
 
+template< typename Real,
+          typename Device,
+          typename Index >
+   template< typename Vector_ >
+typename VectorView< Real, Device, Index >::NonConstReal
+VectorView< Real, Device, Index >::
+operator,( const Vector_& v ) const
+{
+   static_assert( std::is_same< DeviceType, typename Vector_::DeviceType >::value, "Cannot compute product of vectors allocated on different devices." );
+   return Algorithms::VectorOperations< Device >::getScalarProduct( *this, v );
+}
+
 
 template< typename Real,
           typename Device,
diff --git a/src/UnitTests/Containers/VectorTest-7.h b/src/UnitTests/Containers/VectorTest-7.h
index 6a2c8b5904650d6d47afeeb81f9cd55cb1dc7bc4..4843383d7cf1e7d460ead504c2943fc966ed9072 100644
--- a/src/UnitTests/Containers/VectorTest-7.h
+++ b/src/UnitTests/Containers/VectorTest-7.h
@@ -76,6 +76,27 @@ TYPED_TEST( VectorTest, verticalOperations )
    EXPECT_NEAR( lpNorm( u - v, 2.0 ), l2NormDiff, 2.0e-5 );
 }
 
+TYPED_TEST( VectorTest, scalarProduct )
+{
+   using VectorType = typename TestFixture::VectorType;
+   using ViewType = typename TestFixture::ViewType;
+   using RealType = typename VectorType::RealType;
+   const int size = VECTOR_TEST_SIZE;
+
+   VectorType _u( size ), _v( size );
+   ViewType u( _u ), v( _v );
+   RealType aux( 0.0 );
+   for( int i = 0; i < size; i++ )
+   {
+      const RealType x = i;
+      const RealType y = size / 2 - i;
+      u.setElement( i, x );
+      v.setElement( i, y );
+      aux += x * y;
+   }
+   EXPECT_NEAR( ( u, v ), aux, 1.0e-5 );
+}
+
 #endif // HAVE_GTEST