diff --git a/graphcode b/graphcode
index a83dbb1..2b5a198 160000
--- a/graphcode
+++ b/graphcode
@@ -1 +1 @@
-Subproject commit a83dbb1d3a0ad497c0f405a0c74c4aef4cf7b96d
+Subproject commit 2b5a1980cf5162ea980c12bd1911e31850d407fa
diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h
new file mode 100644
index 0000000..b4402a5
--- /dev/null
+++ b/include/DeviceAllocator.h
@@ -0,0 +1,252 @@
+/*
+  @copyright Russell Standish 2026
+  @author Russell Standish
+  This file is part of EcoLab
+
+  Open source licensed under the MIT license. See LICENSE for details.
+*/
+
+#ifndef DEVICE_ALLOCATOR_H
+#define DEVICE_ALLOCATOR_H
+#include "sycl.h"
+#include "graphcode.h"
+
+namespace ecolab
+{
+  /// minimum allocation size = 2^minOrder, maximum object size = 2^(maxOrder-1)
+  constexpr unsigned minOrder=4, maxOrder=20;
+  /// memory allocated to each order, total memory allocated on device=poolSize*(maxOrder-minOrder)
+  constexpr unsigned poolSize=32*1024*1024;
+
+  struct FatalErrorFlag
+  {
+    bool flag;
+  };
+  
+  inline __attribute__((noinline)) bool& fatalErrorFlag() {
+    return sycl::ext::oneapi::group_local_memory<FatalErrorFlag>(syclGroup(),false)->flag;
+  }
+  
+  // 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
+  {
+    static_assert((size&(size-1))==0,"size must be power of two");
+    constexpr static unsigned mask=size-1;
+
+    struct Slot
+    {
+      unsigned seq;
+      unsigned value;
+    };
+
+    Slot slots[size];
+    unsigned head=size, tail=0;
+
+    using Atomic=sycl::atomic_ref<unsigned,sycl::memory_order::relaxed,sycl::memory_scope::device>;
+
+  public:
+    Queue()
+    {
+      for (unsigned i=0; i<size; ++i) {
+        slots[i].value=i;
+        slots[i].seq=i+1;
+      }
+    }
+
+    void enqueue(unsigned 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)
+        {
+          if (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()
+    {
+      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)
+        {
+          if (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;
+          }
+        }
+        else if (diff<0)
+        {
+          return ~0U; // signal buffer empty, don't wait
+        }
+      }
+    }
+  };
+
+  template <unsigned order> class DeviceAllocator;
+  /// empty allocator to terminate template recursion
+  template <> class DeviceAllocator<maxOrder> {
+  public:
+    void* allocate(size_t sz) {
+      if (groupLeader())
+        sycl::ext::oneapi::experimental::printf("failed to allocate %zu bytes\n",sz);
+#ifdef __SYCL_DEVICE_ONLY__
+      fatalErrorFlag()=true;
+#else
+      throw std::bad_alloc();
+#endif
+      return nullptr;
+    }
+    void deallocate(void* p, size_t) {sycl::ext::oneapi::experimental::printf("%p leaked on device\n",p);}
+  };
+
+  template <unsigned order=minOrder> class DeviceAllocator
+  {
+    constexpr static unsigned pageSize=1<<order;
+    constexpr static unsigned numPages=poolSize/pageSize;
+    Queue<numPages> queue;
+    char memory[poolSize];
+    DeviceAllocator<order+1> nextAllocator; // next size up allocator
+  public:
+    // all members of group get the same pointer
+    void* allocate(size_t size) {
+      if (size==0) return nullptr;
+      if (size<=pageSize) {
+        unsigned offs;
+        if (groupLeader()) offs=queue.dequeue();
+#ifdef __SYCL_DEVICE_ONLY__
+        offs=sycl::group_broadcast(syclGroup(),offs);
+#endif
+        if (offs!=~0U)
+          return memory+(offs<<order);
+      }
+      return nextAllocator.allocate(size);
+    }
+    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);
+        return;
+      }
+      nextAllocator.deallocate(p,size);
+    }
+    
+    bool inAllocator(void* p) const {return p>=this && p<this+1;}
+  };
+  
+  inline DeviceAllocator<>& deviceAllocator() {
+    static DeviceType<DeviceAllocator<>> deviceAllocator;
+    return *deviceAllocator;
+  }
+
+  /// Allocator wrapping the DeviceAllocator singleton
+  template <class T>
+  struct GlobalDeviceAllocator
+  {
+    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;
+
+    DeviceAllocator<>* allocator=&deviceAllocator();
+    
+    GlobalDeviceAllocator() = default; // note: default constructor must be called on host
+    template <class U>
+    GlobalDeviceAllocator(const GlobalDeviceAllocator<U>& other):
+      allocator(other.allocator) {}
+    
+    T* allocate(size_t n)
+    {return reinterpret_cast<T*>(allocator->allocate(n*sizeof(T)));}
+    void deallocate(T* p, size_t n){allocator->deallocate(p,n*sizeof(T));}
+    template<class U> struct rebind {using other=GlobalDeviceAllocator<U>;};
+  };
+
+  template <class T>
+  struct HostSharedAllocator: public graphcode::Allocator<T>
+  {
+    HostSharedAllocator(): graphcode::Allocator<T>(syclQ(), sycl::usm::alloc::shared) {}
+    template<class U> struct rebind {using other=HostSharedAllocator<U>;};
+  };
+
+  constexpr static unsigned LocalAllocatorSize=30*1024; // 32KiB = half typical local storage
+
+  struct LocalAllocatorBuffer
+  {
+    unsigned next;
+    char buffer[LocalAllocatorSize];
+  };
+
+  /*
+    group_local_memory is a weird beast. The address returned is tied to the line of code in which it instantiated, so we need to specify noinline to prevent it from being inlined, and inline to ensure single definition
+  */
+  inline __attribute__((noinline)) LocalAllocatorBuffer& localAllocatorBuffer()
+   {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*,size_t) {} // cleaned up when group exits
+    template<class U> struct rebind {using other=LocalAllocator<U>;};
+  };
+#endif
+   
+}
+#endif
diff --git a/include/arrays.h b/include/arrays.h
index 3dd2ac9..8ad4bc3 100644
--- a/include/arrays.h
+++ b/include/arrays.h
@@ -1477,26 +1477,25 @@ namespace ecolab
       
       friend class WhereContext;
 
-      ///allocate \a n variables of type \a T 
+      ///allocate \a n variables of type \a T
+#ifdef __SYCL_DEVICE_ONLY__
+      __attribute__((always_inline))
+#endif
       array_data<T> *alloc(std::size_t n)
       {
+        if (n==0) return nullptr;
         T *p; 
         array_data<T> *r;
-        //p = (char*)std::malloc((n-array_data<T>::debug_display) * sizeof(T) + sizeof(array_data<T>) + 16);
         // over allocate to allow for alignment and metadata
         auto allocation=n + (sizeof(array_data<T>) + 16)/sizeof(T)+1-array_data<T>::debug_display;
         p = m_allocator.allocate(allocation);
       
         if (!p)
           {
-#ifdef __SYCL_DEVICE_ONLY__
-            printf("failed to allocate %d bytes in array\n",sizeof(T)*n);
-#endif
+            if (groupLeader())
+              printf("failed to allocate %zu bytes in array\n",sizeof(T)*n);
             return nullptr; // SYCL allocator returns nullptr if not initialised
           }
-// #ifdef __SYCL_DEVICE_ONLY__
-//            syclPrintf("succeeded in allocating %d bytes in array\n",sizeof(T)*n);
-//#endif
        
 #ifdef __ICC
         // we need to align data onto 16 byte boundaries
@@ -1517,47 +1516,49 @@ namespace ecolab
       void free(array_data<T> *p)
       {
         assert(p);
-        //std::free(p->allocated_pointer);
         m_allocator.deallocate(p->allocated_pointer,p->allocation);
       }
 
       void set_size(size_t s) {dt = alloc(s);}
 
-#if defined(SYCL_LANGUAGE_VERSION)
-      using AtomicUnsignedRef=sycl::atomic_ref
-        <unsigned,sycl::memory_order::relaxed,sycl::memory_scope::system>;
-#else
-      using AtomicUnsignedRef=unsigned&;
-#endif
-      
-      AtomicUnsignedRef ref() // access reference counter
+      unsigned& ref() // access reference counter
       {
         assert(dt);
-        return AtomicUnsignedRef(dt->cnt);
+        return dt->cnt;
       }
 
+      // increment reference counter
+      void incrRef()
+      {
+        if (dt)
+          {
+            groupBarrier();
+            if (groupLeader()) ref()++;
+            groupBarrier();
+          }
+      }
+
+      // decrement reference counter
+      void decrRef()
+      {
+        if (dt)
+          {
+            groupBarrier();
+            if (groupLeader()) ref()--;
+            groupBarrier();
+          }
+       }
+      
       void release()
       {
         if (dt)
           {
-//#if defined(SYCL_LANGUAGE_VERSION) && !defined(__SYCL_DEVICE_ONLY__)
-//            // dt pointer may be allocated on device, or in device
-//            // memory, and we may be running on the host, in which
-//            // case just leak the memory, otherwise we'll have a crash
-//            // TODO - call release on device in a single_task for the
-//            // first situation
-//            // TODO in second situation, update ref
-//            // count in a single_task, and pass back value of
-//            // allocated pointer for deallocation on host
-//            if (is_same<A,typename CellBase::CellAllocator<T>>::value ||
-//                sycl::get_pointer_type(dt,syclQ().get_context())==sycl::usm::alloc::device) return;
-//#endif
             if (ref()==1)
               {
                 free(dt);
+                return;
               }
-            else
-              ref()--;
+            decrRef();
           }
       }
 
@@ -1565,15 +1566,16 @@ namespace ecolab
 
       // implements copying data from x to this in a GPU friendly way
       template <class E>
-      void asgV(const A& alloc, size_t size, const E& x)
+      void asgV(size_t size, const E& x)
       {
+        // copy into temporary data, as E may contain references to this
 #ifdef __SYCL_DEVICE_ONLY__
-        GroupLocal<array> tmp(size,alloc);
-        array_ns::asg_v(tmp->dt->dt,size,x);
-        groupBarrier();
-        if (syclGroup().leader()) swap(*tmp);
+        array<T,LocalAllocator<T>> tmp(size);
+        asg_v(tmp.data(),size,x);
+        resize(size);
+        asg_v(data(),size,tmp);
 #else
-        array tmp(size,alloc);
+        array tmp(size);
         asg_v(tmp.data(),size,x);
         swap(tmp);
 #endif
@@ -1583,17 +1585,16 @@ namespace ecolab
       {           // to implement copy-on-write semantics
         if (dt && ref()>1)
           {
-#ifdef __SYCL_DEVICE_ONLY__
-            printf("b4 asgV in copy\n");
-            asgV(m_allocator, size(), dt->dt);
-            printf("after asgV in copy\n");
-#else
             array_data<T>* oldData=dt;
-            bool freeMem = ref()-- == 0;
+            decrRef();
+            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));
-            if (freeMem) free(oldData);
 #endif
+            if (freeMem) free(oldData);
           }
       }
 
@@ -1616,11 +1617,9 @@ namespace ecolab
 
       array(const array& x): m_allocator(x.m_allocator) 
       {
-//#ifdef __SYCL_DEVICE_ONLY__
-//        syclPrintf("creating array on group %u, thread  %u with dt=%x\n",syclGroup().get_group_linear_id(), syclGroup().get_local_linear_id(),x.dt);
-//#endif
         dt=x.dt;
-        if (dt) ref()++;
+        incrRef();
+          
       }
 
       template <class expr>
@@ -1628,7 +1627,8 @@ namespace ecolab
             typename enable_if< is_expression<expr>, void*>::T dummy=0): m_allocator(alloc)
       {
         set_size(e.size());
-        operator=(e);
+        asg_v(data(),e.size(),e);
+        //operator=(e);
       }
 
       ~array() {release();}
@@ -1636,12 +1636,14 @@ namespace ecolab
       const Allocator& allocator() const {return m_allocator;}
       const Allocator& allocator(const Allocator& alloc) {
         if (alloc==m_allocator) return m_allocator;
-        asgV(alloc, size(), data());
+        array tmp(size(),alloc);
+        tmp.asgV(size(), data());
+        swap(tmp);
         return m_allocator;
       }
 
       /// current value of the reference counter
-      unsigned refCnt() const {return dt->cnt;}
+      unsigned refCnt() const {return dt? dt->cnt: 0;}
       
       /// resize array to \a s elements
       void resize(size_t s) {
@@ -1680,7 +1682,7 @@ namespace ecolab
         release();
         m_allocator=x.m_allocator;
         dt=x.dt;
-        if (dt) ref()++;
+        incrRef();
         return *this;
       }
 
@@ -1688,7 +1690,7 @@ namespace ecolab
       enable_if<is_expression<expr>, array&>::T 
       operator=(const expr& x) {
         if ((void*)(&x)==(void*)(this)) return *this;
-        asgV(m_allocator, x.size(), x);
+        asgV(x.size(), x);
         return *this;
       }
       template <class expr> typename
@@ -2399,8 +2401,8 @@ namespace ecolab
     return o;
   }
 
-  template <class T>
-  std::ostream& operator<<(std::ostream& o, const array<T>& x)
+  template <class T, class A>
+  std::ostream& operator<<(std::ostream& o, const array<T,A>& x)
   {return put(o,x);}
 
   /// ostream putter
diff --git a/include/ecolab.h b/include/ecolab.h
index 631c8bf..07e4b3a 100644
--- a/include/ecolab.h
+++ b/include/ecolab.h
@@ -12,16 +12,9 @@
 #ifndef ECOLAB_H
 #define ECOLAB_H
 
+#include "sycl.h"
 #ifdef SYCL_LANGUAGE_VERSION
 #include "DeviceAllocator.h"
-#ifdef __INTEL_LLVM_COMPILER
-template <typename... Args>
-inline sycl::nd_item<1> syclItem() {return sycl::ext::oneapi::this_work_item::get_nd_item<1>();}
-inline sycl::group<1> syclGroup() {return sycl::ext::oneapi::this_work_item::get_work_group<1>();}
-inline sycl::sub_group syclSubGroup() {return sycl::ext::oneapi::this_work_item::get_sub_group();}
-#else
-#error "EcoLab requires OneAPI compiler for some experimental functions"
-#endif
 #endif
 
 #include <stdlib.h>
@@ -42,7 +35,7 @@ namespace classdesc
   class pack_t;
 //class unpack_t;
   typedef pack_t unpack_t;
-  class TCL_obj_t;
+  //  class TCL_obj_t;
 }
 
 #define THROW_PTR_EXCEPTION  //Allows more generously for types containing pointers
@@ -73,10 +66,6 @@ namespace ecolab
 {
   using namespace classdesc;
   
-#if defined(__INTEL_LLVM_COMPILER) && defined(SYCL_LANGUAGE_VERSION)
-  using sycl::ext::oneapi::experimental::printf;
-#endif
-  
   /* these are defined to default values, even if MPI is false */
   /// MPI process ID and number of processes
   unsigned myid();
@@ -108,190 +97,14 @@ namespace ecolab
     ~OnExit() {f();}
   };
   
-#ifdef SYCL_LANGUAGE_VERSION
-  extern bool syclQDestroyed;
-  sycl::queue& syclQ();
-  void* reallocSycl(void*,size_t);
-#endif
-
-  inline void groupBarrier()
-  {
-#ifdef __SYCL_DEVICE_ONLY__
-    sycl::group_barrier(syclGroup());
-#endif
-  }
-  
-  template <class T>
-  struct SyclType: public T
-  {
-    template <class... A> SyclType(A... args): T(std::forward<A>(args)...) {}
-#ifdef SYCL_LANGUAGE_VERSION
-    void* operator new(size_t s) {return reallocSycl(nullptr,s);}
-    void operator delete(void* p) {reallocSycl(p,0);}
-    void* operator new[](size_t s) {return reallocSycl(nullptr,s);}
-    void operator delete[](void* p) {reallocSycl(p,0);}
-#endif
-  };
-
-  template <>
-  struct SyclType<size_t>
-  {
-    size_t data;
-    operator size_t() const {return data;}
-    operator size_t&() {return data;}
-    size_t operator=(size_t x) {return data=x;}
-#ifdef SYCL_LANGUAGE_VERSION
-    void* operator new(size_t s) {return reallocSycl(nullptr,s);}
-    void operator delete(void* p) {reallocSycl(p,0);}
-    void* operator new[](size_t s) {return reallocSycl(nullptr,s);}
-    void operator delete[](void* p) {reallocSycl(p,0);}
-#endif
-  };
-
-  template <class T>
-  struct SyclType<T*>
-  {
-    T* data;
-    operator T*() const {return data;}
-    operator T*&() {return data;}
-    T* operator=(T* x) {return data=x;}
-#ifdef SYCL_LANGUAGE_VERSION
-    void* operator new(size_t s) {return reallocSycl(nullptr,s);}
-    void operator delete(void* p) {reallocSycl(p,0);}
-    void* operator new[](size_t s) {return reallocSycl(nullptr,s);}
-    void operator delete[](void* p) {reallocSycl(p,0);}
-#endif
-  };
-
-
-  
-  template <class M>
-  class DeviceType
-  {
-    SyclType<M>* const model;
-  public:
-    using element_type=M;
-    template <class... A>
-    DeviceType(A... args): model(new SyclType<M>(std::forward<A>(args)...)) {}
-    DeviceType(const DeviceType& x): model(new SyclType<M>(*x.model)) {}
-    DeviceType& operator=(const DeviceType& x) {*model=*x.model; return *this;}
-#if defined(__GNUC__) && !defined(__clang__)
-#pragma GCC diagnostic push
-    // gcc incorrectly infers SyclType is polymorphic, which is quite plainly is not
-#pragma GCC diagnostic ignored "-Wdelete-non-virtual-dtor"
-#endif
-    ~DeviceType() {delete model;}
-#if defined(__GNUC__) && !defined(__clang__)
-#pragma GCC diagnostic pop
-#endif
-    M& operator*() {return *model;}
-    const M& operator*() const {return *model;}
-    M* operator->() {return model;}
-    const M* operator->() const {return model;}
-    operator bool() const {return true;} // always defined
-  };
-
-  template <class T, ecolab::USMAlloc UA>
-  struct SyclQAllocator: public graphcode::Allocator<T>
-  {
-#ifdef SYCL_LANGUAGE_VERSION
-    SyclQAllocator(): graphcode::Allocator<T>(syclQ(), UA) {}
-#endif
-    template<class U> struct rebind {using other=SyclQAllocator;};
-  };
-
   class GlobalDeallocateKernelTag;
   struct CellBase
   {
     size_t m_idx=0; // stash the position within the local node vector here
     size_t idx() const {return m_idx;}
-#ifdef SYCL_LANGUAGE_VERSION
-    using MemAllocator=DeviceAllocator<>;
-    MemAllocator* memAlloc=nullptr;
-    template <class T> class CellAllocator
-    {
-    public:
-      MemAllocator* memAlloc=nullptr;
-      template <class U> friend class Allocator;
-      CellAllocator()=default;
-      CellAllocator(MemAllocator* memAlloc): memAlloc(memAlloc) {}
-      template <class U> CellAllocator(const CellAllocator<U>& x): memAlloc(x.memAlloc) {}
-      T* allocate(size_t sz) {
-#ifdef  __SYCL_DEVICE_ONLY__
-        if (memAlloc)  {
-          auto r=reinterpret_cast<T*>(memAlloc->allocate(sz*sizeof(T)));
-          if (!r) printf("Mem allocation failed\n");
-          return r;
-        }
-        printf("Missing allocator memAlloc=%x\n",memAlloc);
-        return nullptr; // TODO raise an error?? how? We can't throw an exception here
-#else
-        auto r=sycl::malloc_shared<T>(sz,syclQ());
-        return r;
-#endif
-      }
-      void deallocate(T* p,size_t n) {
-        if (!p) return;
-        if (memAlloc && memAlloc->inAllocator(p)) {
-          memAlloc->deallocate(p,n);
-          return;
-        }
-#ifdef  __SYCL_DEVICE_ONLY__
-        printf("leaked %d bytes\n",n*sizeof(T));
-#else        
-        sycl::free(p,syclQ());
-#endif
-      }
-      bool operator==(const CellAllocator& x) const {return memAlloc==x.memAlloc;}
-    };
-    template <class T> CellAllocator<T> allocator() const {
-      return CellAllocator<T>(memAlloc);
-    }
-#else //!SYCL
-    template <class T> using CellAllocator=std::allocator<T>;
-    template <class T> CellAllocator<T> allocator() const {return CellAllocator<T>();}
-#endif
   };
 
-  template <class T>
-  void printAllocator(const char* prefix, const T&) {}
-  
-  template <class T>
-  void printAllocator(const char* prefix, const CellBase::CellAllocator<T>& x)
-  {printf("%s: allocator.desc=%p memAlloc=%p\n",prefix,x.desc,x.memAlloc);}
-  
-  
-  class SyclGraphBase
-  {
-  protected:
-#ifdef SYCL_LANGUAGE_VERSION
-    using MemAllocator=CellBase::MemAllocator;
-    MemAllocator *sharedMemAlloc=nullptr, *deviceMemAlloc=nullptr;
-    SyclGraphBase() {
-      // for now, do not use on-device dynamic allocation
-      sharedMemAlloc=sycl::malloc_shared<MemAllocator>(1, syclQ());
-      new(sharedMemAlloc) MemAllocator;
-      sharedMemAlloc->init(syclQ());
-//      deviceMemAlloc=sycl::malloc_shared<MemAllocator>(1, syclQ());
-//      new(deviceMemAlloc) MemAllocator;
-//      deviceMemAlloc->init(syclQ()); // TODO - run this on device?
-//      sharedMemAlloc->initialize(syclQ(), sycl::usm::alloc::shared, 512ULL * 1024ULL * 1024ULL); // TODO - expose this Python?
-//      deviceMemAlloc->initialize(syclQ(), sycl::usm::alloc::device, 512ULL * 1024ULL * 1024ULL); // TODO - expose this Python?
-    }
-    ~SyclGraphBase() {
-      sharedMemAlloc->~MemAllocator();
-      sycl::free(sharedMemAlloc,syclQ());
-//      deviceMemAlloc->~MemAllocator();
-//      sycl::free(deviceMemAlloc,syclQ());
-    }
-    // deleted because we're managing a resource here
-    SyclGraphBase(const SyclGraphBase&)=delete;
-    void operator=(const SyclGraphBase&)=delete;
-#endif      
-  };
-
-  template <class Cell> struct EcolabGraph:
-    public Exclude<SyclGraphBase>, public graphcode::Graph<Cell>
+  template <class Cell> struct EcolabGraph: public graphcode::Graph<Cell>
   {
 #ifdef SYCL_LANGUAGE_VERSION
     EcolabGraph(): graphcode::Graph<Cell>
@@ -310,8 +123,7 @@ namespace ecolab
 #endif
       for (size_t idx=0; idx<sz; ++idx) {
         auto& cell=*(*this)[idx]->template as<Cell>();
-        cell.m_idx=idx;
-        f(cell);
+        f(cell,idx);
       }
     }
     
@@ -328,8 +140,7 @@ namespace ecolab
           auto idx=i.get_global_linear_id();
           if (idx<this->size()) {
             auto& cell=*(*this)[idx]->template as<Cell>();
-            cell.m_idx=idx;
-            f(cell);
+            f(cell,idx);
           }
         });
       });
@@ -345,29 +156,63 @@ namespace ecolab
     template <class F>
     void groupedForAll(F f) {
 #ifdef SYCL_LANGUAGE_VERSION
-      // TODO - pass in workGroupSize as an optional parameter??
-      static size_t workGroupSize=32;//syclQ().get_device().get_info<sycl::info::device::max_work_group_size>();
-      syclQ().submit([&](auto& h) {
-        h.parallel_for(sycl::nd_range<1>(this->size()*workGroupSize, workGroupSize), [=,this](auto i) {
-          auto idx=i.get_group_linear_id();
-          if (idx<this->size()) {
-            auto& cell=*(*this)[idx]->template as<Cell>();
-            cell.m_idx=idx;
-            f(cell);
-          }
+      auto dev=syclQ().get_device();
+      // 1. Max threads per single work-group
+      size_t max_wg_size = dev.get_info<sycl::info::device::max_work_group_size>();
+
+      // 2. Hardware SIMD width (Sub-group size, usually 8, 16, or 32 on Intel)
+      std::vector<size_t> sg_sizes = dev.get_info<sycl::info::device::sub_group_sizes>();
+      size_t native_sg_size = sg_sizes.back(); // Usually 16 or 32 is best for compute
+
+      // 3. Total physical SLM available per hardware compute unit (DSS)
+      size_t max_slm_size = dev.get_info<sycl::info::device::local_mem_size>(); 
+
+      // 4. Maximum number of compute units (Execution units / DSS count)
+      uint32_t max_compute_units = dev.get_info<sycl::info::device::max_compute_units>();
+
+      size_t workGroupSize=native_sg_size;//256;
+//      if (workGroupSize > max_wg_size) 
+//        workGroupSize = max_wg_size; // Fallback for limited devices
+//      
+//      // Ensure it aligns perfectly with hardware SIMD lanes
+//      workGroupSize = (workGroupSize / native_sg_size) * native_sg_size;
+
+      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;
+
+      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;
+
+      DeviceType<int> fatalError(0);
+      for (size_t cellStart=0; cellStart<this->size(); cellStart+=num_work_groups)
+        syclQ().submit([&](auto& h) {
+          h.parallel_for(sycl::nd_range<1>(num_work_groups*workGroupSize, workGroupSize), [=,this,fatalError=&*fatalError](auto i) {
+            if (syclGroup().leader()) {
+              localAllocatorBuffer().next = 0;
+            } // reset local memory allocation
+            sycl::group_barrier(syclGroup()); 
+
+            auto idx = cellStart+i.get_group_linear_id();
+            auto stride = i.get_group_range(0);
+            if (idx < this->size()) {
+              auto& cell=*(*this)[idx]->template as<Cell>();
+              f(cell,idx);
+            }
+            // flag fatal error to throw afterwards.
+            if (fatalErrorFlag())
+              sycl::atomic_ref<int,sycl::memory_order::relaxed,sycl::memory_scope::device>(*fatalError).fetch_or(1);
+          });
         });
-      });
+      syclQ().wait_and_throw();
+      if (*fatalError)
+        throw std::runtime_error("Local Allocator Exhausted");
 #else
       hostForAll(f);
 #endif
     }
 
-    void syncThreads() {
-#ifdef SYCL_LANGUAGE_VERSION
-      syclQ().wait();
-#endif
-    }
-    
     template <class T, class F> T max(T x, F f) {
 #ifdef SYCL_LANGUAGE_VERSION
       DeviceType<T> r(x);
@@ -384,88 +229,13 @@ namespace ecolab
 #endif
       for (size_t idx=0; idx<sz; ++idx) {
         auto& cell=*(*this)[idx]->template as<Cell>();
-        cell.m_idx=idx;
         x=std::max(x,f(cell));
       }
       return x;
 #endif
     }
+    void syncThreads() {ecolab::syncThreads();}
   };
-
-  template <class T>
-  class GroupLocal
-  {
-#ifdef __SYCL_DEVICE_ONLY__
-    using BufferType=char[sizeof(T)];
-    using LocalBufferType=decltype(sycl::ext::oneapi::group_local_memory_for_overwrite<BufferType>(syclGroup()));
-    LocalBufferType buffer=sycl::ext::oneapi::group_local_memory_for_overwrite<BufferType>(syclGroup());
-  public:
-    template <class... Args>
-    GroupLocal(Args&&... args) {
-      if (syclGroup().leader())
-        new (*buffer) T(std::forward<Args>(args)...);
-      sycl::group_barrier(syclGroup());
-    }
-    ~GroupLocal() {
-      sycl::group_barrier(syclGroup());
-      if (syclGroup().leader())
-        (**this).~T();
-    }
-    T& operator*() {return reinterpret_cast<T&>(**buffer);}
-#else
-    T buffer;
-  public:
-    template <class... Args>
-    GroupLocal(Args&&... args): buffer(std::forward<Args>(args)...) {}
-    T& operator*() {return buffer;}
-#endif
-    T* operator->() {return &**this;}
-  };
-
-#ifdef __SYCL_DEVICE_ONLY__
-  template <class T, size_t N, size_t M>
-  T* getGroupBuffer()
-  {
-    if constexpr (N<=M)
-      return *sycl::ext::oneapi::group_local_memory_for_overwrite<T[1<<N]>(syclGroup());
-    assert(false && "Group buffer limit exceeded");
-    return nullptr;
-  }  
-  
-  /// Return a buffer in group local memory, of at least \a n Ts.
-  /// @tparam T element type of buffer
-  /// @tparam N maximum buffer size supported is 2^N.
-  /// This method can only be called in kernel code
-  /// Increasing N hurts performance, or may exceed system limits.
-  /// Individual elements are default initialised
-  /// Buffer is destroyed once all threads in group have completed (not at scope exit)
-  template <class T, size_t N>
-  T* groupBuffer(size_t n) {
-    // work out power of two >= n
-    switch (sizeof(size_t)*8-sycl::clz(n-1))
-      {
-      case 0: return getGroupBuffer<T,0,N>();
-      case 1: return getGroupBuffer<T,1,N>();
-      case 2: return getGroupBuffer<T,2,N>();
-      case 3: return getGroupBuffer<T,3,N>();
-      case 4: return getGroupBuffer<T,4,N>();
-      case 5: return getGroupBuffer<T,5,N>();
-      case 6: return getGroupBuffer<T,6,N>();
-      case 7: return getGroupBuffer<T,7,N>();
-      case 8: return getGroupBuffer<T,8,N>();
-      case 9: return getGroupBuffer<T,9,N>();
-      case 10: return getGroupBuffer<T,10,N>();
-      case 11: return getGroupBuffer<T,11,N>();
-      case 12: return getGroupBuffer<T,12,N>();
-      case 13: return getGroupBuffer<T,13,N>();
-      case 14: return getGroupBuffer<T,14,N>();
-      case 15: return getGroupBuffer<T,15,N>();
-      case 16: return getGroupBuffer<T,16,N>();
-      default: assert(false && "Group buffer limit exceeded"); return nullptr;
-      }
-  }
-#endif
-
 }
 
 namespace classdesc
diff --git a/include/non-sycl.h b/include/non-sycl.h
new file mode 100644
index 0000000..3a8a74c
--- /dev/null
+++ b/include/non-sycl.h
@@ -0,0 +1,34 @@
+/*
+  @copyright Russell Standish 2026
+  @author Russell Standish
+  This file is part of EcoLab
+
+  Open source licensed under the MIT license. See LICENSE for details.
+*/
+#ifndef ECOLAB_NON_SYCL_H
+#define ECOLAB_NON_SYCL_H
+namespace ecolab
+{
+  template <class T, ecolab::USMAlloc U> struct SyclQAllocator: public std::allocator<T> {};
+  
+  inline void syncThreads() {}
+
+  /// Non SYCL version of SyclType
+  template <class T>
+  struct SyclType
+  {
+    T data;
+    template <class... A> SyclType(A... args): data(std::forward<A>(args)...) {}
+    operator T() const {return data;}
+    operator T&() {return data;}
+  };
+  
+  template <class T> struct GroupLocal: public std::unique_ptr<T>
+  {
+  };
+  
+  template <class E, ecolab::USMAlloc UA>
+  struct SyclRandomEngine: public E {};
+
+}
+#endif
diff --git a/include/sparse_mat.h b/include/sparse_mat.h
index dde6b0b..abe66aa 100644
--- a/include/sparse_mat.h
+++ b/include/sparse_mat.h
@@ -27,22 +27,16 @@ namespace ecolab
     array_ns::array<F, A<F>> diag, val;
     /*row and columns of offdiagonal elements*/
     array_ns::array<unsigned, A<unsigned>> row, col;
-    sparse_mat(int s=0, int o=0, const A<F>& falloc={},
-               const A<unsigned>& ialloc={}): rowsz(s), colsz(s),
-                                              diag(s,falloc), val(o,falloc),
-                                              row(o,ialloc), col(o,ialloc) {} 
-    sparse_mat(const sparse_mat& x) 
+    sparse_mat(int s=0, int o=0): rowsz(s), colsz(s),
+                                              diag(s), val(o),
+                                              row(o), col(o) {} 
+    template <template<class> class A1>
+      sparse_mat(const sparse_mat<F,A1>& x) 
     {
       rowsz=x.rowsz; colsz=x.colsz; 
       diag=x.diag; val=x.val; 
       row=x.row; col=x.col;
     }
-    void setAllocators(const A<unsigned>& ialloc, const A<F>& falloc) {
-      diag.allocator(falloc);
-      val.allocator(falloc);
-      row.allocator(ialloc);
-      col.allocator(ialloc);
-    }
 
     /*matrix multiplication*/
     template <class E> typename
@@ -97,7 +91,7 @@ namespace ecolab
     assert(sigma>=0);
     assert(conn<diag.size());
 
-    array<int,A<int>> rsize(diag.size(),row.allocator());
+    array<unsigned,A<unsigned>> rsize(diag.size(),row.allocator());
     array<F,A<F>> tmp(diag.size(),diag.allocator());
     fillgrand(tmp);
     tmp=tmp*sigma+(double)conn;
diff --git a/include/sycl.h b/include/sycl.h
new file mode 100644
index 0000000..654958e
--- /dev/null
+++ b/include/sycl.h
@@ -0,0 +1,159 @@
+/*
+  @copyright Russell Standish 2026
+  @author Russell Standish
+  This file is part of EcoLab
+
+  Open source licensed under the MIT license. See LICENSE for details.
+*/
+
+#ifndef ECOLAB_SYCL_H
+#define ECOLAB_SYCL_H
+#include "usmAlloc.h" // TODO - merge into this file.
+#include "graphcode.h"
+#include <memory>
+#include <utility>
+#ifdef SYCL_LANGUAGE_VERSION
+
+#include <sycl/sycl.hpp>
+
+// we use some experimental features available in the OneAPI compiler
+#ifndef __INTEL_LLVM_COMPILER
+#error "EcoLab requires OneAPI compiler for some experimental functions"
+#endif
+
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+namespace ecolab
+{
+  using sycl::ext::oneapi::experimental::printf;
+
+  inline sycl::nd_item<1> syclItem() {return sycl::ext::oneapi::this_work_item::get_nd_item<1>();}
+  inline sycl::group<1> syclGroup() {return sycl::ext::oneapi::this_work_item::get_work_group<1>();}
+  inline sycl::sub_group syclSubGroup() {return sycl::ext::oneapi::this_work_item::get_sub_group();}
+
+  extern bool syclQDestroyed;
+  sycl::queue& syclQ();
+  void* reallocSycl(void*,size_t);
+
+  inline void syncThreads() {syclQ().wait_and_throw();}
+
+  /// SyclType is a pointer type, that when allocated via new is allocated from USM
+  template <class T>
+  struct SyclType: public T
+  {
+    T data;
+    template <class... A> SyclType(A... args): data(std::forward<A>(args)...) {}
+    operator T() const {return data;}
+    operator T&() {return data;}
+    void* operator new(size_t s) {return reallocSycl(nullptr,s);}
+    void operator delete(void* p) {reallocSycl(p,0);}
+    void* operator new[](size_t s) {return reallocSycl(nullptr,s);}
+    void operator delete[](void* p) {reallocSycl(p,0);}
+  };
+
+  template <class T, ecolab::USMAlloc UA>
+  struct SyclQAllocator: public graphcode::Allocator<T>
+  {
+    SyclQAllocator(): graphcode::Allocator<T>(syclQ(), UA) {}
+    template<class U> struct rebind {using other=SyclQAllocator<U,UA>;};
+  };
+  
+  // random numbers
+  template <class E, ecolab::USMAlloc UA>
+  struct SyclRandomEngine
+  {
+    std::vector<E,SyclQAllocator<E,UA>> rngs;
+    static constexpr unsigned incr=1000000;
+    SyclRandomEngine(size_t numRngs): rngs(numRngs) {}
+    void seed() {for (auto& i: rngs) i.seed();}
+    void seed(unsigned s) {
+#ifdef __SYCL_DEVICE_ONLY__
+      rngs[syclItem().get_global_linear_id() % rngs.size()].
+        seed(s+syclItem().get_global_linear_id()*incr);
+#else
+      size_t numRngs=rngs.size();
+#pragma omp parallel for
+      for (unsigned j=0; j<numRngs; ++j)
+        rngs[j].seed(s+j*incr);
+#endif
+    }
+    unsigned operator()() {
+#ifdef __SYCL_DEVICE_ONLY__
+      return rngs[syclItem().get_global_linear_id() % rngs.size()]();
+#elif defined(_OPENMP)
+      return rngs[omp_get_thread_num()%rngs.size()]();
+#else
+      return rngs[0]();
+#endif
+    }
+    unsigned min() const {return rngs[0].min();}
+    unsigned max() const {return rngs[0].max();}
+  };
+}
+#else // !SYCL
+#include "non-sycl.h"
+#endif
+
+namespace ecolab
+{
+  /// Used to declare an object in unified shared memory, accessible on GPU device
+  /// Smart pointer semantics, but never invalid
+  template <class M>
+  class DeviceType
+  {
+    M* const model;
+  public:
+    using element_type=M;
+    template <class... A>
+#ifdef SYCL_LANGUAGE_VERSION
+    DeviceType(A... args): model(sycl::malloc_shared<M>(1,syclQ()))
+      {new (model) M(std::forward<A>(args)...);}
+#else
+    DeviceType(A... args): model(new M(std::forward<A>(args)...)) {}
+#endif
+    DeviceType& operator=(const DeviceType& x) {*model=*x.model; return *this;}
+    ~DeviceType() {
+#ifdef SYCL_LANGUAGE_VERSION
+      model->~M();
+      sycl::free(model,syclQ());
+#else
+      delete model;
+#endif
+    }
+    M& operator*() {return *model;}
+    const M& operator*() const {return *model;}
+    M* operator->() {return model;}
+    const M* operator->() const {return model;}
+    operator bool() const {return true;} // always defined
+  };
+
+#ifndef __SYCL_DEVICE_ONLY__
+  template <class T> using LocalAllocator=std::allocator<T>;
+#endif
+  
+  inline void groupBarrier() {
+#ifdef __SYCL_DEVICE_ONLY__
+    sycl::group_barrier(syclGroup());
+#endif
+  }
+
+  inline bool groupLeader() {
+#ifdef __SYCL_DEVICE_ONLY__
+    return syclGroup().leader();
+#endif
+    return true;
+  }
+
+
+}
+
+
+#define CLASSDESC_RESTProcess___ecolab__DeviceType_M_
+#define CLASSDESC_json_pack___ecolab__DeviceType_M_
+#define CLASSDESC_json_unpack___ecolab__DeviceType_M_
+#define CLASSDESC_pack___ecolab__DeviceType_M_
+#define CLASSDESC_unpack___ecolab__DeviceType_M_
+#include "sycl.cd"
+#endif
diff --git a/include/usmAlloc.h b/include/usmAlloc.h
index 2b5bb11..f0999e7 100644
--- a/include/usmAlloc.h
+++ b/include/usmAlloc.h
@@ -8,6 +8,10 @@
 
 #ifndef ECOLAB_USMALLOC_H
 #define ECOLAB_USMALLOC_H
+#ifdef SYCL_LANGUAGE_VERSION
+#include <sycl/sycl.hpp>
+#endif
+
 namespace ecolab
 {
   // this type of macro stuff doesn't play nice with classdesc, hence why
diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc
index bbc686a..c90e4c0 100644
--- a/models/ecolab_model.cc
+++ b/models/ecolab_model.cc
@@ -31,8 +31,10 @@ using namespace classdesc;
 
 namespace model
 {
+#ifndef __SYCL_DEVICE_ONLY__
   PanmicticModel panmictic_ecolab;
   CLASSDESC_ADD_GLOBAL(panmictic_ecolab);
+#endif
   DeviceType<SpatialModel> spatial_ecolab;
   CLASSDESC_ADD_GLOBAL(spatial_ecolab);
   CLASSDESC_PYTHON_MODULE(ecolab_model);
@@ -52,14 +54,13 @@ void PanmicticModel::generate(unsigned niter)
 } 
 
 /* Rounding function, randomly round up or down, in the range 0..INT_MAX */
-template <class B>
-int EcolabPoint<B>::ROUND(Float x) 
+int EcolabPoint::ROUND(Float x) 
 {
   Float dum;
   const Float maxInt=Float(std::numeric_limits<int>::max()-1);
-  if (x<0) x=0;
-  if (x>maxInt) x=maxInt;
-  //syclPrintff("ROUND inner x=%g, modf=%g, rand()=%g, rand.max=%g, rand.min=%g\n",x,std::fabs(std::modf(x,&dum)),rand(),rand.max(),rand.min());
+  if (x<=0) return 0;
+  if (x>=maxInt) return maxInt;
+  //printf("ROUND inner x=%g, modf=%g, rand()=%u, rand.max=%u, rand.min=%u\n",x,std::fabs(std::modf(x,&dum)),rand(),rand.max(),rand.min());
   return std::fabs(std::modf(x,&dum))*(rand.max()-rand.min()) > (rand()-rand.min()) ?
     (int)x+1 : (int)x;
 }
@@ -75,7 +76,6 @@ struct RoundArray
   int operator[](size_t i) const //{return point.ROUND(expr[i]);}
   {
     auto r=point.ROUND(expr[i]);
-    //syclPrintff("ROUND: %d, %g=%d\n",i,expr[i],r);
     return r;
   }
 };
@@ -83,81 +83,37 @@ struct RoundArray
 namespace ecolab::array_ns
 {template <class E, class P> struct is_expression<RoundArray<E,P>>: public true_type {};}
 
-template <class B>
 template <class E>
-RoundArray<E,EcolabPoint<B>> EcolabPoint<B>::roundArray(const E& expr)
-{return RoundArray<E,EcolabPoint<B>>(*this,expr);}
-
-template <> void setArray(array<int,std::allocator<int>>& x, const array<int>& y)
-{x=y;}
-template <> array<int> getArray(const array<int,std::allocator<int>>& x) {return x;}
-
-#ifdef SYCL_LANGUAGE_VERSION
-template <>
-array<int> getArray(const array<int,ecolab::CellBase::CellAllocator<int>>& x) 
-{
-  DeviceType<size_t> size;
-  DeviceType<const int*> xData;
-  syclQ().single_task([size=&*size,xData=&*xData,x=&x](){
-    *size=x->size();
-    *xData=x->data();
-  }).wait();
-  array<int> r(*size);
-  syclQ().copy(*xData,r.data(),*size);
-  return r;
-}
-
-template <>
-void setArray(array<int,ecolab::CellBase::CellAllocator<int>>& x, const array<int>& y)
+RoundArray<E,EcolabPoint> EcolabPoint::roundArray(const E& expr)
+{return RoundArray<E,EcolabPoint>(*this,expr);}
+
+//template <> void setArray(array<int,std::allocator<int>>& x, const array<int>& y)
+//{x=y;}
+//template <> array<int> getArray(const array<int,std::allocator<int>>& x) {return x;}
+//
+void EcolabPoint::generate(unsigned niter, const ModelData& model)
 {
-  auto size=y.size();
-  DeviceType<int*> xData;
-  syclQ().single_task([size,x=&x,xData=&*xData](){
-    //array<int,typename B::template CellAllocator<int>> tmp(size,this->template allocator<int>());
-    //m_density.swap(tmp);
-    x->resize(size);
-    *xData=x->data(); //return allocated data pointer to host
-  }).wait();
-  syclQ().copy(y.data(),*xData,size);
-}
-#endif
-
-
-template <class B>
-void EcolabPoint<B>::generate(unsigned niter, const ModelData& model)
-{
-#ifdef __SYCL_DEVICE_ONLY__
-  Float* interactionResult=groupBuffer<Float,SpatialModel::log2MaxNsp>(density.size());
-#else
-  array<Float> interactionResult(density.size());
-#endif
+  array<int,LocalAllocator<int>> lDensity(density), tmp(density.size());
+  
   for (unsigned step=0; step<niter; step++)
     {
-      array_ns::map(density.size(),  [&](size_t i){
-        interactionResult[i]=model.interaction.diag[i]*density[i];
+      array_ns::map(lDensity.size(),  [&](size_t i){
+        Float ir=model.interaction.diag[i]*lDensity[i];
+        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();
-      array_ns::map(model.interaction.row.size(), [&](size_t i){
-#ifdef __SYCL_DEVICE_ONLY__
-        sycl::atomic_ref<Float, sycl::memory_order::relaxed, sycl::memory_scope::work_group>
-#endif
-          (interactionResult[model.interaction.row[i]]) +=
-          model.interaction.val[i]*density[model.interaction.col[i]];
-      });
-      
-      groupBarrier();
-      array_ns::map(density.size(),  [&](size_t i){
-        density[i]=ROUND(density[i] + density[i] * (model.repro_rate[i] + interactionResult[i]));
-      });
-      groupBarrier();
+      lDensity.swap(tmp);
+      groupBarrier(); // synchronises threads on each iteration
     }
+  density=lDensity;
+  assert(all(density>=0));
 //  // sequential/non-GPU version
 //  for (unsigned i=0; i<niter; i++)
 //    {density = roundArray(density + density * (model.repro_rate + model.interaction*density));}
 }
 
-template <class B>
-unsigned EcolabPoint<B>::nsp() const
+unsigned EcolabPoint::nsp() const
 {return sum(density!=0);}
 
 array<unsigned> SpatialModel::nsp() const
@@ -167,8 +123,7 @@ array<unsigned> SpatialModel::nsp() const
   return nsp;
 }
 
-template <class B>
-void EcolabPoint<B>::condense(const array<bool>& mask, size_t mask_true)
+void EcolabPoint::condense(const array<bool>& mask, size_t mask_true)
 {
   density = pack( density, mask, mask_true); 
 }
@@ -191,6 +146,7 @@ void ModelData::condense(const array<bool>& mask, size_t mask_true)   /* remove
   interaction.row = map[pack(interaction.row, mask_off,mask_off_true)];
   interaction.col = map[pack(interaction.col, mask_off,mask_off_true)];
 
+  computeODiagIdx();
   foodweb = interaction;
 }
 
@@ -206,7 +162,7 @@ unsigned PanmicticModel::condense()
 
 unsigned SpatialModel::condense()
 {
-  array<int> total_density(species.size());
+  array<int> total_density(species.size(),0);
   for (auto& i: *this) total_density+=i->as<EcolabCell>()->density; // TODO
 #ifdef MPI_SUPPORT
   array<int> recv(total_density.size());
@@ -241,37 +197,21 @@ void PanmicticModel::mutate()
 
 void SpatialModel::mutate()
 {
-
-  
-  
-  DeviceType<array<Float,graphcode::Allocator<Float>>> mut_scale;
-#ifdef SYCL_LANGUAGE_VERSION
-  using ArrayAlloc=CellBase::CellAllocator<unsigned>;
-  using NewSpAlloc=graphcode::Allocator<array<unsigned,ArrayAlloc>>;
-  vector<array<unsigned,ArrayAlloc>,NewSpAlloc> newSp
-    (size(),NewSpAlloc(syclQ(),sycl::usm::alloc::shared));
-  mut_scale->allocator(graphcode::Allocator<Float>(syclQ(),sycl::usm::alloc::shared));
-#else
-  vector<array<unsigned>> newSp(size());
-#endif
+  DeviceType<array<Float,ModelData::Allocator<Float>>> mut_scale(species.size());
   *mut_scale=sp_sep * repro_rate * mutation * int(tstep - last_mut_tstep);
+  assert(all(*mut_scale<=1));
   last_mut_tstep=tstep;
 
-  groupedForAll([newSp=newSp.data(),mut_scale=&*mut_scale,this](EcolabCell& c) {
-    auto tmp{c.mutate(*mut_scale)};
-#ifdef __SYCL_DEVICE_ONLY__
-    if (syclGroup().leader() && tmp.size())
-#endif
-      newSp[c.idx()]=tmp;
+  vector<EcolabPoint::UnsignedArray,ModelData::Allocator<EcolabPoint::UnsignedArray>> newSp(size());
+
+  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<array<unsigned,graphcode::Allocator<unsigned>>> cell_ids;
-#ifdef SYCL_LANGUAGE_VERSION
+  DeviceType<array<unsigned,ModelData::Allocator<unsigned>>> cell_ids;
   syncThreads();
-#endif
-
-
   
   // TODO - this is a kind of scan - can it be done on device?
   size_t j=0;
@@ -283,10 +223,12 @@ void SpatialModel::mutate()
     }
 
   // deallocate on device
-  groupedForAll([newSp=newSp.data()](EcolabCell& c) {
-    newSp[c.idx()].clear();
+  groupedForAll([newSp=newSp.data()](EcolabCell& c,size_t i) {
+    newSp[i].clear();
+    assert(newSp[i].refCnt()==0);
   });
-  
+   
+
 #ifdef MPI_SUPPORT
   MPIbuf b; b<<new_sp<<(*cell_ids); b.gather(0);
   if (myid()==0)
@@ -302,63 +244,52 @@ void SpatialModel::mutate()
     }
   MPIbuf() << (*cell_ids) << *(ModelData*)this << bcast(0)
            >> (*cell_ids) >> *(ModelData*)this;
+  if (cell_ids->size()==0) return;
 #else
+  if (new_sp.size()==0) return;
   ModelData::mutate(new_sp);
 #endif
-  if (cell_ids->size()==0) return;
+
+  computeODiagIdx();
+  
   // set the new species density to 1 for those created on this cell
   //groupedForAll([cell_ids=&*cell_ids](EcolabCell& c) {
-  hostForAll([cell_ids=&*cell_ids](EcolabCell& c) {
+  hostForAll([cell_ids=&*cell_ids,this](EcolabCell& c,size_t) {
+    assert(all(c.density>=0));
     c.density <<= (*cell_ids)==c.id;
+    assert(c.density.size()==this->species.size());
   });
 }
 
-template <class B>
 template <class E>
-array<unsigned,typename EcolabPoint<B>::template Allocator<unsigned>>
-EcolabPoint<B>::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 {};
 #ifdef __SYCL_DEVICE_ONLY__
+  LocalArray speciations=roundArray(mut_scale * density);
   auto nsp=density.size();
-  auto speciations=groupBuffer<unsigned,SpatialModel::log2MaxNsp>(nsp);
+  //  auto new_sp = gen_index(speciations);
+  LocalArray offsets(nsp+1);
+  sycl::joint_exclusive_scan(syclGroup(),speciations.data(),speciations.data()+nsp,
+                             offsets.data(),sycl::plus<unsigned>());
   groupBarrier();
-  asg_v(speciations, nsp, roundArray(mut_scale * density));
-
-  auto numNewSp=groupBuffer<unsigned,SpatialModel::log2MaxNsp>(syclGroup().get_group_linear_range());
-  auto myId=syclGroup().get_local_linear_id();
-  
-  numNewSp[myId]=0;
-  array_ns::map(nsp, [&](size_t i) {
-    numNewSp[myId]+=speciations[i];
-  });
-
-  auto offsets=groupBuffer<unsigned,SpatialModel::log2MaxNsp>(syclGroup().get_group_linear_range()+1);
-  if (syclGroup().leader())
-    {
-      offsets[0]=0;
-      for (size_t i=1; i<=syclGroup().get_group_linear_range(); ++i)
-        offsets[i]=offsets[i-1]+numNewSp[i-1];
-    }
-  groupBarrier();
-  if (offsets[syclGroup().get_group_linear_range()])
-    array_ns::asg_minus_v(density.data(), nsp, speciations);
-  GroupLocal<array<unsigned,typename EcolabPoint<B>::template Allocator<unsigned>>> new_sp
-    (offsets[syclGroup().get_group_linear_range()], this->template allocator<unsigned>());
+  if (groupLeader())
+    offsets[nsp]=offsets[nsp-1]+speciations[nsp-1];
   groupBarrier();
 
+  if (offsets[nsp]==0) return {};
+  
+  density-=speciations;
 
-  // gen_index
-  auto p=offsets[myId];
-  auto n=new_sp->data();
-  array_ns::map(nsp, [&](size_t i) {
-    for (size_t j=0; j<speciations[i]; ++j) {
-      n[p++]=i;
-    }
+  LocalArray new_sp(offsets[nsp]);
+  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;
   });
-
+  
   groupBarrier();
-  return *new_sp;
+  return new_sp;
 #else
   array<unsigned> speciations=roundArray(mut_scale * density);
   auto new_sp = gen_index(speciations);
@@ -662,7 +593,7 @@ void SpatialModel::setGrid(size_t nx, size_t ny)
 void SpatialModel::generate(unsigned niter)
 {
   if (tstep==0) makeConsistent();
-  groupedForAll([=,this](EcolabCell& c) {
+  groupedForAll([=,this](EcolabCell& c,size_t) {
     c.generate(niter,*this);
   });
   tstep+=niter;
@@ -671,7 +602,8 @@ void SpatialModel::generate(unsigned niter)
 unsigned SpatialModel::migrate()
 {
   /* each cell gets a distinct random salt value */
-  hostForAll([=,this](EcolabCell& c) {c.salt=Float(c.rand()-c.rand.min())/(c.rand.max()-c.rand.min());});
+  hostForAll([=,this](EcolabCell& c,size_t)
+  {c.salt=Float(c.rand()-c.rand.min())/(c.rand.max()-c.rand.min());});
   
   prepareNeighbours();
 
@@ -682,15 +614,16 @@ unsigned SpatialModel::migrate()
   const Float cap=1.0/maxNbrs;
   array<Float> capped_migration = merge(mm>cap,cap,mm);
   
-  hostForAll([&,this](EcolabCell& c) {
+  hostForAll([&,this](EcolabCell& c, size_t i) {
+    assert(all(c.density>=0));
     /* loop over neighbours */
     for (auto& n: c) 
       {
         auto& nbr=*n->as<EcolabCell>();
         Float salt=&c<&nbr? c.salt: nbr.salt;
         array<int> m=capped_migration * (nbr.density-c.density);
-        delta[c.idx()]+=m;//*(1 + salt * (abs(m)<cap*array_ns::min(nbr.density,c.density)));
-        assert(all(c.density>=-delta[c.idx()]));
+        delta[i]+=m;//*(1 + salt * (abs(m)<cap*array_ns::min(nbr.density,c.density)));
+        assert(all(c.density>=-delta[i]));
       }
   });
 
@@ -705,14 +638,14 @@ unsigned SpatialModel::migrate()
   for (size_t i=0; i<numCells; ++i)
     {
       auto& c=*(*this)[i]->as<EcolabCell>();
-      c.density+=delta[c.idx()];
+      c.density+=delta[i];
       assert(all(c.density>=0));
-      totalMigration+=sum(abs(delta[c.idx()]));
+      totalMigration+=sum(abs(delta[i]));
 #if !defined(NDEBUG)
 #ifdef _OPENMP
 #pragma omp critical
 #endif
-      ssum+=delta[c.idx()];
+      ssum+=delta[i];
 #endif
     }
   last_mig_tstep=tstep;
@@ -757,19 +690,6 @@ unsigned SpatialModel::migrate()
 
 void ModelData::makeConsistent(size_t nsp)
 {
-#ifdef SYCL_LANGUAGE_VERSION
-  if (sycl::get_pointer_type(species.data(),syclQ().get_context())!=sycl::usm::alloc::shared)
-    {
-      FAlloc falloc(syclQ(),sycl::usm::alloc::shared);
-      species.allocator(graphcode::Allocator<int>(syclQ(),sycl::usm::alloc::shared));
-      create.allocator(falloc);
-      repro_rate.allocator(falloc);
-      mutation.allocator(falloc);
-      migration.allocator(falloc);
-      interaction.setAllocators
-        (graphcode::Allocator<unsigned>(syclQ(),sycl::usm::alloc::shared),falloc);
-    }
-#endif
   if (!species.size())
     {
       species=pcoord(nsp);
@@ -778,26 +698,14 @@ void ModelData::makeConsistent(size_t nsp)
   if (!create.size()) create.resize(species.size(),0);
   if (!mutation.size()) mutation.resize(species.size(),0);
   if (!migration.size()) migration.resize(species.size(),0);
+  computeODiagIdx();
 }
 
-void SpatialModel::setDensitiesShared()
+void ModelData::computeODiagIdx()
 {
-#ifdef SYCL_LANGUAGE_VERSION
-  groupedForAll([=,this](EcolabCell& c) {
-    c.memAlloc=sharedMemAlloc;
-    c.density.allocator(c.allocator<int>());
-  });
-#endif
-}
-    
-void SpatialModel::setDensitiesDevice()
-{
-#ifdef SYCL_LANGUAGE_VERSION
-  groupedForAll([=,this](EcolabCell& c) {
-    c.memAlloc=deviceMemAlloc;
-    c.density.allocator(c.allocator<int>());
-  });
-#endif
+  oDiagIdx.clear(); oDiagIdx.resize(species.size());
+  for (unsigned i=0; i<interaction.row.size(); ++i)
+    oDiagIdx[interaction.row[i]]<<=i;
 }
 
 void SpatialModel::makeConsistent()
@@ -808,13 +716,6 @@ void SpatialModel::makeConsistent()
 #ifdef MPI_SUPPORT
   MPI_Allreduce(MPI_IN_PLACE,&nsp,1,MPI_UNSIGNED_LONG,MPI_MAX,MPI_COMM_WORLD);
 #endif
-  hostForAll([=,this](EcolabCell& c) {
-#ifdef SYCL_LANGUAGE_VERSION
-        if (!c.memAlloc) c.memAlloc=sharedMemAlloc;
-#endif
-        // not needed, as we're not resizing density on device
-        if (nsp>c.density.size()) c.density.allocator(c.allocator<int>());
-  });
   ModelData::makeConsistent(nsp);
   syncThreads();
 }
diff --git a/models/ecolab_model.h b/models/ecolab_model.h
index a12434b..e6fb9c7 100644
--- a/models/ecolab_model.h
+++ b/models/ecolab_model.h
@@ -38,11 +38,17 @@ struct ConnectionPlot: public Object<ConnectionPlot, ConnectionPlotBase>
 
 struct ModelData
 {
-  using FAlloc=graphcode::Allocator<Float>;
-  using IAlloc=graphcode::Allocator<int>;
-  array<int,IAlloc> species;
-  array<Float,FAlloc> create, repro_rate, mutation, migration;
-  sparse_mat<Float,graphcode::Allocator> interaction;
+#ifdef SYCL_LANGUAGE_VERSION
+  template <class T> using Allocator=HostSharedAllocator<T>;
+#else
+  template <class T> using Allocator=std::allocator<T>;
+#endif
+
+  array<int,Allocator<int>> species;
+  array<Float,Allocator<Float>> create, repro_rate, mutation, migration;
+  sparse_mat<Float,Allocator> interaction;
+  using UnsignedArray=array<unsigned,Allocator<unsigned>>;
+  std::vector<UnsignedArray,Allocator<UnsignedArray>> oDiagIdx;
   sparse_mat_graph foodweb;
   unsigned long long tstep=0, last_mut_tstep=0, last_mig_tstep=0;
   //mutation parameters
@@ -51,6 +57,7 @@ struct ModelData
   bool fixMutation=false, fixMigration=false;
   
   void makeConsistent(size_t nsp);
+  void computeODiagIdx();
   void random_interaction(unsigned conn, double sigma);
   void condense(const array<bool>& mask, size_t mask_true);
   void mutate(const array<int>&); 
@@ -65,33 +72,36 @@ template <class T,class A> void setArray(array<T,A>&, const array<T>&);
 template <class T,class A> array<T> getArray(const array<T,A>&);
 
 /* ecolab cell  */
-template <class CellBase>
-class EcolabPoint: public Exclude<CellBase>
+class EcolabPoint
 {
 public:
+#ifdef SYCL_LANGUAGE_VERSION
+  template <class T> using Allocator=GlobalDeviceAllocator<T>;
+#else
+  template <class T> using Allocator=std::allocator<T>;
+#endif
+  using UnsignedArray=array<unsigned,Allocator<unsigned>>;
+  using LocalArray=array<unsigned,LocalAllocator<unsigned>>;
+
   Float salt;  /* random no. used for migration */
-  template <class T> using Allocator=typename CellBase::template CellAllocator<T>;
-  array<int,Allocator<int>> density{this->template allocator<int>()};
+  array<int,Allocator<int>> density;
   
   void generate(unsigned niter, const ModelData&);
   void condense(const array<bool>& mask, size_t mask_true);
-  template <class E>
-  array<unsigned,Allocator<unsigned>> 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);
   template <class E> RoundArray<E,EcolabPoint> roundArray(const E& expr);
+#ifdef SYCL_LANGUAGE_VERSION
+  Exclude<SyclRandomEngine<std::mt19937,USMAlloc::shared>> rand
+    {syclQ().get_device().get_info<sycl::info::device::max_work_group_size>()}; // TODO make configurable?
+#else
   Exclude<std::mt19937> rand; // random number generator
+#endif
 };
 
-// for the panmictic model, we need to use std::allocator
-struct AllocatorBase
-{
-  template <class T> using CellAllocator=std::allocator<T>;
-  template <class T> CellAllocator<T> allocator() const {return CellAllocator<T>();}
-};
-
-struct PanmicticModel: public ModelData, public EcolabPoint<AllocatorBase>, public ecolab::Model<PanmicticModel>
+struct PanmicticModel: public ModelData, public EcolabPoint, public ecolab::Model<PanmicticModel>
 {
   ConnectionPlot connectionPlot;
   void updateConnectionPlot() {connectionPlot.update(this->density,interaction);}
@@ -109,10 +119,9 @@ struct PanmicticModel: public ModelData, public EcolabPoint<AllocatorBase>, publ
   Float extinctionConnectivity() const;
 };
 
-struct EcolabCell: public EcolabPoint<ecolab::CellBase>, public graphcode::Object<EcolabCell>
+struct EcolabCell: public EcolabPoint, public graphcode::Object<EcolabCell>
 {
-  unsigned id=0; // stash the graphcode node id here - TODO - why can't we call graphcode::Object<EcolabCell>::id()?
-  array<int,EcolabPoint::Allocator<int>> delta{this->template allocator<int>()};
+  unsigned id=0; // stash the graphcode node id here so it is accessible from within the cell, rather than its pointer wrapper
 };
 
 class SpatialModel: public ModelData, public EcolabGraph<EcolabCell>,
@@ -122,19 +131,16 @@ class SpatialModel: public ModelData, public EcolabGraph<EcolabCell>,
   CLASSDESC_ACCESS(SpatialModel);
   size_t maxNbrs=0;
 public:
-  static constexpr size_t log2MaxNsp=10;
+  static constexpr size_t log2MaxNsp=11;
   // function valid for x∈(-numX,∞], y∈(-numY,∞]
   size_t makeId(size_t x, size_t y) const {return (x+numX)%numX + numX*((y+numY)%numY);}
   void setGrid(size_t nx, size_t ny);
   EcolabCell& cell(size_t x, size_t y) {
     return *objects[makeId(x,y)];
   }
-  /// on GPUs, set the 
-  void setDensitiesShared();
-  void setDensitiesDevice();
   array<unsigned> nsp() const;
   void makeConsistent();
-  void seed(unsigned x) {forAll([=](EcolabCell& cell){cell.rand.seed(x);});}
+  void seed(unsigned x) {groupedForAll([=](EcolabCell& cell,size_t){cell.rand.seed(x);});}
   void generate(unsigned niter);
   void generate() {generate(1);}
   /// returns number of extinctions
diff --git a/models/spatial_ecolab.py b/models/spatial_ecolab.py
index 0c59c24..365428d 100644
--- a/models/spatial_ecolab.py
+++ b/models/spatial_ecolab.py
@@ -13,14 +13,14 @@
 array_urand.seed(10+myid())
 
 # initial number of species
-nsp=30
+nsp=100
 
 ecolab.repro_min(-0.1)
 ecolab.repro_max(0.1)
 ecolab.odiag_min(-1e-5)
 ecolab.odiag_max(1e-5)
-#ecolab.mut_max(1e-4)
-ecolab.mut_max(1e-3)
+ecolab.mut_max(1e-5)
+#ecolab.mut_max(1e-3)
 ecolab.sp_sep(0.1)
 
 def randomList(num, min, max):
@@ -28,8 +28,10 @@ def randomList(num, min, max):
 
 ecolab.species(range(nsp))
 
-numX=8
-numY=8
+numX=12
+numY=12
+#numX=2
+#numY=2
 ecolab.setGrid(numX,numY)
 ecolab.partitionObjects()
 
@@ -59,7 +61,6 @@ def randomList(num, min, max):
 extinctions=0
 migrations=0
 def stepImpl():
-    #ecolab.setDensitiesDevice()
     ecolab.generate(100)
     ecolab.mutate()
 
@@ -73,8 +74,6 @@ def stepImpl():
     migrations+=ecolab.migrate()
     extinctions+=ecolab.condense()
     #ecolab.syncThreads()
-    #print(ecolab.nsp()())
-    #ecolab.setDensitiesShared()
     #ecolab.gather()
 
 print(ecolab.nsp()())
@@ -83,15 +82,15 @@ def stepImpl():
 print(ecolab.nsp()())
 
 from timeit import timeit
-print(timeit('stepImpl()', globals=globals(), number=10))
+print(timeit('stepImpl()', globals=globals(), number=100))
                 
 def step():
     global extinctions,migrations
     extinctions=0
     migrations=0
-    for i in range(epoch//10000):
+    for i in range(epoch//1000000):
         stepImpl()
-    print('migrations=',migrations,' extinctions=',extinctions)
+    #print('migrations=',migrations,' extinctions=',extinctions)
     if myid()==0:
         nsp=len(ecolab.species)
         statusBar.configure(text=f't={ecolab.tstep()} nsp:{nsp}')
diff --git a/src/ecolab.cc b/src/ecolab.cc
index 0bd7786..afd5feb 100644
--- a/src/ecolab.cc
+++ b/src/ecolab.cc
@@ -49,7 +49,9 @@ namespace
   
   struct SyclQ: public queue
   {
-    SyclQ(): queue(default_selector_v, errHandler) {}
+    // in_order version for debugging purposes
+    SyclQ(): queue(default_selector_v, errHandler, sycl::property::queue::in_order{}) {}
+    //SyclQ(): queue(default_selector_v, errHandler) {}
     ~SyclQ() {syclQDestroyed=true;}
   };
 }