Clean up around BsplineAllocator.hpp and MultiSpline.hpp

src/QMCWaveFunctions/BsplineFactory/SplineC2C.cpp

@@ -32,8 +32,8 @@ inline void SplineC2C<ST>::set_spline(SingleSplineType* spline_r,
                                       int ispline,
                                       int level)
 {
-  SplineInst->copy_spline(spline_r, 2 * ispline);
-  SplineInst->copy_spline(spline_i, 2 * ispline + 1);
+  copy_spline<double, ST>(*spline_r, *SplineInst->getSplinePtr(), 2 * ispline);
+  copy_spline<double, ST>(*spline_i, *SplineInst->getSplinePtr(), 2 * ispline + 1);
 }
 
 template<typename ST>

src/QMCWaveFunctions/BsplineFactory/SplineC2C.h

@@ -134,8 +134,8 @@ class SplineC2C : public BsplineSet
     gatherv(comm, SplineInst->getSplinePtr(), SplineInst->getSplinePtr()->z_stride, offset);
   }
 
-  template<typename GT, typename BCT>
-  void create_spline(GT& xyz_g, BCT& xyz_bc)
+  template<typename BCT>
+  void create_spline(const Ugrid xyz_g[3], const BCT& xyz_bc)
   {
     resize_kpoints();
     SplineInst = std::make_shared<MultiBspline<ST>>();

src/QMCWaveFunctions/BsplineFactory/SplineC2COMPTarget.cpp

@@ -30,8 +30,8 @@ inline void SplineC2COMPTarget<ST>::set_spline(SingleSplineType* spline_r,
                                                int ispline,
                                                int level)
 {
-  SplineInst->copy_spline(spline_r, 2 * ispline);
-  SplineInst->copy_spline(spline_i, 2 * ispline + 1);
+  copy_spline<double, ST>(*spline_r, *SplineInst->getSplinePtr(), 2 * ispline);
+  copy_spline<double, ST>(*spline_i, *SplineInst->getSplinePtr(), 2 * ispline + 1);
 }
 
 template<typename ST>

src/QMCWaveFunctions/BsplineFactory/SplineC2COMPTarget.h

@@ -169,8 +169,8 @@ class SplineC2COMPTarget : public BsplineSet
     gatherv(comm, SplineInst->getSplinePtr(), SplineInst->getSplinePtr()->z_stride, offset);
   }
 
-  template<typename GT, typename BCT>
-  void create_spline(GT& xyz_g, BCT& xyz_bc)
+  template<typename BCT>
+  void create_spline(const Ugrid xyz_g[3], const BCT& xyz_bc)
   {
     resize_kpoints();
     SplineInst = std::make_shared<MultiBspline<ST, OffloadAllocator<ST>>>();

src/QMCWaveFunctions/BsplineFactory/SplineC2R.cpp

@@ -33,8 +33,8 @@ inline void SplineC2R<ST>::set_spline(SingleSplineType* spline_r,
                                       int ispline,
                                       int level)
 {
-  SplineInst->copy_spline(spline_r, 2 * ispline);
-  SplineInst->copy_spline(spline_i, 2 * ispline + 1);
+  copy_spline<double, ST>(*spline_r, *SplineInst->getSplinePtr(), 2 * ispline);
+  copy_spline<double, ST>(*spline_i, *SplineInst->getSplinePtr(), 2 * ispline + 1);
 }
 
 template<typename ST>

src/QMCWaveFunctions/BsplineFactory/SplineC2R.h

@@ -120,8 +120,8 @@ class SplineC2R : public BsplineSet
     gatherv(comm, SplineInst->getSplinePtr(), SplineInst->getSplinePtr()->z_stride, offset);
   }
 
-  template<typename GT, typename BCT>
-  void create_spline(GT& xyz_g, BCT& xyz_bc)
+  template<typename BCT>
+  void create_spline(const Ugrid xyz_g[3], const BCT& xyz_bc)
   {
     resize_kpoints();
     SplineInst = std::make_shared<MultiBspline<ST>>();

src/QMCWaveFunctions/BsplineFactory/SplineC2ROMPTarget.cpp

@@ -29,8 +29,8 @@ inline void SplineC2ROMPTarget<ST>::set_spline(SingleSplineType* spline_r,
                                                int ispline,
                                                int level)
 {
-  SplineInst->copy_spline(spline_r, 2 * ispline);
-  SplineInst->copy_spline(spline_i, 2 * ispline + 1);
+  copy_spline<double, ST>(*spline_r, *SplineInst->getSplinePtr(), 2 * ispline);
+  copy_spline<double, ST>(*spline_i, *SplineInst->getSplinePtr(), 2 * ispline + 1);
 }
 
 template<typename ST>

src/QMCWaveFunctions/BsplineFactory/SplineC2ROMPTarget.h

@@ -175,8 +175,8 @@ class SplineC2ROMPTarget : public BsplineSet
     gatherv(comm, SplineInst->getSplinePtr(), SplineInst->getSplinePtr()->z_stride, offset);
   }
 
-  template<typename GT, typename BCT>
-  void create_spline(GT& xyz_g, BCT& xyz_bc)
+  template<typename BCT>
+  void create_spline(const Ugrid xyz_g[3], const BCT& xyz_bc)
   {
     resize_kpoints();
     SplineInst = std::make_shared<MultiBspline<ST, OffloadAllocator<ST>>>();

src/QMCWaveFunctions/BsplineFactory/SplineR2R.cpp

@@ -32,7 +32,7 @@ inline void SplineR2R<ST>::set_spline(SingleSplineType* spline_r,
                                       int ispline,
                                       int level)
 {
-  SplineInst->copy_spline(spline_r, ispline);
+  copy_spline<double, ST>(*spline_r, *SplineInst->getSplinePtr(), ispline);
 }
 
 template<typename ST>

src/QMCWaveFunctions/BsplineFactory/SplineR2R.h

@@ -130,8 +130,8 @@ class SplineR2R : public BsplineSet
     gatherv(comm, SplineInst->getSplinePtr(), SplineInst->getSplinePtr()->z_stride, offset);
   }
 
-  template<typename GT, typename BCT>
-  void create_spline(GT& xyz_g, BCT& xyz_bc)
+  template<typename BCT>
+  void create_spline(const Ugrid xyz_g[3], const BCT& xyz_bc)
   {
     GGt        = dot(transpose(PrimLattice.G), PrimLattice.G);
     SplineInst = std::make_shared<MultiBspline<ST>>();

src/Sandbox/einspline_spo.hpp

@@ -173,7 +173,7 @@ struct einspline_spo
         einsplines[i]->create(grid, BC, nSplinesPerBlock);
         if (init_random)
           for (int j = 0; j < nSplinesPerBlock; ++j)
-            einsplines[i]->copy_spline(aspline, j);
+            copy_spline<double, T>(*aspline, *einsplines[i]->getSplinePtr(), j);
       }
       mAllocator.destroy(aspline);
     }

src/spline2/BsplineAllocator.hpp

@@ -255,34 +255,5 @@ typename BsplineAllocator<T, ALLOC>::SingleSplineType* BsplineAllocator<T, ALLOC
   return spline;
 }
 
-template<typename T, typename ALLOC>
-template<typename UBT, typename MBT>
-void BsplineAllocator<T, ALLOC>::copy(UBT* single, MBT* multi, int i, const int* offset, const int* N)
-{
-  using out_type        = typename bspline_type<MBT>::value_type;
-  using in_type         = typename bspline_type<UBT>::value_type;
-  intptr_t x_stride_in  = single->x_stride;
-  intptr_t y_stride_in  = single->y_stride;
-  intptr_t x_stride_out = multi->x_stride;
-  intptr_t y_stride_out = multi->y_stride;
-  intptr_t z_stride_out = multi->z_stride;
-  intptr_t offset0      = static_cast<intptr_t>(offset[0]);
-  intptr_t offset1      = static_cast<intptr_t>(offset[1]);
-  intptr_t offset2      = static_cast<intptr_t>(offset[2]);
-  const intptr_t istart = static_cast<intptr_t>(i);
-  const intptr_t n0 = N[0], n1 = N[1], n2 = N[2];
-  for (intptr_t ix = 0; ix < n0; ++ix)
-    for (intptr_t iy = 0; iy < n1; ++iy)
-    {
-      out_type* restrict out = multi->coefs + ix * x_stride_out + iy * y_stride_out + istart;
-      const in_type* restrict in =
-          single->coefs + (ix + offset0) * x_stride_in + (iy + offset1) * y_stride_in + offset2;
-      for (intptr_t iz = 0; iz < n2; ++iz)
-      {
-        out[iz * z_stride_out] = static_cast<out_type>(in[iz]);
-      }
-    }
-}
-
 } // namespace qmcplusplus
 #endif

src/spline2/MultiBspline.hpp

@@ -40,6 +40,8 @@ class MultiBspline
   using SplineType = typename bspline_traits<T, 3>::SplineType;
   ///define the real type
   using real_type = typename bspline_traits<T, 3>::real_type;
+  ///define the boundary condition type
+  using BoundaryCondition = typename bspline_traits<T, 3>::BCType;
   ///actual einspline multi-bspline object
   SplineType* spline_m;
   ///use allocator
@@ -59,21 +61,20 @@ class MultiBspline
   SplineType* getSplinePtr() { return spline_m; }
 
   /** create the einspline as used in the builder
-   * @tparam GT grid type
    * @tparam BCT boundary type
    * @param bc num_splines number of splines
    *
    * num_splines must be padded to the aligned size. The caller must be aware of padding and pad all result arrays.
    */
-  template<typename GT, typename BCT>
-  void create(GT& grid, BCT& bc, int num_splines)
+  template<typename BCT>
+  inline void create(const Ugrid grid[3], const BCT& bc, int num_splines)
   {
     static_assert(std::is_same<T, typename ALLOC::value_type>::value, "MultiBspline and ALLOC data types must agree!");
     if (getAlignedSize<T, ALLOC::alignment>(num_splines) != num_splines)
       throw std::runtime_error("When creating the data space of MultiBspline, num_splines must be padded!\n");
     if (spline_m == nullptr)
     {
-      typename bspline_traits<T, 3>::BCType xBC, yBC, zBC;
+      BoundaryCondition xBC, yBC, zBC;
       xBC.lCode = bc[0].lCode;
       yBC.lCode = bc[1].lCode;
       zBC.lCode = bc[2].lCode;
@@ -101,26 +102,41 @@ class MultiBspline
   int num_splines() const { return (spline_m == nullptr) ? 0 : spline_m->num_splines; }
 
   size_t sizeInByte() const { return (spline_m == nullptr) ? 0 : spline_m->coefs_size * sizeof(T); }
+};
+
 
-  /** copy a single spline to the big table
-   * @tparam SingleSpline single spline type
-   * @param aSpline UBspline_3d_(d,s)
-   * @param int index of aSpline
+/** copy a single spline to multi spline
+   * @tparam SSDT single spline coefficient data type
+   * @tparam MSDT multi spline coefficient data type
+   * @param single UBspline_3d_(d,s)
+   * @param multi multi_UBspline_3d_(d,s)
+   * @param int index of single in multi
    */
-  template<typename SingleSpline>
-  void copy_spline(SingleSpline* aSpline, int i)
-  {
-    if (spline_m == nullptr)
-      throw std::runtime_error("The internal storage of MultiBspline must be created first!\n");
-    if (aSpline->x_grid.num != spline_m->x_grid.num || aSpline->y_grid.num != spline_m->y_grid.num ||
-        aSpline->z_grid.num != spline_m->z_grid.num)
-      throw std::runtime_error("Cannot copy a single spline to MultiSpline with a different grid!\n");
+template<typename SSDT, typename MSDT>
+void copy_spline(const typename bspline_traits<SSDT, 3>::SingleSplineType& single,
+                 typename bspline_traits<MSDT, 3>::SplineType& multi,
+                 int i)
+{
+  if (single.x_grid.num != multi.x_grid.num || single.y_grid.num != multi.y_grid.num ||
+      single.z_grid.num != multi.z_grid.num)
+    throw std::runtime_error("Cannot copy a single spline to MultiSpline with a different grid!\n");
 
-    const int BaseOffset[3] = {0, 0, 0};
-    const int BaseN[3]      = {spline_m->x_grid.num + 3, spline_m->y_grid.num + 3, spline_m->z_grid.num + 3};
-    myAllocator.copy(aSpline, spline_m, i, BaseOffset, BaseN);
-  }
-};
+  intptr_t x_stride_in  = single.x_stride;
+  intptr_t y_stride_in  = single.y_stride;
+  intptr_t x_stride_out = multi.x_stride;
+  intptr_t y_stride_out = multi.y_stride;
+  intptr_t z_stride_out = multi.z_stride;
+  const intptr_t istart = static_cast<intptr_t>(i);
+  const intptr_t n0 = multi.x_grid.num + 3, n1 = multi.y_grid.num + 3, n2 = multi.z_grid.num + 3;
+  for (intptr_t ix = 0; ix < n0; ++ix)
+    for (intptr_t iy = 0; iy < n1; ++iy)
+    {
+      auto* __restrict__ out      = multi.coefs + ix * x_stride_out + iy * y_stride_out + istart;
+      const auto* __restrict__ in = single.coefs + ix * x_stride_in + iy * y_stride_in;
+      for (intptr_t iz = 0; iz < n2; ++iz)
+        out[iz * z_stride_out] = in[iz];
+    }
+}
 
 } // namespace qmcplusplus
 

src/spline2/tests/test_multi_spline.cpp

@@ -166,7 +166,7 @@ struct test_splines : public test_splines_base<T, GRID_SIZE, NUM_SPLINES>
     UBspline_3d_d* aspline = mAllocator.allocateUBspline(grid[0], grid[1], grid[2], bc[0], bc[1], bc[2], data.data());
 
     for (int i = 0; i < num_splines; i++)
-      bs.copy_spline(aspline, i);
+      copy_spline<double, T>(*aspline, *bs.getSplinePtr(), i);
 
     mAllocator.destroy(aspline);
 
@@ -223,7 +223,7 @@ struct test_splines<T, 5, 1> : public test_splines_base<T, 5, 1>
     UBspline_3d_d* aspline = mAllocator.allocateUBspline(grid[0], grid[1], grid[2], bc[0], bc[1], bc[2], data.data());
 
     for (int i = 0; i < num_splines; i++)
-      bs.copy_spline(aspline, i);
+      copy_spline<double, T>(*aspline, *bs.getSplinePtr(), i);
 
     mAllocator.destroy(aspline);