generate step kind of working under SYCL.

include/arrays.h

@@ -1561,6 +1561,7 @@ namespace ecolab
         // over allocate to allow for alignment and metadata
         auto allocation=n + (sizeof(array_data<T>) + 16)/sizeof(T)+1-array_data<T>::debug_display;
         p = m_allocator.allocate(allocation);
+        if (!p) return nullptr; // SYCL allocator returns nullptr if not initialised
 #ifdef __ICC
         // we need to align data onto 16 byte boundaries
         size_t d = (size_t)(reinterpret_cast<array_data<T>*>(p)->dt);
@@ -1637,7 +1638,7 @@ namespace ecolab
       array(const array& x): m_allocator(x.m_allocator) 
       {
         dt=x.dt;
-        ref()++;
+        if (dt) ref()++;
       }
 
       template <class expr>
@@ -1659,7 +1660,7 @@ namespace ecolab
             release();
             dt = alloc(s);
           } 
-        dt->sz=s;
+        if (dt) dt->sz=s; // in case s is smaller
       } 
 
       void clear() {resize(0);}
@@ -1688,7 +1689,7 @@ namespace ecolab
         if (m_allocator==x.m_allocator) { // shared data optimisation
           release();
           dt=x.dt;
-          ref()++;
+          if (dt) ref()++;
           return *this;
         } 
         array tmp(x.size(),m_allocator);
@@ -1835,7 +1836,7 @@ namespace ecolab
 
       /// returns a writeable pointer to data without copy-on-write semantics
       /// dangerous, but needed to run array expressions on GPU
-      T* dataNoCow() {return dt? dt->dt: 0;}
+      //      T* dataNoCow() {return dt? dt->dt: 0;}
 
       typedef T *iterator;
       typedef const T *const_iterator;

include/ecolab.h

@@ -195,7 +195,7 @@ namespace ecolab
       static size_t workGroupSize=syclQ().get_device().get_info<sycl::info::device::max_work_group_size>();
       size_t range=this->size()/workGroupSize;
       if (range*workGroupSize < this->size()) ++range;
-      syclQ().parallel_for(sycl::nd_range<1>(range, workGroupSize), [=,this](auto i) {
+      syclQ().parallel_for(sycl::nd_range<1>(range*workGroupSize, workGroupSize), [=,this](auto i) {
         auto idx=i.get_global_linear_id();
         if (idx<this->size()) {
           auto& cell=*(*this)[idx]->template as<Cell>();

include/sparse_mat.h

@@ -44,8 +44,8 @@ namespace ecolab
       array_ns::is_expression<E>,
       array_ns::array
       <
-        double,
-        typename array_ns::MakeAllocator<double,typename E::Allocator>::type
+        F,
+        typename array_ns::MakeAllocator<F,typename E::Allocator>::type
         >
       >::T
     operator*(const E& x) const

models/ecolab_model.cc

@@ -543,12 +543,13 @@ void SpatialModel::migrate()
 void SpatialModel::makeConsistent()
 {
   // all cells must have same number of species. Pack out with zero density if necessary
-  unsigned long nsp=0;
+  unsigned long nsp=species.size();
   for (auto& i: *this) nsp=max(nsp,i->as<EcolabCell>()->density.size());
 #ifdef MPI_SUPPORT
   MPI_Allreduce(MPI_IN_PLACE,&nsp,1,MPI_UNSIGNED_LONG,MPI_MAX,MPI_COMM_WORLD);
 #endif
-  for (auto& i: *this)
-    i->as<EcolabCell>()->density<<=array<int>(nsp-i->as<EcolabCell>()->density.size(),0);
+  forAll([=,this](EcolabCell& c) {
+    c.density<<=array<int,EcolabCell::Allocator<int>>(nsp-c.density.size(),0,c.density.allocator());
+  });
   ModelData::makeConsistent(nsp);
 }

models/spatial_ecolab.py

@@ -1,6 +1,6 @@
-from ecolab_model import spatial_ecolab
+from ecolab_model import spatial_ecolab as ecolab
 # spatial_ecolab is a smart ptr to possible Device accessible memery, so must be dereferenced
-ecolab=spatial_ecolab()
+#ecolab=spatial_ecolab()
 
 from random import random, seed as randomSeed
 
@@ -52,10 +52,10 @@ def randomList(num, min, max):
 
 def step():
     ecolab.generate(100)
-    ecolab.mutate()
-    ecolab.migrate()
-    ecolab.condense()
-    ecolab.gather()
+#    ecolab.mutate()
+#    ecolab.migrate()
+#    ecolab.condense()
+#    ecolab.gather()
     if myid()==0:
         nsp=len(ecolab.species)
         statusBar.configure(text=f't={ecolab.tstep()} nsp:{nsp}')