generate step now running, although all cells behave identitically...

highperformancecoder · Dec 10, 2024 · 5b2c222 · 5b2c222
1 parent 77b55bc
commit 5b2c222
Show file tree

Hide file tree

Showing 8 changed files with 109 additions and 72 deletions.
diff --git a/include/arrays.h b/include/arrays.h
@@ -1629,6 +1629,12 @@ namespace ecolab
       ~array() {release();}
 
       const Allocator& allocator() const {return m_allocator;}
+      const Allocator& allocator(const Allocator& alloc) {
+        array tmp(size(),alloc);
+        asg_v(tmp.data(),size(),data());
+        swap(tmp);
+        return m_allocator;
+      }
 
       /// resize array to \a s elements
       void resize(size_t s) {
@@ -1648,6 +1654,7 @@ namespace ecolab
 
       void swap(array& x) {
         std::swap(dt, x.dt);
+        std::swap(m_allocator,x.m_allocator);
       }
 
       T& operator[](size_t i) {assert(i<size()); copy(); return data()[i];}
@@ -1679,8 +1686,9 @@ namespace ecolab
       operator=(const expr& x) {
         if ((void*)(&x)==(void*)(this)) return *this;
         // since expression x may contain a reference to this, assign to a temporary
-        resize(x.size());
-        array_ns::asg_v(data(),size(),x);
+        array tmp(x.size(),m_allocator);
+        array_ns::asg_v(tmp.data(),tmp.size(),x);
+        swap(tmp);
         return *this;
       }
       template <class expr> typename
@@ -2453,10 +2461,17 @@ namespace ecolab
 
     /// fill array with uniformly random numbers from [0,1)
     template <class F> void fillrand(array<F>& x);
+    template <class F, class A> void fillrand(array<F, A>& x)
+    {array<F> tmp(x.size()); fillrand(tmp); x=tmp;}
     /// fill array with gaussian random numbers from \f$N(0,1)\propto\exp(-x^2/2)\f$
     template <class F> void fillgrand(array<F>& x);
+    template <class F, class A> void fillgrand(array<F,A>& x)
+    {array<F> tmp(x.size()); fillgrand(tmp); x=tmp;}
     /// fill array with exponential random numbers \f$x\leftarrow-\ln\xi,\,\xi\in [0,1)\f$
     template <class F> void fillprand(array<F>& x);
+    template <class F, class A> void fillprand(array<F,A>& x)
+    {array<F> tmp(x.size()); fillprand(tmp); x=tmp;}
+
     /// fill with uniform numbers drawn from [0...\a max] without replacement
     void fill_unique_rand(array<int>& x, unsigned max);
 

diff --git a/include/ecolab.h b/include/ecolab.h
@@ -20,6 +20,14 @@
 #include "device/Ouroboros_impl.dp.hpp"
 #include "InstanceDefinitions.dp.hpp"
 #include "device/MemoryInitialization.dp.hpp"
+#if 1 //def __INTEL_LLVM_COMPLER
+#warning "Engaging experimental printf"
+template <typename... Args>
+void syclPrintf(Args... args) {sycl::ext::oneapi::experimental::printf(args...);}
+sycl::nd_item<1> syclItem() {return sycl::ext::oneapi::this_work_item::get_nd_item<1>();}
+#else
+#define syclPrintf(...)
+#endif
 #endif
 
 #include <stdlib.h>
@@ -137,22 +145,19 @@ namespace ecolab
     {
     public:
       Ouro::SyclDesc<1>*const* desc=nullptr;
-      MemAllocator*const* memAlloc;
-      const sycl::stream* out=nullptr;
+      MemAllocator*const* memAlloc=nullptr;
       template <class U> friend class Allocator;
       CellAllocator()=default;
-      CellAllocator(Ouro::SyclDesc<1>* const& desc, MemAllocator*const& memAlloc, const sycl::stream* out):
-        desc(&desc), memAlloc(&memAlloc), out(out) {
+      CellAllocator(Ouro::SyclDesc<1>* const& desc, MemAllocator*const& memAlloc):
+        desc(&desc), memAlloc(&memAlloc) {
       }
       template <class U> CellAllocator(const CellAllocator<U>& x):
-        desc(x.desc) {}
+        desc(x.desc), memAlloc(x.memAlloc) {}
       T* allocate(size_t sz) {
         if (memAlloc && *memAlloc && desc && *desc)  {
           auto r=reinterpret_cast<T*>((*memAlloc)->malloc(**desc,sz*sizeof(T)));
-          if (out) (*out)<<"allocated "<<sz*sizeof(T)<<" bytes, "<<r<<sycl::endl;
           return r;
         }
-        if (out) (*out)<<"failed allocation "<<sz*sizeof(T)<<sycl::endl;
         return nullptr; // TODO raise an error??
       }
       void deallocate(T* p,size_t) {
@@ -162,14 +167,22 @@ namespace ecolab
       bool operator==(const CellAllocator& x) const {return desc==x.desc && memAlloc==x.memAlloc;}
     };
     template <class T> CellAllocator<T> allocator() const {
-      return CellAllocator<T>(desc,memAlloc,out);
+      return CellAllocator<T>(desc,memAlloc);
     }
 #else
     template <class T> using CellAllocator=std::allocator<T>;
     template <class T> CellAllocator<T> allocator() const {return CellAllocator<T>();}
 #endif
   };
 
+  template <class T>
+  void printAllocator(const char* prefix, const T&) {}
+
+  template <class T>
+  void printAllocator(const char* prefix, const CellBase::CellAllocator<T>& x)
+  {syclPrintf("%s: allocator.desc=%p memAlloc=%p\n",prefix,x.desc,x.memAlloc);}
+
+
   class SyclGraphBase
   {
   protected:
@@ -215,7 +228,6 @@ namespace ecolab
         h.parallel_for(sycl::nd_range<1>(range*workGroupSize, workGroupSize), [=,this](auto i) {
           auto idx=i.get_global_linear_id();
           if (idx<this->size()) {
-            out<<"idx="<<idx<<"(*this)[idx]"<<(*this)[idx].payload->get()<<" as "<<(*this)[idx]->template as<Cell>()<<sycl::endl;
             auto& cell=*(*this)[idx]->template as<Cell>();
             Ouro::SyclDesc<> desc(i,{});
             cell.desc=&desc;

diff --git a/include/graph.h b/include/graph.h
@@ -408,8 +408,8 @@ bool GraphAdaptor<BiDirectionalGraph>::directed() const {return false;}
 class sparse_mat_graph: public ConcreteGraph<DiGraph>
 {
 public:
-  template <class F>
-  const sparse_mat_graph& operator=(const sparse_mat<F>& mat) {
+  template <class F, template<class T> class A>
+  const sparse_mat_graph& operator=(const sparse_mat<F,A>& mat) {
     clear();
     std::map<std::pair<size_t, size_t>, double > weights;
     for (size_t i=0; i<mat.row.size(); ++i) {

diff --git a/include/sparse_mat.h b/include/sparse_mat.h
@@ -14,6 +14,8 @@
 
 #include <arrays.h>
 
+using sycl::ext::oneapi::experimental::printf;
+
 namespace ecolab
 {
   /// sparse matrix class
@@ -37,6 +39,13 @@ namespace ecolab
       diag=x.diag; val=x.val; 
       row=x.row; col=x.col;
     }
+    void setAllocators(const A<unsigned>& ialloc, const A<F>& falloc) {
+      diag.allocator(falloc);
+      val.allocator(falloc);
+      row.allocator(ialloc);
+      col.allocator(ialloc);
+    }
+
     /*matrix multiplication*/
     template <class E> typename
     array_ns::enable_if
@@ -52,15 +61,18 @@ namespace ecolab
     {
       auto alloc=array_ns::makeAllocator<F>(x.allocator());
       array_ns::array<F,decltype(alloc)> r(alloc);
-      assert(row.size()==col.size() && row.size()==val.size());      
+      assert(row.size()==col.size() && row.size()==val.size());
       if (diag.size()>0)
 	{
 	  assert(diag.size()==x.size());
-	  r = diag*x; 
+	  r = diag*x;
+          assert(r.size()==x.size());
 	}
       else
 	r.resize(rowsz,0);
+      assert(r.size()>max(row) && x.size()>max(col));
       r[row]+=val*x[col];
+
       return r;
     }
     /* initialise sparse mat in a random configuration */

diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc
@@ -64,21 +64,11 @@ int EcolabPoint<B>::ROUND(Float x)
   Float dum;
   if (x<0) x=0;
   if (x>std::numeric_limits<int>::max()-1) x=std::numeric_limits<int>::max()-1;
+  //syclPrintf("ROUND inner x=%g, modf=%g, rand()=%g, rand.max=%g, rand.min=%g\n",x,std::fabs(std::modf(x,&dum)),rand(),rand.max(),rand.min());
   return std::fabs(std::modf(x,&dum))*(rand.max()-rand.min()) > (rand()-rand.min()) ?
     (int)x+1 : (int)x;
 }
 
-//inline int ROUND(float x) {return ROUND(double(x));}
-
-//  template <class E>
-//  inline array<int> EcolabPoint::RoundArray(const E& x)
-//  {
-//    array<int> r(x.size());
-//    for (size_t i=0; i<x.size(); i++)
-//      r[i]=ROUND(x[i]);
-//    return r;
-//  }
-
 template <class E, class P>
 struct RoundArray
 {
@@ -87,21 +77,27 @@ struct RoundArray
   RoundArray(P& point, const E& expr): expr(expr), point(point) {}
   using value_type=int;
   size_t size() const {return expr.size();}
-  int operator[](size_t i) const {return point.ROUND(expr[i]);}
+  int operator[](size_t i) const //{return point.ROUND(expr[i]);}
+  {
+    auto r=point.ROUND(expr[i]);
+    //syclPrintf("ROUND: %d, %g=%d\n",i,expr[i],r);
+    return r;
+  }
 };
 
 namespace ecolab::array_ns
 {template <class E, class P> struct is_expression<RoundArray<E,P>>: public true_type {};}
 
-template <class E, class P>
-RoundArray<E,P> roundArray(P& point, const E& expr)
-{return RoundArray<E,P>(point,expr);}
+template <class B>
+template <class E>
+RoundArray<E,EcolabPoint<B>> EcolabPoint<B>::roundArray(const E& expr)
+{return RoundArray<E,EcolabPoint<B>>(*this,expr);}
 
 template <class B>
 void EcolabPoint<B>::generate(unsigned niter, const ModelData& model)
-{ 
+{
   for (unsigned i=0; i<niter; i++)
-    {density += roundArray(*this, density * (model.repro_rate + model.interaction*density));}
+    {density = roundArray(density + density * (model.repro_rate + model.interaction*density));}
 }
 
 template <class B>
@@ -224,7 +220,7 @@ array<int> EcolabPoint<B>::mutate(const array<double>& mut_scale)
 {
   array<int> speciations;
   /* calculate the number of mutants each species produces */
-  speciations = roundArray(*this, mut_scale * density); 
+  speciations = roundArray(mut_scale * density); 
 
   /* generate index list of old species that mutate to the new */
   array<int> new_sp = gen_index(speciations); 
@@ -485,7 +481,9 @@ void SpatialModel::setGrid(size_t nx, size_t ny)
 void SpatialModel::generate(unsigned niter)
 {
   if (tstep==0) makeConsistent();
-  forAll([=,this](EcolabCell& c) {c.generate(niter,*this);});
+  forAll([=,this](EcolabCell& c) {
+    c.generate(niter,*this);
+  });
   tstep+=niter;
 }
 
@@ -540,6 +538,27 @@ void SpatialModel::migrate()
 
 }
 
+void ModelData::makeConsistent(size_t nsp)
+{
+#ifdef SYCL_LANGUAGE_VERSION
+  FAlloc falloc(syclQ(),sycl::usm::alloc::shared);
+  species.allocator(graphcode::Allocator<int>(syclQ(),sycl::usm::alloc::shared));
+  create.allocator(falloc);
+  repro_rate.allocator(falloc);
+  mutation.allocator(falloc);
+  migration.allocator(falloc);
+  interaction.setAllocators
+    (graphcode::Allocator<unsigned>(syclQ(),sycl::usm::alloc::shared),falloc);
+#endif
+
+  if (!species.size())
+    species=pcoord(nsp);
+
+  if (!create.size()) create.resize(species.size(),0);
+  if (!mutation.size()) mutation.resize(species.size(),0);
+  if (!migration.size()) migration.resize(species.size(),0);
+}
+
 
 void SpatialModel::makeConsistent()
 {
@@ -549,23 +568,13 @@ void SpatialModel::makeConsistent()
 #ifdef MPI_SUPPORT
   MPI_Allreduce(MPI_IN_PLACE,&nsp,1,MPI_UNSIGNED_LONG,MPI_MAX,MPI_COMM_WORLD);
 #endif
-  cout<<nsp<<endl;
   forAll([=,this](EcolabCell& c) {
     if (nsp>c.density.size())
       {
         array<int,EcolabCell::CellAllocator<int>> tmp(nsp,0,c.allocator<int>());
         asg_v(tmp.data(),c.density.size(),c.density);
-        //(*c.out)<<"tmp.size:"<<tmp.size()<<" "<<tmp.data()<<sycl::endl;
         c.density.swap(tmp);
-        //(*c.out)<<"c.density.size:"<<c.density.size()<<" "<<&c<<sycl::endl;
       }
   });
-  for (auto& i: *this)
-    {
-      auto& c=*i->as<EcolabCell>();
-      //cout<<i.id()<<" "<<c.density.size()<<" "<<&c<<" usm type:"<<
-      //  sycl::get_pointer_type(&c,syclQ().get_context())<<" "<<sycl::get_pointer_type(i.payload,syclQ().get_context())<<endl;
-
-    }
   ModelData::makeConsistent(nsp);
 }
diff --git a/models/ecolab_model.h b/models/ecolab_model.h
@@ -23,11 +23,11 @@ using Float=double;
 
 struct ConnectionPlot: public Object<ConnectionPlot, CairoSurface>
 {
-  array<int> density;
-  sparse_mat<Float> connections;
+  array<int,graphcode::Allocator<int>> density;
+  sparse_mat<Float, graphcode::Allocator> connections;
   bool redraw(int x0, int y0, int width, int height) override;
   void requestRedraw() const {if (surface) surface->requestRedraw();}
-  void update(const array<int>& d, const sparse_mat<Float>& c) {
+  void update(const array<int,graphcode::Allocator<int>>& d, const sparse_mat<Float, graphcode::Allocator>& c) {
     density=d;
     connections=c;
     requestRedraw();
@@ -36,33 +36,25 @@ struct ConnectionPlot: public Object<ConnectionPlot, CairoSurface>
 
 struct ModelData
 {
-  array<int> species;
-  array<Float> create;
-  array<Float> repro_rate, mutation, migration;
-  sparse_mat<Float> interaction;
+  using FAlloc=graphcode::Allocator<Float>;
+  array<int,graphcode::Allocator<int>> species;
+  array<Float,FAlloc> create;
+  array<Float,FAlloc> repro_rate, mutation, migration;
+  sparse_mat<Float,graphcode::Allocator> interaction;
   sparse_mat_graph foodweb;
   unsigned long long tstep=0, last_mut_tstep=0, last_mig_tstep=0;
   //mutation parameters
   float sp_sep=0.1, mut_max=0, repro_min=0, repro_max=1, odiag_min=0, odiag_max=1;
 
-  void makeConsistent(size_t nsp)
-    {
-      if (!species.size())
-        species=pcoord(nsp);
-
-      if (!create.size()) create.resize(species.size(),0);
-      if (!mutation.size()) mutation.resize(species.size(),0);
-      if (!migration.size()) migration.resize(species.size(),0);
-    }
-
+  void makeConsistent(size_t nsp);
   void random_interaction(unsigned conn, double sigma);
-
   void condense(const array<bool>& mask, size_t mask_true);
   void mutate(const array<int>&); 
-
   double complexity() {return ::complexity(foodweb);}
 };
 
+template <class E, class P> struct RoundArray; 
+
 /* ecolab cell  */
 template <class CellBase>
 struct EcolabPoint: public Exclude<CellBase>
@@ -74,8 +66,9 @@ struct EcolabPoint: public Exclude<CellBase>
   array<int> mutate(const array<double>&);
   unsigned nsp() const; ///< number of living species in this cell
   /// Rounding function, randomly round up or down, in the range 0..INT_MAX
-  int ROUND(Float x); 
-  template <class E> array<int> RoundArray(const E& x);
+  int ROUND(Float x);
+  template <class E> RoundArray<E,EcolabPoint> roundArray(const E& expr);
+  //template <class E> array<int> RoundArray(const E& x);
   Exclude<std::mt19937> rand; // random number generator
 };