SYCL related tweaks. Spatial model now runnign on GPU, albeit no fast…

…er than on CPU.

include/Makefile

@@ -108,6 +108,8 @@ else
 ifdef DPCPP
 CPLUSPLUS=icpx
 SYCL+=-fsycl
+# -std=c++20 used to suppress warning messages
+FLAGS+=-std=c++20
 else
 CPLUSPLUS=g++
 endif

include/cairo_base.h

@@ -51,6 +51,7 @@ namespace ecolab
     /// Container object for managing cairo_surface_t lifetimes
     class Surface
     {
+      CLASSDESC_ACCESS(Surface);
     private:
       cairo_surface_t* m_surface;
       cairo_t* m_cairo;
@@ -200,7 +201,7 @@ namespace ecolab
 
       double width() const override {return imageBlock.width;}
       double height() const override {return imageBlock.height;}
-      void clear() {
+      void clear() override {
         if (imageBlock.transparency)
           // make it transparent
           memset(imageData.data(),0,imageData.size());
@@ -228,7 +229,7 @@ namespace ecolab
                   width, height, imageBlock.pitch));
       }
 
-      void resize(size_t width, size_t height)
+      void resize(size_t width, size_t height) override
       {
 //#if TK_MAJOR_VERSION == 8 && TK_MINOR_VERSION < 5
 //        Tk_PhotoSetSize(photo, width, height);
@@ -239,7 +240,7 @@ namespace ecolab
       }
 
       /// copy the cairo output into the TkPhoto
-      void blit(unsigned x, unsigned y, unsigned width, unsigned height) 
+      void blit(unsigned x, unsigned y, unsigned width, unsigned height) override
       {
 //        imageBlock.pixelPtr=&imageData[x*imageBlock.pixelSize + y*imageBlock.pitch];
 //        if (width != unsigned(imageBlock.width) || 
@@ -255,13 +256,13 @@ namespace ecolab
 //          compositing? TK_PHOTO_COMPOSITE_OVERLAY:TK_PHOTO_COMPOSITE_SET);
 //#endif
       }
-      void blit() {blit(0,0,imageBlock.width, imageBlock.height);}
-      void requestRedraw() {blit();}
+      void blit() override {blit(0,0,imageBlock.width, imageBlock.height);}
+      void requestRedraw() override {blit();}
     };
 
     class CairoImage
     {
-
+      CLASSDESC_ACCESS(CairoImage);
     protected:
       void resize(size_t width, size_t height);
 

include/ecolab.h

@@ -11,6 +11,11 @@
 */
 #ifndef ECOLAB_H
 #define ECOLAB_H
+
+#ifdef SYCL_LANGUAGE_VERSION
+#define REALLOC reallocSycl
+#endif
+
 #include <stdlib.h>
 #include "pythonBuffer.h"
 #include "graphcode.h"
@@ -89,9 +94,37 @@ namespace ecolab
 
 #ifdef SYCL_LANGUAGE_VERSION
   sycl::queue& syclQ();
+  void* reallocSycl(void*,size_t);
 #endif
 
-  template <class Cell> struct EcolabGraph: public graphcode::Graph<Cell>
+  template <class T>
+  struct SyclType: public T
+  {
+#ifdef SYCL_LANGUAGE_VERSION
+    void* operator new(size_t s) {return reallocSycl(nullptr,s);}
+    void operator delete(void* p) {reallocSycl(p,0);}
+    void* operator new[](size_t s) {return reallocSycl(nullptr,s);}
+    void operator delete[](void* p) {reallocSycl(p,0);}
+#endif
+  };
+
+  template <class M>
+  struct DeviceType
+  {
+    using element_type=M;
+    SyclType<M>* const model=new SyclType<M>;
+    DeviceType()=default;
+    DeviceType(const DeviceType& x) {*this=x;}
+    DeviceType& operator=(const DeviceType& x) {*model=x.*model; return *this;}
+    ~DeviceType() {delete model;}
+    M& operator*() {return *model;}
+    const M& operator*() const {return *model;}
+    M* operator->() {return model;}
+    const M* operator->() const {return model;}
+    operator bool() const {return true;} // always defined
+  };
+
+  template <class Model, class Cell> struct EcolabGraph: public graphcode::Graph<Cell>
   {
     /// apply a functional to all local cells of this processor in parallel
     /// @param f 
@@ -110,6 +143,11 @@ namespace ecolab
   };
 }
 
+namespace classdesc
+{
+  template <class M> struct is_smart_ptr<ecolab::DeviceType<M>>: public true_type {}; 
+}
+
 #ifdef MPI_SUPPORT
 #include <classdescMP.h>
 namespace ecolab

include/graph.h

@@ -174,16 +174,16 @@ namespace ecolab
   {
     GraphAdaptor_const_iterator_base(const typename G::const_iterator& i): 
       G::const_iterator(i) {}
-    Edge operator*() const {return G::const_iterator::operator*();}
-    const Graph::const_iterator_base& operator++() {
+    Edge operator*() const override {return G::const_iterator::operator*();}
+    const Graph::const_iterator_base& operator++() override {
       G::const_iterator::operator++(); return *this;}
     bool operator==(const GraphAdaptor_const_iterator_base& x) const {
       return static_cast<const typename G::const_iterator&>(*this)==x;
     }
-    bool operator==(const Graph::const_iterator_base& x) const {
+    bool operator==(const Graph::const_iterator_base& x) const override {
       const GraphAdaptor_const_iterator_base* xp=
         dynamic_cast<const GraphAdaptor_const_iterator_base*>(&x);
-      return xp && operator==(*xp);
+      return xp && static_cast<const typename G::const_iterator&>(*this)==static_cast<const typename G::const_iterator&>(*xp);
     }
   };
 

include/pack.h

@@ -13,5 +13,8 @@
 #include "pack_base.h"
 //#include "ref.h"
 #include "pack_stl.h"
+#ifndef SYCL_LANGUAGE_VERSION
+// SYCL cannot call recursive functions, so serialisation of graph structures is impossible
 #include "pack_graph.h"
 #endif
+#endif

include/plot.h

@@ -206,9 +206,9 @@ namespace ecolab
     //void plot(TCL_args args); 
 
     //serialisation support (surface is not auto-serialisable)
-    void pack(classdesc::pack_t& p) const {p<<m_image<<x<<y<<minx<<maxx<<miny<<maxy;}
+    void pack(classdesc::pack_t& p) const override {p<<m_image<<x<<y<<minx<<maxx<<miny<<maxy;}
 
-    void unpack(classdesc::pack_t& p) {
+    void unpack(classdesc::pack_t& p) override {
         p>>m_image>>x>>y>>minx>>maxx>>miny>>maxy;
         Image(m_image);
     }

include/random.h

@@ -18,8 +18,8 @@
 
 #include "classdesc_access.h"
 #include "pack_base.h"
-#include "pack_stl.h"
-#include "pack_graph.h"
+//#include "pack_stl.h"
+//#include "pack_graph.h"
 
 #include <stdlib.h>
 #include <math.h>

lib/GUI.py

@@ -4,13 +4,19 @@
 registerParallel()
 from ecolab import Parallel, myid
 from objectBrowser import Browser
+from math import inf
+from time import perf_counter
 
 # list of windows to update each simulation step
 windows=[]
 runner=Tk()
 runner.title('EcoLab')
-statusBar=ttk.Label(runner,text='Not started')
-statusBar.pack()
+statusFrame=ttk.Frame(runner)
+statusFrame.pack()
+statusBar=ttk.Label(statusFrame,text='Not started')
+statusBar.pack(side='left')
+performance=ttk.Label(statusFrame,text='perf: 0')
+performance.pack(side='left')
 
 class Simulator:
     running=False
@@ -21,7 +27,10 @@ def __init__(self,step):
     def __call__(self):
         self.running=True
         while self.running:
+            start=perf_counter()
             self.doStep()
+            timePerStep=perf_counter()-start
+            performance.configure(text='perf: %f'%(1.0/timePerStep))
             for i in windows:
                 win=i()
                 if win: win.update()

models/ecolab_model.cc

@@ -39,7 +39,7 @@ namespace model
 {
   PanmicticModel panmictic_ecolab;
   CLASSDESC_ADD_GLOBAL(panmictic_ecolab);
-  SpatialModel spatial_ecolab;
+  DeviceType<SpatialModel> spatial_ecolab;
   CLASSDESC_ADD_GLOBAL(spatial_ecolab);
   CLASSDESC_PYTHON_MODULE(ecolab_model);
 }
@@ -236,14 +236,14 @@ void do_row_or_col(array<double>& tmp, double range, double minval, double gdist
   if (r>0)
     for (unsigned j=0; j<ntrue; j++)
       {
-	pos = (int)((tmp.size()-1) * ((float) rand()/RAND_MAX) +.5);
+	pos = (int)((tmp.size()-1) * ((double) rand()/RAND_MAX) +.5);
 	if (tmp[pos]==0.0)
-	  tmp[pos] = range *((float) rand()/RAND_MAX) + minval;
+	  tmp[pos] = range *((double) rand()/RAND_MAX) + minval;
       }
   else
     for (unsigned j=0; j<ntrue; j++)
       {
-	pos = (int)((tmp.size()-1) * ((float) rand()/RAND_MAX) +.5);
+	pos = (int)((tmp.size()-1) * ((double) rand()/RAND_MAX) +.5);
 	tmp[pos]=0;
       }
 
@@ -513,7 +513,7 @@ void SpatialModel::migrate()
   MPI_Reduce(&mig,&m,1,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD);
   mig=m;
 #endif
-  if (myid==0) assert(sum(ssum==0)==int(ssum.size()));
+  if (myid()==0) assert(sum(ssum==0)==int(ssum.size()));
 #endif
 
 }

models/ecolab_model.h

@@ -87,7 +87,8 @@ struct PanmicticModel: public ModelData, public EcolabPoint, public ecolab::Mode
 
 struct EcolabCell: public EcolabPoint, public graphcode::Object<EcolabCell> {};
 
-class SpatialModel: public ModelData, public EcolabGraph<EcolabCell>, public ecolab::Model<SpatialModel>
+class SpatialModel: public ModelData, public EcolabGraph<SpatialModel, EcolabCell>,
+                    public ecolab::Model<SpatialModel>
 {
   size_t numX=1, numY=1;
   CLASSDESC_ACCESS(SpatialModel);

models/spatial_ecolab.py

@@ -1,7 +1,10 @@
-from ecolab_model import spatial_ecolab as ecolab
+from ecolab_model import spatial_ecolab
+# spatial_ecolab is a smart ptr to possible Device accessible memery, so must be dereferenced
+ecolab=spatial_ecolab()
+
 from random import random, seed as randomSeed
 
-from ecolab import array_urand, myid
+from ecolab import array_urand, myid, device
 
 # we want initialisation to be identical across all processes
 randomSeed(1)
@@ -45,8 +48,10 @@ def randomList(num, min, max):
 from plot import plot
 from GUI import gui, statusBar, windows
 
+print(device())
+
 def step():
-    ecolab.generate()
+    ecolab.generate(100)
     ecolab.mutate()
     ecolab.migrate()
     ecolab.condense()

src/ecolab.cc

@@ -30,15 +30,34 @@ queue& ecolab::syclQ() {
   static queue q{default_selector_v};
   return q;
 }
-void* operator new(size_t s) {
-  cout<<s<<" bytes allocated shared"<<std::endl;
-  if (auto r=malloc_shared(s,syclQ()))
-    return r;
-  throw std::bad_alloc();
+
+void* ecolab::reallocSycl(void* pp,size_t s)
+{
+  size_t* p=reinterpret_cast<size_t*>(pp);
+  size_t prevSz=0;
+  if (p)
+    prevSz=*--p;
+
+  if (!s && p)
+    {
+      free(p);
+      return nullptr;
+    }
+
+  if (prevSz>=s) return pp; // nothing further to do
+  auto r=reinterpret_cast<size_t*>(malloc_shared(s+sizeof(size_t),syclQ()));
+  if (r)
+    {
+      // stash allocation size at beginning of allocated region
+      *r++=s;
+      if (p)
+        {
+          memcpy(r,pp,prevSz);
+          free(p,syclQ());
+        }
+    }
+  return r;
 }
-void operator delete(void* p) {return free(p,syclQ());}
-void* operator new[](size_t s) {return operator new(s);}
-void operator delete[](void* p) {return free(p,syclQ());}
 #endif  
 
 namespace ecolab
@@ -174,7 +193,6 @@ namespace ecolab
     static ParallelType parallelType;
     PyModule_AddObject(pythonModule,"Parallel",reinterpret_cast<PyObject*>(&parallelType));
   }
-
   CLASSDESC_ADD_FUNCTION(registerParallel);
 
 
@@ -185,6 +203,17 @@ namespace ecolab
   CLASSDESC_ADD_FUNCTION(myid);
   CLASSDESC_ADD_FUNCTION(nprocs);
   CLASSDESC_DECLARE_TYPE(Plot);
+
+  // device SYCL kernels are running on
+  std::string device() {
+#ifdef SYCL_LANGUAGE_VERSION
+    return syclQ().get_device().get_info<info::device::name>();
+#else
+    return "CPU";
+#endif
+  }
+  CLASSDESC_ADD_FUNCTION(device);
+
 }
 
 CLASSDESC_PYTHON_MODULE(ecolab);