[GPU] Interpolators

src/PYB11/Utilities/QuadraticInterpolator.py

@@ -96,4 +96,3 @@ def lowerBound(self,
     xmin = PYB11property(doc="Minimum x coordinate for table")
     xmax = PYB11property(doc="Maximum x coordinate for table")
     xstep = PYB11property(doc="delta x between tabulated values")
-    coeffs = PYB11property(doc="the fitting coefficients")

src/Utilities/QuadraticInterpolator.cc

@@ -7,32 +7,20 @@
 // Created by JMO, Fri Dec  4 14:28:08 PST 2020
 //----------------------------------------------------------------------------//
 #include "QuadraticInterpolator.hh"
+#include "umpire/Umpire.hpp"
+#include <algorithm>
 
-namespace Spheral {
+#include <Eigen/Dense>
 
-//------------------------------------------------------------------------------
-// Default constructor
-//------------------------------------------------------------------------------
-QuadraticInterpolator::QuadraticInterpolator():
-  mN1(),
-  mXmin(),
-  mXmax(),
-  mXstep(),
-  mcoeffs() {
-}
+namespace Spheral {
 
 //------------------------------------------------------------------------------
 // Constructor with sampled values
 //------------------------------------------------------------------------------
 QuadraticInterpolator::QuadraticInterpolator(double xmin,
                                              double xmax,
-                                             const std::vector<double>& yvals):
-  mN1(),
-  mXmin(),
-  mXmax(),
-  mXstep(),
-  mcoeffs() {
-  this->initialize(xmin, xmax, yvals);
+                                             const std::vector<double>& yvals) {
+  initialize(xmin, xmax, yvals);
 }
 
 //------------------------------------------------------------------------------
@@ -47,11 +35,12 @@ QuadraticInterpolator::initialize(double xmin,
   VERIFY2(n % 2 == 1, "QuadraticInterpolator::initialize requires an odd number of tabulated values");
   VERIFY2(xmax > xmin, "QuadraticInterpolator::initialize requires a positive domain: [" << xmin << " " << xmax << "]");
 
-  mN1 = (n - 1u)/2u - 1u;  // Maximum index into arrays
-  mXmin = xmin;
-  mXmax = xmax;
-  mXstep = (xmax - xmin)/(mN1 + 1u);
-  mcoeffs.resize(3*(mN1 + 1u));
+  double N1 = (n - 1u)/2u - 1u;  // Maximum index into arrays
+  double xstep = (xmax - xmin)/(N1 + 1u);
+  size_t N = 3*(N1 + 1u);
+  auto& rm = umpire::ResourceManager::getInstance();
+  auto host_allocator = rm.getAllocator("HOST");
+  double* vals = static_cast<double*>(host_allocator.allocate(N*sizeof(double)));
 
   typedef Eigen::Matrix<double, 3, 3, Eigen::RowMajor> EMatrix;
   typedef Eigen::Matrix<double, 3, 1> EVector;
@@ -62,37 +51,42 @@ QuadraticInterpolator::initialize(double xmin,
   double x0, x1, x2;
   EMatrix A;
   EVector X, B;
-  for (auto i0 = 0u; i0 <= mN1; ++i0) {
-    x0 = xmin + i0*mXstep;
-    x1 = x0 + 0.5*mXstep;
-    x2 = x0 + mXstep;
+  for (auto i0 = 0u; i0 <= N1; ++i0) {
+    x0 = xmin + i0*xstep;
+    x1 = x0 + 0.5*xstep;
+    x2 = x0 + xstep;
     A << 1.0, x0, x0*x0,
          1.0, x1, x1*x1,
          1.0, x2, x2*x2;
     B << yvals[2u*i0], yvals[2u*i0 + 1u], yvals[2u*i0 + 2u];
     X = A.inverse()*B;
-    mcoeffs[3*i0     ] = X(0);
-    mcoeffs[3*i0 + 1u] = X(1);
-    mcoeffs[3*i0 + 2u] = X(2);
+    vals[3*i0     ] = X(0);
+    vals[3*i0 + 1u] = X(1);
+    vals[3*i0 + 2u] = X(2);
   }
+#ifdef SPHERAL_GPU_ENABLED
+  auto allocator = rm.getAllocator("DEVICE");
+  mDeviceCoeffs = static_cast<double*>(allocator.allocate(N*sizeof(double)));
+  rm.copy(mDeviceCoeffs, vals, N*sizeof(double));
+#endif
+  mcoeffs = vals;
+  mN1 = N1;
+  mXmin = xmin;
+  mXmax = xmax;
+  mXstep = xstep;
 }
 
 //------------------------------------------------------------------------------
 // Destructor
 //------------------------------------------------------------------------------
 QuadraticInterpolator::~QuadraticInterpolator() {
-}
-
-//------------------------------------------------------------------------------
-// Equivalence
-//------------------------------------------------------------------------------
-bool
-QuadraticInterpolator::
-operator==(const QuadraticInterpolator& rhs) const {
-  return ((mN1 == rhs.mN1) and
-          (mXmin == rhs.mXmin) and
-          (mXmax == rhs.mXmax) and
-          (mcoeffs == rhs.mcoeffs));
+  auto& rm = umpire::ResourceManager::getInstance();
+#ifdef SPHERAL_GPU_ENABLED
+  auto allocator = rm.getAllocator("DEVICE");
+  allocator.deallocate(mDeviceCoeffs);
+#endif
+  auto host_allocator = rm.getAllocator("HOST");
+  host_allocator.deallocate(mcoeffs);
 }
 
 }

src/Utilities/QuadraticInterpolator.hh

@@ -9,57 +9,93 @@
 #ifndef __Spheral_QuadraticInterpolator__
 #define __Spheral_QuadraticInterpolator__
 
+#include "chai/ExecutionSpaces.hpp"
+#include "config.hh"
+
 #include <cstddef>
 #include <vector>
 
 namespace Spheral {
 
-class QuadraticInterpolator {
+class QuadraticInterpolatorBase {
 public:
   //--------------------------- Public Interface ---------------------------//
   // Constructors, destructors
-  template<typename Func>
-  QuadraticInterpolator(double xmin, double xmax, size_t n, const Func& F);
-  QuadraticInterpolator(double xmin, double xmax, const std::vector<double>& yvals);
-  QuadraticInterpolator();
-  ~QuadraticInterpolator();
-
-  // Initialize after construction, either with a function or tabulated values
-  template<typename Func>
-  void initialize(double xmin, double xmax, size_t n, const Func& f);
-  void initialize(double xmin, double xmax, const std::vector<double>& yvals);
+  SPHERAL_HOST_DEVICE QuadraticInterpolatorBase() = default;
 
   // Comparisons
-  bool operator==(const QuadraticInterpolator& rhs) const;
+  SPHERAL_HOST_DEVICE bool operator==(const QuadraticInterpolatorBase& rhs) const;
 
   // Interpolate for the y value
-  double operator()(const double x) const;
-  double prime(const double x) const;    // First derivative
-  double prime2(const double x) const;   // Second derivative
+  SPHERAL_HOST_DEVICE double operator()(const double x) const;
+  SPHERAL_HOST_DEVICE double prime(const double x) const;    // First derivative
+  SPHERAL_HOST_DEVICE double prime2(const double x) const;   // Second derivative
 
   // Same as above, but use a pre-computed table position (from lowerBound)
-  double operator()(const double x, const size_t i0) const;
-  double prime(const double x, const size_t i0) const;    // First derivative
-  double prime2(const double x, const size_t i0) const;   // Second derivative
+  SPHERAL_HOST_DEVICE double operator()(const double x, const size_t i0) const;
+  SPHERAL_HOST_DEVICE double prime(const double x, const size_t i0) const;    // First derivative
+  SPHERAL_HOST_DEVICE double prime2(const double x, const size_t i0) const;   // Second derivative
 
   // Return the lower bound index in the table for the given x coordinate
-  size_t lowerBound(const double x) const;
+  SPHERAL_HOST_DEVICE size_t lowerBound(const double x) const;
 
   // Allow read access the internal data representation
-  size_t size() const;                        // The size of the tabulated coefficient arrays
-  double xmin() const;                        // Minimum x coordinate for table              
-  double xmax() const;                        // Maximum x coordinate for table              
-  double xstep() const;                       // delta x between tabulated values            
-  const std::vector<double>& coeffs() const;  // the fitting coefficients
-
-private:
+  SPHERAL_HOST_DEVICE size_t size() const;                        // The size of the tabulated coefficient arrays
+  SPHERAL_HOST_DEVICE double xmin() const;                        // Minimum x coordinate for table
+  SPHERAL_HOST_DEVICE double xmax() const;                        // Maximum x coordinate for table
+  SPHERAL_HOST_DEVICE double xstep() const;                       // delta x between tabulated values
+
+protected:
+  SPHERAL_HOST QuadraticInterpolatorBase(size_t N1,
+                                         double xmin,
+                                         double xmax,
+                                         double xstep,
+                                         double* vals) :
+    mN1(N1),
+    mXmin(xmin),
+    mXmax(xmax),
+    mXstep(xstep),
+    mcoeffs(vals) { }
   //--------------------------- Private Interface --------------------------//
   // Member data
-  size_t mN1;
-  double mXmin, mXmax, mXstep;
-  std::vector<double> mcoeffs;
+  size_t mN1 = 0u;
+  double mXmin = 0.;
+  double mXmax = 0.;
+  double mXstep = 0.;
+  double* mcoeffs = nullptr;
+};
+
+class QuadraticInterpolator : public QuadraticInterpolatorBase {
+public:
+  template<typename Func>
+  QuadraticInterpolator(double xmin, double xmax, size_t n, const Func& F);
+  QuadraticInterpolator(double xmin, double xmax, const std::vector<double>& yvals);
+  ~QuadraticInterpolator();
+
+  // Initialize after construction, either with a function or tabulated values
+  template<typename Func>
+  void initialize(double xmin, double xmax, size_t n, const Func& f);
+  void initialize(double xmin, double xmax, const std::vector<double>& yvals);
+
+  template<typename QIView>
+  QIView view(chai::ExecutionSpace space) {
+    if (space == chai::CPU) {
+      return QIView(mN1, mXmin, mXmax, mXstep, mcoeffs);
+    } else {
+      return QIView(mN1, mXmin, mXmax, mXstep, mDeviceCoeffs);
+    }
+  }
+private:
+  double* mDeviceCoeffs = nullptr;
 };
 
+// For use on device
+class QIView : public QuadraticInterpolatorBase {
+public:
+  SPHERAL_HOST_DEVICE QIView() = default;
+  SPHERAL_HOST QIView(size_t N1, double xmin, double xmax, double xstep, double* vals);
+  SPHERAL_HOST_DEVICE ~QIView() { }
+};
 }
 
 #include "QuadraticInterpolatorInline.hh"