From 167fa2c575c8a9316e6b3a389346110b2bbad722 Mon Sep 17 00:00:00 2001
From: Russell Standish <hpcoder@hpcoders.com.au>
Date: Mon, 27 Jul 2026 16:48:05 +1000
Subject: [PATCH 01/19] Temporary checin to ask copilot a question

---
 include/DeviceAllocator.h | 20 +++++++++++++++-----
 include/arrays.h          | 32 +++++++++++++++++++++++++-------
 include/non-sycl.h        |  6 +++---
 include/sycl.h            |  4 ----
 models/ecolab_model.cc    | 20 +++++++++++++-------
 models/ecolab_model.h     |  2 +-
 models/spatial_ecolab.py  | 38 ++++++++++++++++++++------------------
 7 files changed, 77 insertions(+), 45 deletions(-)

diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h
index fb3a2c2..0761f74 100644
--- a/include/DeviceAllocator.h
+++ b/include/DeviceAllocator.h
@@ -74,6 +74,7 @@ namespace ecolab
             slot.value=x;
             Atomic publish(slot.seq);
             publish.store(pos+1,sycl::memory_order::release);
+            printf("dealloc: %zu\n",size-head+tail);
             return;
           }
       }
@@ -95,6 +96,7 @@ namespace ecolab
             unsigned v=slot.value;
             Atomic release(slot.seq);
             release.store(pos+size,sycl::memory_order::release);
+            printf("alloc: %zu\n",size-head+tail);
             return v;
           }
         if (diff<0)
@@ -132,9 +134,14 @@ namespace ecolab
     DeviceAllocator<order+2> nextAllocator; // next size up allocator
   public:
     void init() {
-      for (int pagesLeftToInit=numPages; pagesLeftToInit>0; pagesLeftToInit-=workGroupSize) 
-        syclQ().parallel_for(std::min(workGroupSize,unsigned(pagesLeftToInit)),
-                             [this](size_t) {queue.init();});
+      auto chunkOWork=syclQ().get_device().
+        get_info<sycl::info::device::max_compute_units>()*workGroupSize;
+      for (int pagesLeftToInit=numPages; pagesLeftToInit>0; )
+        {
+          syclQ().parallel_for(std::min(chunkOWork,unsigned(pagesLeftToInit)),
+                               [this](size_t) {queue.init();});
+          pagesLeftToInit-=chunkOWork;
+        }
       nextAllocator.init();
     }
     // all members of group get the same pointer
@@ -155,8 +162,9 @@ namespace ecolab
       if (!p) return;
       if (p>=memory && p<memory+poolSize) {
         groupBarrier();
-        if (groupLeader())
+        if (groupLeader()) {
           queue.enqueue((reinterpret_cast<char*>(p)-memory)>>order);
+        }
         return;
       }
       nextAllocator.deallocate(p,size);
@@ -253,11 +261,13 @@ namespace ecolab
       char* alloc=b.buffer+offs;
       return reinterpret_cast<T*>(alloc);
     }
-    void deallocate(T*,size_t) {} // cleaned up when group exits
+    void deallocate(T*p,size_t) {if (groupLeader()) printf("local dealloc %p\n",p);} // cleaned up when group exits
     template<class U> struct rebind {using other=LocalAllocator<U>;};
     // allocator is stateless
     bool operator==(const LocalAllocator&) const {return true;}
   };
+#else
+  template <class T> using LocalAllocator=std::allocator<T>;
 #endif
    
 }
diff --git a/include/arrays.h b/include/arrays.h
index 4ef1813..e925a9d 100644
--- a/include/arrays.h
+++ b/include/arrays.h
@@ -1523,10 +1523,15 @@ namespace ecolab
 
       void set_size(size_t s) {dt = alloc(s);}
 
-      unsigned& ref() // access reference counter
+#ifdef __SYCL_DEVICE_ONLY__
+      using Ref=sycl::atomic_ref<unsigned,sycl::memory_order::relaxed,sycl::memory_scope::work_group>;
+#else
+      using Ref=unsigned&;
+#endif
+      Ref ref() // access reference counter
       {
         assert(dt);
-        return dt->cnt;
+        return Ref(dt->cnt);
       }
 
       // increment reference counter
@@ -1649,15 +1654,14 @@ namespace ecolab
       
       /// resize array to \a s elements. Id \a copy is true, then ensure data is retained
       void resize(size_t s, bool copy) {
+        if (s==size()) return;
+        groupBarrier();
         if (!dt || s>dt->sz || ref()>1)
           {
-//            array tmp(*this);
-//            release();
-//            dt = alloc(s);
-//            if (dt && copy) asg_v(dt->dt,std::min(s,tmp.size()),tmp.data());
             array tmp(s,m_allocator);
-            if (copy) asg_v(tmp.dt->dt,std::min(s,tmp.size()),dt->dt);
+            if (dt && copy) asg_v(tmp.dt->dt,std::min(s,tmp.size()),dt->dt);
             swap(tmp);
+            if (groupLeader()) printf("resize: deallocating %p, refCnt=%u\n",tmp.dt,tmp.refCnt());
           } 
         if (dt) dt->sz=s; // in case s is smaller
       } 
@@ -1672,10 +1676,12 @@ namespace ecolab
 
       void swap(array& x) {
 #ifndef __SYCL_DEVICE_ONLY__
+        printf("host swapping %p & %p\n",dt, x.dt);
         std::swap(dt, x.dt);
         std::swap(m_allocator,x.m_allocator);
 #else
         if (onDevice && x.onDevice) {
+          if (groupLeader()) printf("on device swapping %p & %p\n",dt, x.dt);
           std::swap(dt, x.dt);
           std::swap(m_allocator,x.m_allocator);
         } else {
@@ -1684,10 +1690,12 @@ namespace ecolab
           auto lhs=dt, rhs=x.dt;
           auto lalloc=m_allocator, ralloc=x.m_allocator;
           groupBarrier();
+          if (groupLeader()) printf("swapping %p & %p\n",lhs,rhs);
           dt=rhs;
           x.dt=lhs;
           m_allocator=ralloc;
           x.m_allocator=lalloc;
+          groupBarrier();
         }
 #endif
       }
@@ -1785,7 +1793,17 @@ namespace ecolab
       typename enable_if<is_expression<E>,array&>::T
       operator<<=(const E& x) {
         auto origSize=size();
+#ifdef __SYCL_DEVICE_ONLY__
+        if (groupLeader()) printf("<<= resize from %zu to %zu\n",origSize,origSize+x.size());
+#else
+        printf("<<= host resize from %zu to %zu\n",origSize,origSize+x.size());
+#endif
         resize(origSize+x.size());
+#ifdef __SYCL_DEVICE_ONLY__
+        if (groupLeader()) printf("<<= after resize from %zu to %zu\n",origSize,origSize+x.size());
+#else
+        printf("<<= host after resize from %zu to %zu\n",origSize,origSize+x.size());
+#endif
         asg_v(data()+origSize,x.size(),x);
         return *this;
       }
diff --git a/include/non-sycl.h b/include/non-sycl.h
index 3a8a74c..de48dff 100644
--- a/include/non-sycl.h
+++ b/include/non-sycl.h
@@ -23,9 +23,9 @@ namespace ecolab
     operator T&() {return data;}
   };
   
-  template <class T> struct GroupLocal: public std::unique_ptr<T>
-  {
-  };
+//  template <class T> struct GroupLocal: public std::unique_ptr<T>
+//  {
+//  };
   
   template <class E, ecolab::USMAlloc UA>
   struct SyclRandomEngine: public E {};
diff --git a/include/sycl.h b/include/sycl.h
index a1569dc..0f282f6 100644
--- a/include/sycl.h
+++ b/include/sycl.h
@@ -129,10 +129,6 @@ namespace ecolab
     operator bool() const {return true;} // always defined
   };
 
-#ifndef __SYCL_DEVICE_ONLY__
-  template <class T> using LocalAllocator=std::allocator<T>;
-#endif
-  
   inline void groupBarrier() {
 #ifdef __SYCL_DEVICE_ONLY__
     sycl::group_barrier(syclGroup());
diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc
index 8a87bbf..0c82285 100644
--- a/models/ecolab_model.cc
+++ b/models/ecolab_model.cc
@@ -223,7 +223,8 @@ void SpatialModel::mutate()
 
   groupedForAll([newSp=newSp.data(),mut_scale=&*mut_scale,this](EcolabCell& c,size_t i) {
     assert(all(c.density>=0));
-    newSp[i]=c.mutate(*mut_scale);
+    EcolabPoint::UnsignedArray tmp(c.mutate(*mut_scale),c.density.allocator());
+    newSp[i]=tmp;
   });
 
   array<unsigned> new_sp;
@@ -240,7 +241,7 @@ void SpatialModel::mutate()
     }
 
   // deallocate on device
-  groupedForAll([newSp=newSp.data()](EcolabCell& c,size_t i) {
+  hostForAll([newSp=newSp.data()](EcolabCell& c,size_t i) {
     newSp[i].clear();
     assert(newSp[i].refCnt()==0);
   });
@@ -264,19 +265,24 @@ void SpatialModel::mutate()
 #else
   ModelData::mutate(new_sp);
 #endif
-  if (new_sp.size()==0) return;
+  //if (new_sp.size()==0) return;
 
-  computeODiagIdx();
-  
+  //  computeODiagIdx();
+  mut_scale->clear();
+  newSp.clear();
+
+  (*cell_ids)<<=0;
+  cout<<"b4 append:"<<endl;
   // set the new species density to 1 for those created on this cell
   groupedForAll([cell_ids=&*cell_ids](EcolabCell& c,size_t) {
     //hostForAll([cell_ids=&*cell_ids,this](EcolabCell& c,size_t) {
     c.density <<= (*cell_ids)==c.id;
   });
+  cout<<"after append:"<<endl;
 }
 
 template <class E>
-EcolabPoint::LocalArray EcolabPoint::mutate(const E& mut_scale)
+EcolabPoint::UnsignedArray EcolabPoint::mutate(const E& mut_scale)
 {
   /* calculate the number of mutants each species produces */
   if (density.size()==0) return {};
@@ -296,7 +302,7 @@ EcolabPoint::LocalArray EcolabPoint::mutate(const E& mut_scale)
   
   density-=speciations;
 
-  LocalArray new_sp(offsets[nsp]);
+  UnsignedArray new_sp(offsets[nsp], density.allocator());
   array_ns::map(nsp, [offsets=offsets.data(),new_sp=new_sp.data()](size_t i) {
     for (auto j=offsets[i]; j<offsets[i+1]; ++j)
       new_sp[j]=i;
diff --git a/models/ecolab_model.h b/models/ecolab_model.h
index 288fbbe..8a95e63 100644
--- a/models/ecolab_model.h
+++ b/models/ecolab_model.h
@@ -89,7 +89,7 @@ class EcolabPoint
   
   void generate(unsigned niter, const ModelData&);
   void condense(const ModelData::BoolArray& mask, size_t mask_true);
-  template <class E> LocalArray mutate(const E&);
+  template <class E> UnsignedArray mutate(const E&);
   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);
diff --git a/models/spatial_ecolab.py b/models/spatial_ecolab.py
index c42a5e3..82b9638 100644
--- a/models/spatial_ecolab.py
+++ b/models/spatial_ecolab.py
@@ -21,7 +21,7 @@
 ecolab.repro_max(0.1)
 ecolab.odiag_min(-1e-5)
 ecolab.odiag_max(1e-5)
-ecolab.mut_max(1e-5)
+ecolab.mut_max(1000)
 #ecolab.mut_max(1e-3)
 ecolab.sp_sep(0.1)
 
@@ -32,8 +32,8 @@ def randomList(num, min, max):
 
 #numX=12
 #numY=12
-numX=2
-numY=2
+numX=1
+numY=1
 ecolab.setGrid(numX,numY)
 ecolab.partitionObjects()
 
@@ -64,7 +64,7 @@ def randomList(num, min, max):
 extinctions=0
 migrations=0
 def stepImpl():
-    ecolab.generate(100)
+    #ecolab.generate(100)
     ecolab.mutate()
 
     epochTs=ecolab.tstep()%epoch
@@ -75,34 +75,36 @@ def stepImpl():
 
     global extinctions, migrations
     #migrations+=ecolab.migrate()
-    extinctions+=ecolab.condense()
+    #extinctions+=ecolab.condense()
     #print(ecolab.nsp()())
     #ecolab.syncThreads()
     #ecolab.gather()
 
-print(ecolab.nsp()())
+#print(ecolab.nsp()())
 ecolab.makeConsistent()
 ecolab.syncThreads()
-print(ecolab.nsp()())
+#print(ecolab.nsp()())
 
 from timeit import timeit
-print(timeit('stepImpl()', globals=globals(), number=100))
-                
+print('-----')
+print(timeit('ecolab.mutate()', globals=globals(), number=1))
+print('-----')
+               
 def step():
     global extinctions,migrations
     extinctions=0
     migrations=0
     for i in range(epoch//1000000):
         stepImpl()
-    #print('migrations=',migrations,' extinctions=',extinctions)
-    if myid()==0:
-        nsp=len(ecolab.species)
-        statusBar.configure(text=f't={ecolab.tstep()} nsp:{nsp}')
-        plot('No. species',ecolab.tstep(),nsp,200*(ecolab.tstep()%epoch<0.5*epoch))
-        #plot('No. species',ecolab.tstep(),nsp)
-        plot('No. species by cell',ecolab.tstep(),ecolab.nsp()())
-        plot('Extinctions',ecolab.tstep(),extinctions)
-        plot('Migration',ecolab.tstep(),migrations)
+#    #print('migrations=',migrations,' extinctions=',extinctions)
+#    if myid()==0:
+#        nsp=len(ecolab.species)
+#        statusBar.configure(text=f't={ecolab.tstep()} nsp:{nsp}')
+#        plot('No. species',ecolab.tstep(),nsp,200*(ecolab.tstep()%epoch<0.5*epoch))
+#        #plot('No. species',ecolab.tstep(),nsp)
+#        plot('No. species by cell',ecolab.tstep(),ecolab.nsp()())
+#        plot('Extinctions',ecolab.tstep(),extinctions)
+#        plot('Migration',ecolab.tstep(),migrations)
 #        for i in range(numX):
 #            for j in range(numY):
 #                plot(f'Density({i},{j})',ecolab.tstep(),ecolab.cell(i,j).density(), pens=ecolab.species())

From f704151605466d52ab33d1d1934f096c8add43f1 Mon Sep 17 00:00:00 2001
From: Russell Standish <hpcoder@hpcoders.com.au>
Date: Tue, 28 Jul 2026 09:20:57 +1000
Subject: [PATCH 02/19] Ask Copilot another question.

---
 include/DeviceAllocator.h | 19 +++++++++++++++----
 include/arrays.h          |  3 ++-
 models/ecolab_model.cc    | 17 ++++++++++-------
 3 files changed, 27 insertions(+), 12 deletions(-)

diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h
index 0761f74..1e243a7 100644
--- a/include/DeviceAllocator.h
+++ b/include/DeviceAllocator.h
@@ -74,7 +74,7 @@ namespace ecolab
             slot.value=x;
             Atomic publish(slot.seq);
             publish.store(pos+1,sycl::memory_order::release);
-            printf("dealloc: %zu\n",size-head+tail);
+            printf("dealloc: %zu, pageSize=%zu\n",size-pos+Atomic(tail),poolSize/size);
             return;
           }
       }
@@ -96,7 +96,7 @@ namespace ecolab
             unsigned v=slot.value;
             Atomic release(slot.seq);
             release.store(pos+size,sycl::memory_order::release);
-            printf("alloc: %zu\n",size-head+tail);
+            printf("alloc: %zu, pageSize=%zu\n",size-Atomic(head)+pos,poolSize/size);
             return v;
           }
         if (diff<0)
@@ -149,6 +149,7 @@ namespace ecolab
       if (size==0) return nullptr;
       if (size<=pageSize) {
         unsigned offs;
+        groupBarrier();
         if (groupLeader()) offs=queue.dequeue();
 #ifdef __SYCL_DEVICE_ONLY__
         offs=sycl::group_broadcast(syclGroup(),offs);
@@ -174,10 +175,15 @@ namespace ecolab
   };
   
   inline DeviceAllocator<>& deviceAllocator() {
+//#ifdef __SYCL_DEVICE_ONLY__
+//    printf("deviceAllocator() illegally called on device\n",0);
+//    return *reinterpret_cast<DeviceAllocator<>*>(0);
+//#else
     static DeviceType<DeviceAllocator<>> deviceAllocator;
     static int dummy=
       (deviceAllocator->init(), syclQ().wait_and_throw(), 0);      
     return *deviceAllocator;
+    //#endif
   }
 
   /// Allocator wrapping the DeviceAllocator singleton
@@ -190,9 +196,14 @@ namespace ecolab
     using difference_type=std::ptrdiff_t;
     using propagate_on_container_move_assignment=std::true_type;
 
-    DeviceAllocator<>* allocator=&deviceAllocator();
+    DeviceAllocator<>* allocator;
     
-    GlobalDeviceAllocator() = default; // note: default constructor must be called on host
+#ifdef __SYCL_DEVICE_ONLY__
+    GlobalDeviceAllocator() = delete;
+#else
+    GlobalDeviceAllocator() // note: default constructor must be called on host
+    {allocator=&deviceAllocator();}
+#endif
     template <class U>
     GlobalDeviceAllocator(const GlobalDeviceAllocator<U>& other):
       allocator(other.allocator) {}
diff --git a/include/arrays.h b/include/arrays.h
index e925a9d..f670681 100644
--- a/include/arrays.h
+++ b/include/arrays.h
@@ -1610,7 +1610,8 @@ namespace ecolab
       typedef size_t size_type; 
       using Allocator=A;
 
-      array(const Allocator& alloc={}): m_allocator(alloc) {}
+      array()=default;
+      array(const Allocator& alloc): m_allocator(alloc) {}
       explicit array(size_t s, const Allocator& alloc=Allocator()): m_allocator(alloc)
       {
         set_size(s);
diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc
index 0c82285..1b0fd52 100644
--- a/models/ecolab_model.cc
+++ b/models/ecolab_model.cc
@@ -219,16 +219,17 @@ void SpatialModel::mutate()
   assert(all(*mut_scale<=1));
   last_mut_tstep=tstep;
 
-  vector<EcolabPoint::UnsignedArray,ModelData::Allocator<EcolabPoint::UnsignedArray>> newSp(size());
-
+  auto deviceAllocator=cell(0,0).density.allocator();
+  vector<EcolabPoint::UnsignedArray,ModelData::Allocator<EcolabPoint::UnsignedArray>>
+    newSp(size(),{deviceAllocator});
+  
   groupedForAll([newSp=newSp.data(),mut_scale=&*mut_scale,this](EcolabCell& c,size_t i) {
     assert(all(c.density>=0));
-    EcolabPoint::UnsignedArray tmp(c.mutate(*mut_scale),c.density.allocator());
-    newSp[i]=tmp;
+    newSp[i]=c.mutate(*mut_scale);
   });
 
   array<unsigned> new_sp;
-  DeviceType<EcolabPoint::UnsignedArray> cell_ids;
+  DeviceType<EcolabPoint::UnsignedArray> cell_ids(deviceAllocator);
   syncThreads();
   
   // TODO - this is a kind of scan - can it be done on device?
@@ -285,7 +286,7 @@ template <class E>
 EcolabPoint::UnsignedArray EcolabPoint::mutate(const E& mut_scale)
 {
   /* calculate the number of mutants each species produces */
-  if (density.size()==0) return {};
+  if (density.size()==0) return {density.allocator()};
 #ifdef __SYCL_DEVICE_ONLY__
   LocalArray speciations=roundArray(mut_scale * density);
   auto nsp=density.size();
@@ -298,7 +299,7 @@ EcolabPoint::UnsignedArray EcolabPoint::mutate(const E& mut_scale)
     offsets[nsp]=offsets[nsp-1]+speciations[nsp-1];
   groupBarrier();
 
-  if (offsets[nsp]==0) return {};
+  if (offsets[nsp]==0) return {density.allocator()};
   
   density-=speciations;
 
@@ -583,6 +584,7 @@ bool ConnectionPlot::redraw(int x0, int y0, int width, int height)
   return true;
 }
 
+#ifndef __SYCL_DEVICE_ONLY__
 void SpatialModel::setGrid(size_t nx, size_t ny)
 {
   numX=nx; numY=ny;
@@ -609,6 +611,7 @@ void SpatialModel::setGrid(size_t nx, size_t ny)
   for (auto& i: objects)
     maxNbrs=std::max(maxNbrs, i->neighbours.size());
 }
+#endif
 
 void SpatialModel::generate(unsigned niter)
 {

From 8616caa3f79be42af95a8e46510ef0e0f066275f Mon Sep 17 00:00:00 2001
From: Russell Standish <hpcoder@hpcoders.com.au>
Date: Tue, 28 Jul 2026 12:24:52 +1000
Subject: [PATCH 03/19] Commit latest for copilot question.

---
 include/DeviceAllocator.h |  6 +++---
 include/arrays.h          | 20 ++++++++++++++------
 include/sycl.h            |  1 +
 3 files changed, 18 insertions(+), 9 deletions(-)

diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h
index 1e243a7..b27c484 100644
--- a/include/DeviceAllocator.h
+++ b/include/DeviceAllocator.h
@@ -74,7 +74,7 @@ namespace ecolab
             slot.value=x;
             Atomic publish(slot.seq);
             publish.store(pos+1,sycl::memory_order::release);
-            printf("dealloc: %zu, pageSize=%zu\n",size-pos+Atomic(tail),poolSize/size);
+            printf("dealloc: %zu, pageSize=%zu\n",size-headAtomic+Atomic(tail),poolSize/size);
             return;
           }
       }
@@ -96,7 +96,7 @@ namespace ecolab
             unsigned v=slot.value;
             Atomic release(slot.seq);
             release.store(pos+size,sycl::memory_order::release);
-            printf("alloc: %zu, pageSize=%zu\n",size-Atomic(head)+pos,poolSize/size);
+            printf("alloc: %zu, pageSize=%zu\n",size-Atomic(head)+tailAtomic,poolSize/size);
             return v;
           }
         if (diff<0)
@@ -148,7 +148,7 @@ namespace ecolab
     void* allocate(size_t size) {
       if (size==0) return nullptr;
       if (size<=pageSize) {
-        unsigned offs;
+        unsigned offs=~0U;
         groupBarrier();
         if (groupLeader()) offs=queue.dequeue();
 #ifdef __SYCL_DEVICE_ONLY__
diff --git a/include/arrays.h b/include/arrays.h
index f670681..1c33e7b 100644
--- a/include/arrays.h
+++ b/include/arrays.h
@@ -1524,7 +1524,7 @@ namespace ecolab
       void set_size(size_t s) {dt = alloc(s);}
 
 #ifdef __SYCL_DEVICE_ONLY__
-      using Ref=sycl::atomic_ref<unsigned,sycl::memory_order::relaxed,sycl::memory_scope::work_group>;
+      using Ref=sycl::atomic_ref<unsigned,sycl::memory_order::acq_rel,sycl::memory_scope::work_group>;
 #else
       using Ref=unsigned&;
 #endif
@@ -1558,11 +1558,17 @@ namespace ecolab
       
       void release()
       {
+        groupBarrier();
         if (dt)
           {
-            if (ref()==1)
+            bool freeMem=ref()==1;
+#ifdef __SYCL_DEVICE_ONLY__
+            freeMem=sycl::group_broadcast(syclGroup(),freeMem);
+#endif
+            if (freeMem)
               {
                 free(dt);
+                dt=nullptr;
                 return;
               }
             decrRef();
@@ -1594,14 +1600,14 @@ namespace ecolab
           {
             array_data<T>* oldData=dt;
             decrRef();
-            bool freeMem=ref()==0;
+            //bool freeMem=ref()==0;
             dt=alloc(size());
 #ifdef __SYCL_DEVICE_ONLY__
             asg_v(dt->dt,size(),oldData->dt);
 #else
             memcpy(dt->dt,oldData->dt,size()*sizeof(T));
 #endif
-            if (freeMem) free(oldData);
+            //if (freeMem) free(oldData);
           }
       }
 
@@ -1664,7 +1670,9 @@ namespace ecolab
             swap(tmp);
             if (groupLeader()) printf("resize: deallocating %p, refCnt=%u\n",tmp.dt,tmp.refCnt());
           } 
-        if (dt) dt->sz=s; // in case s is smaller
+        groupBarrier();
+        if (groupLeader() && dt) dt->sz=s; // in case s is smaller
+        groupBarrier();
       } 
 
       // note using default argument for copy above breaks classdesc::has_resize.
@@ -1718,8 +1726,8 @@ namespace ecolab
           release();
           if (groupLeader()||onDevice) {
             dt=x.dt;
-            incrRef();
           }
+          incrRef();
         } else
           asgV(x.size(), x);
         return *this;
diff --git a/include/sycl.h b/include/sycl.h
index 0f282f6..9b6c42f 100644
--- a/include/sycl.h
+++ b/include/sycl.h
@@ -138,6 +138,7 @@ namespace ecolab
   inline bool groupLeader() {
 #ifdef __SYCL_DEVICE_ONLY__
     return syclGroup().leader();
+    //return syclGroup().get_local_linear_id()==0;
 #endif
     return true;
   }

From 3994067b85a23993945b17a3e087acb0e9db7f76 Mon Sep 17 00:00:00 2001
From: Russell Standish <hpcoder@hpcoders.com.au>
Date: Tue, 28 Jul 2026 12:51:32 +1000
Subject: [PATCH 04/19] Sync to repo for copilot question

---
 include/arrays.h | 23 +++++++++++++----------
 1 file changed, 13 insertions(+), 10 deletions(-)

diff --git a/include/arrays.h b/include/arrays.h
index 1c33e7b..52f49c1 100644
--- a/include/arrays.h
+++ b/include/arrays.h
@@ -1570,8 +1570,9 @@ namespace ecolab
                 free(dt);
                 dt=nullptr;
                 return;
-              }
-            decrRef();
+              } else
+              decrRef();
+            groupBarrier();
           }
       }
 
@@ -1689,23 +1690,25 @@ namespace ecolab
         std::swap(dt, x.dt);
         std::swap(m_allocator,x.m_allocator);
 #else
-        if (onDevice && x.onDevice) {
-          if (groupLeader()) printf("on device swapping %p & %p\n",dt, x.dt);
-          std::swap(dt, x.dt);
-          std::swap(m_allocator,x.m_allocator);
-        } else {
+//        if (onDevice && x.onDevice) {
+//          if (groupLeader()) printf("on device swapping %p & %p\n",dt, x.dt);
+//          std::swap(dt, x.dt);
+//          std::swap(m_allocator,x.m_allocator);
+//        } else {
           // assumption here is these array may be per thread, or
           // maybe shared by all threads in a group, hence std::swap as above won't work
           auto lhs=dt, rhs=x.dt;
           auto lalloc=m_allocator, ralloc=x.m_allocator;
-          groupBarrier();
+          //groupBarrier();
+          lhs=sycl::group_broadcast(syclGroup(),lhs);
+          rhs=sycl::group_broadcast(syclGroup(),rhs);
           if (groupLeader()) printf("swapping %p & %p\n",lhs,rhs);
           dt=rhs;
           x.dt=lhs;
           m_allocator=ralloc;
           x.m_allocator=lalloc;
           groupBarrier();
-        }
+          //        }
 #endif
       }
     
@@ -1724,7 +1727,7 @@ namespace ecolab
         if (x.dt==dt) return *this;
         if (m_allocator==x.m_allocator) {
           release();
-          if (groupLeader()||onDevice) {
+          /*if (groupLeader()||onDevice)*/ {
             dt=x.dt;
           }
           incrRef();

From 880ea47d0aba0d8c237c30544eac8b3bf0cb0510 Mon Sep 17 00:00:00 2001
From: Russell Standish <hpcoder@hpcoders.com.au>
Date: Tue, 28 Jul 2026 12:59:43 +1000
Subject: [PATCH 05/19] Again for copilot conversation.

---
 include/arrays.h | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/include/arrays.h b/include/arrays.h
index 52f49c1..d76a0ce 100644
--- a/include/arrays.h
+++ b/include/arrays.h
@@ -1569,7 +1569,6 @@ namespace ecolab
               {
                 free(dt);
                 dt=nullptr;
-                return;
               } else
               decrRef();
             groupBarrier();
@@ -1602,11 +1601,12 @@ namespace ecolab
             array_data<T>* oldData=dt;
             decrRef();
             //bool freeMem=ref()==0;
-            dt=alloc(size());
+            auto sz=size();
+            dt=alloc(sz);
 #ifdef __SYCL_DEVICE_ONLY__
-            asg_v(dt->dt,size(),oldData->dt);
+            asg_v(dt->dt,sz,oldData->dt);
 #else
-            memcpy(dt->dt,oldData->dt,size()*sizeof(T));
+            memcpy(dt->dt,oldData->dt,sz*sizeof(T));
 #endif
             //if (freeMem) free(oldData);
           }

From cdf4ef9159d8075bd2df3442e22a405bdf8188ae Mon Sep 17 00:00:00 2001
From: Russell Standish <hpcoder@hpcoders.com.au>
Date: Tue, 28 Jul 2026 13:13:52 +1000
Subject: [PATCH 06/19] More copilot questions: After disabling cow protocol.

---
 include/arrays.h | 30 ++++++++++++++++--------------
 1 file changed, 16 insertions(+), 14 deletions(-)

diff --git a/include/arrays.h b/include/arrays.h
index d76a0ce..8f568a4 100644
--- a/include/arrays.h
+++ b/include/arrays.h
@@ -1561,7 +1561,7 @@ namespace ecolab
         groupBarrier();
         if (dt)
           {
-            bool freeMem=ref()==1;
+            bool freeMem=true; //ref()==1;
 #ifdef __SYCL_DEVICE_ONLY__
             freeMem=sycl::group_broadcast(syclGroup(),freeMem);
 #endif
@@ -1596,7 +1596,7 @@ namespace ecolab
       
       void copy() //any nonconst method needs to call this
       {           // to implement copy-on-write semantics
-        if (dt && ref()>1)
+        if (dt /*&& ref()>1*/)
           {
             array_data<T>* oldData=dt;
             decrRef();
@@ -1608,7 +1608,8 @@ namespace ecolab
 #else
             memcpy(dt->dt,oldData->dt,sz*sizeof(T));
 #endif
-            //if (freeMem) free(oldData);
+            //if (freeMem)
+            free(oldData);
           }
       }
 
@@ -1632,9 +1633,10 @@ namespace ecolab
 
       array(const array& x): m_allocator(x.m_allocator) 
       {
-        dt=x.dt;
-        incrRef();
-          
+        //dt=x.dt;
+        //incrRef();
+        set_size(x.size());
+        asg_v(data(),x.size(),x);
       }
 
       template <class expr>
@@ -1664,7 +1666,7 @@ namespace ecolab
       void resize(size_t s, bool copy) {
         if (s==size()) return;
         groupBarrier();
-        if (!dt || s>dt->sz || ref()>1)
+        if (!dt || s>dt->sz /*|| ref()>1*/)
           {
             array tmp(s,m_allocator);
             if (dt && copy) asg_v(tmp.dt->dt,std::min(s,tmp.size()),dt->dt);
@@ -1725,13 +1727,13 @@ namespace ecolab
 
       array& operator=(const array& x) {
         if (x.dt==dt) return *this;
-        if (m_allocator==x.m_allocator) {
-          release();
-          /*if (groupLeader()||onDevice)*/ {
-            dt=x.dt;
-          }
-          incrRef();
-        } else
+//        if (m_allocator==x.m_allocator) {
+//          release();
+//          /*if (groupLeader()||onDevice)*/ {
+//            dt=x.dt;
+//          }
+//          incrRef();
+//        } else
           asgV(x.size(), x);
         return *this;
       }

From 2e96857b771110d26defa431437eff5f91ef1908 Mon Sep 17 00:00:00 2001
From: Russell Standish <hpcoder@hpcoders.com.au>
Date: Tue, 28 Jul 2026 17:39:10 +1000
Subject: [PATCH 07/19] For copilot questioning

---
 include/DeviceAllocator.h | 92 +++++++++++++++++++--------------------
 include/arrays.h          | 20 ++++-----
 include/sycl.h            |  8 ++++
 models/ecolab_model.cc    | 16 ++++---
 models/ecolab_model.h     |  2 +-
 5 files changed, 73 insertions(+), 65 deletions(-)

diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h
index b27c484..f2ce4f6 100644
--- a/include/DeviceAllocator.h
+++ b/include/DeviceAllocator.h
@@ -136,12 +136,8 @@ namespace ecolab
     void init() {
       auto chunkOWork=syclQ().get_device().
         get_info<sycl::info::device::max_compute_units>()*workGroupSize;
-      for (int pagesLeftToInit=numPages; pagesLeftToInit>0; )
-        {
-          syclQ().parallel_for(std::min(chunkOWork,unsigned(pagesLeftToInit)),
-                               [this](size_t) {queue.init();});
-          pagesLeftToInit-=chunkOWork;
-        }
+      syclQ().parallel_for(std::min(chunkOWork,unsigned(numPages)),
+                           [this](size_t) {queue.init();});
       nextAllocator.init();
     }
     // all members of group get the same pointer
@@ -150,9 +146,9 @@ namespace ecolab
       if (size<=pageSize) {
         unsigned offs=~0U;
         groupBarrier();
-        if (groupLeader()) offs=queue.dequeue();
+        if (localThreadId()==0) offs=queue.dequeue();
 #ifdef __SYCL_DEVICE_ONLY__
-        offs=sycl::group_broadcast(syclGroup(),offs);
+        offs=sycl::group_broadcast(syclGroup(),offs,0);
 #endif
         if (offs!=~0U)
           return memory+(offs<<order);
@@ -240,46 +236,46 @@ namespace ecolab
    {return *sycl::ext::oneapi::group_local_memory_for_overwrite<LocalAllocatorBuffer>(syclGroup());}
 
   
-#ifdef __SYCL_DEVICE_ONLY__
-  
-  /**
-     A Local allocator allocates memory from device local memory,
-     which is shared between threads of a work group, and has the same
-     lifetime as the kernel
-  */
-  template <class T>
-  class LocalAllocator
-  {
-  public:
-    using value_type=T;
-    using pointer=T*;
-    using reference=T&;
-    using difference_type=std::ptrdiff_t;
-    using propagate_on_container_move_assignment=std::true_type;
-
-    // no need for destructor, as Impl has nothing to tear down
-    T* allocate(size_t n) {
-      auto& b=localAllocatorBuffer();
-      unsigned offs=b.next;
-      if (offs+n*sizeof(T)>LocalAllocatorSize)
-        {
-          fatalErrorFlag()=true;
-          return nullptr;
-        }
-      sycl::group_barrier(syclGroup());
-      if (syclGroup().leader()) b.next+=n*sizeof(T);
-      sycl::group_barrier(syclGroup());
-      char* alloc=b.buffer+offs;
-      return reinterpret_cast<T*>(alloc);
-    }
-    void deallocate(T*p,size_t) {if (groupLeader()) printf("local dealloc %p\n",p);} // cleaned up when group exits
-    template<class U> struct rebind {using other=LocalAllocator<U>;};
-    // allocator is stateless
-    bool operator==(const LocalAllocator&) const {return true;}
-  };
-#else
-  template <class T> using LocalAllocator=std::allocator<T>;
-#endif
+//#ifdef __SYCL_DEVICE_ONLY__
+//  
+//  /**
+//     A Local allocator allocates memory from device local memory,
+//     which is shared between threads of a work group, and has the same
+//     lifetime as the kernel
+//  */
+//  template <class T>
+//  class LocalAllocator
+//  {
+//  public:
+//    using value_type=T;
+//    using pointer=T*;
+//    using reference=T&;
+//    using difference_type=std::ptrdiff_t;
+//    using propagate_on_container_move_assignment=std::true_type;
+//
+//    // no need for destructor, as Impl has nothing to tear down
+//    T* allocate(size_t n) {
+//      auto& b=localAllocatorBuffer();
+//      unsigned offs=b.next;
+//      if (offs+n*sizeof(T)>LocalAllocatorSize)
+//        {
+//          fatalErrorFlag()=true;
+//          return nullptr;
+//        }
+//      sycl::group_barrier(syclGroup());
+//      if (syclGroup().leader()) b.next+=n*sizeof(T);
+//      sycl::group_barrier(syclGroup());
+//      char* alloc=b.buffer+offs;
+//      return reinterpret_cast<T*>(alloc);
+//    }
+//    void deallocate(T*p,size_t) {if (groupLeader()) printf("local dealloc %p\n",p);} // cleaned up when group exits
+//    template<class U> struct rebind {using other=LocalAllocator<U>;};
+//    // allocator is stateless
+//    bool operator==(const LocalAllocator&) const {return true;}
+//  };
+//#else
+//  template <class T> using LocalAllocator=std::allocator<T>;
+//#endif
    
 }
 #endif
diff --git a/include/arrays.h b/include/arrays.h
index 8f568a4..f179bfe 100644
--- a/include/arrays.h
+++ b/include/arrays.h
@@ -1582,16 +1582,16 @@ namespace ecolab
       void asgV(size_t size, const E& x)
       {
         // copy into temporary data, as E may contain references to this
-#ifdef __SYCL_DEVICE_ONLY__
-        array<T,LocalAllocator<T>> tmp(size);
-        asg_v(tmp.data(),size,x);
-        resize(size, false);
-        asg_v(data(),size,tmp);
-#else
+//#ifdef __SYCL_DEVICE_ONLY__
+//        array tmp(size,m_allocator);
+//        asg_v(tmp.data(),size,x);
+//        resize(size, false);
+//        asg_v(data(),size,tmp);
+//#else
         array tmp(size,m_allocator);
         asg_v(tmp.data(),size,x);
         swap(tmp);
-#endif
+        //#endif
       }
       
       void copy() //any nonconst method needs to call this
@@ -1701,9 +1701,9 @@ namespace ecolab
           // maybe shared by all threads in a group, hence std::swap as above won't work
           auto lhs=dt, rhs=x.dt;
           auto lalloc=m_allocator, ralloc=x.m_allocator;
-          //groupBarrier();
-          lhs=sycl::group_broadcast(syclGroup(),lhs);
-          rhs=sycl::group_broadcast(syclGroup(),rhs);
+          groupBarrier();
+          //lhs=sycl::group_broadcast(syclGroup(),lhs);
+          //rhs=sycl::group_broadcast(syclGroup(),rhs);
           if (groupLeader()) printf("swapping %p & %p\n",lhs,rhs);
           dt=rhs;
           x.dt=lhs;
diff --git a/include/sycl.h b/include/sycl.h
index 9b6c42f..6b7140e 100644
--- a/include/sycl.h
+++ b/include/sycl.h
@@ -143,6 +143,14 @@ namespace ecolab
     return true;
   }
 
+  inline size_t localThreadId() {
+#ifdef __SYCL_DEVICE_ONLY__
+    return syclGroup().get_local_linear_id();
+#else
+    return 0;
+#endif
+  }
+  
   inline bool onDevice() {
 #ifdef __SYCL_DEVICE_ONLY__
     return true;
diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc
index 1b0fd52..4706bb5 100644
--- a/models/ecolab_model.cc
+++ b/models/ecolab_model.cc
@@ -89,7 +89,8 @@ RoundArray<E,EcolabPoint> EcolabPoint::roundArray(const E& expr)
 
 void EcolabPoint::generate(unsigned niter, const ModelData& model)
 {
-  array<int,LocalAllocator<int>> lDensity(density), tmp(density.size());
+  //array<int,LocalAllocator<int>> lDensity(density), tmp(density.size());
+  array<int,Allocator<int>> lDensity(density), tmp(density.size(), density.allocator());
   //auto& lDensity=density;
   //array<int,GlobalDeviceAllocator<int>> tmp(density.size(), density.allocator());
   
@@ -127,7 +128,8 @@ void EcolabPoint::condense(const ModelData::BoolArray& mask, size_t mask_true)
     density.clear();
     return;
   }
-  LocalArray tmp(mask_true);
+  //LocalArray tmp(mask_true);
+  UnsignedArray tmp(mask_true,density.allocator());
   if (groupLeader())
     for (size_t i=0, j=0; i<density.size(); ++i)
       if (mask[i] && j<mask_true) 
@@ -272,10 +274,10 @@ void SpatialModel::mutate()
   mut_scale->clear();
   newSp.clear();
 
-  (*cell_ids)<<=0;
+  //(*cell_ids)<<=0;
   cout<<"b4 append:"<<endl;
   // set the new species density to 1 for those created on this cell
-  groupedForAll([cell_ids=&*cell_ids](EcolabCell& c,size_t) {
+  hostForAll([cell_ids=&*cell_ids](EcolabCell& c,size_t) {
     //hostForAll([cell_ids=&*cell_ids,this](EcolabCell& c,size_t) {
     c.density <<= (*cell_ids)==c.id;
   });
@@ -288,10 +290,12 @@ EcolabPoint::UnsignedArray EcolabPoint::mutate(const E& mut_scale)
   /* calculate the number of mutants each species produces */
   if (density.size()==0) return {density.allocator()};
 #ifdef __SYCL_DEVICE_ONLY__
-  LocalArray speciations=roundArray(mut_scale * density);
+  //LocalArray speciations=roundArray(mut_scale * density);
+  UnsignedArray speciations(roundArray(mut_scale * density), density.allocator());
   auto nsp=density.size();
   //  auto new_sp = gen_index(speciations);
-  LocalArray offsets(nsp+1);
+  //LocalArray offsets(nsp+1);
+  UnsignedArray offsets(nsp+1,density.allocator());
   sycl::joint_exclusive_scan(syclGroup(),speciations.data(),speciations.data()+nsp,
                              offsets.data(),sycl::plus<unsigned>());
   groupBarrier();
diff --git a/models/ecolab_model.h b/models/ecolab_model.h
index 8a95e63..6cd2fc8 100644
--- a/models/ecolab_model.h
+++ b/models/ecolab_model.h
@@ -82,7 +82,7 @@ class EcolabPoint
   template <class T> using Allocator=std::allocator<T>;
 #endif
   using UnsignedArray=array<unsigned,Allocator<unsigned>>;
-  using LocalArray=array<unsigned,LocalAllocator<unsigned>>;
+  //using LocalArray=array<unsigned,LocalAllocator<unsigned>>;
 
   Float salt;  /* random no. used for migration */
   array<int,Allocator<int>> density;

From 5d5a79baf12c130099335df9106d84925b9cd980 Mon Sep 17 00:00:00 2001
From: Russell Standish <hpcoder@hpcoders.com.au>
Date: Wed, 29 Jul 2026 14:02:37 +1000
Subject: [PATCH 08/19] Comitting now to put a pin it it. Big reveal is that
 assorted array methods cannot be used in an envirionment without full threads
 in a group participation.

---
 include/DeviceAllocator.h | 87 ++++++++++++++++++++-------------------
 include/arrays.h          | 59 +++++++++-----------------
 include/ecolab.h          |  1 +
 models/ecolab_model.cc    | 66 +++++++++++++++--------------
 models/ecolab_model.h     |  2 +-
 models/spatial_ecolab.py  | 26 ++++++------
 6 files changed, 117 insertions(+), 124 deletions(-)

diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h
index f2ce4f6..d7b685b 100644
--- a/include/DeviceAllocator.h
+++ b/include/DeviceAllocator.h
@@ -60,6 +60,8 @@ namespace ecolab
     
     void enqueue(unsigned x)
     {
+//      Atomic t(tail);
+//      slots[--t].value=x;
       while (true)
       {
         Atomic headAtomic(head);
@@ -74,7 +76,6 @@ namespace ecolab
             slot.value=x;
             Atomic publish(slot.seq);
             publish.store(pos+1,sycl::memory_order::release);
-            printf("dealloc: %zu, pageSize=%zu\n",size-headAtomic+Atomic(tail),poolSize/size);
             return;
           }
       }
@@ -82,6 +83,9 @@ namespace ecolab
 
     unsigned dequeue()
     {
+//      Atomic t(tail);
+//      if (t>=size) return ~0;
+//      return slots[t++].value;
       while (true)
       {
         Atomic tailAtomic(tail);
@@ -96,7 +100,6 @@ namespace ecolab
             unsigned v=slot.value;
             Atomic release(slot.seq);
             release.store(pos+size,sycl::memory_order::release);
-            printf("alloc: %zu, pageSize=%zu\n",size-Atomic(head)+tailAtomic,poolSize/size);
             return v;
           }
         if (diff<0)
@@ -236,46 +239,46 @@ namespace ecolab
    {return *sycl::ext::oneapi::group_local_memory_for_overwrite<LocalAllocatorBuffer>(syclGroup());}
 
   
-//#ifdef __SYCL_DEVICE_ONLY__
-//  
-//  /**
-//     A Local allocator allocates memory from device local memory,
-//     which is shared between threads of a work group, and has the same
-//     lifetime as the kernel
-//  */
-//  template <class T>
-//  class LocalAllocator
-//  {
-//  public:
-//    using value_type=T;
-//    using pointer=T*;
-//    using reference=T&;
-//    using difference_type=std::ptrdiff_t;
-//    using propagate_on_container_move_assignment=std::true_type;
-//
-//    // no need for destructor, as Impl has nothing to tear down
-//    T* allocate(size_t n) {
-//      auto& b=localAllocatorBuffer();
-//      unsigned offs=b.next;
-//      if (offs+n*sizeof(T)>LocalAllocatorSize)
-//        {
-//          fatalErrorFlag()=true;
-//          return nullptr;
-//        }
-//      sycl::group_barrier(syclGroup());
-//      if (syclGroup().leader()) b.next+=n*sizeof(T);
-//      sycl::group_barrier(syclGroup());
-//      char* alloc=b.buffer+offs;
-//      return reinterpret_cast<T*>(alloc);
-//    }
-//    void deallocate(T*p,size_t) {if (groupLeader()) printf("local dealloc %p\n",p);} // cleaned up when group exits
-//    template<class U> struct rebind {using other=LocalAllocator<U>;};
-//    // allocator is stateless
-//    bool operator==(const LocalAllocator&) const {return true;}
-//  };
-//#else
-//  template <class T> using LocalAllocator=std::allocator<T>;
-//#endif
+#ifdef __SYCL_DEVICE_ONLY__
+  
+  /**
+     A Local allocator allocates memory from device local memory,
+     which is shared between threads of a work group, and has the same
+     lifetime as the kernel
+  */
+  template <class T>
+  class LocalAllocator
+  {
+  public:
+    using value_type=T;
+    using pointer=T*;
+    using reference=T&;
+    using difference_type=std::ptrdiff_t;
+    using propagate_on_container_move_assignment=std::true_type;
+
+    // no need for destructor, as Impl has nothing to tear down
+    T* allocate(size_t n) {
+      auto& b=localAllocatorBuffer();
+      unsigned offs=b.next;
+      if (offs+n*sizeof(T)>LocalAllocatorSize)
+        {
+          fatalErrorFlag()=true;
+          return nullptr;
+        }
+      sycl::group_barrier(syclGroup());
+      if (syclGroup().leader()) b.next+=n*sizeof(T);
+      sycl::group_barrier(syclGroup());
+      char* alloc=b.buffer+offs;
+      return reinterpret_cast<T*>(alloc);
+    }
+    void deallocate(T*p,size_t) {} // cleaned up when group exits
+    template<class U> struct rebind {using other=LocalAllocator<U>;};
+    // allocator is stateless
+    bool operator==(const LocalAllocator&) const {return true;}
+  };
+#else
+  template <class T> using LocalAllocator=std::allocator<T>;
+#endif
    
 }
 #endif
diff --git a/include/arrays.h b/include/arrays.h
index f179bfe..d95edb4 100644
--- a/include/arrays.h
+++ b/include/arrays.h
@@ -1561,11 +1561,12 @@ namespace ecolab
         groupBarrier();
         if (dt)
           {
-            bool freeMem=true; //ref()==1;
+            unsigned refCnt=ref();
 #ifdef __SYCL_DEVICE_ONLY__
-            freeMem=sycl::group_broadcast(syclGroup(),freeMem);
+            // ensure absolute consistency of the reference count
+            refCnt=sycl::reduce_over_group(syclGroup(),refCnt,sycl::minimum<unsigned>());
 #endif
-            if (freeMem)
+            if (refCnt==1)
               {
                 free(dt);
                 dt=nullptr;
@@ -1596,11 +1597,10 @@ namespace ecolab
       
       void copy() //any nonconst method needs to call this
       {           // to implement copy-on-write semantics
-        if (dt /*&& ref()>1*/)
+        if (dt && ref()>1)
           {
             array_data<T>* oldData=dt;
             decrRef();
-            //bool freeMem=ref()==0;
             auto sz=size();
             dt=alloc(sz);
 #ifdef __SYCL_DEVICE_ONLY__
@@ -1608,7 +1608,6 @@ namespace ecolab
 #else
             memcpy(dt->dt,oldData->dt,sz*sizeof(T));
 #endif
-            //if (freeMem)
             free(oldData);
           }
       }
@@ -1633,10 +1632,8 @@ namespace ecolab
 
       array(const array& x): m_allocator(x.m_allocator) 
       {
-        //dt=x.dt;
-        //incrRef();
-        set_size(x.size());
-        asg_v(data(),x.size(),x);
+        dt=x.dt;
+        incrRef();
       }
 
       template <class expr>
@@ -1669,9 +1666,8 @@ namespace ecolab
         if (!dt || s>dt->sz /*|| ref()>1*/)
           {
             array tmp(s,m_allocator);
-            if (dt && copy) asg_v(tmp.dt->dt,std::min(s,tmp.size()),dt->dt);
+            if (dt && copy) asg_v(tmp.dt->dt,std::min(s,dt->sz),dt->dt);
             swap(tmp);
-            if (groupLeader()) printf("resize: deallocating %p, refCnt=%u\n",tmp.dt,tmp.refCnt());
           } 
         groupBarrier();
         if (groupLeader() && dt) dt->sz=s; // in case s is smaller
@@ -1688,29 +1684,24 @@ namespace ecolab
 
       void swap(array& x) {
 #ifndef __SYCL_DEVICE_ONLY__
-        printf("host swapping %p & %p\n",dt, x.dt);
         std::swap(dt, x.dt);
         std::swap(m_allocator,x.m_allocator);
 #else
-//        if (onDevice && x.onDevice) {
-//          if (groupLeader()) printf("on device swapping %p & %p\n",dt, x.dt);
-//          std::swap(dt, x.dt);
-//          std::swap(m_allocator,x.m_allocator);
-//        } else {
+        if (onDevice && x.onDevice) {
+          std::swap(dt, x.dt);
+          std::swap(m_allocator,x.m_allocator);
+        } else {
           // assumption here is these array may be per thread, or
           // maybe shared by all threads in a group, hence std::swap as above won't work
           auto lhs=dt, rhs=x.dt;
           auto lalloc=m_allocator, ralloc=x.m_allocator;
           groupBarrier();
-          //lhs=sycl::group_broadcast(syclGroup(),lhs);
-          //rhs=sycl::group_broadcast(syclGroup(),rhs);
-          if (groupLeader()) printf("swapping %p & %p\n",lhs,rhs);
           dt=rhs;
           x.dt=lhs;
           m_allocator=ralloc;
           x.m_allocator=lalloc;
           groupBarrier();
-          //        }
+        }
 #endif
       }
     
@@ -1727,13 +1718,13 @@ namespace ecolab
 
       array& operator=(const array& x) {
         if (x.dt==dt) return *this;
-//        if (m_allocator==x.m_allocator) {
-//          release();
-//          /*if (groupLeader()||onDevice)*/ {
-//            dt=x.dt;
-//          }
-//          incrRef();
-//        } else
+        if (m_allocator==x.m_allocator) {
+          release();
+          if (groupLeader()||onDevice) {
+            dt=x.dt;
+          }
+          incrRef();
+        } else
           asgV(x.size(), x);
         return *this;
       }
@@ -1807,17 +1798,7 @@ namespace ecolab
       typename enable_if<is_expression<E>,array&>::T
       operator<<=(const E& x) {
         auto origSize=size();
-#ifdef __SYCL_DEVICE_ONLY__
-        if (groupLeader()) printf("<<= resize from %zu to %zu\n",origSize,origSize+x.size());
-#else
-        printf("<<= host resize from %zu to %zu\n",origSize,origSize+x.size());
-#endif
         resize(origSize+x.size());
-#ifdef __SYCL_DEVICE_ONLY__
-        if (groupLeader()) printf("<<= after resize from %zu to %zu\n",origSize,origSize+x.size());
-#else
-        printf("<<= host after resize from %zu to %zu\n",origSize,origSize+x.size());
-#endif
         asg_v(data()+origSize,x.size(),x);
         return *this;
       }
diff --git a/include/ecolab.h b/include/ecolab.h
index aca28ae..673d354 100644
--- a/include/ecolab.h
+++ b/include/ecolab.h
@@ -147,6 +147,7 @@ namespace ecolab
           }
         });
       });
+      syclQ().wait_and_throw();
 #else
       hostForAll(f);
 #endif
diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc
index 4706bb5..ac24fab 100644
--- a/models/ecolab_model.cc
+++ b/models/ecolab_model.cc
@@ -89,8 +89,8 @@ RoundArray<E,EcolabPoint> EcolabPoint::roundArray(const E& expr)
 
 void EcolabPoint::generate(unsigned niter, const ModelData& model)
 {
-  //array<int,LocalAllocator<int>> lDensity(density), tmp(density.size());
-  array<int,Allocator<int>> lDensity(density), tmp(density.size(), density.allocator());
+  array<int,LocalAllocator<int>> lDensity(density), tmp(density.size());
+  //array<int,Allocator<int>> lDensity(density), tmp(density.size(), density.allocator());
   //auto& lDensity=density;
   //array<int,GlobalDeviceAllocator<int>> tmp(density.size(), density.allocator());
   
@@ -128,8 +128,8 @@ void EcolabPoint::condense(const ModelData::BoolArray& mask, size_t mask_true)
     density.clear();
     return;
   }
-  //LocalArray tmp(mask_true);
-  UnsignedArray tmp(mask_true,density.allocator());
+  LocalArray tmp(mask_true);
+  //UnsignedArray tmp(mask_true,density.allocator());
   if (groupLeader())
     for (size_t i=0, j=0; i<density.size(); ++i)
       if (mask[i] && j<mask_true) 
@@ -221,18 +221,19 @@ void SpatialModel::mutate()
   assert(all(*mut_scale<=1));
   last_mut_tstep=tstep;
 
+  // this bit of merde is because this line of code needs to compile
+  // in kernel code, even though it runs on the host.
   auto deviceAllocator=cell(0,0).density.allocator();
   vector<EcolabPoint::UnsignedArray,ModelData::Allocator<EcolabPoint::UnsignedArray>>
-    newSp(size(),{deviceAllocator});
-  
+    newSp(size(),EcolabPoint::UnsignedArray(deviceAllocator));
+ 
   groupedForAll([newSp=newSp.data(),mut_scale=&*mut_scale,this](EcolabCell& c,size_t i) {
     assert(all(c.density>=0));
     newSp[i]=c.mutate(*mut_scale);
   });
 
   array<unsigned> new_sp;
-  DeviceType<EcolabPoint::UnsignedArray> cell_ids(deviceAllocator);
-  syncThreads();
+  DeviceType<array<unsigned,ModelData::Allocator<unsigned>>> cell_ids;
   
   // TODO - this is a kind of scan - can it be done on device?
   size_t j=0;
@@ -244,12 +245,14 @@ void SpatialModel::mutate()
     }
 
   // deallocate on device
-  hostForAll([newSp=newSp.data()](EcolabCell& c,size_t i) {
+  //cout<<"deallocate"<<endl;
+  groupedForAll([newSp=newSp.data()](EcolabCell& c,size_t i) {
     newSp[i].clear();
     assert(newSp[i].refCnt()==0);
   });
    
 
+  //cout<<"ModelData::mutate"<<endl;
 #ifdef MPI_SUPPORT
   MPIbuf b; b<<new_sp<<(*cell_ids); b.gather(0);
   if (myid()==0)
@@ -268,20 +271,18 @@ void SpatialModel::mutate()
 #else
   ModelData::mutate(new_sp);
 #endif
-  //if (new_sp.size()==0) return;
+  if (new_sp.size()==0) return;
 
-  //  computeODiagIdx();
+  //cout<<"computeODiagIdx"<<endl;
+  computeODiagIdx();
   mut_scale->clear();
   newSp.clear();
 
-  //(*cell_ids)<<=0;
-  cout<<"b4 append:"<<endl;
+  //cout<<"set 1"<<endl;
   // set the new species density to 1 for those created on this cell
-  hostForAll([cell_ids=&*cell_ids](EcolabCell& c,size_t) {
-    //hostForAll([cell_ids=&*cell_ids,this](EcolabCell& c,size_t) {
+  groupedForAll([cell_ids=&*cell_ids](EcolabCell& c,size_t) {
     c.density <<= (*cell_ids)==c.id;
   });
-  cout<<"after append:"<<endl;
 }
 
 template <class E>
@@ -290,26 +291,31 @@ EcolabPoint::UnsignedArray EcolabPoint::mutate(const E& mut_scale)
   /* calculate the number of mutants each species produces */
   if (density.size()==0) return {density.allocator()};
 #ifdef __SYCL_DEVICE_ONLY__
-  //LocalArray speciations=roundArray(mut_scale * density);
-  UnsignedArray speciations(roundArray(mut_scale * density), density.allocator());
+  LocalArray speciations=roundArray(mut_scale * density);
+  //UnsignedArray speciations(roundArray(mut_scale * density), density.allocator());
   auto nsp=density.size();
-  //  auto new_sp = gen_index(speciations);
-  //LocalArray offsets(nsp+1);
-  UnsignedArray offsets(nsp+1,density.allocator());
-  sycl::joint_exclusive_scan(syclGroup(),speciations.data(),speciations.data()+nsp,
-                             offsets.data(),sycl::plus<unsigned>());
+//  //  auto new_sp = gen_index(speciations);
+  LocalArray offsets(nsp+1);
+  //UnsignedArray offsets(nsp+1,density.allocator());
+  unsigned* offs_p=offsets.data();
+  const unsigned* sp_p=speciations.data();
+  sycl::joint_exclusive_scan(syclGroup(),sp_p,sp_p+nsp,offs_p,sycl::plus<unsigned>());
   groupBarrier();
-  if (groupLeader())
-    offsets[nsp]=offsets[nsp-1]+speciations[nsp-1];
+  if (groupLeader()) {
+    // do not do array operations: data(), operator[] when not in full
+    // group scope, because of COW semantics
+    offs_p[nsp]=offs_p[nsp-1]+sp_p[nsp-1];
+  }
   groupBarrier();
-
-  if (offsets[nsp]==0) return {density.allocator()};
+  auto numSpeciations=offs_p[nsp];
+  
+  if (numSpeciations==0) return {density.allocator()};
   
   density-=speciations;
 
-  UnsignedArray new_sp(offsets[nsp], density.allocator());
-  array_ns::map(nsp, [offsets=offsets.data(),new_sp=new_sp.data()](size_t i) {
-    for (auto j=offsets[i]; j<offsets[i+1]; ++j)
+  UnsignedArray new_sp(numSpeciations, density.allocator());
+  array_ns::map(nsp, [offs_p,new_sp=new_sp.data()](size_t i) {
+    for (auto j=offs_p[i]; j<offs_p[i+1]; ++j)
       new_sp[j]=i;
   });
   
diff --git a/models/ecolab_model.h b/models/ecolab_model.h
index 6cd2fc8..8a95e63 100644
--- a/models/ecolab_model.h
+++ b/models/ecolab_model.h
@@ -82,7 +82,7 @@ class EcolabPoint
   template <class T> using Allocator=std::allocator<T>;
 #endif
   using UnsignedArray=array<unsigned,Allocator<unsigned>>;
-  //using LocalArray=array<unsigned,LocalAllocator<unsigned>>;
+  using LocalArray=array<unsigned,LocalAllocator<unsigned>>;
 
   Float salt;  /* random no. used for migration */
   array<int,Allocator<int>> density;
diff --git a/models/spatial_ecolab.py b/models/spatial_ecolab.py
index 82b9638..c3f9e71 100644
--- a/models/spatial_ecolab.py
+++ b/models/spatial_ecolab.py
@@ -21,8 +21,7 @@
 ecolab.repro_max(0.1)
 ecolab.odiag_min(-1e-5)
 ecolab.odiag_max(1e-5)
-ecolab.mut_max(1000)
-#ecolab.mut_max(1e-3)
+ecolab.mut_max(1e-3)
 ecolab.sp_sep(0.1)
 
 def randomList(num, min, max):
@@ -30,10 +29,10 @@ def randomList(num, min, max):
 
 ecolab.species(range(nsp))
 
-#numX=12
-#numY=12
-numX=1
-numY=1
+numX=12
+numY=12
+numX=2
+numY=2
 ecolab.setGrid(numX,numY)
 ecolab.partitionObjects()
 
@@ -64,7 +63,9 @@ def randomList(num, min, max):
 extinctions=0
 migrations=0
 def stepImpl():
-    #ecolab.generate(100)
+    #print("b4 generate")
+    ecolab.generate(100)
+    #print("b4 mutate")
     ecolab.mutate()
 
     epochTs=ecolab.tstep()%epoch
@@ -75,7 +76,8 @@ def stepImpl():
 
     global extinctions, migrations
     #migrations+=ecolab.migrate()
-    #extinctions+=ecolab.condense()
+    #print("b4 condense")
+    extinctions+=ecolab.condense()
     #print(ecolab.nsp()())
     #ecolab.syncThreads()
     #ecolab.gather()
@@ -87,7 +89,7 @@ def stepImpl():
 
 from timeit import timeit
 print('-----')
-print(timeit('ecolab.mutate()', globals=globals(), number=1))
+print(timeit('stepImpl()', globals=globals(), number=1))
 print('-----')
                
 def step():
@@ -97,9 +99,9 @@ def step():
     for i in range(epoch//1000000):
         stepImpl()
 #    #print('migrations=',migrations,' extinctions=',extinctions)
-#    if myid()==0:
-#        nsp=len(ecolab.species)
-#        statusBar.configure(text=f't={ecolab.tstep()} nsp:{nsp}')
+    if myid()==0:
+        nsp=len(ecolab.species)
+        statusBar.configure(text=f't={ecolab.tstep()} nsp:{nsp}')
 #        plot('No. species',ecolab.tstep(),nsp,200*(ecolab.tstep()%epoch<0.5*epoch))
 #        #plot('No. species',ecolab.tstep(),nsp)
 #        plot('No. species by cell',ecolab.tstep(),ecolab.nsp()())

From dc57275c8a5a623863276c44eb4e4d5eadcd23b2 Mon Sep 17 00:00:00 2001
From: Russell Standish <hpcoder@hpcoders.com.au>
Date: Wed, 29 Jul 2026 15:15:05 +1000
Subject: [PATCH 09/19] Switching to dell laptop to perform debugging.

---
 include/DeviceAllocator.h |  2 +-
 models/ecolab_model.cc    | 33 ++++++++++++++++++++++++---------
 models/spatial_ecolab.py  |  4 ++--
 3 files changed, 27 insertions(+), 12 deletions(-)

diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h
index d7b685b..263e544 100644
--- a/include/DeviceAllocator.h
+++ b/include/DeviceAllocator.h
@@ -224,7 +224,7 @@ namespace ecolab
     bool operator==(const HostSharedAllocator&) const {return true;}
   };
 
-  constexpr static unsigned LocalAllocatorSize=30*1024; // 32KiB = half typical local storage
+  constexpr static unsigned LocalAllocatorSize=8*1024; // 32KiB = half typical local storage
 
   struct LocalAllocatorBuffer
   {
diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc
index ac24fab..9eb5335 100644
--- a/models/ecolab_model.cc
+++ b/models/ecolab_model.cc
@@ -89,17 +89,20 @@ RoundArray<E,EcolabPoint> EcolabPoint::roundArray(const E& expr)
 
 void EcolabPoint::generate(unsigned niter, const ModelData& model)
 {
-  array<int,LocalAllocator<int>> lDensity(density), tmp(density.size());
+  //array<int,LocalAllocator<int>> lDensity(density), tmp(density.size());
   //array<int,Allocator<int>> lDensity(density), tmp(density.size(), density.allocator());
-  //auto& lDensity=density;
-  //array<int,GlobalDeviceAllocator<int>> tmp(density.size(), density.allocator());
+  auto& lDensity=density;
+  array<int,GlobalDeviceAllocator<int>> tmp(density.size(), density.allocator());
   
   for (unsigned step=0; step<niter; step++)
     {
       array_ns::map(lDensity.size(),  [&](size_t i){
         Float ir=model.interaction.diag[i]*lDensity[i];
-        for (auto& j: model.oDiagIdx[i])
+        //for (auto j: model.oDiagIdx[i])
+        for (size_t k=0; k<model.oDiagIdx[i].size(); ++k) {
+          unsigned j=model.oDiagIdx[i][k];
           ir+=model.interaction.val[j]*lDensity[model.interaction.col[j]];
+        }
         tmp[i]=ROUND(lDensity[i] + lDensity[i] * (model.repro_rate[i] + ir));
       });
       groupBarrier(); // synchronises threads on each iteration
@@ -282,6 +285,7 @@ void SpatialModel::mutate()
   // set the new species density to 1 for those created on this cell
   groupedForAll([cell_ids=&*cell_ids](EcolabCell& c,size_t) {
     c.density <<= (*cell_ids)==c.id;
+    assert(all(c.density>=0));
   });
 }
 
@@ -290,16 +294,27 @@ EcolabPoint::UnsignedArray EcolabPoint::mutate(const E& mut_scale)
 {
   /* calculate the number of mutants each species produces */
   if (density.size()==0) return {density.allocator()};
-#ifdef __SYCL_DEVICE_ONLY__
-  LocalArray speciations=roundArray(mut_scale * density);
-  //UnsignedArray speciations(roundArray(mut_scale * density), density.allocator());
+  //#ifdef __SYCL_DEVICE_ONLY__
+#if 1
+  //LocalArray speciations=roundArray(mut_scale * density);
+  UnsignedArray speciations(roundArray(mut_scale * density), density.allocator());
   auto nsp=density.size();
 //  //  auto new_sp = gen_index(speciations);
-  LocalArray offsets(nsp+1);
-  //UnsignedArray offsets(nsp+1,density.allocator());
+  //LocalArray offsets(nsp+1);
+  UnsignedArray offsets(nsp+1,density.allocator());
   unsigned* offs_p=offsets.data();
   const unsigned* sp_p=speciations.data();
+  //#ifdef __SYCL_DEVICE_ONLY__
+#if 0
   sycl::joint_exclusive_scan(syclGroup(),sp_p,sp_p+nsp,offs_p,sycl::plus<unsigned>());
+#else
+  if (groupLeader())
+    {
+      // sequential version for diagnostic purposes
+      offs_p[0]=0;
+      for (size_t i=0; i<nsp; ++i) offs_p[i+1]=offs_p[i]+sp_p[i];
+    }
+#endif
   groupBarrier();
   if (groupLeader()) {
     // do not do array operations: data(), operator[] when not in full
diff --git a/models/spatial_ecolab.py b/models/spatial_ecolab.py
index c3f9e71..fca068b 100644
--- a/models/spatial_ecolab.py
+++ b/models/spatial_ecolab.py
@@ -31,8 +31,8 @@ def randomList(num, min, max):
 
 numX=12
 numY=12
-numX=2
-numY=2
+#numX=2
+#numY=2
 ecolab.setGrid(numX,numY)
 ecolab.partitionObjects()
 

From 2de48bda874c9e7ddd9fc2ea9b3ccd61bcee1527 Mon Sep 17 00:00:00 2001
From: Russell Standish <hpcoder@hpcoders.com.au>
Date: Thu, 30 Jul 2026 12:15:31 +1000
Subject: [PATCH 10/19] Implementing DeviceAllocator using a stack rather than
 a queue, but using a cas loop to avoid race condition of empty stack test.

---
 include/DeviceAllocator.h | 91 +++++++++++++++++++++------------------
 models/ecolab_model.cc    |  4 +-
 models/spatial_ecolab.py  |  4 +-
 3 files changed, 54 insertions(+), 45 deletions(-)

diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h
index 263e544..c867add 100644
--- a/include/DeviceAllocator.h
+++ b/include/DeviceAllocator.h
@@ -47,7 +47,7 @@ namespace ecolab
     Slot slots[size];
     unsigned head=size, tail=0;
 
-    using Atomic=sycl::atomic_ref<unsigned,sycl::memory_order::relaxed,sycl::memory_scope::device>;
+    using Atomic=sycl::atomic_ref<unsigned,sycl::memory_order::seq_cst,sycl::memory_scope::device>;
 
   public:
     void init() {
@@ -60,53 +60,60 @@ namespace ecolab
     
     void enqueue(unsigned x)
     {
-//      Atomic t(tail);
-//      slots[--t].value=x;
-      while (true)
-      {
-        Atomic headAtomic(head);
-        unsigned pos=headAtomic.load();
-        Slot& slot=slots[pos & mask];
-        Atomic seqAtomic(slot.seq);
-        unsigned seq=seqAtomic.load(sycl::memory_order::acquire);
-        int diff=int(seq)-int(pos);
-
-        if (diff==0 && headAtomic.compare_exchange_strong(pos,pos+1))
-          {
-            slot.value=x;
-            Atomic publish(slot.seq);
-            publish.store(pos+1,sycl::memory_order::release);
-            return;
-          }
-      }
+      Atomic t(tail);
+      slots[--t].value=x;
+//      while (true)
+//      {
+//        Atomic headAtomic(head);
+//        unsigned pos=headAtomic.load();
+//        Slot& slot=slots[pos & mask];
+//        Atomic seqAtomic(slot.seq);
+//        unsigned seq=seqAtomic.load(sycl::memory_order::acquire);
+//        int diff=int(seq)-int(pos);
+//
+//        if (diff==0 && headAtomic.compare_exchange_strong(pos,pos+1))
+//          {
+//            slot.value=x;
+//            Atomic publish(slot.seq);
+//            publish.store(pos+1,sycl::memory_order::release);
+//            return;
+//          }
+//      }
     }
 
     unsigned dequeue()
     {
-//      Atomic t(tail);
-//      if (t>=size) return ~0;
-//      return slots[t++].value;
+      Atomic t(tail);
+      unsigned pos=t;
+      // updating tail in a cas loop avoids the race condition
+      // between the test and increment
       while (true)
-      {
-        Atomic tailAtomic(tail);
-        unsigned pos=tailAtomic.load();
-        Slot& slot=slots[pos & mask];
-        Atomic seqAtomic(slot.seq);
-        unsigned seq=seqAtomic.load(sycl::memory_order::acquire);
-        int diff=int(seq)-int(pos+1);
-
-        if (diff==0 && tailAtomic.compare_exchange_strong(pos,pos+1))
-          {
-            unsigned v=slot.value;
-            Atomic release(slot.seq);
-            release.store(pos+size,sycl::memory_order::release);
-            return v;
-          }
-        if (diff<0)
         {
-          return ~0U; // signal buffer empty, don't wait
+          if (pos>=size) return ~0; // stack empty
+          if (t.compare_exchange_weak(pos,pos+1))
+            return slots[pos].value;
         }
-      }
+//      while (true)
+//      {
+//        Atomic tailAtomic(tail);
+//        unsigned pos=tailAtomic.load();
+//        Slot& slot=slots[pos & mask];
+//        Atomic seqAtomic(slot.seq);
+//        unsigned seq=seqAtomic.load(sycl::memory_order::acquire);
+//        int diff=int(seq)-int(pos+1);
+//
+//        if (diff==0 && tailAtomic.compare_exchange_strong(pos,pos+1))
+//          {
+//            unsigned v=slot.value;
+//            Atomic release(slot.seq);
+//            release.store(pos+size,sycl::memory_order::release);
+//            return v;
+//          }
+//        if (diff<0)
+//        {
+//          return ~0U; // signal buffer empty, don't wait
+//        }
+//      }
     }
   };
 
@@ -152,6 +159,8 @@ namespace ecolab
         if (localThreadId()==0) offs=queue.dequeue();
 #ifdef __SYCL_DEVICE_ONLY__
         offs=sycl::group_broadcast(syclGroup(),offs,0);
+        if (groupLeader()) printf("alloc pageSize=%u offs=%u on group %u\n",
+                                  pageSize,offs,syclGroup().get_group_linear_id());
 #endif
         if (offs!=~0U)
           return memory+(offs<<order);
diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc
index 9eb5335..a565624 100644
--- a/models/ecolab_model.cc
+++ b/models/ecolab_model.cc
@@ -304,8 +304,8 @@ EcolabPoint::UnsignedArray EcolabPoint::mutate(const E& mut_scale)
   UnsignedArray offsets(nsp+1,density.allocator());
   unsigned* offs_p=offsets.data();
   const unsigned* sp_p=speciations.data();
-  //#ifdef __SYCL_DEVICE_ONLY__
-#if 0
+#ifdef __SYCL_DEVICE_ONLY__
+  //#if 0
   sycl::joint_exclusive_scan(syclGroup(),sp_p,sp_p+nsp,offs_p,sycl::plus<unsigned>());
 #else
   if (groupLeader())
diff --git a/models/spatial_ecolab.py b/models/spatial_ecolab.py
index fca068b..c3f9e71 100644
--- a/models/spatial_ecolab.py
+++ b/models/spatial_ecolab.py
@@ -31,8 +31,8 @@ def randomList(num, min, max):
 
 numX=12
 numY=12
-#numX=2
-#numY=2
+numX=2
+numY=2
 ecolab.setGrid(numX,numY)
 ecolab.partitionObjects()
 

From 5a4acc2ffbf82f3b2dcbade1e89571f08aca5aa0 Mon Sep 17 00:00:00 2001
From: Russell Standish <hpcoder@hpcoders.com.au>
Date: Thu, 30 Jul 2026 14:07:53 +1000
Subject: [PATCH 11/19] Revert to MPMC.

---
 include/DeviceAllocator.h | 99 ++++++++++++++++++++-------------------
 1 file changed, 50 insertions(+), 49 deletions(-)

diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h
index c867add..f50cfbd 100644
--- a/include/DeviceAllocator.h
+++ b/include/DeviceAllocator.h
@@ -45,8 +45,8 @@ namespace ecolab
     };
 
     Slot slots[size];
-    unsigned head=size, tail=0;
-
+    unsigned head=size, tail=0, pushing=0;
+    
     using Atomic=sycl::atomic_ref<unsigned,sycl::memory_order::seq_cst,sycl::memory_scope::device>;
 
   public:
@@ -60,60 +60,61 @@ namespace ecolab
     
     void enqueue(unsigned x)
     {
-      Atomic t(tail);
-      slots[--t].value=x;
-//      while (true)
-//      {
-//        Atomic headAtomic(head);
-//        unsigned pos=headAtomic.load();
-//        Slot& slot=slots[pos & mask];
-//        Atomic seqAtomic(slot.seq);
-//        unsigned seq=seqAtomic.load(sycl::memory_order::acquire);
-//        int diff=int(seq)-int(pos);
-//
-//        if (diff==0 && headAtomic.compare_exchange_strong(pos,pos+1))
-//          {
-//            slot.value=x;
-//            Atomic publish(slot.seq);
-//            publish.store(pos+1,sycl::memory_order::release);
-//            return;
-//          }
-//      }
+//      Atomic t(tail), p(pushing);
+//      slots[--t].value=x;
+      while (true)
+      {
+        Atomic headAtomic(head);
+        unsigned pos=headAtomic.load();
+        Slot& slot=slots[pos & mask];
+        Atomic seqAtomic(slot.seq);
+        unsigned seq=seqAtomic.load(sycl::memory_order::acquire);
+        int diff=int(seq)-int(pos);
+
+        if (diff==0 && headAtomic.compare_exchange_strong(pos,pos+1))
+          {
+            slot.value=x;
+            Atomic publish(slot.seq);
+            publish.store(pos+1,sycl::memory_order::release);
+            return;
+          }
+      }
     }
 
     unsigned dequeue()
     {
-      Atomic t(tail);
-      unsigned pos=t;
-      // updating tail in a cas loop avoids the race condition
-      // between the test and increment
-      while (true)
-        {
-          if (pos>=size) return ~0; // stack empty
-          if (t.compare_exchange_weak(pos,pos+1))
-            return slots[pos].value;
-        }
+//      Atomic t(tail);
+//      unsigned pos=t;
+//      // updating tail in a cas loop avoids the race condition
+//      // between the test and increment
 //      while (true)
-//      {
-//        Atomic tailAtomic(tail);
-//        unsigned pos=tailAtomic.load();
-//        Slot& slot=slots[pos & mask];
-//        Atomic seqAtomic(slot.seq);
-//        unsigned seq=seqAtomic.load(sycl::memory_order::acquire);
-//        int diff=int(seq)-int(pos+1);
-//
-//        if (diff==0 && tailAtomic.compare_exchange_strong(pos,pos+1))
-//          {
-//            unsigned v=slot.value;
-//            Atomic release(slot.seq);
-//            release.store(pos+size,sycl::memory_order::release);
-//            return v;
-//          }
-//        if (diff<0)
 //        {
-//          return ~0U; // signal buffer empty, don't wait
+//          if (p)
+//          if (pos>=size) return ~0; // stack empty
+//          if (t.compare_exchange_weak(pos,pos+1))
+//            return slots[pos].value;
 //        }
-//      }
+      while (true)
+      {
+        Atomic tailAtomic(tail);
+        unsigned pos=tailAtomic.load();
+        Slot& slot=slots[pos & mask];
+        Atomic seqAtomic(slot.seq);
+        unsigned seq=seqAtomic.load(sycl::memory_order::acquire);
+        int diff=int(seq)-int(pos+1);
+
+        if (diff==0 && tailAtomic.compare_exchange_strong(pos,pos+1))
+          {
+            unsigned v=slot.value;
+            Atomic release(slot.seq);
+            release.store(pos+size,sycl::memory_order::release);
+            return v;
+          }
+        if (diff<0)
+        {
+          return ~0U; // signal buffer empty, don't wait
+        }
+      }
     }
   };
 

From f8f2a4055cb44828e1ab7963f44ba46bbebd817d Mon Sep 17 00:00:00 2001
From: Russell Standish <hpcoder@hpcoders.com.au>
Date: Thu, 30 Jul 2026 17:40:51 +1000
Subject: [PATCH 12/19] tweaks from dellGPU

---
 include/arrays.h | 1 +
 include/ecolab.h | 2 +-
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/include/arrays.h b/include/arrays.h
index d95edb4..cb2b617 100644
--- a/include/arrays.h
+++ b/include/arrays.h
@@ -1591,6 +1591,7 @@ namespace ecolab
 //#else
         array tmp(size,m_allocator);
         asg_v(tmp.data(),size,x);
+        groupBarrier();
         swap(tmp);
         //#endif
       }
diff --git a/include/ecolab.h b/include/ecolab.h
index 673d354..63d8147 100644
--- a/include/ecolab.h
+++ b/include/ecolab.h
@@ -205,7 +205,7 @@ namespace ecolab
             }
             // flag fatal error to throw afterwards.
             if (fatalErrorFlag())
-              sycl::atomic_ref<int,sycl::memory_order::relaxed,sycl::memory_scope::device>(*fatalError).fetch_or(1);
+              sycl::atomic_ref<int,sycl::memory_order::seq_cst,sycl::memory_scope::device>(*fatalError).fetch_or(1);
           });
         });
       syclQ().wait_and_throw();

From 3e7fb9e853da807b71b03c2dc8db01f9cf0ed7f4 Mon Sep 17 00:00:00 2001
From: Russell Standish <hpcoder@hpcoders.com.au>
Date: Sun, 2 Aug 2026 18:42:12 +1000
Subject: [PATCH 13/19] Discard stack

---
 include/DeviceAllocator.h | 151 ++++++++++++++++++++++----------------
 include/arrays.h          |  13 ++--
 include/ecolab.h          |   3 +-
 models/ecolab_model.cc    |   4 +-
 models/ecolab_model.h     |   1 +
 models/spatial_ecolab.py  |  12 +--
 6 files changed, 107 insertions(+), 77 deletions(-)

diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h
index f50cfbd..363bfe5 100644
--- a/include/DeviceAllocator.h
+++ b/include/DeviceAllocator.h
@@ -40,81 +40,95 @@ namespace ecolab
 
     struct Slot
     {
-      unsigned seq;
+      uint64_t seq;
       unsigned value;
     };
 
     Slot slots[size];
-    unsigned head=size, tail=0, pushing=0;
+    uint64_t head=0, tail=size, pushing=0;
     
-    using Atomic=sycl::atomic_ref<unsigned,sycl::memory_order::seq_cst,sycl::memory_scope::device>;
+    template <class T> using Atomic=sycl::atomic_ref<T,sycl::memory_order::acq_rel,sycl::memory_scope::device>;
 
   public:
     void init() {
+      tail=0;
       for (unsigned i=syclItem().get_global_linear_id(); i<size;
            i+=syclItem().get_global_range().size()) {
-        slots[i].value=i;
-        slots[i].seq=i+1;
-      }
+          slots[i].value=i;
+          slots[i].seq=i+1;
+        }
     }
     
     void enqueue(unsigned x)
     {
-//      Atomic t(tail), p(pushing);
-//      slots[--t].value=x;
-      while (true)
-      {
-        Atomic headAtomic(head);
-        unsigned pos=headAtomic.load();
-        Slot& slot=slots[pos & mask];
-        Atomic seqAtomic(slot.seq);
-        unsigned seq=seqAtomic.load(sycl::memory_order::acquire);
-        int diff=int(seq)-int(pos);
-
-        if (diff==0 && headAtomic.compare_exchange_strong(pos,pos+1))
-          {
-            slot.value=x;
-            Atomic publish(slot.seq);
-            publish.store(pos+1,sycl::memory_order::release);
-            return;
-          }
-      }
+      Atomic<uint64_t> t(tail), p(pushing);
+      slots[--t].value=x;
+//      for (bool published=false; !published;)
+//      {
+//        if (localThreadId()==0)
+//          {
+//            Atomic<uint64_t> headAtomic(head);
+//            auto pos=headAtomic.load();
+//            Slot& slot=slots[pos & mask];
+//            Atomic<uint64_t> seqAtomic(slot.seq);
+//            auto seq=seqAtomic.load(sycl::memory_order::acquire);
+//            auto diff=int64_t(seq)-int64_t(pos);
+//
+//            if (diff==0 && headAtomic.compare_exchange_weak(pos,pos+1))
+//              {
+//                Atomic<unsigned> valueAtomic(slot.value);
+//                valueAtomic=x;
+//                seqAtomic=pos+1;
+//                published=true;
+//              }
+//          }
+//#ifdef __SYCL_DEVICE_ONLY__
+//        sycl::atomic_fence(sycl::memory_order::seq_cst, sycl::memory_scope::device);
+//        published=sycl::group_broadcast(syclGroup(),published,0);
+//#endif
+//      }
     }
 
     unsigned dequeue()
     {
-//      Atomic t(tail);
-//      unsigned pos=t;
-//      // updating tail in a cas loop avoids the race condition
-//      // between the test and increment
-//      while (true)
-//        {
-//          if (p)
-//          if (pos>=size) return ~0; // stack empty
-//          if (t.compare_exchange_weak(pos,pos+1))
-//            return slots[pos].value;
+      Atomic<uint64_t> t(tail);
+      // updating tail in a cas loop avoids the race condition
+      // between the test and increment
+      uint64_t p=++t;
+      if (p>=size) {t=size; return ~0;} // stack empty
+      return slots[p].value;
+//      unsigned v=~0U-1;
+//      unsigned niter=0;
+//      while (v==~0U-1)
+//      {
+//        if (localThreadId()==0)
+//          {
+//            Atomic<uint64_t> tailAtomic(tail);
+//            auto pos=tailAtomic.load();
+//            Slot& slot=slots[pos & mask];
+//            Atomic<uint64_t> seqAtomic(slot.seq);
+//            auto seq=seqAtomic.load(sycl::memory_order::acquire);
+//            auto diff=int64_t(seq)-int64_t(pos+1);
+//
+//            if (diff==0 && tailAtomic.compare_exchange_weak(pos,pos+1))
+//              {
+//                Atomic<unsigned> valueAtomic(slot.value);
+//                v=valueAtomic;
+//                Atomic<uint64_t> release(slot.seq);
+//                release.store(pos+size,sycl::memory_order::release);
+//              }
+//            if (diff<0) v=~0U; // empty queue
+//          }
+//#ifdef __SYCL_DEVICE_ONLY__
+//        sycl::atomic_fence(sycl::memory_order::seq_cst, sycl::memory_scope::device);
+//        v=sycl::group_broadcast(syclGroup(),v,0);
+//#endif
+//        if (niter++>10) {
+//          break;
+//          if (groupLeader()) printf("enqueue failed\n",0);
 //        }
-      while (true)
-      {
-        Atomic tailAtomic(tail);
-        unsigned pos=tailAtomic.load();
-        Slot& slot=slots[pos & mask];
-        Atomic seqAtomic(slot.seq);
-        unsigned seq=seqAtomic.load(sycl::memory_order::acquire);
-        int diff=int(seq)-int(pos+1);
-
-        if (diff==0 && tailAtomic.compare_exchange_strong(pos,pos+1))
-          {
-            unsigned v=slot.value;
-            Atomic release(slot.seq);
-            release.store(pos+size,sycl::memory_order::release);
-            return v;
-          }
-        if (diff<0)
-        {
-          return ~0U; // signal buffer empty, don't wait
-        }
-      }
+//      }
+//      return v;
     }
   };
 
@@ -134,6 +148,7 @@ namespace ecolab
     }
     void deallocate(void* p, size_t) {sycl::ext::oneapi::experimental::printf("%p leaked on device\n",p);}
     void init() {}
+    void recycleDiscardPile() {}
   };
 
   template <unsigned order=minOrder> class DeviceAllocator
@@ -141,6 +156,7 @@ namespace ecolab
     constexpr static unsigned pageSize=1<<order;
     constexpr static unsigned numPages=poolSize/pageSize;
     Queue<numPages> queue;
+    Queue<numPages> discard;
     char memory[poolSize];
     DeviceAllocator<order+2> nextAllocator; // next size up allocator
   public:
@@ -151,18 +167,25 @@ namespace ecolab
                            [this](size_t) {queue.init();});
       nextAllocator.init();
     }
+    void recycleDiscardPile() {
+      for (auto v=discard.dequeue(); v!=~0U; v=discard.dequeue())
+        queue.enqueue(v);
+      nextAllocator.recycleDiscardPile();
+    }
     // all members of group get the same pointer
     void* allocate(size_t size) {
       if (size==0) return nullptr;
       if (size<=pageSize) {
         unsigned offs=~0U;
-        groupBarrier();
         if (localThreadId()==0) offs=queue.dequeue();
 #ifdef __SYCL_DEVICE_ONLY__
         offs=sycl::group_broadcast(syclGroup(),offs,0);
-        if (groupLeader()) printf("alloc pageSize=%u offs=%u on group %u\n",
-                                  pageSize,offs,syclGroup().get_group_linear_id());
 #endif
+       
+//#ifdef __SYCL_DEVICE_ONLY__
+//        if (groupLeader()) printf("alloc pageSize=%u offs=%u on group %u\n",
+//                                  pageSize,offs,syclGroup().get_group_linear_id());
+//#endif
         if (offs!=~0U)
           return memory+(offs<<order);
       }
@@ -171,10 +194,12 @@ namespace ecolab
     void deallocate(void* p, size_t size) {
       if (!p) return;
       if (p>=memory && p<memory+poolSize) {
-        groupBarrier();
-        if (groupLeader()) {
-          queue.enqueue((reinterpret_cast<char*>(p)-memory)>>order);
-        }
+#ifdef __SYCL_DEVICE_ONLY__
+        if (groupLeader())
+          discard.enqueue((reinterpret_cast<char*>(p)-memory)>>order);
+#else
+        queue.enqueue((reinterpret_cast<char*>(p)-memory)>>order);
+#endif
         return;
       }
       nextAllocator.deallocate(p,size);
diff --git a/include/arrays.h b/include/arrays.h
index cb2b617..8c224f2 100644
--- a/include/arrays.h
+++ b/include/arrays.h
@@ -1609,7 +1609,6 @@ namespace ecolab
 #else
             memcpy(dt->dt,oldData->dt,sz*sizeof(T));
 #endif
-            free(oldData);
           }
       }
 
@@ -1664,13 +1663,17 @@ namespace ecolab
       void resize(size_t s, bool copy) {
         if (s==size()) return;
         groupBarrier();
-        if (!dt || s>dt->sz /*|| ref()>1*/)
+        if (!dt || s>dt->sz || ref()>1)
           {
             array tmp(s,m_allocator);
-            if (dt && copy) asg_v(tmp.dt->dt,std::min(s,dt->sz),dt->dt);
+            if (dt && tmp.dt && copy) asg_v(tmp.dt->dt,std::min(s,dt->sz),dt->dt);
             swap(tmp);
-          } 
-        groupBarrier();
+            groupBarrier();
+          }
+#ifdef __SYCL_DEVICE_ONLY__
+        // assert all pointers are the same
+        assert(sycl::reduce_over_group(syclGroup(),size_t(dt),sycl::minimum<size_t>())==size_t(dt));
+#endif
         if (groupLeader() && dt) dt->sz=s; // in case s is smaller
         groupBarrier();
       } 
diff --git a/include/ecolab.h b/include/ecolab.h
index 63d8147..c2bf93e 100644
--- a/include/ecolab.h
+++ b/include/ecolab.h
@@ -183,7 +183,7 @@ namespace ecolab
       size_t wg_per_compute_unit = max_slm_size / LocalAllocatorSize;
       // To maximize latency hiding, it's often beneficial to double or triple this 
       // so the GPU can switch to a waiting wave while another wave is blocked by a barrier.
-      size_t num_work_groups = max_compute_units * wg_per_compute_unit;
+      size_t num_work_groups = max_compute_units;
 
       num_work_groups=std::min(num_work_groups,this->size());
       //std::cout<<max_slm_size<<" max_slm_size "<<max_compute_units<<" max_compute_units "<<num_work_groups<<" work groups of "<<workGroupSize<<" threads"<<std::endl;
@@ -209,6 +209,7 @@ namespace ecolab
           });
         });
       syclQ().wait_and_throw();
+      deviceAllocator().recycleDiscardPile();
       if (*fatalError)
         throw std::runtime_error("Local Allocator Exhausted");
 #else
diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc
index a565624..5a020a4 100644
--- a/models/ecolab_model.cc
+++ b/models/ecolab_model.cc
@@ -131,8 +131,8 @@ void EcolabPoint::condense(const ModelData::BoolArray& mask, size_t mask_true)
     density.clear();
     return;
   }
-  LocalArray tmp(mask_true);
-  //UnsignedArray tmp(mask_true,density.allocator());
+  //LocalArray tmp(mask_true);
+  UnsignedArray tmp(mask_true,density.allocator());
   if (groupLeader())
     for (size_t i=0, j=0; i<density.size(); ++i)
       if (mask[i] && j<mask_true) 
diff --git a/models/ecolab_model.h b/models/ecolab_model.h
index 8a95e63..bc94f0b 100644
--- a/models/ecolab_model.h
+++ b/models/ecolab_model.h
@@ -14,6 +14,7 @@ using classdesc::Object;
 #include <random>
 #include <vector>
 #include <pack_stl.h>
+#include "DeviceAllocator.h"
 
 #ifdef USE_FLOAT
 using Float=float;
diff --git a/models/spatial_ecolab.py b/models/spatial_ecolab.py
index c3f9e71..e4bdf2d 100644
--- a/models/spatial_ecolab.py
+++ b/models/spatial_ecolab.py
@@ -29,10 +29,10 @@ def randomList(num, min, max):
 
 ecolab.species(range(nsp))
 
-numX=12
-numY=12
-numX=2
-numY=2
+numX=8
+numY=8
+#numX=2
+#numY=2
 ecolab.setGrid(numX,numY)
 ecolab.partitionObjects()
 
@@ -64,7 +64,7 @@ def randomList(num, min, max):
 migrations=0
 def stepImpl():
     #print("b4 generate")
-    ecolab.generate(100)
+    ecolab.generate(1)
     #print("b4 mutate")
     ecolab.mutate()
 
@@ -77,7 +77,7 @@ def stepImpl():
     global extinctions, migrations
     #migrations+=ecolab.migrate()
     #print("b4 condense")
-    extinctions+=ecolab.condense()
+    #extinctions+=ecolab.condense()
     #print(ecolab.nsp()())
     #ecolab.syncThreads()
     #ecolab.gather()

From bdbdec3a0a3a1d32c36bf5f4b614498c9bd140c6 Mon Sep 17 00:00:00 2001
From: Russell Standish <hpcoder@hpcoders.com.au>
Date: Fri, 7 Aug 2026 16:23:41 +1000
Subject: [PATCH 14/19] GPU allocation working stably now. Note use of
 LocalArray in return of EcolabPoint::mutate.

---
 include/DeviceAllocator.h |  2 +-
 models/ecolab_model.cc    | 36 ++++++++++++++++++++++--------------
 models/ecolab_model.h     |  2 +-
 models/spatial_ecolab.py  |  6 +++---
 4 files changed, 27 insertions(+), 19 deletions(-)

diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h
index 363bfe5..ce67c1a 100644
--- a/include/DeviceAllocator.h
+++ b/include/DeviceAllocator.h
@@ -94,7 +94,7 @@ namespace ecolab
       Atomic<uint64_t> t(tail);
       // updating tail in a cas loop avoids the race condition
       // between the test and increment
-      uint64_t p=++t;
+      uint64_t p=t++;
       if (p>=size) {t=size; return ~0;} // stack empty
       return slots[p].value;
 //      unsigned v=~0U-1;
diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc
index 5a020a4..d9decb6 100644
--- a/models/ecolab_model.cc
+++ b/models/ecolab_model.cc
@@ -131,8 +131,8 @@ void EcolabPoint::condense(const ModelData::BoolArray& mask, size_t mask_true)
     density.clear();
     return;
   }
-  //LocalArray tmp(mask_true);
-  UnsignedArray tmp(mask_true,density.allocator());
+  LocalArray tmp(mask_true);
+  //UnsignedArray tmp(mask_true,density.allocator());
   if (groupLeader())
     for (size_t i=0, j=0; i<density.size(); ++i)
       if (mask[i] && j<mask_true) 
@@ -285,23 +285,28 @@ void SpatialModel::mutate()
   // set the new species density to 1 for those created on this cell
   groupedForAll([cell_ids=&*cell_ids](EcolabCell& c,size_t) {
     c.density <<= (*cell_ids)==c.id;
+    //auto oldSize=c.density.size();
+    //groupBarrier();
+    //c.density.resize(c.density.size()+cell_ids->size());
+    //groupBarrier();
+    //array_ns::asg_v(c.density.data()+oldSize, cell_ids->size(), (*cell_ids)==c.id);
     assert(all(c.density>=0));
   });
 }
 
 template <class E>
-EcolabPoint::UnsignedArray EcolabPoint::mutate(const E& mut_scale)
+EcolabPoint::LocalArray EcolabPoint::mutate(const E& mut_scale)
 {
   /* calculate the number of mutants each species produces */
-  if (density.size()==0) return {density.allocator()};
+  if (density.size()==0) return {};//{density.allocator()};
   //#ifdef __SYCL_DEVICE_ONLY__
 #if 1
-  //LocalArray speciations=roundArray(mut_scale * density);
-  UnsignedArray speciations(roundArray(mut_scale * density), density.allocator());
+  LocalArray speciations=roundArray(mut_scale * density);
+  //UnsignedArray speciations(roundArray(mut_scale * density), density.allocator());
   auto nsp=density.size();
 //  //  auto new_sp = gen_index(speciations);
-  //LocalArray offsets(nsp+1);
-  UnsignedArray offsets(nsp+1,density.allocator());
+  LocalArray offsets(nsp+1);
+  //UnsignedArray offsets(nsp+1,density.allocator());
   unsigned* offs_p=offsets.data();
   const unsigned* sp_p=speciations.data();
 #ifdef __SYCL_DEVICE_ONLY__
@@ -316,19 +321,22 @@ EcolabPoint::UnsignedArray EcolabPoint::mutate(const E& mut_scale)
     }
 #endif
   groupBarrier();
-  if (groupLeader()) {
+  unsigned numSpeciations=0;
+  if (localThreadId()==0) {
     // do not do array operations: data(), operator[] when not in full
     // group scope, because of COW semantics
-    offs_p[nsp]=offs_p[nsp-1]+sp_p[nsp-1];
+    numSpeciations=offs_p[nsp]=offs_p[nsp-1]+sp_p[nsp-1];
   }
-  groupBarrier();
-  auto numSpeciations=offs_p[nsp];
+#ifdef __SYCL_DEVICE_ONLY__
+  numSpeciations=sycl::group_broadcast(syclGroup(),numSpeciations,0);
+#endif
   
-  if (numSpeciations==0) return {density.allocator()};
+  if (numSpeciations==0) return {};//{density.allocator()};
   
   density-=speciations;
 
-  UnsignedArray new_sp(numSpeciations, density.allocator());
+  LocalArray new_sp(numSpeciations);
+  //UnsignedArray new_sp(numSpeciations, density.allocator());
   array_ns::map(nsp, [offs_p,new_sp=new_sp.data()](size_t i) {
     for (auto j=offs_p[i]; j<offs_p[i+1]; ++j)
       new_sp[j]=i;
diff --git a/models/ecolab_model.h b/models/ecolab_model.h
index bc94f0b..c7d6750 100644
--- a/models/ecolab_model.h
+++ b/models/ecolab_model.h
@@ -90,7 +90,7 @@ class EcolabPoint
   
   void generate(unsigned niter, const ModelData&);
   void condense(const ModelData::BoolArray& mask, size_t mask_true);
-  template <class E> UnsignedArray mutate(const E&);
+  template <class E> LocalArray mutate(const E&);
   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);
diff --git a/models/spatial_ecolab.py b/models/spatial_ecolab.py
index e4bdf2d..f75b6e3 100644
--- a/models/spatial_ecolab.py
+++ b/models/spatial_ecolab.py
@@ -21,7 +21,7 @@
 ecolab.repro_max(0.1)
 ecolab.odiag_min(-1e-5)
 ecolab.odiag_max(1e-5)
-ecolab.mut_max(1e-3)
+ecolab.mut_max(1e-5)
 ecolab.sp_sep(0.1)
 
 def randomList(num, min, max):
@@ -64,7 +64,7 @@ def randomList(num, min, max):
 migrations=0
 def stepImpl():
     #print("b4 generate")
-    ecolab.generate(1)
+    ecolab.generate(100)
     #print("b4 mutate")
     ecolab.mutate()
 
@@ -77,7 +77,7 @@ def stepImpl():
     global extinctions, migrations
     #migrations+=ecolab.migrate()
     #print("b4 condense")
-    #extinctions+=ecolab.condense()
+    extinctions+=ecolab.condense()
     #print(ecolab.nsp()())
     #ecolab.syncThreads()
     #ecolab.gather()

From da128ee40c23a41445ba110ff22d58d561c97e17 Mon Sep 17 00:00:00 2001
From: Russell Standish <hpcoder@hpcoders.com.au>
Date: Sat, 8 Aug 2026 14:15:58 +1000
Subject: [PATCH 15/19] revised onDevice() test, but swap algorithm still
 doesn't work.

---
 include/arrays.h       | 21 +++++++++++++++------
 models/ecolab_model.cc | 10 +++++-----
 models/ecolab_model.h  |  2 +-
 3 files changed, 21 insertions(+), 12 deletions(-)

diff --git a/include/arrays.h b/include/arrays.h
index 8c224f2..90a2e08 100644
--- a/include/arrays.h
+++ b/include/arrays.h
@@ -1473,9 +1473,6 @@ namespace ecolab
     {
       array_data<T> *dt=nullptr;
       A m_allocator;
-#ifdef SYCL_LANGUAGE_VERSION
-      bool onDevice=ecolab::onDevice(); // true if created in a parallel section
-#endif
       
       friend class WhereContext;
 
@@ -1686,15 +1683,27 @@ namespace ecolab
       template <class V>
       void resizeAndInit(size_t s, const V& val) {resize(s,false); operator=(val);}
 
+      bool onDevice() const {
+#ifdef SYCL_LANGUAGE_VERSION
+      // true if created in a parallel section
+      return sycl::address_space_cast
+        <sycl::access::address_space::private_space,
+         sycl::access::decorated::yes>(this).get();
+#else
+      return false;
+#endif
+       }
+      
       void swap(array& x) {
 #ifndef __SYCL_DEVICE_ONLY__
         std::swap(dt, x.dt);
         std::swap(m_allocator,x.m_allocator);
 #else
-        if (onDevice && x.onDevice) {
+        // this code still doesn't work, even on the revised onDevice test above
+        if (onDevice() && x.onDevice()) {
           std::swap(dt, x.dt);
           std::swap(m_allocator,x.m_allocator);
-        } else {
+          } else {
           // assumption here is these array may be per thread, or
           // maybe shared by all threads in a group, hence std::swap as above won't work
           auto lhs=dt, rhs=x.dt;
@@ -1724,7 +1733,7 @@ namespace ecolab
         if (x.dt==dt) return *this;
         if (m_allocator==x.m_allocator) {
           release();
-          if (groupLeader()||onDevice) {
+          if (groupLeader()||onDevice()) {
             dt=x.dt;
           }
           incrRef();
diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc
index d9decb6..e89f764 100644
--- a/models/ecolab_model.cc
+++ b/models/ecolab_model.cc
@@ -295,10 +295,10 @@ void SpatialModel::mutate()
 }
 
 template <class E>
-EcolabPoint::LocalArray EcolabPoint::mutate(const E& mut_scale)
+EcolabPoint::UnsignedArray EcolabPoint::mutate(const E& mut_scale)
 {
   /* calculate the number of mutants each species produces */
-  if (density.size()==0) return {};//{density.allocator()};
+  if (density.size()==0) return {density.allocator()};
   //#ifdef __SYCL_DEVICE_ONLY__
 #if 1
   LocalArray speciations=roundArray(mut_scale * density);
@@ -331,12 +331,12 @@ EcolabPoint::LocalArray EcolabPoint::mutate(const E& mut_scale)
   numSpeciations=sycl::group_broadcast(syclGroup(),numSpeciations,0);
 #endif
   
-  if (numSpeciations==0) return {};//{density.allocator()};
+  if (numSpeciations==0) return {density.allocator()};
   
   density-=speciations;
 
-  LocalArray new_sp(numSpeciations);
-  //UnsignedArray new_sp(numSpeciations, density.allocator());
+  //LocalArray new_sp(numSpeciations);
+  UnsignedArray new_sp(numSpeciations, density.allocator());
   array_ns::map(nsp, [offs_p,new_sp=new_sp.data()](size_t i) {
     for (auto j=offs_p[i]; j<offs_p[i+1]; ++j)
       new_sp[j]=i;
diff --git a/models/ecolab_model.h b/models/ecolab_model.h
index c7d6750..bc94f0b 100644
--- a/models/ecolab_model.h
+++ b/models/ecolab_model.h
@@ -90,7 +90,7 @@ class EcolabPoint
   
   void generate(unsigned niter, const ModelData&);
   void condense(const ModelData::BoolArray& mask, size_t mask_true);
-  template <class E> LocalArray mutate(const E&);
+  template <class E> UnsignedArray mutate(const E&);
   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);

From 4ba1c06ab4bbc2a27c8970b2f6ad6e88a3cc4ad8 Mon Sep 17 00:00:00 2001
From: Russell Standish <hpcoder@hpcoders.com.au>
Date: Sat, 8 Aug 2026 14:25:09 +1000
Subject: [PATCH 16/19] chore: code cleanup

---
 include/arrays.h       | 26 +-------------------------
 include/sycl.h         |  7 -------
 models/ecolab_model.cc | 31 ++-----------------------------
 3 files changed, 3 insertions(+), 61 deletions(-)

diff --git a/include/arrays.h b/include/arrays.h
index 90a2e08..fc8757a 100644
--- a/include/arrays.h
+++ b/include/arrays.h
@@ -1580,17 +1580,10 @@ namespace ecolab
       void asgV(size_t size, const E& x)
       {
         // copy into temporary data, as E may contain references to this
-//#ifdef __SYCL_DEVICE_ONLY__
-//        array tmp(size,m_allocator);
-//        asg_v(tmp.data(),size,x);
-//        resize(size, false);
-//        asg_v(data(),size,tmp);
-//#else
         array tmp(size,m_allocator);
         asg_v(tmp.data(),size,x);
         groupBarrier();
         swap(tmp);
-        //#endif
       }
       
       void copy() //any nonconst method needs to call this
@@ -1683,27 +1676,11 @@ namespace ecolab
       template <class V>
       void resizeAndInit(size_t s, const V& val) {resize(s,false); operator=(val);}
 
-      bool onDevice() const {
-#ifdef SYCL_LANGUAGE_VERSION
-      // true if created in a parallel section
-      return sycl::address_space_cast
-        <sycl::access::address_space::private_space,
-         sycl::access::decorated::yes>(this).get();
-#else
-      return false;
-#endif
-       }
-      
       void swap(array& x) {
 #ifndef __SYCL_DEVICE_ONLY__
         std::swap(dt, x.dt);
         std::swap(m_allocator,x.m_allocator);
 #else
-        // this code still doesn't work, even on the revised onDevice test above
-        if (onDevice() && x.onDevice()) {
-          std::swap(dt, x.dt);
-          std::swap(m_allocator,x.m_allocator);
-          } else {
           // assumption here is these array may be per thread, or
           // maybe shared by all threads in a group, hence std::swap as above won't work
           auto lhs=dt, rhs=x.dt;
@@ -1714,7 +1691,6 @@ namespace ecolab
           m_allocator=ralloc;
           x.m_allocator=lalloc;
           groupBarrier();
-        }
 #endif
       }
     
@@ -1733,7 +1709,7 @@ namespace ecolab
         if (x.dt==dt) return *this;
         if (m_allocator==x.m_allocator) {
           release();
-          if (groupLeader()||onDevice()) {
+          /*if (groupLeader()||onDevice())*/ {
             dt=x.dt;
           }
           incrRef();
diff --git a/include/sycl.h b/include/sycl.h
index 6b7140e..19773ad 100644
--- a/include/sycl.h
+++ b/include/sycl.h
@@ -150,13 +150,6 @@ namespace ecolab
     return 0;
 #endif
   }
-  
-  inline bool onDevice() {
-#ifdef __SYCL_DEVICE_ONLY__
-    return true;
-#endif
-    return false;
-  }
 }
 
 
diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc
index e89f764..69fb1c0 100644
--- a/models/ecolab_model.cc
+++ b/models/ecolab_model.cc
@@ -73,11 +73,7 @@ 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]);}
-  {
-    auto r=point.ROUND(expr[i]);
-    return r;
-  }
+  int operator[](size_t i) const {return point.ROUND(expr[i]);}
 };
 
 namespace ecolab::array_ns
@@ -132,7 +128,6 @@ void EcolabPoint::condense(const ModelData::BoolArray& mask, size_t mask_true)
     return;
   }
   LocalArray tmp(mask_true);
-  //UnsignedArray tmp(mask_true,density.allocator());
   if (groupLeader())
     for (size_t i=0, j=0; i<density.size(); ++i)
       if (mask[i] && j<mask_true) 
@@ -248,14 +243,12 @@ void SpatialModel::mutate()
     }
 
   // deallocate on device
-  //cout<<"deallocate"<<endl;
   groupedForAll([newSp=newSp.data()](EcolabCell& c,size_t i) {
     newSp[i].clear();
     assert(newSp[i].refCnt()==0);
   });
    
 
-  //cout<<"ModelData::mutate"<<endl;
 #ifdef MPI_SUPPORT
   MPIbuf b; b<<new_sp<<(*cell_ids); b.gather(0);
   if (myid()==0)
@@ -276,20 +269,13 @@ void SpatialModel::mutate()
 #endif
   if (new_sp.size()==0) return;
 
-  //cout<<"computeODiagIdx"<<endl;
   computeODiagIdx();
   mut_scale->clear();
   newSp.clear();
 
-  //cout<<"set 1"<<endl;
   // set the new species density to 1 for those created on this cell
   groupedForAll([cell_ids=&*cell_ids](EcolabCell& c,size_t) {
     c.density <<= (*cell_ids)==c.id;
-    //auto oldSize=c.density.size();
-    //groupBarrier();
-    //c.density.resize(c.density.size()+cell_ids->size());
-    //groupBarrier();
-    //array_ns::asg_v(c.density.data()+oldSize, cell_ids->size(), (*cell_ids)==c.id);
     assert(all(c.density>=0));
   });
 }
@@ -299,8 +285,7 @@ EcolabPoint::UnsignedArray EcolabPoint::mutate(const E& mut_scale)
 {
   /* calculate the number of mutants each species produces */
   if (density.size()==0) return {density.allocator()};
-  //#ifdef __SYCL_DEVICE_ONLY__
-#if 1
+#ifdef __SYCL_DEVICE_ONLY__
   LocalArray speciations=roundArray(mut_scale * density);
   //UnsignedArray speciations(roundArray(mut_scale * density), density.allocator());
   auto nsp=density.size();
@@ -309,17 +294,7 @@ EcolabPoint::UnsignedArray EcolabPoint::mutate(const E& mut_scale)
   //UnsignedArray offsets(nsp+1,density.allocator());
   unsigned* offs_p=offsets.data();
   const unsigned* sp_p=speciations.data();
-#ifdef __SYCL_DEVICE_ONLY__
-  //#if 0
   sycl::joint_exclusive_scan(syclGroup(),sp_p,sp_p+nsp,offs_p,sycl::plus<unsigned>());
-#else
-  if (groupLeader())
-    {
-      // sequential version for diagnostic purposes
-      offs_p[0]=0;
-      for (size_t i=0; i<nsp; ++i) offs_p[i+1]=offs_p[i]+sp_p[i];
-    }
-#endif
   groupBarrier();
   unsigned numSpeciations=0;
   if (localThreadId()==0) {
@@ -327,9 +302,7 @@ EcolabPoint::UnsignedArray EcolabPoint::mutate(const E& mut_scale)
     // group scope, because of COW semantics
     numSpeciations=offs_p[nsp]=offs_p[nsp-1]+sp_p[nsp-1];
   }
-#ifdef __SYCL_DEVICE_ONLY__
   numSpeciations=sycl::group_broadcast(syclGroup(),numSpeciations,0);
-#endif
   
   if (numSpeciations==0) return {density.allocator()};
   

From 9ea9bf978f4e46cf8d0e6bf2ddeaa23ec44b7397 Mon Sep 17 00:00:00 2001
From: Russell Standish <hpcoder@hpcoders.com.au>
Date: Sat, 8 Aug 2026 15:13:13 +1000
Subject: [PATCH 17/19] Fix up conventional build. Clean up device allocators.
 Add profiling output to spatial_ecolab

---
 include/DeviceAllocator.h | 131 ++++++++++----------------------------
 include/non-sycl.h        |   2 +
 models/ecolab_model.cc    |   2 +-
 models/ecolab_model.h     |   2 +
 models/spatial_ecolab.py  |  11 +++-
 5 files changed, 47 insertions(+), 101 deletions(-)

diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h
index ce67c1a..ed69f11 100644
--- a/include/DeviceAllocator.h
+++ b/include/DeviceAllocator.h
@@ -27,108 +27,53 @@ namespace ecolab
   };
   
   inline __attribute__((noinline)) bool& fatalErrorFlag() {
+#ifdef  SYCL_LANGUAGE_VERSION
     return sycl::ext::oneapi::group_local_memory<FatalErrorFlag>(syclGroup(),false)->flag;
+#else
+    static bool flag;
+    return flag;
+#endif
   }
   
   // Bounded MPMC circular buffer queue for SYCL using per-slot sequence numbers.
   // dequeue() returns ~0U when queue appears empty (non-blocking empty signal).
   template <unsigned size>
-  class Queue
+  class Stack
   {
     static_assert((size&(size-1))==0,"size must be power of two");
     constexpr static unsigned mask=size-1;
 
-    struct Slot
-    {
-      uint64_t seq;
-      unsigned value;
-    };
+    unsigned slots[size];
+    unsigned top=size; //empty stack, stack grows down
 
-    Slot slots[size];
-    uint64_t head=0, tail=size, pushing=0;
+    using Atomic=sycl::atomic_ref<unsigned,sycl::memory_order::acq_rel,sycl::memory_scope::device>;
     
-    template <class T> using Atomic=sycl::atomic_ref<T,sycl::memory_order::acq_rel,sycl::memory_scope::device>;
-
   public:
     void init() {
-      tail=0;
+      top=0; // full stack
       for (unsigned i=syclItem().get_global_linear_id(); i<size;
-           i+=syclItem().get_global_range().size()) {
-          slots[i].value=i;
-          slots[i].seq=i+1;
-        }
+           i+=syclItem().get_global_range().size()) 
+          slots[i]=i;
     }
     
-    void enqueue(unsigned x)
+    void push(unsigned x)
     {
-      Atomic<uint64_t> t(tail), p(pushing);
-      slots[--t].value=x;
-//      for (bool published=false; !published;)
-//      {
-//        if (localThreadId()==0)
-//          {
-//            Atomic<uint64_t> headAtomic(head);
-//            auto pos=headAtomic.load();
-//            Slot& slot=slots[pos & mask];
-//            Atomic<uint64_t> seqAtomic(slot.seq);
-//            auto seq=seqAtomic.load(sycl::memory_order::acquire);
-//            auto diff=int64_t(seq)-int64_t(pos);
-//
-//            if (diff==0 && headAtomic.compare_exchange_weak(pos,pos+1))
-//              {
-//                Atomic<unsigned> valueAtomic(slot.value);
-//                valueAtomic=x;
-//                seqAtomic=pos+1;
-//                published=true;
-//              }
-//          }
-//#ifdef __SYCL_DEVICE_ONLY__
-//        sycl::atomic_fence(sycl::memory_order::seq_cst, sycl::memory_scope::device);
-//        published=sycl::group_broadcast(syclGroup(),published,0);
-//#endif
-//      }
+      slots[--Atomic(top)]=x;
     }
 
-    unsigned dequeue()
+    unsigned pop()
     {
-      Atomic<uint64_t> t(tail);
-      // updating tail in a cas loop avoids the race condition
-      // between the test and increment
-      uint64_t p=t++;
+      Atomic t(top);
+      unsigned p=t++;
       if (p>=size) {t=size; return ~0;} // stack empty
-      return slots[p].value;
-//      unsigned v=~0U-1;
-//      unsigned niter=0;
-//      while (v==~0U-1)
-//      {
-//        if (localThreadId()==0)
-//          {
-//            Atomic<uint64_t> tailAtomic(tail);
-//            auto pos=tailAtomic.load();
-//            Slot& slot=slots[pos & mask];
-//            Atomic<uint64_t> seqAtomic(slot.seq);
-//            auto seq=seqAtomic.load(sycl::memory_order::acquire);
-//            auto diff=int64_t(seq)-int64_t(pos+1);
-//
-//            if (diff==0 && tailAtomic.compare_exchange_weak(pos,pos+1))
-//              {
-//                Atomic<unsigned> valueAtomic(slot.value);
-//                v=valueAtomic;
-//                Atomic<uint64_t> release(slot.seq);
-//                release.store(pos+size,sycl::memory_order::release);
-//              }
-//            if (diff<0) v=~0U; // empty queue
-//          }
-//#ifdef __SYCL_DEVICE_ONLY__
-//        sycl::atomic_fence(sycl::memory_order::seq_cst, sycl::memory_scope::device);
-//        v=sycl::group_broadcast(syclGroup(),v,0);
-//#endif
-//        if (niter++>10) {
-//          break;
-//          if (groupLeader()) printf("enqueue failed\n",0);
-//        }
-//      }
-//      return v;
+      return slots[p];
+    }
+
+    // move contents of \a x onto this. Not threadsafe, call from host
+    void appendAndDiscard(Stack& x) {
+      top-=size-x.top;
+      memcpy(slots+top, x.slots+x.top, (size-x.top)*sizeof(slots[0]));
+      x.top=size;
     }
   };
 
@@ -155,8 +100,8 @@ namespace ecolab
   {
     constexpr static unsigned pageSize=1<<order;
     constexpr static unsigned numPages=poolSize/pageSize;
-    Queue<numPages> queue;
-    Queue<numPages> discard;
+    Stack<numPages> queue;
+    Stack<numPages> discard; // discard pile
     char memory[poolSize];
     DeviceAllocator<order+2> nextAllocator; // next size up allocator
   public:
@@ -168,8 +113,7 @@ namespace ecolab
       nextAllocator.init();
     }
     void recycleDiscardPile() {
-      for (auto v=discard.dequeue(); v!=~0U; v=discard.dequeue())
-        queue.enqueue(v);
+      queue.appendAndDiscard(discard);
       nextAllocator.recycleDiscardPile();
     }
     // all members of group get the same pointer
@@ -177,15 +121,10 @@ namespace ecolab
       if (size==0) return nullptr;
       if (size<=pageSize) {
         unsigned offs=~0U;
-        if (localThreadId()==0) offs=queue.dequeue();
+        if (localThreadId()==0) offs=queue.pop();
 #ifdef __SYCL_DEVICE_ONLY__
         offs=sycl::group_broadcast(syclGroup(),offs,0);
 #endif
-       
-//#ifdef __SYCL_DEVICE_ONLY__
-//        if (groupLeader()) printf("alloc pageSize=%u offs=%u on group %u\n",
-//                                  pageSize,offs,syclGroup().get_group_linear_id());
-//#endif
         if (offs!=~0U)
           return memory+(offs<<order);
       }
@@ -196,9 +135,12 @@ namespace ecolab
       if (p>=memory && p<memory+poolSize) {
 #ifdef __SYCL_DEVICE_ONLY__
         if (groupLeader())
-          discard.enqueue((reinterpret_cast<char*>(p)-memory)>>order);
+          // push onto discard pile to avoid race condition
+          discard.push((reinterpret_cast<char*>(p)-memory)>>order);
 #else
-        queue.enqueue((reinterpret_cast<char*>(p)-memory)>>order);
+        // on host, we can push back onto stack. Note this is not
+        // threadsafe, so not to be used with OpenMP.
+        queue.push((reinterpret_cast<char*>(p)-memory)>>order);
 #endif
         return;
       }
@@ -209,15 +151,10 @@ namespace ecolab
   };
   
   inline DeviceAllocator<>& deviceAllocator() {
-//#ifdef __SYCL_DEVICE_ONLY__
-//    printf("deviceAllocator() illegally called on device\n",0);
-//    return *reinterpret_cast<DeviceAllocator<>*>(0);
-//#else
     static DeviceType<DeviceAllocator<>> deviceAllocator;
     static int dummy=
       (deviceAllocator->init(), syclQ().wait_and_throw(), 0);      
     return *deviceAllocator;
-    //#endif
   }
 
   /// Allocator wrapping the DeviceAllocator singleton
diff --git a/include/non-sycl.h b/include/non-sycl.h
index de48dff..d41b7c4 100644
--- a/include/non-sycl.h
+++ b/include/non-sycl.h
@@ -29,6 +29,8 @@ namespace ecolab
   
   template <class E, ecolab::USMAlloc UA>
   struct SyclRandomEngine: public E {};
+  
+  template <class T> using LocalAllocator=std::allocator<T>;
 
 }
 #endif
diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc
index 69fb1c0..1502cd4 100644
--- a/models/ecolab_model.cc
+++ b/models/ecolab_model.cc
@@ -88,7 +88,7 @@ void EcolabPoint::generate(unsigned niter, const ModelData& model)
   //array<int,LocalAllocator<int>> lDensity(density), tmp(density.size());
   //array<int,Allocator<int>> lDensity(density), tmp(density.size(), density.allocator());
   auto& lDensity=density;
-  array<int,GlobalDeviceAllocator<int>> tmp(density.size(), density.allocator());
+  array<int,Allocator<int>> tmp(density.size(), density.allocator());
   
   for (unsigned step=0; step<niter; step++)
     {
diff --git a/models/ecolab_model.h b/models/ecolab_model.h
index bc94f0b..ace8b23 100644
--- a/models/ecolab_model.h
+++ b/models/ecolab_model.h
@@ -14,7 +14,9 @@ using classdesc::Object;
 #include <random>
 #include <vector>
 #include <pack_stl.h>
+#ifdef SYCL_LANGUAGE_VERSION
 #include "DeviceAllocator.h"
+#endif
 
 #ifdef USE_FLOAT
 using Float=float;
diff --git a/models/spatial_ecolab.py b/models/spatial_ecolab.py
index f75b6e3..4720bb0 100644
--- a/models/spatial_ecolab.py
+++ b/models/spatial_ecolab.py
@@ -56,16 +56,20 @@ def randomList(num, min, max):
                   
 from plot import plot
 from GUI import gui, statusBar, windows
+from time import perf_counter
 
 epoch=2000000
 mut_factor=1000
 
 extinctions=0
 migrations=0
+
+out=open("time.dat","w")
+
 def stepImpl():
-    #print("b4 generate")
+    start=perf_counter()
+    nsp=len(ecolab.species)
     ecolab.generate(100)
-    #print("b4 mutate")
     ecolab.mutate()
 
     epochTs=ecolab.tstep()%epoch
@@ -76,8 +80,9 @@ def stepImpl():
 
     global extinctions, migrations
     #migrations+=ecolab.migrate()
-    #print("b4 condense")
     extinctions+=ecolab.condense()
+    timePerStep=perf_counter()-start
+    print(nsp,timePerStep,file=out,flush=True)
     #print(ecolab.nsp()())
     #ecolab.syncThreads()
     #ecolab.gather()

From 288c3281834a86315d31a9316a7859b2a017135c Mon Sep 17 00:00:00 2001
From: Russell Standish <hpcoder@hpcoders.com.au>
Date: Mon, 10 Aug 2026 15:00:57 +1000
Subject: [PATCH 18/19] optimisation:

---
 classdesc                 |  2 +-
 include/DeviceAllocator.h | 20 +++++++++++++-------
 include/ecolab.h          |  2 ++
 models/ecolab_model.cc    | 20 +++++++++-----------
 models/ecolab_model.h     |  2 +-
 models/spatial_ecolab.py  |  2 +-
 6 files changed, 27 insertions(+), 21 deletions(-)

diff --git a/classdesc b/classdesc
index 025ca7a..3cab790 160000
--- a/classdesc
+++ b/classdesc
@@ -1 +1 @@
-Subproject commit 025ca7aac67e2676698f5ef8d7f7c92f210617ca
+Subproject commit 3cab7902a007ec006bcb9d1b470447c87b04e81e
diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h
index ed69f11..0d21e69 100644
--- a/include/DeviceAllocator.h
+++ b/include/DeviceAllocator.h
@@ -47,7 +47,7 @@ namespace ecolab
     unsigned top=size; //empty stack, stack grows down
 
     using Atomic=sycl::atomic_ref<unsigned,sycl::memory_order::acq_rel,sycl::memory_scope::device>;
-    
+    CLASSDESC_ACCESS(Stack);
   public:
     void init() {
       top=0; // full stack
@@ -79,7 +79,7 @@ namespace ecolab
 
   template <unsigned order> class DeviceAllocator;
   /// empty allocator to terminate template recursion
-  template <> class DeviceAllocator<maxOrder> {
+  template <> class DeviceAllocator<ecolab::maxOrder> {
   public:
     void* allocate(size_t sz) {
       if (groupLeader())
@@ -104,6 +104,7 @@ namespace ecolab
     Stack<numPages> discard; // discard pile
     char memory[poolSize];
     DeviceAllocator<order+2> nextAllocator; // next size up allocator
+    CLASSDESC_ACCESS(DeviceAllocator);
   public:
     void init() {
       auto chunkOWork=syclQ().get_device().
@@ -170,7 +171,7 @@ namespace ecolab
     DeviceAllocator<>* allocator;
     
 #ifdef __SYCL_DEVICE_ONLY__
-    GlobalDeviceAllocator() = delete;
+    GlobalDeviceAllocator(): allocator(nullptr) {} // = delete;
 #else
     GlobalDeviceAllocator() // note: default constructor must be called on host
     {allocator=&deviceAllocator();}
@@ -216,10 +217,11 @@ namespace ecolab
   /**
      A Local allocator allocates memory from device local memory,
      which is shared between threads of a work group, and has the same
-     lifetime as the kernel
+     lifetime as the kernel.
+     LocalAllocatorT so we can expose LocalAllocator as a template alias on both host and device branches
   */
   template <class T>
-  class LocalAllocator
+  class LocalAllocatorT
   {
   public:
     using value_type=T;
@@ -244,13 +246,17 @@ namespace ecolab
       return reinterpret_cast<T*>(alloc);
     }
     void deallocate(T*p,size_t) {} // cleaned up when group exits
-    template<class U> struct rebind {using other=LocalAllocator<U>;};
+    template<class U> struct rebind {using other=LocalAllocatorT<U>;};
     // allocator is stateless
-    bool operator==(const LocalAllocator&) const {return true;}
+    bool operator==(const LocalAllocatorT&) const {return true;}
   };
+  template <class T> using LocalAllocator=LocalAllocatorT<T>;
 #else
+  template <class T> class LocalAllocatorT {};
   template <class T> using LocalAllocator=std::allocator<T>;
 #endif
    
 }
+
+#include "DeviceAllocator.cd"
 #endif
diff --git a/include/ecolab.h b/include/ecolab.h
index c2bf93e..d5186a9 100644
--- a/include/ecolab.h
+++ b/include/ecolab.h
@@ -152,6 +152,8 @@ namespace ecolab
       hostForAll(f);
 #endif
     }
+
+    // TODO - we need a const version as well. Implement via a free template function for the this pointer, and ensure constness is passed through to the arguments of the functional
     
     /// apply a functional to all local cells of this processor in
     /// parallel, where each cell is allocated SIMD parallel computer
diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc
index 1502cd4..d1bb426 100644
--- a/models/ecolab_model.cc
+++ b/models/ecolab_model.cc
@@ -85,20 +85,16 @@ RoundArray<E,EcolabPoint> EcolabPoint::roundArray(const E& expr)
 
 void EcolabPoint::generate(unsigned niter, const ModelData& model)
 {
-  //array<int,LocalAllocator<int>> lDensity(density), tmp(density.size());
-  //array<int,Allocator<int>> lDensity(density), tmp(density.size(), density.allocator());
-  auto& lDensity=density;
-  array<int,Allocator<int>> tmp(density.size(), density.allocator());
+  array<int,LocalAllocator<int>> lDensity(density), tmp(density.size());
+  //auto& lDensity=density;
+  //array<int,Allocator<int>> tmp(density.size(), density.allocator());
   
   for (unsigned step=0; step<niter; step++)
     {
       array_ns::map(lDensity.size(),  [&](size_t i){
         Float ir=model.interaction.diag[i]*lDensity[i];
-        //for (auto j: model.oDiagIdx[i])
-        for (size_t k=0; k<model.oDiagIdx[i].size(); ++k) {
-          unsigned j=model.oDiagIdx[i][k];
+        for (auto j: model.oDiagIdx[i])
           ir+=model.interaction.val[j]*lDensity[model.interaction.col[j]];
-        }
         tmp[i]=ROUND(lDensity[i] + lDensity[i] * (model.repro_rate[i] + ir));
       });
       groupBarrier(); // synchronises threads on each iteration
@@ -114,10 +110,12 @@ void EcolabPoint::generate(unsigned niter, const ModelData& model)
 unsigned EcolabPoint::nsp() const
 {return sum(density!=0);}
 
-array<unsigned> SpatialModel::nsp() const
+array<unsigned> SpatialModel::nsp()
 {
-  array<unsigned> nsp;
-  for (auto& i: objects) nsp<<=i->nsp();
+  EcolabPoint::UnsignedArray nsp(size());
+  groupedForAll([nsp=nsp.data()](const EcolabCell& c,size_t i) {
+    nsp[i]=c.nsp();
+  });
   return nsp;
 }
 
diff --git a/models/ecolab_model.h b/models/ecolab_model.h
index ace8b23..0e9716e 100644
--- a/models/ecolab_model.h
+++ b/models/ecolab_model.h
@@ -142,7 +142,7 @@ class SpatialModel: public ModelData, public EcolabGraph<EcolabCell>,
   EcolabCell& cell(size_t x, size_t y) {
     return *objects[makeId(x,y)];
   }
-  array<unsigned> nsp() const;
+  array<unsigned> nsp();
   void makeConsistent();
   void seed(unsigned x) {groupedForAll([=](EcolabCell& cell,size_t){cell.rand.seed(x);});}
   void generate(unsigned niter);
diff --git a/models/spatial_ecolab.py b/models/spatial_ecolab.py
index 4720bb0..661436d 100644
--- a/models/spatial_ecolab.py
+++ b/models/spatial_ecolab.py
@@ -69,7 +69,7 @@ def randomList(num, min, max):
 def stepImpl():
     start=perf_counter()
     nsp=len(ecolab.species)
-    ecolab.generate(100)
+    ecolab.generate(1000)
     ecolab.mutate()
 
     epochTs=ecolab.tstep()%epoch

From ba794750d719bc503c1865af3258a9d0931f36b3 Mon Sep 17 00:00:00 2001
From: Russell Standish <hpcoder@hpcoders.com.au>
Date: Tue, 11 Aug 2026 09:53:09 +1000
Subject: [PATCH 19/19] chore: address code review comments

---
 include/arrays.h | 1 +
 include/sycl.h   | 1 -
 2 files changed, 1 insertion(+), 1 deletion(-)

diff --git a/include/arrays.h b/include/arrays.h
index fc8757a..166956d 100644
--- a/include/arrays.h
+++ b/include/arrays.h
@@ -1594,6 +1594,7 @@ namespace ecolab
             decrRef();
             auto sz=size();
             dt=alloc(sz);
+            if (!dt) return;
 #ifdef __SYCL_DEVICE_ONLY__
             asg_v(dt->dt,sz,oldData->dt);
 #else
diff --git a/include/sycl.h b/include/sycl.h
index 19773ad..5e09ffa 100644
--- a/include/sycl.h
+++ b/include/sycl.h
@@ -138,7 +138,6 @@ namespace ecolab
   inline bool groupLeader() {
 #ifdef __SYCL_DEVICE_ONLY__
     return syclGroup().leader();
-    //return syclGroup().get_local_linear_id()==0;
 #endif
     return true;
   }