diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/abs.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/abs.hpp
index 9c42ba1bb4..f96b31c137 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/abs.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/abs.hpp
@@ -61,7 +61,7 @@ template <typename argT, typename resT> struct AbsFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &x)
+    resT operator()(const argT &x) const
     {
 
         if constexpr (std::is_same_v<argT, bool> ||
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/acos.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/acos.hpp
index 5dcd1268b5..c091005f42 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/acos.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/acos.hpp
@@ -63,7 +63,7 @@ template <typename argT, typename resT> struct AcosFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         if constexpr (is_complex<argT>::value) {
             using realT = typename argT::value_type;
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/acosh.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/acosh.hpp
index 1dca4d36d6..9601f6b3a7 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/acosh.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/acosh.hpp
@@ -63,7 +63,7 @@ template <typename argT, typename resT> struct AcoshFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         if constexpr (is_complex<argT>::value) {
             using realT = typename argT::value_type;
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/add.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/add.hpp
index 0b675c9685..18e946a7d0 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/add.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/add.hpp
@@ -58,14 +58,15 @@ template <typename argT1, typename argT2, typename resT> struct AddFunctor
     using supports_vec = std::negation<
         std::disjunction<tu_ns::is_complex<argT1>, tu_ns::is_complex<argT2>>>;
 
-    resT operator()(const argT1 &in1, const argT2 &in2)
+    resT operator()(const argT1 &in1, const argT2 &in2) const
     {
         return in1 + in2;
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT1, vec_sz> &in1,
-                                       const sycl::vec<argT2, vec_sz> &in2)
+    sycl::vec<resT, vec_sz>
+    operator()(const sycl::vec<argT1, vec_sz> &in1,
+               const sycl::vec<argT2, vec_sz> &in2) const
     {
         auto tmp = in1 + in2;
         if constexpr (std::is_same_v<resT,
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/asin.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/asin.hpp
index a150b8e2e1..f559171b54 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/asin.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/asin.hpp
@@ -63,7 +63,7 @@ template <typename argT, typename resT> struct AsinFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         if constexpr (is_complex<argT>::value) {
             using realT = typename argT::value_type;
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/asinh.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/asinh.hpp
index a42ea02598..729d364abf 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/asinh.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/asinh.hpp
@@ -63,7 +63,7 @@ template <typename argT, typename resT> struct AsinhFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         if constexpr (is_complex<argT>::value) {
             using realT = typename argT::value_type;
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/atan.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/atan.hpp
index cc988bd56d..8f26d0cb42 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/atan.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/atan.hpp
@@ -64,7 +64,7 @@ template <typename argT, typename resT> struct AtanFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         if constexpr (is_complex<argT>::value) {
             using realT = typename argT::value_type;
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/atan2.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/atan2.hpp
index 613e2bff85..fdc9c4e66e 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/atan2.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/atan2.hpp
@@ -55,7 +55,7 @@ template <typename argT1, typename argT2, typename resT> struct Atan2Functor
     using supports_sg_loadstore = std::true_type;
     using supports_vec = std::false_type;
 
-    resT operator()(const argT1 &in1, const argT2 &in2)
+    resT operator()(const argT1 &in1, const argT2 &in2) const
     {
         if (std::isinf(in2) && !std::signbit(in2)) {
             if (std::isfinite(in1)) {
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/atanh.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/atanh.hpp
index 4842f785a2..6eae85f7b7 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/atanh.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/atanh.hpp
@@ -64,7 +64,7 @@ template <typename argT, typename resT> struct AtanhFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         if constexpr (is_complex<argT>::value) {
             using realT = typename argT::value_type;
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/bitwise_and.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/bitwise_and.hpp
index 12dece02bd..24a547eacf 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/bitwise_and.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/bitwise_and.hpp
@@ -57,7 +57,7 @@ struct BitwiseAndFunctor
     using supports_sg_loadstore = typename std::true_type;
     using supports_vec = typename std::true_type;
 
-    resT operator()(const argT1 &in1, const argT2 &in2)
+    resT operator()(const argT1 &in1, const argT2 &in2) const
     {
         using tu_ns::convert_impl;
 
@@ -70,8 +70,9 @@ struct BitwiseAndFunctor
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT1, vec_sz> &in1,
-                                       const sycl::vec<argT2, vec_sz> &in2)
+    sycl::vec<resT, vec_sz>
+    operator()(const sycl::vec<argT1, vec_sz> &in1,
+               const sycl::vec<argT2, vec_sz> &in2) const
     {
 
         if constexpr (std::is_same_v<resT, bool>) {
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/bitwise_invert.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/bitwise_invert.hpp
index 44ac214b27..f258af34e2 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/bitwise_invert.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/bitwise_invert.hpp
@@ -72,7 +72,7 @@ template <typename argT, typename resT> struct BitwiseInvertFunctor
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT, vec_sz> &in)
+    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT, vec_sz> &in) const
     {
         if constexpr (std::is_same_v<argT, bool>) {
             auto res_vec = !in;
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/bitwise_left_shift.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/bitwise_left_shift.hpp
index 96deccc0a3..6e12d80246 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/bitwise_left_shift.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/bitwise_left_shift.hpp
@@ -66,8 +66,9 @@ struct BitwiseLeftShiftFunctor
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT1, vec_sz> &in1,
-                                       const sycl::vec<argT2, vec_sz> &in2)
+    sycl::vec<resT, vec_sz>
+    operator()(const sycl::vec<argT1, vec_sz> &in1,
+               const sycl::vec<argT2, vec_sz> &in2) const
     {
         sycl::vec<resT, vec_sz> res;
 #pragma unroll
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/bitwise_or.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/bitwise_or.hpp
index 47d34ac182..fd1f3b357c 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/bitwise_or.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/bitwise_or.hpp
@@ -56,7 +56,7 @@ template <typename argT1, typename argT2, typename resT> struct BitwiseOrFunctor
     using supports_sg_loadstore = typename std::true_type;
     using supports_vec = typename std::true_type;
 
-    resT operator()(const argT1 &in1, const argT2 &in2)
+    resT operator()(const argT1 &in1, const argT2 &in2) const
     {
         using tu_ns::convert_impl;
 
@@ -69,8 +69,9 @@ template <typename argT1, typename argT2, typename resT> struct BitwiseOrFunctor
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT1, vec_sz> &in1,
-                                       const sycl::vec<argT2, vec_sz> &in2)
+    sycl::vec<resT, vec_sz>
+    operator()(const sycl::vec<argT1, vec_sz> &in1,
+               const sycl::vec<argT2, vec_sz> &in2) const
     {
 
         if constexpr (std::is_same_v<resT, bool>) {
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/bitwise_xor.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/bitwise_xor.hpp
index 2e5568fe8c..9569cd53cd 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/bitwise_xor.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/bitwise_xor.hpp
@@ -57,7 +57,7 @@ struct BitwiseXorFunctor
     using supports_sg_loadstore = typename std::true_type;
     using supports_vec = typename std::true_type;
 
-    resT operator()(const argT1 &in1, const argT2 &in2)
+    resT operator()(const argT1 &in1, const argT2 &in2) const
     {
         if constexpr (std::is_same_v<resT, bool>) {
             // (false != false) -> false, (false != true) -> true
@@ -70,8 +70,9 @@ struct BitwiseXorFunctor
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT1, vec_sz> &in1,
-                                       const sycl::vec<argT2, vec_sz> &in2)
+    sycl::vec<resT, vec_sz>
+    operator()(const sycl::vec<argT1, vec_sz> &in1,
+               const sycl::vec<argT2, vec_sz> &in2) const
     {
 
         if constexpr (std::is_same_v<resT, bool>) {
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/ceil.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/ceil.hpp
index cc6b1cd534..7cabf9b7e3 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/ceil.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/ceil.hpp
@@ -63,7 +63,7 @@ template <typename argT, typename resT> struct CeilFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         if constexpr (std::is_integral_v<argT>) {
             return in;
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/conj.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/conj.hpp
index 24c5a128d0..7000d83ff1 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/conj.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/conj.hpp
@@ -65,7 +65,7 @@ template <typename argT, typename resT> struct ConjFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         if constexpr (is_complex<argT>::value) {
             return std::conj(in);
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/cos.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/cos.hpp
index 313fb565e8..4988cbdce0 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/cos.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/cos.hpp
@@ -63,7 +63,7 @@ template <typename argT, typename resT> struct CosFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         if constexpr (is_complex<argT>::value) {
             using realT = typename argT::value_type;
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/cosh.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/cosh.hpp
index 317cd101af..1fb6bcc0c4 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/cosh.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/cosh.hpp
@@ -63,7 +63,7 @@ template <typename argT, typename resT> struct CoshFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         if constexpr (is_complex<argT>::value) {
             using realT = typename argT::value_type;
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/equal.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/equal.hpp
index 46c4e58f0e..7418bd1de9 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/equal.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/equal.hpp
@@ -60,14 +60,15 @@ template <typename argT1, typename argT2, typename resT> struct EqualFunctor
         std::negation<std::disjunction<tu_ns::is_complex<argT1>,
                                        tu_ns::is_complex<argT2>>>>;
 
-    resT operator()(const argT1 &in1, const argT2 &in2)
+    resT operator()(const argT1 &in1, const argT2 &in2) const
     {
         return (in1 == in2);
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT1, vec_sz> &in1,
-                                       const sycl::vec<argT2, vec_sz> &in2)
+    sycl::vec<resT, vec_sz>
+    operator()(const sycl::vec<argT1, vec_sz> &in1,
+               const sycl::vec<argT2, vec_sz> &in2) const
     {
         auto tmp = (in1 == in2);
         if constexpr (std::is_same_v<resT,
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/exp.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/exp.hpp
index 2f7657f5b4..9a02402381 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/exp.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/exp.hpp
@@ -62,7 +62,7 @@ template <typename argT, typename resT> struct ExpFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         if constexpr (is_complex<argT>::value) {
             using realT = typename argT::value_type;
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/expm1.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/expm1.hpp
index e1c23113c6..b4bfecb22f 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/expm1.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/expm1.hpp
@@ -65,7 +65,7 @@ template <typename argT, typename resT> struct Expm1Functor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         if constexpr (is_complex<argT>::value) {
             using realT = typename argT::value_type;
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/floor.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/floor.hpp
index db6b669a2c..9c16a805b3 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/floor.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/floor.hpp
@@ -63,7 +63,7 @@ template <typename argT, typename resT> struct FloorFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         if constexpr (std::is_integral_v<argT>) {
             return in;
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/floor_divide.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/floor_divide.hpp
index 32e97df58d..e8d4f524fd 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/floor_divide.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/floor_divide.hpp
@@ -55,7 +55,7 @@ struct FloorDivideFunctor
     using supports_sg_loadstore = std::true_type;
     using supports_vec = std::true_type;
 
-    resT operator()(const argT1 &in1, const argT2 &in2)
+    resT operator()(const argT1 &in1, const argT2 &in2) const
     {
         if constexpr (std::is_same_v<argT1, bool> &&
                       std::is_same_v<argT2, bool>) {
@@ -83,8 +83,9 @@ struct FloorDivideFunctor
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT1, vec_sz> &in1,
-                                       const sycl::vec<argT2, vec_sz> &in2)
+    sycl::vec<resT, vec_sz>
+    operator()(const sycl::vec<argT1, vec_sz> &in1,
+               const sycl::vec<argT2, vec_sz> &in2) const
     {
         if constexpr (std::is_same_v<argT1, bool> &&
                       std::is_same_v<argT2, bool>) {
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/greater.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/greater.hpp
index ab2ce9b5c3..23186b7633 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/greater.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/greater.hpp
@@ -62,7 +62,7 @@ template <typename argT1, typename argT2, typename resT> struct GreaterFunctor
         std::negation<std::disjunction<tu_ns::is_complex<argT1>,
                                        tu_ns::is_complex<argT2>>>>;
 
-    resT operator()(const argT1 &in1, const argT2 &in2)
+    resT operator()(const argT1 &in1, const argT2 &in2) const
     {
         if constexpr (tu_ns::is_complex<argT1>::value ||
                       tu_ns::is_complex<argT2>::value)
@@ -77,8 +77,9 @@ template <typename argT1, typename argT2, typename resT> struct GreaterFunctor
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT1, vec_sz> &in1,
-                                       const sycl::vec<argT2, vec_sz> &in2)
+    sycl::vec<resT, vec_sz>
+    operator()(const sycl::vec<argT1, vec_sz> &in1,
+               const sycl::vec<argT2, vec_sz> &in2) const
     {
 
         auto tmp = (in1 > in2);
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/greater_equal.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/greater_equal.hpp
index fa2ece210f..92b7189fc7 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/greater_equal.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/greater_equal.hpp
@@ -63,7 +63,7 @@ struct GreaterEqualFunctor
         std::negation<std::disjunction<tu_ns::is_complex<argT1>,
                                        tu_ns::is_complex<argT2>>>>;
 
-    resT operator()(const argT1 &in1, const argT2 &in2)
+    resT operator()(const argT1 &in1, const argT2 &in2) const
     {
         if constexpr (tu_ns::is_complex<argT1>::value ||
                       tu_ns::is_complex<argT2>::value)
@@ -78,8 +78,9 @@ struct GreaterEqualFunctor
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT1, vec_sz> &in1,
-                                       const sycl::vec<argT2, vec_sz> &in2)
+    sycl::vec<resT, vec_sz>
+    operator()(const sycl::vec<argT1, vec_sz> &in1,
+               const sycl::vec<argT2, vec_sz> &in2) const
     {
 
         auto tmp = (in1 >= in2);
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/hypot.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/hypot.hpp
index c7f50deb75..a12b84e4c6 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/hypot.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/hypot.hpp
@@ -57,14 +57,15 @@ template <typename argT1, typename argT2, typename resT> struct HypotFunctor
     using supports_vec = std::negation<
         std::disjunction<tu_ns::is_complex<argT1>, tu_ns::is_complex<argT2>>>;
 
-    resT operator()(const argT1 &in1, const argT2 &in2)
+    resT operator()(const argT1 &in1, const argT2 &in2) const
     {
         return std::hypot(in1, in2);
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT1, vec_sz> &in1,
-                                       const sycl::vec<argT2, vec_sz> &in2)
+    sycl::vec<resT, vec_sz>
+    operator()(const sycl::vec<argT1, vec_sz> &in1,
+               const sycl::vec<argT2, vec_sz> &in2) const
     {
         auto res = sycl::hypot(in1, in2);
         if constexpr (std::is_same_v<resT,
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/imag.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/imag.hpp
index 6c85de5561..b094355ad0 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/imag.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/imag.hpp
@@ -65,7 +65,7 @@ template <typename argT, typename resT> struct ImagFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         if constexpr (is_complex<argT>::value) {
             return std::imag(in);
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/isfinite.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/isfinite.hpp
index cc91bdbb97..b075521fdb 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/isfinite.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/isfinite.hpp
@@ -86,7 +86,7 @@ template <typename argT, typename resT> struct IsFiniteFunctor
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT, vec_sz> &in)
+    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT, vec_sz> &in) const
     {
         auto const &res_vec = sycl::isfinite(in);
 
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/isinf.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/isinf.hpp
index 8cb93ea456..22a83b1128 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/isinf.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/isinf.hpp
@@ -84,7 +84,7 @@ template <typename argT, typename resT> struct IsInfFunctor
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT, vec_sz> &in)
+    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT, vec_sz> &in) const
     {
         auto const &res_vec = sycl::isinf(in);
 
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/isnan.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/isnan.hpp
index 5d628437ed..da5e08a9ce 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/isnan.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/isnan.hpp
@@ -82,7 +82,7 @@ template <typename argT, typename resT> struct IsNanFunctor
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT, vec_sz> &in)
+    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT, vec_sz> &in) const
     {
         auto const &res_vec = sycl::isnan(in);
 
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/less.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/less.hpp
index 36675a2e19..3562ee08d8 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/less.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/less.hpp
@@ -61,7 +61,7 @@ template <typename argT1, typename argT2, typename resT> struct LessFunctor
         std::negation<std::disjunction<tu_ns::is_complex<argT1>,
                                        tu_ns::is_complex<argT2>>>>;
 
-    resT operator()(const argT1 &in1, const argT2 &in2)
+    resT operator()(const argT1 &in1, const argT2 &in2) const
     {
         if constexpr (tu_ns::is_complex<argT1>::value ||
                       tu_ns::is_complex<argT2>::value)
@@ -76,8 +76,9 @@ template <typename argT1, typename argT2, typename resT> struct LessFunctor
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT1, vec_sz> &in1,
-                                       const sycl::vec<argT2, vec_sz> &in2)
+    sycl::vec<resT, vec_sz>
+    operator()(const sycl::vec<argT1, vec_sz> &in1,
+               const sycl::vec<argT2, vec_sz> &in2) const
     {
 
         auto tmp = (in1 < in2);
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/less_equal.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/less_equal.hpp
index 9b740219fc..a316f8abf3 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/less_equal.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/less_equal.hpp
@@ -62,7 +62,7 @@ template <typename argT1, typename argT2, typename resT> struct LessEqualFunctor
         std::negation<std::disjunction<tu_ns::is_complex<argT1>,
                                        tu_ns::is_complex<argT2>>>>;
 
-    resT operator()(const argT1 &in1, const argT2 &in2)
+    resT operator()(const argT1 &in1, const argT2 &in2) const
     {
         if constexpr (tu_ns::is_complex<argT1>::value ||
                       tu_ns::is_complex<argT2>::value)
@@ -77,8 +77,9 @@ template <typename argT1, typename argT2, typename resT> struct LessEqualFunctor
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT1, vec_sz> &in1,
-                                       const sycl::vec<argT2, vec_sz> &in2)
+    sycl::vec<resT, vec_sz>
+    operator()(const sycl::vec<argT1, vec_sz> &in1,
+               const sycl::vec<argT2, vec_sz> &in2) const
     {
 
         auto tmp = (in1 <= in2);
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/log.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/log.hpp
index f6998d4e0f..9fd366ba39 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/log.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/log.hpp
@@ -63,7 +63,7 @@ template <typename argT, typename resT> struct LogFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         return std::log(in);
     }
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/log10.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/log10.hpp
index 6b1185416f..8cb723f2b9 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/log10.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/log10.hpp
@@ -66,7 +66,7 @@ template <typename argT, typename resT> struct Log10Functor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         if constexpr (is_complex<argT>::value) {
             using realT = typename argT::value_type;
@@ -78,7 +78,7 @@ template <typename argT, typename resT> struct Log10Functor
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT, vec_sz> &in)
+    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT, vec_sz> &in) const
     {
         auto const &res_vec = sycl::log10(in);
         using deducedT = typename std::remove_cv_t<
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/log1p.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/log1p.hpp
index 22f16238db..f227760f46 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/log1p.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/log1p.hpp
@@ -65,7 +65,7 @@ template <typename argT, typename resT> struct Log1pFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         if constexpr (is_complex<argT>::value) {
             // log1p(z) = ln((x + 1) + yI)
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/log2.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/log2.hpp
index 91bd063117..e7ae2911ac 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/log2.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/log2.hpp
@@ -66,7 +66,7 @@ template <typename argT, typename resT> struct Log2Functor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         if constexpr (is_complex<argT>::value) {
             using realT = typename argT::value_type;
@@ -78,7 +78,7 @@ template <typename argT, typename resT> struct Log2Functor
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT, vec_sz> &in)
+    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT, vec_sz> &in) const
     {
         auto const &res_vec = sycl::log2(in);
         using deducedT = typename std::remove_cv_t<
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/logaddexp.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/logaddexp.hpp
index b718a5f991..0721b7582b 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/logaddexp.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/logaddexp.hpp
@@ -59,14 +59,15 @@ template <typename argT1, typename argT2, typename resT> struct LogAddExpFunctor
     using supports_sg_loadstore = std::true_type;
     using supports_vec = std::true_type;
 
-    resT operator()(const argT1 &in1, const argT2 &in2)
+    resT operator()(const argT1 &in1, const argT2 &in2) const
     {
         return impl<resT>(in1, in2);
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT1, vec_sz> &in1,
-                                       const sycl::vec<argT2, vec_sz> &in2)
+    sycl::vec<resT, vec_sz>
+    operator()(const sycl::vec<argT1, vec_sz> &in1,
+               const sycl::vec<argT2, vec_sz> &in2) const
     {
         sycl::vec<resT, vec_sz> res;
         auto diff = in1 - in2; // take advantange of faster vec arithmetic
@@ -86,7 +87,7 @@ template <typename argT1, typename argT2, typename resT> struct LogAddExpFunctor
     }
 
 private:
-    template <typename T> T impl(T const &in1, T const &in2)
+    template <typename T> T impl(T const &in1, T const &in2) const
     {
         if (in1 == in2) { // handle signed infinities
             const T log2 = std::log(T(2));
@@ -106,7 +107,7 @@ template <typename argT1, typename argT2, typename resT> struct LogAddExpFunctor
         }
     }
 
-    template <typename T> T impl_finite(T const &in)
+    template <typename T> T impl_finite(T const &in) const
     {
         return (in > 0) ? (in + std::log1p(std::exp(-in)))
                         : std::log1p(std::exp(in));
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/logical_and.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/logical_and.hpp
index 1b58c9792d..7e5091c04c 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/logical_and.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/logical_and.hpp
@@ -62,7 +62,7 @@ struct LogicalAndFunctor
         std::negation<std::disjunction<tu_ns::is_complex<argT1>,
                                        tu_ns::is_complex<argT2>>>>;
 
-    resT operator()(const argT1 &in1, const argT2 &in2)
+    resT operator()(const argT1 &in1, const argT2 &in2) const
     {
         using tu_ns::convert_impl;
 
@@ -71,8 +71,9 @@ struct LogicalAndFunctor
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT1, vec_sz> &in1,
-                                       const sycl::vec<argT2, vec_sz> &in2)
+    sycl::vec<resT, vec_sz>
+    operator()(const sycl::vec<argT1, vec_sz> &in1,
+               const sycl::vec<argT2, vec_sz> &in2) const
     {
 
         auto tmp = (in1 && in2);
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/logical_or.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/logical_or.hpp
index 7f2f0b8f7a..88e17d37cd 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/logical_or.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/logical_or.hpp
@@ -61,7 +61,7 @@ template <typename argT1, typename argT2, typename resT> struct LogicalOrFunctor
         std::negation<std::disjunction<tu_ns::is_complex<argT1>,
                                        tu_ns::is_complex<argT2>>>>;
 
-    resT operator()(const argT1 &in1, const argT2 &in2)
+    resT operator()(const argT1 &in1, const argT2 &in2) const
     {
         using tu_ns::convert_impl;
 
@@ -70,8 +70,9 @@ template <typename argT1, typename argT2, typename resT> struct LogicalOrFunctor
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT1, vec_sz> &in1,
-                                       const sycl::vec<argT2, vec_sz> &in2)
+    sycl::vec<resT, vec_sz>
+    operator()(const sycl::vec<argT1, vec_sz> &in1,
+               const sycl::vec<argT2, vec_sz> &in2) const
     {
 
         auto tmp = (in1 || in2);
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/logical_xor.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/logical_xor.hpp
index 0db9003970..5ba2ca6bdb 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/logical_xor.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/logical_xor.hpp
@@ -62,7 +62,7 @@ struct LogicalXorFunctor
         std::negation<std::disjunction<tu_ns::is_complex<argT1>,
                                        tu_ns::is_complex<argT2>>>>;
 
-    resT operator()(const argT1 &in1, const argT2 &in2)
+    resT operator()(const argT1 &in1, const argT2 &in2) const
     {
         using tu_ns::convert_impl;
 
@@ -71,8 +71,9 @@ struct LogicalXorFunctor
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT1, vec_sz> &in1,
-                                       const sycl::vec<argT2, vec_sz> &in2)
+    sycl::vec<resT, vec_sz>
+    operator()(const sycl::vec<argT1, vec_sz> &in1,
+               const sycl::vec<argT2, vec_sz> &in2) const
     {
         using tu_ns::vec_cast;
         auto tmp1 = vec_cast<bool, argT1, vec_sz>(in1);
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/maximum.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/maximum.hpp
index 60afc2fc15..57fcd46102 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/maximum.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/maximum.hpp
@@ -60,7 +60,7 @@ template <typename argT1, typename argT2, typename resT> struct MaximumFunctor
         std::negation<std::disjunction<tu_ns::is_complex<argT1>,
                                        tu_ns::is_complex<argT2>>>>;
 
-    resT operator()(const argT1 &in1, const argT2 &in2)
+    resT operator()(const argT1 &in1, const argT2 &in2) const
     {
         if constexpr (tu_ns::is_complex<argT1>::value ||
                       tu_ns::is_complex<argT2>::value)
@@ -77,8 +77,9 @@ template <typename argT1, typename argT2, typename resT> struct MaximumFunctor
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT1, vec_sz> &in1,
-                                       const sycl::vec<argT2, vec_sz> &in2)
+    sycl::vec<resT, vec_sz>
+    operator()(const sycl::vec<argT1, vec_sz> &in1,
+               const sycl::vec<argT2, vec_sz> &in2) const
     {
         sycl::vec<resT, vec_sz> res;
 #pragma unroll
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/minimum.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/minimum.hpp
index d19c829a6c..0ae1d056d0 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/minimum.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/minimum.hpp
@@ -60,7 +60,7 @@ template <typename argT1, typename argT2, typename resT> struct MinimumFunctor
         std::negation<std::disjunction<tu_ns::is_complex<argT1>,
                                        tu_ns::is_complex<argT2>>>>;
 
-    resT operator()(const argT1 &in1, const argT2 &in2)
+    resT operator()(const argT1 &in1, const argT2 &in2) const
     {
         if constexpr (tu_ns::is_complex<argT1>::value ||
                       tu_ns::is_complex<argT2>::value)
@@ -77,8 +77,9 @@ template <typename argT1, typename argT2, typename resT> struct MinimumFunctor
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT1, vec_sz> &in1,
-                                       const sycl::vec<argT2, vec_sz> &in2)
+    sycl::vec<resT, vec_sz>
+    operator()(const sycl::vec<argT1, vec_sz> &in1,
+               const sycl::vec<argT2, vec_sz> &in2) const
     {
         sycl::vec<resT, vec_sz> res;
 #pragma unroll
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/multiply.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/multiply.hpp
index c35a096254..aa416d196d 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/multiply.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/multiply.hpp
@@ -58,14 +58,15 @@ template <typename argT1, typename argT2, typename resT> struct MultiplyFunctor
     using supports_vec = std::negation<
         std::disjunction<tu_ns::is_complex<argT1>, tu_ns::is_complex<argT2>>>;
 
-    resT operator()(const argT1 &in1, const argT2 &in2)
+    resT operator()(const argT1 &in1, const argT2 &in2) const
     {
         return in1 * in2;
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT1, vec_sz> &in1,
-                                       const sycl::vec<argT2, vec_sz> &in2)
+    sycl::vec<resT, vec_sz>
+    operator()(const sycl::vec<argT1, vec_sz> &in1,
+               const sycl::vec<argT2, vec_sz> &in2) const
     {
         auto tmp = in1 * in2;
         if constexpr (std::is_same_v<resT,
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/negative.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/negative.hpp
index 830794680c..4cfc231d7a 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/negative.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/negative.hpp
@@ -61,7 +61,7 @@ template <typename argT, typename resT> struct NegativeFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &x)
+    resT operator()(const argT &x) const
     {
         return -x;
     }
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/not_equal.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/not_equal.hpp
index f01068ab49..7d408e7f5c 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/not_equal.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/not_equal.hpp
@@ -60,7 +60,7 @@ template <typename argT1, typename argT2, typename resT> struct NotEqualFunctor
         std::negation<std::disjunction<tu_ns::is_complex<argT1>,
                                        tu_ns::is_complex<argT2>>>>;
 
-    resT operator()(const argT1 &in1, const argT2 &in2)
+    resT operator()(const argT1 &in1, const argT2 &in2) const
     {
         if constexpr (std::is_same_v<argT1, std::complex<float>> &&
                       std::is_same_v<argT2, float>)
@@ -78,8 +78,9 @@ template <typename argT1, typename argT2, typename resT> struct NotEqualFunctor
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT1, vec_sz> &in1,
-                                       const sycl::vec<argT2, vec_sz> &in2)
+    sycl::vec<resT, vec_sz>
+    operator()(const sycl::vec<argT1, vec_sz> &in1,
+               const sycl::vec<argT2, vec_sz> &in2) const
     {
         auto tmp = (in1 != in2);
         if constexpr (std::is_same_v<resT,
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/positive.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/positive.hpp
index d5fab29315..2d81a0cc48 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/positive.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/positive.hpp
@@ -62,13 +62,13 @@ template <typename argT, typename resT> struct PositiveFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &x)
+    resT operator()(const argT &x) const
     {
         return x;
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT, vec_sz> &in)
+    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT, vec_sz> &in) const
     {
         auto const &res_vec = in;
         using deducedT = typename std::remove_cv_t<
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/pow.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/pow.hpp
index 62cd826c75..e3b44f3191 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/pow.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/pow.hpp
@@ -61,7 +61,7 @@ template <typename argT1, typename argT2, typename resT> struct PowFunctor
                                        std::is_integral<argT1>,
                                        std::is_integral<argT2>>>;
 
-    resT operator()(argT1 in1, argT2 in2)
+    resT operator()(argT1 in1, argT2 in2) const
     {
         if constexpr (std::is_integral_v<argT1> || std::is_integral_v<argT2>) {
             if constexpr (std::is_signed_v<argT2>) {
@@ -89,8 +89,9 @@ template <typename argT1, typename argT2, typename resT> struct PowFunctor
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT1, vec_sz> &in1,
-                                       const sycl::vec<argT2, vec_sz> &in2)
+    sycl::vec<resT, vec_sz>
+    operator()(const sycl::vec<argT1, vec_sz> &in1,
+               const sycl::vec<argT2, vec_sz> &in2) const
     {
         auto res = sycl::pow(in1, in2);
         if constexpr (std::is_same_v<resT,
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/proj.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/proj.hpp
index 6fbc3ae012..b8ad1042ff 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/proj.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/proj.hpp
@@ -65,17 +65,28 @@ template <typename argT, typename resT> struct ProjFunctor
     // do both argTy and resTy support sugroup store/load operation
     using supports_sg_loadstore = typename std::false_type;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         using realT = typename argT::value_type;
         const realT x = std::real(in);
         const realT y = std::imag(in);
 
-        if (std::isinf(x) || std::isinf(y)) {
-            const realT res_im = std::copysign(realT(0), y);
-            return resT{std::numeric_limits<realT>::infinity(), res_im};
+        if (std::isinf(x)) {
+            return value_at_infinity(y);
         }
-        return in;
+        else if (std::isinf(y)) {
+            return value_at_infinity(y);
+        }
+        else {
+            return in;
+        }
+    }
+
+private:
+    template <typename T> std::complex<T> value_at_infinity(const T &y) const
+    {
+        const T res_im = std::copysign(T(0), y);
+        return std::complex<T>{std::numeric_limits<T>::infinity(), res_im};
     }
 };
 
@@ -114,6 +125,30 @@ sycl::event proj_contig_impl(sycl::queue exec_q,
         exec_q, nelems, arg_p, res_p, depends);
 }
 
+template <typename argTy>
+sycl::event
+proj_workaround_contig_impl(sycl::queue exec_q,
+                            size_t nelems,
+                            const char *arg_p,
+                            char *res_p,
+                            const std::vector<sycl::event> &depends = {})
+{
+    using resTy = typename ProjOutputType<argTy>::value_type;
+
+    const argTy *arg_tp = reinterpret_cast<const argTy *>(arg_p);
+    resTy *res_tp = reinterpret_cast<resTy *>(res_p);
+
+    sycl::event e = exec_q.submit([&](sycl::handler &cgh) {
+        cgh.depends_on(depends);
+        cgh.parallel_for({nelems}, [=](sycl::id<1> id) {
+            size_t i = id[0];
+            res_tp[i] = ProjFunctor<argTy, resTy>{}(arg_tp[i]);
+        });
+    });
+
+    return e;
+}
+
 template <typename fnT, typename T> struct ProjContigFactory
 {
     fnT get()
@@ -124,8 +159,14 @@ template <typename fnT, typename T> struct ProjContigFactory
             return fn;
         }
         else {
-            fnT fn = proj_contig_impl<T>;
-            return fn;
+            if constexpr (std::is_same_v<T, std::complex<double>>) {
+                fnT fn = proj_workaround_contig_impl<T>;
+                return fn;
+            }
+            else {
+                fnT fn = proj_contig_impl<T>;
+                return fn;
+            }
         }
     }
 };
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/real.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/real.hpp
index 43e256913f..3faf0ce553 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/real.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/real.hpp
@@ -65,7 +65,7 @@ template <typename argT, typename resT> struct RealFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         if constexpr (is_complex<argT>::value) {
             return std::real(in);
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/remainder.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/remainder.hpp
index e09edf9ec4..54bc42c7c6 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/remainder.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/remainder.hpp
@@ -56,7 +56,7 @@ template <typename argT1, typename argT2, typename resT> struct RemainderFunctor
     using supports_sg_loadstore = std::true_type;
     using supports_vec = std::true_type;
 
-    resT operator()(const argT1 &in1, const argT2 &in2)
+    resT operator()(const argT1 &in1, const argT2 &in2) const
     {
         if constexpr (std::is_integral_v<argT1> || std::is_integral_v<argT2>) {
             if (in2 == argT2(0)) {
@@ -88,8 +88,9 @@ template <typename argT1, typename argT2, typename resT> struct RemainderFunctor
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT1, vec_sz> &in1,
-                                       const sycl::vec<argT2, vec_sz> &in2)
+    sycl::vec<resT, vec_sz>
+    operator()(const sycl::vec<argT1, vec_sz> &in1,
+               const sycl::vec<argT2, vec_sz> &in2) const
     {
         if constexpr (std::is_integral_v<argT1> || std::is_integral_v<argT2>) {
             sycl::vec<resT, vec_sz> rem;
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/round.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/round.hpp
index ac61829a1d..8fceacfcef 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/round.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/round.hpp
@@ -63,7 +63,7 @@ template <typename argT, typename resT> struct RoundFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
 
         if constexpr (std::is_integral_v<argT>) {
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sign.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sign.hpp
index 8f0488746e..b59b0c49ad 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sign.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sign.hpp
@@ -61,7 +61,7 @@ template <typename argT, typename resT> struct SignFunctor
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
     using supports_sg_loadstore = std::false_type;
 
-    resT operator()(const argT &x)
+    resT operator()(const argT &x) const
     {
         if constexpr (std::is_integral_v<argT>) {
             if constexpr (std::is_unsigned_v<argT>) {
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/signbit.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/signbit.hpp
index 016e9cb0b2..98b9df0716 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/signbit.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/signbit.hpp
@@ -65,7 +65,7 @@ template <typename argT, typename resT> struct SignbitFunctor
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT, vec_sz> &in)
+    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT, vec_sz> &in) const
     {
         auto const &res_vec = sycl::signbit(in);
 
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sin.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sin.hpp
index 74cab65ec4..70ebf053ff 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sin.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sin.hpp
@@ -62,7 +62,7 @@ template <typename argT, typename resT> struct SinFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         if constexpr (is_complex<argT>::value) {
             using realT = typename argT::value_type;
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sinh.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sinh.hpp
index f535f3e946..9ce2925246 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sinh.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sinh.hpp
@@ -63,7 +63,7 @@ template <typename argT, typename resT> struct SinhFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         if constexpr (is_complex<argT>::value) {
             using realT = typename argT::value_type;
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sqrt.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sqrt.hpp
index 9b9fa95061..40fda0e147 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sqrt.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/sqrt.hpp
@@ -66,7 +66,7 @@ template <typename argT, typename resT> struct SqrtFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         if constexpr (is_complex<argT>::value) {
             // #ifdef _WINDOWS
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/square.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/square.hpp
index 29d096ae88..517cf2dc26 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/square.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/square.hpp
@@ -66,13 +66,13 @@ template <typename argT, typename resT> struct SquareFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         return in * in;
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT, vec_sz> &in)
+    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT, vec_sz> &in) const
     {
         auto const &res_vec = in * in;
         using deducedT = typename std::remove_cv_t<
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/subtract.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/subtract.hpp
index f286a6d39a..652308b2c4 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/subtract.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/subtract.hpp
@@ -57,14 +57,15 @@ template <typename argT1, typename argT2, typename resT> struct SubtractFunctor
     using supports_vec = std::negation<
         std::disjunction<tu_ns::is_complex<argT1>, tu_ns::is_complex<argT2>>>;
 
-    resT operator()(const argT1 &in1, const argT2 &in2)
+    resT operator()(const argT1 &in1, const argT2 &in2) const
     {
         return in1 - in2;
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT1, vec_sz> &in1,
-                                       const sycl::vec<argT2, vec_sz> &in2)
+    sycl::vec<resT, vec_sz>
+    operator()(const sycl::vec<argT1, vec_sz> &in1,
+               const sycl::vec<argT2, vec_sz> &in2) const
     {
         auto tmp = in1 - in2;
         if constexpr (std::is_same_v<resT,
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/tan.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/tan.hpp
index 3f11b0e721..25bcea68e3 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/tan.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/tan.hpp
@@ -64,7 +64,7 @@ template <typename argT, typename resT> struct TanFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         if constexpr (is_complex<argT>::value) {
 
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/tanh.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/tanh.hpp
index 4ee18eccb4..41f5b5bd5e 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/tanh.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/tanh.hpp
@@ -65,7 +65,7 @@ template <typename argT, typename resT> struct TanhFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         if constexpr (is_complex<argT>::value) {
             using realT = typename argT::value_type;
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/true_divide.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/true_divide.hpp
index 110eaf3a4f..89e5f5c67c 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/true_divide.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/true_divide.hpp
@@ -58,14 +58,15 @@ struct TrueDivideFunctor
     using supports_vec = std::negation<
         std::disjunction<tu_ns::is_complex<argT1>, tu_ns::is_complex<argT2>>>;
 
-    resT operator()(const argT1 &in1, const argT2 &in2)
+    resT operator()(const argT1 &in1, const argT2 &in2) const
     {
         return in1 / in2;
     }
 
     template <int vec_sz>
-    sycl::vec<resT, vec_sz> operator()(const sycl::vec<argT1, vec_sz> &in1,
-                                       const sycl::vec<argT2, vec_sz> &in2)
+    sycl::vec<resT, vec_sz>
+    operator()(const sycl::vec<argT1, vec_sz> &in1,
+               const sycl::vec<argT2, vec_sz> &in2) const
     {
         auto tmp = in1 / in2;
         if constexpr (std::is_same_v<resT,
diff --git a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/trunc.hpp b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/trunc.hpp
index c49bd2e69f..5fb4074ec9 100644
--- a/dpctl/tensor/libtensor/include/kernels/elementwise_functions/trunc.hpp
+++ b/dpctl/tensor/libtensor/include/kernels/elementwise_functions/trunc.hpp
@@ -63,7 +63,7 @@ template <typename argT, typename resT> struct TruncFunctor
     using supports_sg_loadstore = typename std::negation<
         std::disjunction<is_complex<resT>, is_complex<argT>>>;
 
-    resT operator()(const argT &in)
+    resT operator()(const argT &in) const
     {
         if constexpr (std::is_integral_v<argT>) {
             return in;
diff --git a/dpctl/tests/conftest.py b/dpctl/tests/conftest.py
index 06e92c1151..c4666fe5ad 100644
--- a/dpctl/tests/conftest.py
+++ b/dpctl/tests/conftest.py
@@ -53,16 +53,15 @@ def pytest_addoption(parser):
         "--runcomplex",
         action="store_true",
         default=False,
-        help="run broken complex tests on Windows",
+        help="run broken complex tests",
     )
 
 
 def pytest_collection_modifyitems(config, items):
     if config.getoption("--runcomplex"):
         return
-    skip_complex = pytest.mark.skipif(
-        os.name == "nt",
-        reason="need --runcomplex option to run on Windows",
+    skip_complex = pytest.mark.skip(
+        reason="need --runcomplex option to run",
     )
     for item in items:
         if "broken_complex" in item.keywords:
diff --git a/dpctl/tests/elementwise/test_exp.py b/dpctl/tests/elementwise/test_exp.py
index 96314ad46c..4886c0cb78 100644
--- a/dpctl/tests/elementwise/test_exp.py
+++ b/dpctl/tests/elementwise/test_exp.py
@@ -198,6 +198,7 @@ def test_exp_complex_strided(dtype):
             )
 
 
+@pytest.mark.broken_complex
 @pytest.mark.parametrize("dtype", ["c8", "c16"])
 def test_exp_complex_special_cases(dtype):
     q = get_queue_or_skip()
diff --git a/dpctl/tests/elementwise/test_sqrt.py b/dpctl/tests/elementwise/test_sqrt.py
index 7e705f0721..862f64ccbd 100644
--- a/dpctl/tests/elementwise/test_sqrt.py
+++ b/dpctl/tests/elementwise/test_sqrt.py
@@ -157,6 +157,7 @@ def test_sqrt_real_fp_special_values(dtype):
     assert dpt.allclose(r, expected, atol=tol, rtol=tol, equal_nan=True)
 
 
+@pytest.mark.broken_complex
 @pytest.mark.parametrize("dtype", _complex_fp_dtypes)
 def test_sqrt_complex_fp_special_values(dtype):
     q = get_queue_or_skip()
diff --git a/pyproject.toml b/pyproject.toml
index b5a8d33c32..5882128dab 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -34,7 +34,7 @@ omit = [
 
 [tool.pytest.ini.options]
 markers = [
-    "broken_complex: mark a test that is skipped on Windows due to complex implementation",
+    "broken_complex: mark a test that is skipped due to complex implementation issues in DPC++ compiler",
 ]
 minversion = "6.0"
 norecursedirs= [