Update SCTL and add example for sctl::ParticleFMM wrapper

MakeVariables.in

@@ -40,5 +40,5 @@ LDLIBS_PVFMM = @LDFLAGS@ $(PVFMM_LIBS) $(FFTW_LIBS_PVFMM) $(BLAS_LAPACK_LIB_PVFM
 LDFLAGS_PVFMM = $(LDLIBS_PVFMM) # Deprecated, use LDLIBS_PVFMM instead
 
 # Add SCTL
-CXXFLAGS_PVFMM += -DSCTL_HAVE_BLAS -DSCTL_HAVE_LAPACK -DSCTL_HAVE_FFTW -DSCTL_QUAD_T=__float128
+CXXFLAGS_PVFMM += -DSCTL_HAVE_PVFMM -DSCTL_HAVE_BLAS -DSCTL_HAVE_LAPACK -DSCTL_HAVE_FFTW -DSCTL_PROFILE=5 -DSCTL_VERBOSE
 CXXFLAGS_PVFMM += -I$(TOP_SRCDIR_PVFMM)/SCTL/include

examples/Makefile

@@ -25,6 +25,7 @@ INCDIR = ./include
 TARGET_BIN = \
        $(BINDIR)/example1 \
        $(BINDIR)/example2 \
+       $(BINDIR)/example-sctl \
        $(BINDIR)/fmm_pts \
        $(BINDIR)/fmm_cheb
 

examples/src/example-sctl.cpp

@@ -0,0 +1,78 @@
+#include "sctl.hpp"
+
+using namespace sctl;
+
+template <class Real> void test_particle_fmm(const Comm& comm) {
+  constexpr Integer DIM = 3;
+
+  Stokes3D_FSxU kernel_m2l;
+  Stokes3D_FxU kernel_sl;
+  Stokes3D_DxU kernel_dl;
+  srand48(comm.Rank());
+
+  // Create target and source vectors.
+  const Long N = 50000/comm.Size();
+  Vector<Real> trg_coord(N*DIM);
+  Vector<Real>  sl_coord(N*DIM);
+  Vector<Real>  dl_coord(N*DIM);
+  Vector<Real>  dl_norml(N*DIM);
+  for (auto& a : trg_coord) a = drand48();
+  for (auto& a :  sl_coord) a = drand48();
+  for (auto& a :  dl_coord) a = drand48();
+  for (auto& a :  dl_norml) a = drand48();
+  Long n_sl  =  sl_coord.Dim()/DIM;
+  Long n_dl  =  dl_coord.Dim()/DIM;
+
+  // Set source charges.
+  Vector<Real> sl_den(n_sl*kernel_sl.SrcDim());
+  Vector<Real> dl_den(n_dl*kernel_dl.SrcDim());
+  for (auto& a : sl_den) a = drand48() - 0.5;
+  for (auto& a : dl_den) a = drand48() - 0.5;
+
+  ParticleFMM<Real,DIM> fmm(comm);
+  fmm.SetAccuracy(10);
+  fmm.SetKernels(kernel_m2l, kernel_m2l, kernel_sl);
+  fmm.AddTrg("Potential", kernel_m2l, kernel_sl);
+  fmm.AddSrc("SingleLayer", kernel_sl, kernel_sl);
+  fmm.AddSrc("DoubleLayer", kernel_dl, kernel_dl);
+  fmm.SetKernelS2T("SingleLayer", "Potential",kernel_sl);
+  fmm.SetKernelS2T("DoubleLayer", "Potential",kernel_dl);
+
+  fmm.SetTrgCoord("Potential", trg_coord);
+  fmm.SetSrcCoord("SingleLayer", sl_coord);
+  fmm.SetSrcCoord("DoubleLayer", dl_coord, dl_norml);
+
+  fmm.SetSrcDensity("SingleLayer", sl_den);
+  fmm.SetSrcDensity("DoubleLayer", dl_den);
+
+  Vector<Real> Ufmm, Uref;
+  fmm.Eval(Ufmm, "Potential"); // Warm-up run
+  Ufmm = 0;
+
+  Profile::Enable(true);
+  Profile::Tic("FMM-Eval", &comm);
+  fmm.Eval(Ufmm, "Potential");
+  Profile::Toc();
+  Profile::Tic("Direct", &comm);
+  fmm.EvalDirect(Uref, "Potential");
+  Profile::Toc();
+  Profile::print(&comm);
+
+  Vector<Real> Uerr = Uref - Ufmm;
+  { // Print error
+    StaticArray<Real,2> loc_err{0,0}, glb_err{0,0};
+    for (const auto& a : Uerr) loc_err[0] = std::max<Real>(loc_err[0], fabs(a));
+    for (const auto& a : Uref) loc_err[1] = std::max<Real>(loc_err[1], fabs(a));
+    comm.Allreduce<Real>(loc_err, glb_err, 2, Comm::CommOp::MAX);
+    if (!comm.Rank()) std::cout<<"Maximum relative error: "<<glb_err[0]/glb_err[1]<<'\n';
+  }
+}
+
+int main(int argc, char** argv) {
+  sctl::Comm::MPI_Init(&argc, &argv);
+
+  test_particle_fmm<double>(sctl::Comm::World());
+
+  sctl::Comm::MPI_Finalize();
+  return 0;
+}

include/cheb_utils.hpp

@@ -9,13 +9,14 @@
 
 #include <pvfmm_common.hpp>
 #include <vector.hpp>
-#include <kernel.hpp>
 
 #ifndef _PVFMM_CHEB_UTILS_HPP_
 #define _PVFMM_CHEB_UTILS_HPP_
 
 namespace pvfmm{
 
+template <class Real> class Kernel;
+
 /**
  * \brief Returns the sum of the absolute value of coeffecients of the highest
  * order polynomial as an estimate of error.

include/cheb_utils.txx

@@ -17,6 +17,7 @@
 #include <mem_mgr.hpp>
 #include <matrix.hpp>
 #include <profile.hpp>
+#include <kernel.hpp>
 
 namespace pvfmm{
 

include/kernel.hpp

@@ -17,6 +17,10 @@
 #ifndef _PVFMM_FMM_KERNEL_HPP_
 #define _PVFMM_FMM_KERNEL_HPP_
 
+namespace sctl {
+template <class Real, long N> class Vec;
+};
+
 namespace pvfmm{
 
 template <class T>
@@ -47,7 +51,7 @@ struct Kernel{
   /**
    * \brief Constructor.
    */
-  Kernel(Ker_t poten, Ker_t dbl_poten, const char* name, int dim_, std::pair<int,int> k_dim,
+  Kernel(Ker_t poten = nullptr, Ker_t dbl_poten = nullptr, const char* name = "", int dim_ = 0, std::pair<int,int> k_dim = std::make_pair<int,int>(0,0),
          size_t dev_poten=(size_t)NULL, size_t dev_dbl_poten=(size_t)NULL);
 
   /**
@@ -94,11 +98,6 @@ struct Kernel{
   mutable const Kernel<T>* k_l2l;
   mutable const Kernel<T>* k_l2t;
   mutable VolPoten vol_poten;
-
-  private:
-
-  Kernel();
-
 };
 
 template<typename T, void (*A)(T*, int, T*, int, T*, int, T*, mem::MemoryManager* mem_mgr),
@@ -164,6 +163,20 @@ Kernel<T> BuildKernel(const char* name, int dim, std::pair<int,int> k_dim,
   return K;
 }
 
+template <class uKernel> class GenericKernel {
+  template <class VecType, int D, int K0, int K1> static constexpr int get_DIM  (void (*uKer)(VecType (&u)[K1], const VecType (&r)[D], const VecType (&f)[K0], const void* ctx_ptr)) { return D; }
+  template <class VecType, int D, int K0, int K1> static constexpr int get_KDIM0(void (*uKer)(VecType (&u)[K1], const VecType (&r)[D], const VecType (&f)[K0], const void* ctx_ptr)) { return K0; }
+  template <class VecType, int D, int K0, int K1> static constexpr int get_KDIM1(void (*uKer)(VecType (&u)[K1], const VecType (&r)[D], const VecType (&f)[K0], const void* ctx_ptr)) { return K1; }
+
+  static constexpr int DIM   = get_DIM  (uKernel::template uKerEval<sctl::Vec<double,1>,0>);
+  static constexpr int KDIM0 = get_KDIM0(uKernel::template uKerEval<sctl::Vec<double,1>,0>);
+  static constexpr int KDIM1 = get_KDIM1(uKernel::template uKerEval<sctl::Vec<double,1>,0>);
+
+  public:
+
+  template <class Real, int digits = -1> static void Eval(Real* r_src, int src_cnt, Real* v_src, int dof, Real* r_trg, int trg_cnt, Real* v_trg, mem::MemoryManager* mem_mgr);
+};
+
 }//end namespace
 
 #ifdef __INTEL_OFFLOAD