The kernel is now

k(x,y) = exp(-| x - y |^2 / (\sigma))

include/scalfmm/matrix_kernels/gaussian.hpp

@@ -13,10 +13,11 @@ namespace scalfmm::matrix_kernels
     * @brief This structure corresponds to the Gaussian kernel.
     *
     *   The kernel is defined as \f$k(x,y): R^{km} -> R^{kn}\f$ with \f$ kn = km = 1\f$
-     *                            \f$k(x,y) = exp(-| x - y |/(\sigma^2))\f$
+     *                            \f$k(x,y) = exp(-| x - y |^2 / (\sigma))\f$
     *
     * - The kernel is not homogeneous.
     * - The kernel is not symmetric.
+    * - The kernel is smooth.
     *
     * @GlobalValueType The value type of the coefficient \f$\sigma\f$.
     */
@@ -28,8 +29,10 @@ namespace scalfmm::matrix_kernels
         // Mandatory constants
         static constexpr auto homogeneity_tag{homogeneity::non_homogenous};   // Specify the homogeneity of the kernel.
         static constexpr auto symmetry_tag{symmetry::non_symmetric};          // Specify the symmetry of the kernel.
-        static constexpr std::size_t km{1};                                   // The number of inputs of the kernel.
-        static constexpr std::size_t kn{1};                                   // The number of outputs of the kernel.
+        static constexpr bool is_smooth = true;   // Specify the kernel is smooth and can be evaluated at x == y.
+
+        static constexpr std::size_t km{1};   // The number of inputs of the kernel.
+        static constexpr std::size_t kn{1};   // The number of outputs of the kernel.
 
         // Mandatory type
         template<typename ValueType>
@@ -38,8 +41,6 @@ namespace scalfmm::matrix_kernels
         using vector_type = std::array<ValueType, kn>;   // Vector type that is used in the kernel.
 
         // Optional constant
-        static constexpr bool is_smooth = true;   // Specify the kernel is smooth and can be evaluated at x == y.
-
         value_type m_coeff{value_type(1.)};   // parameter/coefficient of the kernel.
 
         /**
@@ -136,7 +137,7 @@ namespace scalfmm::matrix_kernels
                               matrix_type<std::decay_t<ValueType1>>>
         {
             using decayed_type = std::decay_t<ValueType1>;
-            auto l2 = decayed_type(-1.0) / (m_coeff * m_coeff);
+            auto l2 = decayed_type(-1.0) / (m_coeff);
             decayed_type r2 = (((xs.at(Is) - ys.at(Is)) * (xs.at(Is) - ys.at(Is))) + ...);
             return matrix_type<decayed_type>{xsimd::exp(r2 * l2)};
         }