action: template type for U and Omega, to allow differentiation vs. them

src/action.hpp

@@ -55,10 +55,13 @@ Eigen::Matrix<typename XType::value_type,4,4> nambu_diag_mul(const LRType& L, co
     return Y;
 }
 
-template <typename ConstScalar, typename Shape0, typename Shape, typename Shape2, typename Scalar = typename std::remove_const<ConstScalar>::type>
+template <typename ConstScalar, typename ConstScalar2,
+          typename Shape0, typename Shape, typename Shape2,
+          typename Scalar = typename std::remove_const<ConstScalar>::type,
+          typename Scalar2 = typename std::remove_const<ConstScalar2>::type>
 inline Scalar S_2(const array::ArrayView<const Mask, Shape0>& mask,
                   const array::ArrayView<ConstScalar, Shape>& Qmat,
-                  const array::ArrayView<const Complex, Shape2>& U,
+                  const array::ArrayView<ConstScalar2, Shape2>& U,
                   const double& Lx,
                   const double& Ly,
                   const Complex& D)
@@ -109,11 +112,14 @@ inline Scalar S_2(const array::ArrayView<const Mask, Shape0>& mask,
     return (D*dd[0]*dd[1]*2.0/4.0) * S;
 }
 
-template <typename ConstScalar, typename Shape0, typename Shape, typename Shape2, typename Scalar = typename std::remove_const<ConstScalar>::type>
+template <typename ConstScalar, typename ConstScalar2,
+          typename Shape0, typename Shape, typename Shape2,
+          typename Scalar = typename std::remove_const<ConstScalar>::type,
+          typename Scalar2 = typename std::remove_const<ConstScalar>::type>
 inline Scalar S_Omega(const array::ArrayView<const Mask, Shape0>& mask,
                       const array::ArrayView<ConstScalar, Shape>& Qmat,
                       const Scalar& omega,
-                      const array::ArrayView<const Complex, Shape2>& Omega,
+                      const array::ArrayView<ConstScalar2, Shape2>& Omega,
                       const double& Lx,
                       const double& Ly)
 {
@@ -140,10 +146,11 @@ inline Scalar S_Omega(const array::ArrayView<const Mask, Shape0>& mask,
     return (dd[0]*dd[1]) * S;
 }
 
-template <typename ConstScalar, typename Shape0, typename Shape, typename Shape2, typename Scalar = typename std::remove_const<ConstScalar>::type>
+template <typename ConstScalar, typename ConstScalar2, typename Shape0, typename Shape, typename Shape2,
+          typename Scalar = typename std::remove_const<ConstScalar>::type, typename Scalar2 = typename std::remove_const<ConstScalar2>::type>
 inline Scalar S_H(const array::ArrayView<const Mask, Shape0>& mask,
                   const array::ArrayView<ConstScalar, Shape>& Qmat,
-                  const array::ArrayView<const Complex, Shape2>& U,
+                  const array::ArrayView<ConstScalar2, Shape2>& U,
                   const double& Lx,
                   const double& Ly,
                   const Complex& alpha)
@@ -241,12 +248,16 @@ Complex flat_value(const Complex& value)
 }
 CPPAD_DISCRETE_FUNCTION(Complex, flat_value);
 
-template <typename ConstScalar, typename Shape0, typename Shape1, typename Shape2, typename Shape3, typename Scalar = typename std::remove_const<ConstScalar>::type>
+template <typename ConstScalar, typename ConstScalar2, typename ConstScalar3,
+          typename Shape0, typename Shape1, typename Shape2, typename Shape3,
+          typename Scalar = typename std::remove_const<ConstScalar>::type,
+          typename Scalar2 = typename std::remove_const<ConstScalar2>::type,
+          typename Scalar3 = typename std::remove_const<ConstScalar3>::type>
 inline Scalar S(const array::ArrayView<const Mask, Shape0>& mask,
                 const array::ArrayView<ConstScalar, Shape1>& Q,
                 const Scalar& omega,
-                const array::ArrayView<const Complex, Shape2>& U,
-                const array::ArrayView<const Complex, Shape3>& Omega,
+                const array::ArrayView<ConstScalar2, Shape2>& U,
+                const array::ArrayView<ConstScalar3, Shape3>& Omega,
                 const double& Lx,
                 const double& Ly,
                 const Complex& alpha,
@@ -292,43 +303,5 @@ inline Scalar S(const array::ArrayView<const Mask, Shape0>& mask,
     return S;
 }
 
-
-template <typename Scalar, typename Shape>
-inline Scalar laplacian(array::ArrayView<Scalar, Shape> Q)
-{
-    const int N_i[4] = {0, 1, -1, 0};
-    const int N_j[4] = {1, 0, 0, -1};
-    Scalar S = Scalar(0);
-
-    size_t nx = Q.dim(0);
-    size_t ny = Q.dim(1);
-
-    for (size_t i = 0; i < nx; ++i) {
-        for (size_t j = 0; j < ny; ++j) {
-            /* Sum over neighbors */
-            for (size_t k = 0; k < 4; ++k) {
-                int di = N_i[k];
-                int dj = N_j[k];
-                
-                if (i == 0 && di < 0)
-                    continue;
-                if (i == nx-1 && di > 0)
-                    continue;
-                if (j == 0 && dj < 0)
-                    continue;
-                if (j == ny-1 && dj > 0)
-                    continue;
-
-                size_t i2 = i + di;
-                size_t j2 = j + dj;
-
-                S += (Q(i,j) - Q(i2,j2)) * (Q(i,j) - Q(i2,j2));
-            }
-        }
-    }
-
-    return S;
-}
-
 #endif