diff --git a/src/algebra/BiCGStab.hpp b/src/algebra/BiCGStab.hpp
index ff34b89e148d844215cf7fffc21c2ad78bca5ec9..dc62e370e6ec304cb2e4fd44d220218d9bd94283 100644
--- a/src/algebra/BiCGStab.hpp
+++ b/src/algebra/BiCGStab.hpp
@@ -27,7 +27,7 @@ struct BiCGStab
 
     r_k_1 = b - A * x;
 
-    double residu = std::sqrt((r_k_1, r_k_1));   // Norm(r_k_1);
+    double residu = std::sqrt(dot(r_k_1, r_k_1));   // Norm(r_k_1);
 
     if (residu != 0) {
       double resid0 = residu;
@@ -53,23 +53,23 @@ struct BiCGStab
 
         Ap_k = A * p_k;
 
-        const double alpha_k = (r_k_1, rTilda_0) / (Ap_k, rTilda_0);
+        const double alpha_k = dot(r_k_1, rTilda_0) / dot(Ap_k, rTilda_0);
 
         s_k  = r_k_1 - alpha_k * Ap_k;
         As_k = A * s_k;
 
-        const double w_k = (As_k, s_k) / (As_k, As_k);
+        const double w_k = dot(As_k, s_k) / dot(As_k, As_k);
 
         x += alpha_k * p_k + w_k * s_k;
         r_k = s_k - w_k * As_k;
 
-        const double beta_k = (r_k, rTilda_0) / (r_k_1, rTilda_0) * (alpha_k / w_k);
+        const double beta_k = dot(r_k, rTilda_0) / dot(r_k_1, rTilda_0) * (alpha_k / w_k);
 
         p_k -= w_k * Ap_k;
         p_k *= beta_k;
         p_k += r_k;
 
-        if ((residu = std::sqrt((r_k, r_k))) / resid0 < epsilon) {
+        if ((residu = std::sqrt(dot(r_k, r_k))) / resid0 < epsilon) {
           break;
         }
 
diff --git a/src/algebra/CG.hpp b/src/algebra/CG.hpp
index 538c77571dac4d551d7476ea79bf855e8026fb53..7e653d8f95cddab4d5d1b7eb36481fd7690ef170 100644
--- a/src/algebra/CG.hpp
+++ b/src/algebra/CG.hpp
@@ -28,7 +28,7 @@ struct CG
     if (verbose) {
       h = A * x;
       h -= f;
-      std::cout << "- initial " << rang::style::bold << "real" << rang::style::reset << " residu :   " << (h, h)
+      std::cout << "- initial " << rang::style::bold << "real" << rang::style::reset << " residu :   " << dot(h, h)
                 << '\n';
     }
 
@@ -48,14 +48,14 @@ struct CG
 
         g = copy(cg);   // TODO: precond: g = g/C
 
-        gcg = (g, cg);
+        gcg = dot(g, cg);
 
         h = copy(g);
       }
 
       b = A * h;
 
-      double hAh = (h, b);
+      double hAh = dot(h, b);
 
       if (hAh == 0) {
         hAh = 1.;
@@ -69,7 +69,7 @@ struct CG
       g -= ro * b;
 
       double gamma = gcg;
-      gcg          = (g, cg);
+      gcg          = dot(g, cg);
 
       if ((i == 1) && (gcg != 0)) {
         relativeEpsilon = epsilon * gcg;
diff --git a/src/algebra/TinyVector.hpp b/src/algebra/TinyVector.hpp
index c3146a2be9932353388d34d15862056f952631fb..f24503aade42275e8b96f6958535fc8439da356c 100644
--- a/src/algebra/TinyVector.hpp
+++ b/src/algebra/TinyVector.hpp
@@ -64,11 +64,11 @@ class [[nodiscard]] TinyVector
   }
 
   PUGS_INLINE
-  constexpr T operator,(const TinyVector& v) const
+  constexpr friend T dot(const TinyVector& u, const TinyVector& v)
   {
-    T t = m_values[0] * v.m_values[0];
+    T t = u.m_values[0] * v.m_values[0];
     for (size_t i = 1; i < N; ++i) {
-      t += m_values[i] * v.m_values[i];
+      t += u.m_values[i] * v.m_values[i];
     }
     return t;
   }
@@ -239,7 +239,7 @@ l2Norm(const TinyVector<N, T>& x)
 {
   static_assert(std::is_arithmetic<T>(), "Cannot compute L2 norm for non-arithmetic types");
   static_assert(std::is_floating_point<T>::value, "L2 norm is defined for floating point types only");
-  return std::sqrt((x, x));
+  return std::sqrt(dot(x, x));
 }
 
 // Cross product is only defined for K^3 vectors
diff --git a/src/algebra/Vector.hpp b/src/algebra/Vector.hpp
index efe64b41a8940be9fadf7066dd5dc27bff90a633..089efa633161f4d3e83d16a5fc666ad5f9fc1696 100644
--- a/src/algebra/Vector.hpp
+++ b/src/algebra/Vector.hpp
@@ -33,18 +33,18 @@ class Vector   // LCOV_EXCL_LINE
     return x_copy;
   }
 
-  friend Vector operator*(const DataType& a, const Vector& x)
+  friend Vector
+  operator*(const DataType& a, const Vector& x)
   {
     Vector<std::remove_const_t<DataType>> y = copy(x);
     return y *= a;
   }
 
   template <typename DataType2>
-  PUGS_INLINE
-  auto
-  operator,(const Vector<DataType2>& y)
+  PUGS_INLINE friend auto
+  dot(const Vector& x, const Vector<DataType2>& y)
   {
-    Assert(this->size() == y.size());
+    Assert(x.size() == y.size());
     // Quite ugly, TODO: fix me in C++20
     auto promoted = [] {
       DataType a{0};
@@ -55,8 +55,8 @@ class Vector   // LCOV_EXCL_LINE
     decltype(promoted) sum = 0;
 
     // Does not use parallel_for to preserve sum order
-    for (index_type i = 0; i < this->size(); ++i) {
-      sum += m_values[i] * y.m_values[i];
+    for (index_type i = 0; i < x.size(); ++i) {
+      sum += x.m_values[i] * y.m_values[i];
     }
 
     return sum;
@@ -130,7 +130,8 @@ class Vector   // LCOV_EXCL_LINE
   }
 
   PUGS_INLINE
-  DataType& operator[](index_type i) const noexcept(NO_ASSERT)
+  DataType&
+  operator[](index_type i) const noexcept(NO_ASSERT)
   {
     return m_values[i];
   }
diff --git a/src/language/modules/MathModule.cpp b/src/language/modules/MathModule.cpp
index fd38710a5b2232e12b54d69217880188ab292fc9..69f6890fbc5e8c3e3515901acd5b5102100e8c35 100644
--- a/src/language/modules/MathModule.cpp
+++ b/src/language/modules/MathModule.cpp
@@ -72,16 +72,16 @@ MathModule::MathModule()
 
   this->_addBuiltinFunction("dot",
                             std::make_shared<BuiltinFunctionEmbedder<double(const TinyVector<1>, const TinyVector<1>)>>(
-                              [](const TinyVector<1> x, const TinyVector<1> y) -> double { return (x, y); }));
+                              [](const TinyVector<1> x, const TinyVector<1> y) -> double { return dot(x, y); }));
 
   this->_addBuiltinFunction("dot",
                             std::make_shared<BuiltinFunctionEmbedder<double(const TinyVector<2>, const TinyVector<2>)>>(
-                              [](const TinyVector<2> x, const TinyVector<2> y) -> double { return (x, y); }));
+                              [](const TinyVector<2> x, const TinyVector<2> y) -> double { return dot(x, y); }));
 
   this
     ->_addBuiltinFunction("dot",
                           std::make_shared<BuiltinFunctionEmbedder<double(const TinyVector<3>&, const TinyVector<3>&)>>(
-                            [](const TinyVector<3>& x, const TinyVector<3>& y) -> double { return (x, y); }));
+                            [](const TinyVector<3>& x, const TinyVector<3>& y) -> double { return dot(x, y); }));
 }
 
 void
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index 9866c9ad8389c509565ba26c8e94915c3f597931..3a83ef5c3f49d068fbfd8e62d4ff002df50a8845 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -1233,12 +1233,12 @@ GmshReader::__proceedData()
   dimension3_mask[12]                 = 1;
   dimension3_mask[13]                 = 1;
 
-  std::cout << "- dimension 0 entities: " << (dimension0_mask, elementNumber) << '\n';
-  std::cout << "- dimension 1 entities: " << (dimension1_mask, elementNumber) << '\n';
-  std::cout << "- dimension 2 entities: " << (dimension2_mask, elementNumber) << '\n';
-  std::cout << "- dimension 3 entities: " << (dimension3_mask, elementNumber) << '\n';
-  if ((dimension3_mask, elementNumber) > 0) {
-    const size_t nb_cells = (dimension3_mask, elementNumber);
+  std::cout << "- dimension 0 entities: " << dot(dimension0_mask, elementNumber) << '\n';
+  std::cout << "- dimension 1 entities: " << dot(dimension1_mask, elementNumber) << '\n';
+  std::cout << "- dimension 2 entities: " << dot(dimension2_mask, elementNumber) << '\n';
+  std::cout << "- dimension 3 entities: " << dot(dimension3_mask, elementNumber) << '\n';
+  if (dot(dimension3_mask, elementNumber) > 0) {
+    const size_t nb_cells = dot(dimension3_mask, elementNumber);
 
     GmshConnectivityBuilder<3> connectivity_builder(m_mesh_data, nb_cells);
 
@@ -1257,8 +1257,8 @@ GmshReader::__proceedData()
     }
     m_mesh = std::make_shared<MeshType>(p_connectivity, xr);
 
-  } else if ((dimension2_mask, elementNumber) > 0) {
-    const size_t nb_cells = (dimension2_mask, elementNumber);
+  } else if (dot(dimension2_mask, elementNumber) > 0) {
+    const size_t nb_cells = dot(dimension2_mask, elementNumber);
 
     GmshConnectivityBuilder<2> connectivity_builder(m_mesh_data, nb_cells);
 
@@ -1276,8 +1276,8 @@ GmshReader::__proceedData()
     }
     m_mesh = std::make_shared<MeshType>(p_connectivity, xr);
 
-  } else if ((dimension1_mask, elementNumber) > 0) {
-    const size_t nb_cells = (dimension1_mask, elementNumber);
+  } else if (dot(dimension1_mask, elementNumber) > 0) {
+    const size_t nb_cells = dot(dimension1_mask, elementNumber);
 
     GmshConnectivityBuilder<1> connectivity_builder(m_mesh_data, nb_cells);
 
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index f891f5ccb40e92cfbd0aa9e4cf3a600f352a8790..c15bb1a66ddacc34fe392a08d133e6f683b0e757 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -265,7 +265,7 @@ class MeshData : public IMeshData
                 const Rd xjxr0 = xr0 - xj;
                 const Rd xjxr1 = xr1 - xj;
 
-                const double Vjlr = (crossProduct(xjxr0, xjxr1), xjxl);
+                const double Vjlr = dot(crossProduct(xjxr0, xjxr1), xjxl);
 
                 sum += (sign * Vjlr) * (xl_plus_xj + xr0 + xr1);
               }
@@ -297,7 +297,7 @@ class MeshData : public IMeshData
         const auto& cell_nodes = cell_to_node_matrix[j];
 
         for (size_t R = 0; R < cell_nodes.size(); ++R) {
-          sum_cjr_xr += (xr[cell_nodes[R]], Cjr(j, R));
+          sum_cjr_xr += dot(xr[cell_nodes[R]], Cjr(j, R));
         }
         Vj[j] = inv_Dimension * sum_cjr_xr;
       });
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index a49db8b4b358f52e7dd199a64a55a43b1123de2e..92d354fc52ddfd3e35ab69a8d892edb32f6432de 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -157,7 +157,7 @@ MeshFlatNodeBoundary<2>::_checkBoundaryIsFlat(TinyVector<2, double> normal,
   parallel_for(
     m_node_list.size(), PUGS_LAMBDA(size_t r) {
       const R2& x = xr[m_node_list[r]];
-      if ((x - origin, normal) > 1E-13 * length) {
+      if (dot(x - origin, normal) > 1E-13 * length) {
         throw NormalError("this FlatBoundary is not flat!");
       }
     });
@@ -330,13 +330,13 @@ MeshFlatNodeBoundary<3>::_getNormal(const std::shared_ptr<const MeshType>& mesh)
   const R3 w = zmax - zmin;
 
   const R3 uv        = crossProduct(u, v);
-  const double uv_l2 = (uv, uv);
+  const double uv_l2 = dot(uv, uv);
 
   R3 normal        = uv;
   double normal_l2 = uv_l2;
 
   const R3 uw        = crossProduct(u, w);
-  const double uw_l2 = (uw, uw);
+  const double uw_l2 = dot(uw, uw);
 
   if (uw_l2 > uv_l2) {
     normal    = uw;
@@ -344,7 +344,7 @@ MeshFlatNodeBoundary<3>::_getNormal(const std::shared_ptr<const MeshType>& mesh)
   }
 
   const R3 vw        = crossProduct(v, w);
-  const double vw_l2 = (vw, vw);
+  const double vw_l2 = dot(vw, vw);
 
   if (vw_l2 > normal_l2) {
     normal    = vw;
@@ -385,7 +385,7 @@ MeshFlatNodeBoundary<1>::_getOutgoingNormal(const std::shared_ptr<const MeshType
     const auto& j0_nodes = cell_to_node_matrix[j0];
 
     for (size_t r = 0; r < j0_nodes.size(); ++r) {
-      const double height = (xr[j0_nodes[r]] - xr[r0], normal);
+      const double height = dot(xr[j0_nodes[r]] - xr[r0], normal);
       if (std::abs(height) > std::abs(max_height)) {
         max_height = height;
       }
@@ -429,7 +429,7 @@ MeshFlatNodeBoundary<2>::_getOutgoingNormal(const std::shared_ptr<const MeshType
     const CellId j0      = node_to_cell_matrix[r0][0];
     const auto& j0_nodes = cell_to_node_matrix[j0];
     for (size_t r = 0; r < j0_nodes.size(); ++r) {
-      const double height = (xr[j0_nodes[r]] - xr[r0], normal);
+      const double height = dot(xr[j0_nodes[r]] - xr[r0], normal);
       if (std::abs(height) > std::abs(max_height)) {
         max_height = height;
       }
@@ -474,7 +474,7 @@ MeshFlatNodeBoundary<3>::_getOutgoingNormal(const std::shared_ptr<const MeshType
     const auto& j0_nodes = cell_to_node_matrix[j0];
 
     for (size_t r = 0; r < j0_nodes.size(); ++r) {
-      const double height = (xr[j0_nodes[r]] - xr[r0], normal);
+      const double height = dot(xr[j0_nodes[r]] - xr[r0], normal);
       if (std::abs(height) > std::abs(max_height)) {
         max_height = height;
       }
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index 6bac5a6b3d920f9f3b5654e826bf9b974a9fcd0f..6973037571b23f861cec5c70ed28fbe67c036fcf 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -453,7 +453,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
         for (size_t R = 0; R < cell_nodes.size(); ++R) {
           const NodeId r = cell_nodes[R];
           momentum_fluxes += m_Fjr(j, R);
-          energy_fluxes += (m_Fjr(j, R), m_ur[r]);
+          energy_fluxes += dot(m_Fjr(j, R), m_ur[r]);
         }
         const double dt_over_Mj = dt / (rho[j] * Vj[j]);
         new_u[j] -= dt_over_Mj * momentum_fluxes;
diff --git a/src/scheme/PerfectGas.cpp b/src/scheme/PerfectGas.cpp
index b371dea3c9e4aac8ee867dd0f938b72bcd02791e..56cca334dce784737bf1e47ad537fbb7c74f1f89 100644
--- a/src/scheme/PerfectGas.cpp
+++ b/src/scheme/PerfectGas.cpp
@@ -167,7 +167,7 @@ totalEnergyFromEpsilonU(const DiscreteFunctionP0<Dimension, DataType>& epsilon,
 
   parallel_for(
     mesh->numberOfCells(),
-    PUGS_LAMBDA(CellId cell_id) { E[cell_id] = epsilon[cell_id] + 0.5 * (u[cell_id], u[cell_id]); });
+    PUGS_LAMBDA(CellId cell_id) { E[cell_id] = epsilon[cell_id] + 0.5 * dot(u[cell_id], u[cell_id]); });
 
   return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(epsilon);
 }
@@ -212,7 +212,7 @@ epsilonFromTotalEnergyU(const DiscreteFunctionP0<Dimension, DataType>& E,
 
   parallel_for(
     mesh->numberOfCells(),
-    PUGS_LAMBDA(CellId cell_id) { epsilon[cell_id] = E[cell_id] - 0.5 * (u[cell_id], u[cell_id]); });
+    PUGS_LAMBDA(CellId cell_id) { epsilon[cell_id] = E[cell_id] - 0.5 * dot(u[cell_id], u[cell_id]); });
 
   return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(epsilon);
 }
diff --git a/tests/test_BiCGStab.cpp b/tests/test_BiCGStab.cpp
index ac6b547a82d532e0edf2e5c57d68c41ae97afeee..e1be8069885b80a3e26ede233dac3941b7525521 100644
--- a/tests/test_BiCGStab.cpp
+++ b/tests/test_BiCGStab.cpp
@@ -48,7 +48,7 @@ TEST_CASE("BiCGStab", "[algebra]")
 
     BiCGStab{A, x, b, 1e-12, 10, false};
     Vector error = x - x_exact;
-    REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x, x)));
+    REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x, x)));
   }
 
   SECTION("no preconditionner non-converged")
@@ -91,6 +91,6 @@ TEST_CASE("BiCGStab", "[algebra]")
 
     BiCGStab{A, x, b, 1e-12, 1, true};
     Vector error = x - x_exact;
-    REQUIRE(std::sqrt((error, error)) > 1E-5 * std::sqrt((x, x)));
+    REQUIRE(std::sqrt(dot(error, error)) > 1E-5 * std::sqrt(dot(x, x)));
   }
 }
diff --git a/tests/test_BuiltinFunctionProcessor.cpp b/tests/test_BuiltinFunctionProcessor.cpp
index c36399857d5a9b2235111d1ca105559001edf69d..a1fb3fe82b982bc50d2b5247c1158941da24ce0d 100644
--- a/tests/test_BuiltinFunctionProcessor.cpp
+++ b/tests/test_BuiltinFunctionProcessor.cpp
@@ -297,7 +297,7 @@ let x:R^2, x = (-2, 3);
 let y:R^2, y = (4, 3);
 let s:R, s = dot(x,y);
 )";
-      CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "s", (TinyVector<2>{-2, 3}, TinyVector<2>{4, 3}));
+      CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "s", dot(TinyVector<2>{-2, 3}, TinyVector<2>{4, 3}));
     }
 
     {   // dot
@@ -308,7 +308,7 @@ let x:R^3, x = (-2, 3, 4);
 let y:R^3, y = (4, 3, 5);
 let s:R, s = dot(x,y);
 )";
-      CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "s", (TinyVector<3>{-2, 3, 4}, TinyVector<3>{4, 3, 5}));
+      CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "s", dot(TinyVector<3>{-2, 3, 4}, TinyVector<3>{4, 3, 5}));
     }
 
     MathModule math_module;
diff --git a/tests/test_CG.cpp b/tests/test_CG.cpp
index 00751efb3ccbdf5fe86da3a242c0b6dad92b097b..05919858f3ada681eb371883c6fcf6fd7d15a724 100644
--- a/tests/test_CG.cpp
+++ b/tests/test_CG.cpp
@@ -48,7 +48,7 @@ TEST_CASE("CG", "[algebra]")
 
     CG{A, x, b, 1e-12, 10, true};
     Vector error = x - x_exact;
-    REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x, x)));
+    REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x, x)));
   }
 
   SECTION("singular matrix with RHS in its image")
@@ -64,7 +64,7 @@ TEST_CASE("CG", "[algebra]")
     x = 0;
 
     CG{A, x, b, 1e-12, 10};
-    REQUIRE(std::sqrt((x, x)) == 0);
+    REQUIRE(std::sqrt(dot(x, x)) == 0);
   }
 
   SECTION("no preconditionner non-converged")
@@ -107,6 +107,6 @@ TEST_CASE("CG", "[algebra]")
 
     CG{A, x, b, 1e-12, 1, false};
     Vector error = x - x_exact;
-    REQUIRE(std::sqrt((error, error)) > 1E-10 * std::sqrt((x, x)));
+    REQUIRE(std::sqrt(dot(error, error)) > 1E-10 * std::sqrt(dot(x, x)));
   }
 }
diff --git a/tests/test_LinearSolver.cpp b/tests/test_LinearSolver.cpp
index 9ced9e9424a7d1c80a41cd9dd5a807476adea9d7..5504ae390fd259eaa700a911b158bf6a6e46a625 100644
--- a/tests/test_LinearSolver.cpp
+++ b/tests/test_LinearSolver.cpp
@@ -194,7 +194,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
           solver.solveLocalSystem(A, x, b);
           Vector error = x - x_exact;
-          REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+          REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
         }
       }
 
@@ -221,7 +221,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
             solver.solveLocalSystem(A, x, b);
             Vector error = x - x_exact;
-            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
 
           SECTION("CG Diagonal")
@@ -235,7 +235,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
             solver.solveLocalSystem(A, x, b);
             Vector error = x - x_exact;
-            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
 
           SECTION("CG ICholeski")
@@ -249,7 +249,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
             solver.solveLocalSystem(A, x, b);
             Vector error = x - x_exact;
-            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
 
           SECTION("CG AMG")
@@ -263,7 +263,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
             solver.solveLocalSystem(A, x, b);
             Vector error = x - x_exact;
-            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
         }
 
@@ -281,7 +281,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
           solver.solveLocalSystem(A, x, b);
           Vector error = x - x_exact;
-          REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+          REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
         }
 
 #else    // PUGS_HAS_PETSC
@@ -349,7 +349,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
           solver.solveLocalSystem(A, x, b);
           Vector error = x - x_exact;
-          REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+          REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
         }
       }
 
@@ -376,7 +376,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
             solver.solveLocalSystem(A, x, b);
             Vector error = x - x_exact;
-            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
 
           SECTION("BICGStab Diagonal")
@@ -390,7 +390,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
             solver.solveLocalSystem(A, x, b);
             Vector error = x - x_exact;
-            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
 
           SECTION("BICGStab ILU")
@@ -404,7 +404,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
             solver.solveLocalSystem(A, x, b);
             Vector error = x - x_exact;
-            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
         }
 
@@ -426,7 +426,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
             solver.solveLocalSystem(A, x, b);
             Vector error = x - x_exact;
-            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
 
           SECTION("BICGStab2 Diagonal")
@@ -440,7 +440,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
             solver.solveLocalSystem(A, x, b);
             Vector error = x - x_exact;
-            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
         }
 
@@ -462,7 +462,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
             solver.solveLocalSystem(A, x, b);
             Vector error = x - x_exact;
-            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
 
           SECTION("GMRES Diagonal")
@@ -476,7 +476,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
             solver.solveLocalSystem(A, x, b);
             Vector error = x - x_exact;
-            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
 
           SECTION("GMRES ILU")
@@ -490,7 +490,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
             solver.solveLocalSystem(A, x, b);
             Vector error = x - x_exact;
-            REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+            REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
           }
         }
 
@@ -508,7 +508,7 @@ TEST_CASE("LinearSolver", "[algebra]")
 
           solver.solveLocalSystem(A, x, b);
           Vector error = x - x_exact;
-          REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x_exact, x_exact)));
+          REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
         }
 
 #else    // PUGS_HAS_PETSC
diff --git a/tests/test_MathModule.cpp b/tests/test_MathModule.cpp
index fe87699887d6d587094bcd33999d403a20ce10b0..7c2a9bc5438367b7f3cf70d14d3b8acc58a29ae7 100644
--- a/tests/test_MathModule.cpp
+++ b/tests/test_MathModule.cpp
@@ -341,8 +341,8 @@ TEST_CASE("MathModule", "[language]")
       IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
       DataVariant result_variant                  = function_embedder.apply({arg0_variant, arg1_variant});
 
-      auto result = (arg0, arg1);
-      REQUIRE(std::get<decltype(result)>(result_variant) == result);
+      const double result = dot(arg0, arg1);
+      REQUIRE(std::get<double>(result_variant) == result);
     }
 
     SECTION("dot:R^2*R^2")
@@ -359,8 +359,8 @@ TEST_CASE("MathModule", "[language]")
       IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
       DataVariant result_variant                  = function_embedder.apply({arg0_variant, arg1_variant});
 
-      auto result = (arg0, arg1);
-      REQUIRE(std::get<decltype(result)>(result_variant) == result);
+      const double result = dot(arg0, arg1);
+      REQUIRE(std::get<double>(result_variant) == result);
     }
 
     SECTION("dot:R^3*R^3")
@@ -377,8 +377,8 @@ TEST_CASE("MathModule", "[language]")
       IBuiltinFunctionEmbedder& function_embedder = *i_function->second;
       DataVariant result_variant                  = function_embedder.apply({arg0_variant, arg1_variant});
 
-      auto result = (arg0, arg1);
-      REQUIRE(std::get<decltype(result)>(result_variant) == result);
+      const double result = dot(arg0, arg1);
+      REQUIRE(std::get<double>(result_variant) == result);
     }
   }
 }
diff --git a/tests/test_TinyVector.cpp b/tests/test_TinyVector.cpp
index 2775588af687aa274dbba5374af53a85cacce36a..499536f9c702d5c3a1f2fcf8c4da7b1ecb461038 100644
--- a/tests/test_TinyVector.cpp
+++ b/tests/test_TinyVector.cpp
@@ -26,7 +26,7 @@ TEST_CASE("TinyVector", "[algebra]")
   REQUIRE(Catch::Detail::stringify(v) == "(3,2,4)");
 
   TinyVector<3, int> w(1, 2, 6);
-  REQUIRE((v, w) == 31);
+  REQUIRE(dot(v, w) == 31);
 
   w = 2 * v;
   REQUIRE(w == TinyVector<3, int>(6, 4, 8));
diff --git a/tests/test_Vector.cpp b/tests/test_Vector.cpp
index d02481859e1c41c49509f973f09006ba50373cc3..89fedcec868e026cfdf569eb8667b57b730a6e16 100644
--- a/tests/test_Vector.cpp
+++ b/tests/test_Vector.cpp
@@ -145,8 +145,8 @@ TEST_CASE("Vector", "[algebra]")
     y[3] = 2;
     y[4] = 8;
 
-    const int dot = (x, y);
-    REQUIRE(dot == 49);
+    const int s = dot(x, y);
+    REQUIRE(s == 49);
   }
 
   SECTION("self scalar division")