Towards device-enabled project function

Had to copy a bunch of headers to remove the SimpleTensor/FixedTensor
and adjust for potential compilation in CUDA context.
Does not yet compile because we cannot co_await from within
a function that is not a coroutine. Need to come up with a way to
switch between coroutines inside a task.

Signed-off-by: Joseph Schuchart <[email protected]>

examples/CMakeLists.txt

@@ -58,6 +58,11 @@ if (TARGET tiledarray)
     endif(TTG_HAVE_HIP)
 endif()
 
+add_ttg_executable(tensortest madness/mra-device/tests/tensor_test.cc)
+#if (TTG_HAVE_CUDA)
+  add_ttg_executable(mradevice-cuda madness/mra-device/mrattg-device.cc madness/mra-device/gl.cc RUNTIMES "parsec")
+ #endif (TTG_HAVE_CUDA)
+
 if (TARGET MADworld)
   add_ttg_executable(madness-1d madness/madness-1d/madness-1d.cc RUNTIMES "mad")
   if (TARGET blaspp) #(CBLAS_FOUND AND MKL_FOUND)
@@ -66,7 +71,6 @@ if (TARGET MADworld)
     add_ttg_executable(mrattg-streaming madness/mrattg_streaming.cc mragl.cc mratwoscale.cc mradomain.h mrafunctiondata.h mrafunctionfunctor.h mrafunctionnode.h mragl.h mrahash.h mrakey.h mramisc.h mramxm.h mrarange.h mrasimpletensor.h mratwoscale.h mratypes.h LINK_LIBRARIES blaspp MADworld)
   endif ()
 
-  add_ttg_executable(mradevice madness/mra-device/tests/tensor_test.cc)
 endif (TARGET MADworld)
 
 add_ttg_executable(wavefront-wf wavefront/wavefront-wf.cc SINGLERANKONLY)

examples/madness/mra-device/functiondata.h

@@ -0,0 +1,149 @@
+#ifndef MADFUNCTIONDATA_H_INCL
+#define MADFUNCTIONDATA_H_INCL
+
+#include "../../mratypes.h"
+#include "../../mradomain.h"
+#include "../../mratwoscale.h"
+#include "tensor.h"
+#include "gl.h"
+
+namespace mra {
+
+    /// Convenient co-location of frequently used data
+    template <typename T, Dimension NDIM>
+    class FunctionData {
+        std::size_t K;
+        Tensor<T,2> phi; // phi(mu,i) = phi(x[mu],i) --- value of scaling functions at quadrature points on level 0
+        Tensor<T,2> phibar; // phibar(mu,i) = w[mu]*phi(x[mu],i)
+        Tensor<T,2> HG; // Two scale filter applied from left to scaling function coeffs
+        Tensor<T,2> HGT; // Two scale filter applied from right to scaling function coeffs
+        Tensor<T,2> rm, r0, rp; // blocks of the ABGV central derivative operator
+        std::unique_ptr<T[]> x, w; // Quadrature points and weights on level 0
+
+        void make_abgv_diff_operator() {
+            double iphase = 1.0;
+            auto r0_view = r0.current_view();
+            auto rm_view = rm.current_view();
+            auto rp_view = rp.current_view();
+            for (auto i: range(K)) {
+                double jphase = 1.0;
+                for (auto j : range(K)) {
+                    double gammaij = std::sqrt(double((2*i+1)*(2*j+1)));
+                    double Kij;
+                    if (((i-j)>0) && (((i-j)%2)==1))
+                        Kij = 2.0;
+                    else
+                        Kij = 0.0;
+
+                    r0_view(i,j) = T(0.5*(1.0 - iphase*jphase - 2.0*Kij)*gammaij);
+                    rm_view(i,j) = T(0.5*jphase*gammaij);
+                    rp_view(i,j) = T(-0.5*iphase*gammaij);
+                }
+            }
+        }
+
+        /// Set phi(mu,i) to be phi(x[mu],i)
+        void make_phi() {
+            /* retrieve x, w from constant memory */
+            const T *x, *w;
+            detail::GLget(&x, &w, K);
+            T* p = new T[K];
+            auto phi_view = phi.current_view();
+            for (size_t mu : range(K)) {
+                detail::legendre_scaling_functions(x[mu], K, &p[0]);
+                for (size_t i : range(K)) {
+                    phi_view(mu,i) = p[i];
+                }
+            }
+            delete[] p;
+        }
+
+        /// Set phibar(mu,i) to be w[mu]*phi(x[mu],i)
+        void make_phibar() {
+            /* retrieve x, w from constant memory */
+            const T *x, *w;
+            detail::GLget(&x, &w, K);
+            T *p = new T[K];
+            auto phibar_view = phibar.current_view();
+            for (size_t mu : range(K)) {
+                detail::legendre_scaling_functions(x[mu], K, &p[0]);
+                for (size_t i : range(K)) {
+                    phibar_view(mu,i) = w[mu]*p[i];
+                }
+            }
+            delete[] p;
+            // FixedTensor<T,K,2> phi, r;
+            // make_phi<T,K>(phi);
+            // mTxmq(K, K, K, r.ptr(), phi.ptr(), phibar.ptr());
+            // std::cout << r << std::endl; // should be identify matrix
+        }
+
+    public:
+
+        FunctionData(std::size_t K)
+        : K(K)
+        , phi(K, K)
+        , phibar(K, K)
+        , HG(2*K, 2*K)
+        , HGT(2*K, 2*K)
+        , rm(K, K)
+        , r0(K, K)
+        , rp(K, K)
+        {
+            make_phi();
+            make_phibar();
+            twoscale_get(K, HG.data());
+            auto HG_view  = HG.current_view();
+            auto HGT_view = HGT.current_view();
+            for (size_t i : range(2*K)) {
+                for (size_t j : range(2*K)) {
+                    HGT_view(j,i) = HG_view(i,j);
+                }
+            }
+            make_abgv_diff_operator();
+        }
+
+        const auto& get_phi() const {return phi;}
+        const auto& get_phibar() const {return phibar;}
+        const auto& get_hg() const {return HG;}
+        const auto& get_hgT() const {return HGT;}
+        const auto& get_rm() const {return rm;}
+        const auto& get_r0() const {return r0;}
+        const auto& get_rp() const {return rp;}
+
+        /// Set X(d,mu) to be the mu'th quadrature point in dimension d for the box described by key
+        void make_quadrature_pts(const Key<NDIM>& key, TensorView<T,2>& X) {
+            const Level n = key.level();
+            const std::array<Translation,NDIM>& l = key.translation();
+            const T h = std::pow(T(0.5),T(n));
+            /* retrieve x[] from constant memory */
+            const T *x, *w;
+            detail::GLget(&x, &w, K);
+#ifdef __CUDA_ARCH__
+            if (blockIdx.z == 0) {
+                for (int d = blockIdx.y; d < x.Dim(0); d += blockDim.y) {
+                    T lo, hi; std::tie(lo,hi) = Domain<NDIM>::get(d);
+                    T width = h*Domain<NDIM>::get_width(d);
+                    for (int i = blockIdx.x; i < X.dim(1); i += blockDim.y) {
+                        X(d,i) = lo + width*(l[d] + x[i]);
+                    }
+                }
+            }
+            /* wait for all to complete */
+            __syncthreads();
+#else  // __CUDA_ARCH__
+            for (Dimension d : range(NDIM)) {
+                T lo, hi; std::tie(lo,hi) = Domain<NDIM>::get(d);
+                T width = h*Domain<NDIM>::get_width(d);
+                for (size_t i :  X.dim(1)) {
+                    X(d,i) = lo + width*(l[d] + x[i]);
+                }
+            }
+#endif // __CUDA_ARCH__
+        }
+    };
+
+}
+
+
+#endif

examples/madness/mra-device/functionfunctor.h

@@ -0,0 +1,114 @@
+#ifndef MAD_FUNCTION_FUNCTOR_H_INCL
+#define MAD_FUNCTION_FUNCTOR_H_INCL
+
+#include "tensor.h"
+
+#include "../../mratypes.h"
+
+namespace mra {
+
+
+    /// Function functor is going away ... don't use!
+    /// \code
+    ///template <typename T, Dimension NDIM>
+    /// class FunctionFunctor {
+    ///public:
+    /// T operator()(const Coordinate<T,NDIM>& r) const;
+    /// template <size_t K> void operator()(const SimpleTensor<T,NDIM,K>& x, FixedTensor<T,K,NDIM>& values) const;
+    /// Level initial_level() const;
+    /// bool is_negligible(const std::pair<Coordinate<T,NDIM>,Coordinate<T,NDIM>>& box, T thresh) const;
+    /// special point interface to be added
+    ///  }
+    /// \endcode
+
+    /// Adapts a simple callable to the API needed for evaluation --- implement your own for full vectorization
+    template <typename T, Dimension NDIM>
+    class FunctionFunctor {
+        std::function<T(const Coordinate<T,NDIM>&)> f;
+
+    public:
+        static const Level default_initial_level = 3; //< needs to become user configurable
+
+        template <typename functionT>
+        FunctionFunctor(functionT f) : f(f) {}
+
+        /// Evaluate at a single point
+        T operator()(const Coordinate<T,NDIM>& r) const {return f(r);}
+
+    };
+
+    /**
+     * TODO: adapt eval_cube to CUDA
+     */
+
+    /// Evaluate at points formed by tensor product of npt points in each dimension
+    template <typename functorT, typename T> void eval_cube(const functorT& f, const TensorView<T,1>& x, TensorView<T,1>& values) {
+        for (size_t i=0; i<x.dim(0); i++) values(i) = f(Coordinate<T,1>{x(0,i)});
+    }
+
+    /// Evaluate at points formed by tensor product of npt points in each dimension
+    template <typename functorT, typename T> void eval_cube(const functorT& f, const TensorView<T,2>& x, TensorView<T,2>& values) {
+        for (size_t i=0; i<x.dim(0); i++) {
+            for (size_t j=0; j<x.dim(1); j++) {
+                values(i,j) = f(Coordinate<T,2>{x(0,i),x(1,j)});
+            }
+        }
+    }
+
+    /// Evaluate at points formed by tensor product of K points in each dimension
+    template <typename functorT, typename T> void eval_cube(const functorT& f, const TensorView<T,3>& x, TensorView<T,3>& values) {
+        for (size_t i=0; i<x.dim(0); i++) {
+            for (size_t j=0; j<x.dim(1); j++) {
+                for (size_t k=0; k<x.dim(2); k++) {
+                    values(i,j,k) = f(Coordinate<T,3>{x(0,i),x(1,j),x(2,k)});
+                }
+            }
+        }
+    }
+
+    /// Evaluate at points formed by tensor product of K points in each dimension using vectorized form
+    template <typename functorT, typename T> void eval_cube_vec(const functorT& f, const TensorView<T,1>& x, T* values) {
+        for (size_t i=0; i<x.dim(0); i++) {
+            values[i] = f(x(0,i));
+        }
+    }
+
+    /// Evaluate at points formed by tensor product of K points in each dimension using vectorized form
+    template <typename functorT, typename T, size_t K2NDIM> void eval_cube_vec(const functorT& f, const TensorView<T,2>& x, T* values) {
+        for (size_t i=0; i<x.dim(0); i++) {
+            values[i] = f(x(0,i),x(1,i));
+        }
+    }
+
+    /// Evaluate at points formed by tensor product of K points in each dimension using vectorized form
+    template <typename functorT, typename T, size_t K2NDIM> void eval_cube_vec(const functorT& f, const TensorView<T,3>& x, T* values) {
+        for (size_t i=0; i<K2NDIM; i++) {
+            values[i] = f(x(0,i),x(1,i),x(2,i));
+        }
+    }
+
+    namespace detail {
+        template <class functorT> using initial_level_t =
+            decltype(std::declval<const functorT>().initial_level());
+        template <class functorT> using supports_initial_level =
+            ::mra::is_detected<initial_level_t,functorT>;
+
+        template <class functorT, class pairT> using is_negligible_t =
+            decltype(std::declval<const functorT>().is_negligible(std::declval<pairT>(),std::declval<double>()));
+        template <class functorT, class pairT> using supports_is_negligible =
+            ::mra::is_detected<is_negligible_t,functorT,pairT>;
+    }
+
+    template <typename functionT> Level initial_level(const functionT& f) {
+        if constexpr (detail::supports_initial_level<functionT>()) return f.initial_level();
+        else return 2; // <<<<<<<<<<<<<<< needs updating to make user configurable
+    }
+
+    template <typename functionT, typename T, Dimension NDIM>
+    bool is_negligible(const functionT& f, const std::pair<Coordinate<T,NDIM>,Coordinate<T,NDIM>>& box, T thresh) {
+        using pairT = std::pair<Coordinate<T,NDIM>,Coordinate<T,NDIM>>;
+        if constexpr (detail::supports_is_negligible<functionT,pairT>()) return f.is_negligible(box, thresh);
+        else return false;
+    }
+}
+#endif