From 3ecbf9fdc9297ddde53ebf513adc92a1ec6a54a9 Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Mon, 6 Jul 2026 12:28:42 +1000 Subject: [PATCH 01/28] feat: GlobalDeviceAllocator for SYCL. --- include/DeviceAllocator.h | 126 +++++++++++++++++++++++ include/ecolab.h | 141 +++---------------------- include/sparse_mat.h | 18 ++-- include/sycl.h | 116 +++++++++++++++++++++ models/ecolab_model.cc | 211 ++++++++++++++++++-------------------- models/ecolab_model.h | 50 ++++----- models/spatial_ecolab.py | 4 +- src/ecolab.cc | 4 +- 8 files changed, 393 insertions(+), 277 deletions(-) create mode 100644 include/DeviceAllocator.h create mode 100644 include/sycl.h diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h new file mode 100644 index 0000000..6c46396 --- /dev/null +++ b/include/DeviceAllocator.h @@ -0,0 +1,126 @@ +/* + @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; + + // Lock free circular buffer queue for SYCL + template + class Queue + { + static_assert((size&(size-1))==0,"size must be power of two"); + constexpr static unsigned mask=size-1; + unsigned queue[size]; + unsigned head=0, tail=0; + using Atomic=sycl::atomic_ref; + + public: + void enqueue(unsigned x) + { + Atomic h(head); + queue[h++ & mask]=x; + } + unsigned dequeue() + { + Atomic h(head), t(tail); + if (h==t) return ~0U; // signal buffer empty, don't wait + return queue[t++ & mask]; + } + }; + + template class DeviceAllocator; + /// empty allocator to terminate template recursion + template <> class DeviceAllocator { + public: + void init(sycl::queue& q) {} + void* allocate(size_t sz) {return nullptr;} + void deallocate(void* p, size_t) {sycl::ext::oneapi::experimental::printf("%x leaked on device\n",p);} + }; + + template class DeviceAllocator + { + constexpr static unsigned pageSize=1< queue; + char memory[poolSize]; + DeviceAllocator nextAllocator; // next size up allocator + public: + void init(sycl::queue& q) { + q.parallel_for(numPages, [this](auto i) { + queue.enqueue(i<=memory && p(p)-memory)>>order); + return; + } + nextAllocator.deallocate(p,size); + } + + bool inAllocator(void* p) const {return p>=this && p& deviceAllocator() { + static DeviceType> deviceAllocator; + static int initialised=(deviceAllocator->init(syclQ()),0); + return *deviceAllocator; + } + + /// Allocator wrapping the DeviceAllocator singleton + template + 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 + GlobalDeviceAllocator(const GlobalDeviceAllocator& other): + allocator(other.allocator) {} + + T* allocate(size_t n) + {return reinterpret_cast(allocator->allocate(n*sizeof(T)));} + void deallocate(T* p, size_t n){allocator->deallocate(p,n*sizeof(T));} + template struct rebind {using other=GlobalDeviceAllocator;}; + }; + + template + struct HostSharedAllocator: public graphcode::Allocator + { + HostSharedAllocator(): graphcode::Allocator(syclQ(), sycl::usm::alloc::shared) {} + template struct rebind {using other=HostSharedAllocator;}; + }; + +} +#endif diff --git a/include/ecolab.h b/include/ecolab.h index 631c8bf..2a95682 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 -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 @@ -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,98 +97,6 @@ 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 - struct SyclType: public T - { - template SyclType(A... args): T(std::forward(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 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 - struct SyclType - { - 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 DeviceType - { - SyclType* const model; - public: - using element_type=M; - template - DeviceType(A... args): model(new SyclType(std::forward(args)...)) {} - DeviceType(const DeviceType& x): model(new SyclType(*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 - struct SyclQAllocator: public graphcode::Allocator - { -#ifdef SYCL_LANGUAGE_VERSION - SyclQAllocator(): graphcode::Allocator(syclQ(), UA) {} -#endif - template struct rebind {using other=SyclQAllocator;}; - }; - class GlobalDeallocateKernelTag; struct CellBase { @@ -217,7 +114,7 @@ namespace ecolab CellAllocator(MemAllocator* memAlloc): memAlloc(memAlloc) {} template CellAllocator(const CellAllocator& x): memAlloc(x.memAlloc) {} T* allocate(size_t sz) { -#ifdef __SYCL_DEVICE_ONLY__ + //#ifdef __SYCL_DEVICE_ONLY__ if (memAlloc) { auto r=reinterpret_cast(memAlloc->allocate(sz*sizeof(T))); if (!r) printf("Mem allocation failed\n"); @@ -225,10 +122,10 @@ namespace ecolab } 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(sz,syclQ()); - return r; -#endif +//#else +// auto r=sycl::malloc_shared(sz,syclQ()); +// return r; +//#endif } void deallocate(T* p,size_t n) { if (!p) return; @@ -236,11 +133,11 @@ namespace ecolab memAlloc->deallocate(p,n); return; } -#ifdef __SYCL_DEVICE_ONLY__ + //#ifdef __SYCL_DEVICE_ONLY__ printf("leaked %d bytes\n",n*sizeof(T)); -#else - sycl::free(p,syclQ()); -#endif +//#else +// sycl::free(p,syclQ()); +//#endif } bool operator==(const CellAllocator& x) const {return memAlloc==x.memAlloc;} }; @@ -310,8 +207,7 @@ namespace ecolab #endif for (size_t idx=0; idxtemplate as(); - cell.m_idx=idx; - f(cell); + f(cell,idx); } } @@ -328,8 +224,7 @@ namespace ecolab auto idx=i.get_global_linear_id(); if (idxsize()) { auto& cell=*(*this)[idx]->template as(); - cell.m_idx=idx; - f(cell); + f(cell,idx); } }); }); @@ -352,8 +247,7 @@ namespace ecolab auto idx=i.get_group_linear_id(); if (idxsize()) { auto& cell=*(*this)[idx]->template as(); - cell.m_idx=idx; - f(cell); + f(cell,idx); } }); }); @@ -362,12 +256,6 @@ namespace ecolab #endif } - void syncThreads() { -#ifdef SYCL_LANGUAGE_VERSION - syclQ().wait(); -#endif - } - template T max(T x, F f) { #ifdef SYCL_LANGUAGE_VERSION DeviceType r(x); @@ -390,6 +278,7 @@ namespace ecolab return x; #endif } + void syncThreads() {ecolab::syncThreads();} }; template 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> diag, val; /*row and columns of offdiagonal elements*/ array_ns::array> row, col; - sparse_mat(int s=0, int o=0, const A& falloc={}, - const A& 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 class A1> + sparse_mat(const sparse_mat& x) { rowsz=x.rowsz; colsz=x.colsz; diag=x.diag; val=x.val; row=x.row; col=x.col; } - void setAllocators(const A& ialloc, const A& falloc) { - diag.allocator(falloc); - val.allocator(falloc); - row.allocator(ialloc); - col.allocator(ialloc); - } /*matrix multiplication*/ template typename @@ -97,7 +91,7 @@ namespace ecolab assert(sigma>=0); assert(conn> rsize(diag.size(),row.allocator()); + array> rsize(diag.size(),row.allocator()); array> 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..a22eb90 --- /dev/null +++ b/include/sycl.h @@ -0,0 +1,116 @@ +/* + @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 +#ifdef SYCL_LANGUAGE_VERSION + +#include +#include "usmAlloc.h" // TODO - merge into this file. +#include "graphcode.h" + +// we use some experimental features available in the OneAPI compiler +#ifndef __INTEL_LLVM_COMPILER +#error "EcoLab requires OneAPI compiler for some experimental functions" +#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 groupBarrier() {sycl::group_barrier(syclGroup());} + + inline void syncThreads() {syclQ().wait_and_throw();} + + template + struct SyclType: public T + { + template SyclType(A... args): T(std::forward(args)...) {} + 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 <> + struct SyclType + { + size_t data; + operator size_t() const {return data;} + operator size_t&() {return data;} + size_t operator=(size_t x) {return data=x;} + 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 + struct SyclType + { + T* data; + operator T*() const {return data;} + operator T*&() {return data;} + T* operator=(T* x) {return data=x;} + 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 DeviceType + { + SyclType* const model; + public: + using element_type=M; + template + DeviceType(A... args): model(new SyclType(std::forward(args)...)) {} + DeviceType(const DeviceType& x): model(new SyclType(*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 + struct SyclQAllocator: public graphcode::Allocator + { + SyclQAllocator(): graphcode::Allocator(syclQ(), UA) {} + template struct rebind {using other=SyclQAllocator;}; + }; + +} +#else // !SYCL +namespace ecolab +{ + inline void groupBarrier() {} + inline void syncThreads() {} +} +#endif +#include "sycl.cd" +#endif diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc index bbc686a..dc78f1d 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 spatial_ecolab; CLASSDESC_ADD_GLOBAL(spatial_ecolab); CLASSDESC_PYTHON_MODULE(ecolab_model); @@ -52,8 +54,7 @@ void PanmicticModel::generate(unsigned niter) } /* Rounding function, randomly round up or down, in the range 0..INT_MAX */ -template -int EcolabPoint::ROUND(Float x) +int EcolabPoint::ROUND(Float x) { Float dum; const Float maxInt=Float(std::numeric_limits::max()-1); @@ -83,10 +84,9 @@ struct RoundArray namespace ecolab::array_ns {template struct is_expression>: public true_type {};} -template template -RoundArray> EcolabPoint::roundArray(const E& expr) -{return RoundArray>(*this,expr);} +RoundArray EcolabPoint::roundArray(const E& expr) +{return RoundArray(*this,expr);} template <> void setArray(array>& x, const array& y) {x=y;} @@ -123,41 +123,39 @@ void setArray(array>& x, const array -void EcolabPoint::generate(unsigned niter, const ModelData& model) +void EcolabPoint::generate(unsigned niter, const ModelData& model) { -#ifdef __SYCL_DEVICE_ONLY__ - Float* interactionResult=groupBuffer(density.size()); -#else - array interactionResult(density.size()); -#endif - for (unsigned step=0; step -#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(); - } -// // sequential/non-GPU version -// for (unsigned i=0; i(density.size()); +//#else +// array interactionResult(density.size()); +//#endif +// for (unsigned step=0; step +//#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(); +// } + // sequential/non-GPU version + for (unsigned i=0; i -unsigned EcolabPoint::nsp() const +unsigned EcolabPoint::nsp() const {return sum(density!=0);} array SpatialModel::nsp() const @@ -167,8 +165,7 @@ array SpatialModel::nsp() const return nsp; } -template -void EcolabPoint::condense(const array& mask, size_t mask_true) +void EcolabPoint::condense(const array& mask, size_t mask_true) { density = pack( density, mask, mask_true); } @@ -241,32 +238,22 @@ void PanmicticModel::mutate() void SpatialModel::mutate() { - - - - DeviceType>> mut_scale; -#ifdef SYCL_LANGUAGE_VERSION - using ArrayAlloc=CellBase::CellAllocator; - using NewSpAlloc=graphcode::Allocator>; - vector,NewSpAlloc> newSp - (size(),NewSpAlloc(syclQ(),sycl::usm::alloc::shared)); - mut_scale->allocator(graphcode::Allocator(syclQ(),sycl::usm::alloc::shared)); -#else - vector> newSp(size()); -#endif + DeviceType>> mut_scale(species.size()); + using Array=array>; + vector> newSp(size()); *mut_scale=sp_sep * repro_rate * mutation * int(tstep - last_mut_tstep); last_mut_tstep=tstep; - - groupedForAll([newSp=newSp.data(),mut_scale=&*mut_scale,this](EcolabCell& c) { + + groupedForAll([newSp=newSp.data(),mut_scale=&*mut_scale,this](EcolabCell& c,size_t i) { auto tmp{c.mutate(*mut_scale)}; #ifdef __SYCL_DEVICE_ONLY__ if (syclGroup().leader() && tmp.size()) #endif - newSp[c.idx()]=tmp; + newSp[i]=tmp; }); array new_sp; - DeviceType>> cell_ids; + DeviceType>> cell_ids; #ifdef SYCL_LANGUAGE_VERSION syncThreads(); #endif @@ -283,10 +270,13 @@ 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(); }); + if (new_sp.size()==0) return; + + cout<size()==0) return; // 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](EcolabCell& c,size_t) { c.density <<= (*cell_ids)==c.id; }); } -template template -array::template Allocator> -EcolabPoint::mutate(const E& mut_scale) +EcolabPoint::UnsignedArray EcolabPoint::mutate(const E& mut_scale) { /* calculate the number of mutants each species produces */ #ifdef __SYCL_DEVICE_ONLY__ @@ -343,8 +331,8 @@ EcolabPoint::mutate(const E& mut_scale) groupBarrier(); if (offsets[syclGroup().get_group_linear_range()]) array_ns::asg_minus_v(density.data(), nsp, speciations); - GroupLocal::template Allocator>> new_sp - (offsets[syclGroup().get_group_linear_range()], this->template allocator()); + GroupLocal new_sp(offsets[syclGroup().get_group_linear_range()], + density.allocator()); groupBarrier(); @@ -662,7 +650,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 +659,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 +671,15 @@ unsigned SpatialModel::migrate() const Float cap=1.0/maxNbrs; array capped_migration = merge(mm>cap,cap,mm); - hostForAll([&,this](EcolabCell& c) { + hostForAll([&,this](EcolabCell& c, size_t i) { /* loop over neighbours */ for (auto& n: c) { auto& nbr=*n->as(); Float salt=&c<&nbr? c.salt: nbr.salt; array m=capped_migration * (nbr.density-c.density); - delta[c.idx()]+=m;//*(1 + salt * (abs(m)=-delta[c.idx()])); + delta[i]+=m;//*(1 + salt * (abs(m)=-delta[i])); } }); @@ -705,14 +694,14 @@ unsigned SpatialModel::migrate() for (size_t i=0; ias(); - 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 +746,19 @@ 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(syclQ(),sycl::usm::alloc::shared)); - create.allocator(falloc); - repro_rate.allocator(falloc); - mutation.allocator(falloc); - migration.allocator(falloc); - interaction.setAllocators - (graphcode::Allocator(syclQ(),sycl::usm::alloc::shared),falloc); - } -#endif +//#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(syclQ(),sycl::usm::alloc::shared)); +// create.allocator(falloc); +// repro_rate.allocator(falloc); +// mutation.allocator(falloc); +// migration.allocator(falloc); +// interaction.setAllocators +// (graphcode::Allocator(syclQ(),sycl::usm::alloc::shared),falloc); +// } +//#endif if (!species.size()) { species=pcoord(nsp); @@ -780,25 +769,25 @@ void ModelData::makeConsistent(size_t nsp) if (!migration.size()) migration.resize(species.size(),0); } -void SpatialModel::setDensitiesShared() -{ -#ifdef SYCL_LANGUAGE_VERSION - groupedForAll([=,this](EcolabCell& c) { - c.memAlloc=sharedMemAlloc; - c.density.allocator(c.allocator()); - }); -#endif -} - -void SpatialModel::setDensitiesDevice() -{ -#ifdef SYCL_LANGUAGE_VERSION - groupedForAll([=,this](EcolabCell& c) { - c.memAlloc=deviceMemAlloc; - c.density.allocator(c.allocator()); - }); -#endif -} +//void SpatialModel::setDensitiesShared() +//{ +//#ifdef SYCL_LANGUAGE_VERSION +// groupedForAll([=,this](EcolabCell& c) { +// c.memAlloc=sharedMemAlloc; +// c.density.allocator(c.allocator()); +// }); +//#endif +//} +// +//void SpatialModel::setDensitiesDevice() +//{ +//#ifdef SYCL_LANGUAGE_VERSION +// groupedForAll([=,this](EcolabCell& c) { +// c.memAlloc=deviceMemAlloc; +// c.density.allocator(c.allocator()); +// }); +//#endif +//} void SpatialModel::makeConsistent() { @@ -808,13 +797,13 @@ 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()); - }); +// 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()); +// }); ModelData::makeConsistent(nsp); syncThreads(); } diff --git a/models/ecolab_model.h b/models/ecolab_model.h index a12434b..a02e8bf 100644 --- a/models/ecolab_model.h +++ b/models/ecolab_model.h @@ -38,11 +38,15 @@ struct ConnectionPlot: public Object struct ModelData { - using FAlloc=graphcode::Allocator; - using IAlloc=graphcode::Allocator; - array species; - array create, repro_rate, mutation, migration; - sparse_mat interaction; +#ifdef SYCL_LANGUAGE_VERSION + template using Allocator=HostSharedAllocator; +#else + template using Allocator=std::allocator; +#endif + + array> species; + array> create, repro_rate, mutation, migration; + sparse_mat interaction; sparse_mat_graph foodweb; unsigned long long tstep=0, last_mut_tstep=0, last_mig_tstep=0; //mutation parameters @@ -65,18 +69,21 @@ template void setArray(array&, const array&); template array getArray(const array&); /* ecolab cell */ -template -class EcolabPoint: public Exclude +class EcolabPoint { public: Float salt; /* random no. used for migration */ - template using Allocator=typename CellBase::template CellAllocator; - array> density{this->template allocator()}; +#ifdef SYCL_LANGUAGE_VERSION + template using Allocator=GlobalDeviceAllocator; +#else + template using Allocator=std::allocator; +#endif + array> density; void generate(unsigned niter, const ModelData&); void condense(const array& mask, size_t mask_true); - template - array> mutate(const E&); + using UnsignedArray=array>; + template 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); @@ -84,14 +91,7 @@ class EcolabPoint: public Exclude Exclude rand; // random number generator }; -// for the panmictic model, we need to use std::allocator -struct AllocatorBase -{ - template using CellAllocator=std::allocator; - template CellAllocator allocator() const {return CellAllocator();} -}; - -struct PanmicticModel: public ModelData, public EcolabPoint, public ecolab::Model +struct PanmicticModel: public ModelData, public EcolabPoint, public ecolab::Model { ConnectionPlot connectionPlot; void updateConnectionPlot() {connectionPlot.update(this->density,interaction);} @@ -109,10 +109,10 @@ struct PanmicticModel: public ModelData, public EcolabPoint, publ Float extinctionConnectivity() const; }; -struct EcolabCell: public EcolabPoint, public graphcode::Object +struct EcolabCell: public EcolabPoint, public graphcode::Object { - unsigned id=0; // stash the graphcode node id here - TODO - why can't we call graphcode::Object::id()? - array> delta{this->template allocator()}; + unsigned id=0; // stash the graphcode node id here so it is accessible from within the cell, rather than its pointer wrapper +// array> delta; }; class SpatialModel: public ModelData, public EcolabGraph, @@ -130,11 +130,11 @@ class SpatialModel: public ModelData, public EcolabGraph, return *objects[makeId(x,y)]; } /// on GPUs, set the - void setDensitiesShared(); - void setDensitiesDevice(); +// void setDensitiesShared(); +// void setDensitiesDevice(); array nsp() const; void makeConsistent(); - void seed(unsigned x) {forAll([=](EcolabCell& cell){cell.rand.seed(x);});} + void seed(unsigned x) {forAll([=](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..6ddc4ff 100644 --- a/models/spatial_ecolab.py +++ b/models/spatial_ecolab.py @@ -28,8 +28,8 @@ def randomList(num, min, max): ecolab.species(range(nsp)) -numX=8 -numY=8 +numX=12 +numY=12 ecolab.setGrid(numX,numY) ecolab.partitionObjects() 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;} }; } From 97109af97ed05fee1cd5d8ecfd0ea8af2e56ac17 Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Mon, 6 Jul 2026 12:33:24 +1000 Subject: [PATCH 02/28] bugfix: revert switch to standard seqential code in generate. --- models/ecolab_model.cc | 56 +++++++++++++++++++++--------------------- 1 file changed, 28 insertions(+), 28 deletions(-) diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc index dc78f1d..562100d 100644 --- a/models/ecolab_model.cc +++ b/models/ecolab_model.cc @@ -125,34 +125,34 @@ void setArray(array>& x, const array(density.size()); -//#else -// array interactionResult(density.size()); -//#endif -// for (unsigned step=0; step -//#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(); -// } - // sequential/non-GPU version - for (unsigned i=0; i(density.size()); +#else + array interactionResult(density.size()); +#endif + for (unsigned step=0; step +#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(); + } +// // sequential/non-GPU version +// for (unsigned i=0; i Date: Tue, 7 Jul 2026 13:34:43 +1000 Subject: [PATCH 03/28] bugfix: Fix SYCL version of EcolabPoint::mutate bugfix: Thread aware random number generator for SYCL chod: split out sycl support to separate header. --- include/DeviceAllocator.h | 2 +- include/arrays.h | 5 +- include/ecolab.h | 78 +-------------- include/non-sycl.h | 35 +++++++ include/sycl.h | 203 +++++++++++++++++++++++++++++--------- include/usmAlloc.h | 4 + models/ecolab_model.cc | 59 +++++------ models/ecolab_model.h | 7 +- models/spatial_ecolab.py | 10 +- 9 files changed, 234 insertions(+), 169 deletions(-) create mode 100644 include/non-sycl.h diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h index 6c46396..810b9bb 100644 --- a/include/DeviceAllocator.h +++ b/include/DeviceAllocator.h @@ -76,7 +76,7 @@ namespace ecolab } void deallocate(void* p, size_t size) { if (!p) return; - if (p>=memory && p=memory && p(p)-memory)>>order); return; } diff --git a/include/arrays.h b/include/arrays.h index 3dd2ac9..bd5037c 100644 --- a/include/arrays.h +++ b/include/arrays.h @@ -1480,6 +1480,7 @@ namespace ecolab ///allocate \a n variables of type \a T array_data *alloc(std::size_t n) { + if (n==0) return nullptr; T *p; array_data *r; //p = (char*)std::malloc((n-array_data::debug_display) * sizeof(T) + sizeof(array_data) + 16); @@ -2399,8 +2400,8 @@ namespace ecolab return o; } - template - std::ostream& operator<<(std::ostream& o, const array& x) + template + std::ostream& operator<<(std::ostream& o, const array& x) {return put(o,x);} /// ostream putter diff --git a/include/ecolab.h b/include/ecolab.h index 2a95682..6632aa1 100644 --- a/include/ecolab.h +++ b/include/ecolab.h @@ -250,7 +250,7 @@ namespace ecolab f(cell,idx); } }); - }); + }).wait_and_throw(); #else hostForAll(f); #endif @@ -272,7 +272,6 @@ namespace ecolab #endif for (size_t idx=0; idxtemplate as(); - cell.m_idx=idx; x=std::max(x,f(cell)); } return x; @@ -280,81 +279,6 @@ namespace ecolab } void syncThreads() {ecolab::syncThreads();} }; - - template - class GroupLocal - { -#ifdef __SYCL_DEVICE_ONLY__ - using BufferType=char[sizeof(T)]; - using LocalBufferType=decltype(sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup())); - LocalBufferType buffer=sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup()); - public: - template - GroupLocal(Args&&... args) { - if (syclGroup().leader()) - new (*buffer) T(std::forward(args)...); - sycl::group_barrier(syclGroup()); - } - ~GroupLocal() { - sycl::group_barrier(syclGroup()); - if (syclGroup().leader()) - (**this).~T(); - } - T& operator*() {return reinterpret_cast(**buffer);} -#else - T buffer; - public: - template - GroupLocal(Args&&... args): buffer(std::forward(args)...) {} - T& operator*() {return buffer;} -#endif - T* operator->() {return &**this;} - }; - -#ifdef __SYCL_DEVICE_ONLY__ - template - T* getGroupBuffer() - { - if constexpr (N<=M) - return *sycl::ext::oneapi::group_local_memory_for_overwrite(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 - T* groupBuffer(size_t n) { - // work out power of two >= n - switch (sizeof(size_t)*8-sycl::clz(n-1)) - { - case 0: return getGroupBuffer(); - case 1: return getGroupBuffer(); - case 2: return getGroupBuffer(); - case 3: return getGroupBuffer(); - case 4: return getGroupBuffer(); - case 5: return getGroupBuffer(); - case 6: return getGroupBuffer(); - case 7: return getGroupBuffer(); - case 8: return getGroupBuffer(); - case 9: return getGroupBuffer(); - case 10: return getGroupBuffer(); - case 11: return getGroupBuffer(); - case 12: return getGroupBuffer(); - case 13: return getGroupBuffer(); - case 14: return getGroupBuffer(); - case 15: return getGroupBuffer(); - case 16: return getGroupBuffer(); - 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..0a0eb89 --- /dev/null +++ b/include/non-sycl.h @@ -0,0 +1,35 @@ +/* + @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 struct SyclQAllocator: public std::allocator {}; + + inline void groupBarrier() {} + inline void syncThreads() {} + + /// Non SYCL version of SyclType + template + struct SyclType + { + T data; + template SyclType(A... args): data(std::forward(args)...) {} + operator T() const {return data;} + operator T&() {return data;} + }; + + template struct GroupLocal: public std::unique_ptr + { + }; + + template + struct SyclRandomEngine: public E {}; + +} +#endif diff --git a/include/sycl.h b/include/sycl.h index a22eb90..d743938 100644 --- a/include/sycl.h +++ b/include/sycl.h @@ -8,11 +8,13 @@ #ifndef ECOLAB_SYCL_H #define ECOLAB_SYCL_H +#include "usmAlloc.h" // TODO - merge into this file. +#include "graphcode.h" +#include +#include #ifdef SYCL_LANGUAGE_VERSION #include -#include "usmAlloc.h" // TODO - merge into this file. -#include "graphcode.h" // we use some experimental features available in the OneAPI compiler #ifndef __INTEL_LLVM_COMPILER @@ -34,62 +36,172 @@ namespace ecolab inline void groupBarrier() {sycl::group_barrier(syclGroup());} inline void syncThreads() {syclQ().wait_and_throw();} - + + /// SyclType is a pointer type, that when allocated via new is allocated from USM template struct SyclType: public T { - template SyclType(A... args): T(std::forward(args)...) {} + T data; + template SyclType(A... args): data(std::forward(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 <> - struct SyclType + template + struct SyclQAllocator: public graphcode::Allocator { - size_t data; - operator size_t() const {return data;} - operator size_t&() {return data;} - size_t operator=(size_t x) {return data=x;} - 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);} + SyclQAllocator(): graphcode::Allocator(syclQ(), UA) {} + template struct rebind {using other=SyclQAllocator;}; }; - + template - struct SyclType + class GroupLocal { - T* data; - operator T*() const {return data;} - operator T*&() {return data;} - T* operator=(T* x) {return data=x;} - 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);} +#ifdef __SYCL_DEVICE_ONLY__ + using BufferType=char[sizeof(T)]; + using LocalBufferType=decltype(sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup())); + LocalBufferType buffer=sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup()); + public: + template + GroupLocal(Args&&... args) { + if (syclGroup().leader()) + new (*buffer) T(std::forward(args)...); + sycl::group_barrier(syclGroup()); + } + ~GroupLocal() { + sycl::group_barrier(syclGroup()); + if (syclGroup().leader()) + (**this).~T(); + } + T& operator*() {return reinterpret_cast(**buffer);} +#else + T buffer; + public: + template + GroupLocal(Args&&... args): buffer(std::forward(args)...) {} + T& operator*() {return buffer;} +#endif + T* operator->() {return &**this;} }; + +#ifdef __SYCL_DEVICE_ONLY__ + template + T* getGroupBuffer() + { + if constexpr (N<=M) + return *sycl::ext::oneapi::group_local_memory_for_overwrite(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 + T* groupBuffer(size_t n) { + // work out power of two >= n + switch (sizeof(size_t)*8-sycl::clz(n-1)) + { + case 0: return getGroupBuffer(); + case 1: return getGroupBuffer(); + case 2: return getGroupBuffer(); + case 3: return getGroupBuffer(); + case 4: return getGroupBuffer(); + case 5: return getGroupBuffer(); + case 6: return getGroupBuffer(); + case 7: return getGroupBuffer(); + case 8: return getGroupBuffer(); + case 9: return getGroupBuffer(); + case 10: return getGroupBuffer(); + case 11: return getGroupBuffer(); + case 12: return getGroupBuffer(); + case 13: return getGroupBuffer(); + case 14: return getGroupBuffer(); + case 15: return getGroupBuffer(); + case 16: return getGroupBuffer(); + default: assert(false && "Group buffer limit exceeded"); return nullptr; + } + } +#endif + + // random numbers + template + struct SyclRandomEngine + { + std::vector> 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 class DeviceType { - SyclType* const model; + M* const model; public: using element_type=M; template - DeviceType(A... args): model(new SyclType(std::forward(args)...)) {} - DeviceType(const DeviceType& x): model(new SyclType(*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" +#ifdef SYCL_LANGUAGE_VERSION + DeviceType(A... args): model(sycl::malloc_shared(1,syclQ())) + {new (model) M(std::forward(args)...);} +#else + DeviceType(A... args): model(new M(std::forward(args)...)) {} #endif - ~DeviceType() {delete model;} -#if defined(__GNUC__) && !defined(__clang__) -#pragma GCC diagnostic pop + 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() { +#ifdef SYCL_LANGUAGE_VERSION + model->~M(); + sycl::free(model,syclQ()); +#else + delete model; #endif + } +//#if defined(__GNUC__) && !defined(__clang__) +//#pragma GCC diagnostic pop +//#endif M& operator*() {return *model;} const M& operator*() const {return *model;} M* operator->() {return model;} @@ -97,20 +209,13 @@ namespace ecolab operator bool() const {return true;} // always defined }; - template - struct SyclQAllocator: public graphcode::Allocator - { - SyclQAllocator(): graphcode::Allocator(syclQ(), UA) {} - template struct rebind {using other=SyclQAllocator;}; - }; - } -#else // !SYCL -namespace ecolab -{ - inline void groupBarrier() {} - inline void syncThreads() {} -} -#endif + + +#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 +#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 562100d..45e39ac 100644 --- a/models/ecolab_model.cc +++ b/models/ecolab_model.cc @@ -60,7 +60,7 @@ int EcolabPoint::ROUND(Float x) const Float maxInt=Float(std::numeric_limits::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()); + //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; } @@ -76,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; } }; @@ -239,11 +238,11 @@ void PanmicticModel::mutate() void SpatialModel::mutate() { DeviceType>> mut_scale(species.size()); - using Array=array>; - vector> newSp(size()); *mut_scale=sp_sep * repro_rate * mutation * int(tstep - last_mut_tstep); last_mut_tstep=tstep; - + + vector> newSp(size()); + groupedForAll([newSp=newSp.data(),mut_scale=&*mut_scale,this](EcolabCell& c,size_t i) { auto tmp{c.mutate(*mut_scale)}; #ifdef __SYCL_DEVICE_ONLY__ @@ -254,11 +253,7 @@ void SpatialModel::mutate() array new_sp; DeviceType>> cell_ids; -#ifdef SYCL_LANGUAGE_VERSION syncThreads(); -#endif - - // TODO - this is a kind of scan - can it be done on device? size_t j=0; @@ -271,12 +266,17 @@ void SpatialModel::mutate() // deallocate on device groupedForAll([newSp=newSp.data()](EcolabCell& c,size_t i) { - newSp[i].clear(); + #ifdef __SYCL_DEVICE_ONLY__ + if (syclGroup().leader()) +#endif + { + newSp[i].clear(); + assert(newSp[i].refCnt()==0); + } }); - + if (new_sp.size()==0) return; - cout<(syclGroup().get_group_linear_range()); + auto numWorkItems=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(syclGroup().get_group_linear_range()+1); + auto offsets=groupBuffer(nsp+1); + sycl::joint_exclusive_scan(syclGroup(),speciations,speciations+nsp,offsets,sycl::plus()); + groupBarrier(); 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]; - } + offsets[nsp]=offsets[nsp-1]+speciations[nsp-1]; groupBarrier(); - if (offsets[syclGroup().get_group_linear_range()]) + if (offsets[nsp]) array_ns::asg_minus_v(density.data(), nsp, speciations); - GroupLocal new_sp(offsets[syclGroup().get_group_linear_range()], - density.allocator()); - groupBarrier(); + GroupLocal new_sp(offsets[nsp], density.allocator()); - // gen_index - auto p=offsets[myId]; - auto n=new_sp->data(); + //gen_index array_ns::map(nsp, [&](size_t i) { - for (size_t j=0; jdata()+offsets[i]; + for (auto j=0; j speciations=roundArray(mut_scale * density); auto new_sp = gen_index(speciations); diff --git a/models/ecolab_model.h b/models/ecolab_model.h index a02e8bf..029f333 100644 --- a/models/ecolab_model.h +++ b/models/ecolab_model.h @@ -88,7 +88,12 @@ class EcolabPoint /// Rounding function, randomly round up or down, in the range 0..INT_MAX int ROUND(Float x); template RoundArray roundArray(const E& expr); +#ifdef SYCL_LANGUAGE_VERSION + Exclude> rand + {syclQ().get_device().get_info()}; // TODO make configurable? +#else Exclude rand; // random number generator +#endif }; struct PanmicticModel: public ModelData, public EcolabPoint, public ecolab::Model @@ -122,7 +127,7 @@ class SpatialModel: public ModelData, public EcolabGraph, 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); diff --git a/models/spatial_ecolab.py b/models/spatial_ecolab.py index 6ddc4ff..d1de78d 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): @@ -30,6 +30,8 @@ def randomList(num, min, max): numX=12 numY=12 +#numX=2 +#numY=2 ecolab.setGrid(numX,numY) ecolab.partitionObjects() @@ -89,7 +91,7 @@ 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) if myid()==0: From 9d683113fcc44a7fb721e0bf2ffa3e7be584f5a6 Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Tue, 7 Jul 2026 17:32:27 +1000 Subject: [PATCH 04/28] This commit crashes NUC: Build using DPCPP=1 DEBUG=1 The run using "env ONEAPI_DEVICE_SELECTOR=*:cpu python3 spatial_ecolab.py" --- include/DeviceAllocator.h | 49 ++++++++++++- include/arrays.h | 26 ++++--- include/sycl.h | 146 +++++++++++++++++++------------------- models/ecolab_model.cc | 104 ++++++++++++++++++++------- models/ecolab_model.h | 3 +- 5 files changed, 214 insertions(+), 114 deletions(-) diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h index 810b9bb..b8ee1fa 100644 --- a/include/DeviceAllocator.h +++ b/include/DeviceAllocator.h @@ -121,6 +121,53 @@ namespace ecolab HostSharedAllocator(): graphcode::Allocator(syclQ(), sycl::usm::alloc::shared) {} template struct rebind {using other=HostSharedAllocator;}; }; - + +#ifdef __SYCL_DEVICE_ONLY__ + template + class LocalAllocatorImpl + { + static_assert((size&(size-1))==0); // assert power of 2 + unsigned next=0; + char buffer[size-1]; + public: + void* allocate(size_t n) { + unsigned offs=next; + sycl::group_barrier(syclGroup()); + if (syclGroup().leader()) next+=n; + return buffer+offs; + } + }; + + template + class LocalAllocator + { + constexpr static unsigned size=32*1024; + using BufferType=char[sizeof(LocalAllocatorImpl)]; + using LocalBufferType=decltype(sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup())); + // pointer is same for all instantiations in this group + LocalBufferType buffer=sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup()); + LocalAllocatorImpl& impl() const + {return reinterpret_cast&>(**buffer);} + 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; + + LocalAllocator() { + if (syclGroup().leader()) + new (*buffer) LocalAllocatorImpl(); + sycl::group_barrier(syclGroup()); + } + // no need for destructor, as Impl has nothing to tear down + T* allocate(size_t n) {return reinterpret_cast(impl().allocate(n*sizeof(T)));} + void deallocate(T*,size_t) {} // cleaned up when group exits + template struct rebind {using other=LocalAllocator;}; + }; +#else + template using LocalAllocator=std::allocator; +#endif + } #endif diff --git a/include/arrays.h b/include/arrays.h index bd5037c..7285c9c 100644 --- a/include/arrays.h +++ b/include/arrays.h @@ -1568,33 +1568,31 @@ namespace ecolab template void asgV(const A& alloc, size_t size, const E& x) { -#ifdef __SYCL_DEVICE_ONLY__ - GroupLocal tmp(size,alloc); - array_ns::asg_v(tmp->dt->dt,size,x); - groupBarrier(); - if (syclGroup().leader()) swap(*tmp); -#else +//#ifdef __SYCL_DEVICE_ONLY__ +// GroupLocal tmp(size,alloc); +// array_ns::asg_v(tmp->dt->dt,size,x); +// groupBarrier(); +// if (syclGroup().leader()) swap(*tmp); +//#else array tmp(size,alloc); asg_v(tmp.data(),size,x); swap(tmp); -#endif + //#endif } void copy() //any nonconst method needs to call this { // 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* oldData=dt; 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); } } @@ -1642,7 +1640,7 @@ namespace ecolab } /// 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) { diff --git a/include/sycl.h b/include/sycl.h index d743938..28607ea 100644 --- a/include/sycl.h +++ b/include/sycl.h @@ -58,79 +58,79 @@ namespace ecolab template struct rebind {using other=SyclQAllocator;}; }; - template - class GroupLocal - { -#ifdef __SYCL_DEVICE_ONLY__ - using BufferType=char[sizeof(T)]; - using LocalBufferType=decltype(sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup())); - LocalBufferType buffer=sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup()); - public: - template - GroupLocal(Args&&... args) { - if (syclGroup().leader()) - new (*buffer) T(std::forward(args)...); - sycl::group_barrier(syclGroup()); - } - ~GroupLocal() { - sycl::group_barrier(syclGroup()); - if (syclGroup().leader()) - (**this).~T(); - } - T& operator*() {return reinterpret_cast(**buffer);} -#else - T buffer; - public: - template - GroupLocal(Args&&... args): buffer(std::forward(args)...) {} - T& operator*() {return buffer;} -#endif - T* operator->() {return &**this;} - }; - -#ifdef __SYCL_DEVICE_ONLY__ - template - T* getGroupBuffer() - { - if constexpr (N<=M) - return *sycl::ext::oneapi::group_local_memory_for_overwrite(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 - T* groupBuffer(size_t n) { - // work out power of two >= n - switch (sizeof(size_t)*8-sycl::clz(n-1)) - { - case 0: return getGroupBuffer(); - case 1: return getGroupBuffer(); - case 2: return getGroupBuffer(); - case 3: return getGroupBuffer(); - case 4: return getGroupBuffer(); - case 5: return getGroupBuffer(); - case 6: return getGroupBuffer(); - case 7: return getGroupBuffer(); - case 8: return getGroupBuffer(); - case 9: return getGroupBuffer(); - case 10: return getGroupBuffer(); - case 11: return getGroupBuffer(); - case 12: return getGroupBuffer(); - case 13: return getGroupBuffer(); - case 14: return getGroupBuffer(); - case 15: return getGroupBuffer(); - case 16: return getGroupBuffer(); - default: assert(false && "Group buffer limit exceeded"); return nullptr; - } - } -#endif +// template +// class GroupLocal +// { +//#ifdef __SYCL_DEVICE_ONLY__ +// using BufferType=char[sizeof(T)]; +// using LocalBufferType=decltype(sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup())); +// LocalBufferType buffer=sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup()); +// public: +// template +// GroupLocal(Args&&... args) { +// if (syclGroup().leader()) +// new (*buffer) T(std::forward(args)...); +// sycl::group_barrier(syclGroup()); +// } +// ~GroupLocal() { +// sycl::group_barrier(syclGroup()); +// if (syclGroup().leader()) +// (**this).~T(); +// } +// T& operator*() {return reinterpret_cast(**buffer);} +//#else +// T buffer; +// public: +// template +// GroupLocal(Args&&... args): buffer(std::forward(args)...) {} +// T& operator*() {return buffer;} +//#endif +// T* operator->() {return &**this;} +// }; +// +//#ifdef __SYCL_DEVICE_ONLY__ +// template +// T* getGroupBuffer() +// { +// if constexpr (N<=M) +// return *sycl::ext::oneapi::group_local_memory_for_overwrite(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 +// T* groupBuffer(size_t n) { +// // work out power of two >= n +// switch (sizeof(size_t)*8-sycl::clz(n-1)) +// { +// case 0: return getGroupBuffer(); +// case 1: return getGroupBuffer(); +// case 2: return getGroupBuffer(); +// case 3: return getGroupBuffer(); +// case 4: return getGroupBuffer(); +// case 5: return getGroupBuffer(); +// case 6: return getGroupBuffer(); +// case 7: return getGroupBuffer(); +// case 8: return getGroupBuffer(); +// case 9: return getGroupBuffer(); +// case 10: return getGroupBuffer(); +// case 11: return getGroupBuffer(); +// case 12: return getGroupBuffer(); +// case 13: return getGroupBuffer(); +// case 14: return getGroupBuffer(); +// case 15: return getGroupBuffer(); +// case 16: return getGroupBuffer(); +// default: assert(false && "Group buffer limit exceeded"); return nullptr; +// } +// } +//#endif // random numbers template diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc index 45e39ac..aac4e8b 100644 --- a/models/ecolab_model.cc +++ b/models/ecolab_model.cc @@ -58,8 +58,8 @@ int EcolabPoint::ROUND(Float x) { Float dum; const Float maxInt=Float(std::numeric_limits::max()-1); - if (x<0) x=0; - if (x>maxInt) x=maxInt; + 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; @@ -124,11 +124,12 @@ void setArray(array>& x, const array(density.size()); -#else - array interactionResult(density.size()); -#endif +//#ifdef __SYCL_DEVICE_ONLY__ +// //Float* interactionResult=groupBuffer(density.size()); +// array interactionResult(density.size()); +//#else + array> interactionResult(density.size()); + //#endif for (unsigned step=0; step=0)); // // sequential/non-GPU version // for (unsigned i=0; i>> 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; vector> newSp(size()); + hostForAll([](EcolabCell& c,size_t) { + assert(all(c.density>=0)); + }); + groupedForAll([newSp=newSp.data(),mut_scale=&*mut_scale,this](EcolabCell& c,size_t i) { + assert(all(c.density>=0)); auto tmp{c.mutate(*mut_scale)}; #ifdef __SYCL_DEVICE_ONLY__ if (syclGroup().leader() && tmp.size()) @@ -251,6 +259,10 @@ void SpatialModel::mutate() newSp[i]=tmp; }); + hostForAll([](EcolabCell& c,size_t) { + assert(all(c.density>=0)); + }); + array new_sp; DeviceType>> cell_ids; syncThreads(); @@ -299,42 +311,83 @@ void SpatialModel::mutate() // 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,size_t) { + assert(all(c.density>=0)); c.density <<= (*cell_ids)==c.id; }); } template -EcolabPoint::UnsignedArray EcolabPoint::mutate(const E& mut_scale) +EcolabPoint::LocalArray EcolabPoint::mutate(const E& mut_scale) { /* calculate the number of mutants each species produces */ #ifdef __SYCL_DEVICE_ONLY__ + array> speciations=roundArray(mut_scale * density); auto nsp=density.size(); - auto speciations=groupBuffer(nsp); - groupBarrier(); - asg_v(speciations, nsp, roundArray(mut_scale * density)); - - auto numWorkItems=syclGroup().get_group_linear_range(); - auto myId=syclGroup().get_local_linear_id(); - - auto offsets=groupBuffer(nsp+1); - sycl::joint_exclusive_scan(syclGroup(),speciations,speciations+nsp,offsets,sycl::plus()); + // auto new_sp = gen_index(speciations); + array> offsets(density.size()+1); + sycl::joint_exclusive_scan(syclGroup(),speciations.data(),speciations.data()+nsp, + offsets.data(),sycl::plus()); groupBarrier(); if (syclGroup().leader()) offsets[nsp]=offsets[nsp-1]+speciations[nsp-1]; - groupBarrier(); - if (offsets[nsp]) - array_ns::asg_minus_v(density.data(), nsp, speciations); - GroupLocal new_sp(offsets[nsp], density.allocator()); + if (offsets[nsp]==0) return {}; + + density-=speciations; - //gen_index + vector> new_sp{offsets[nsp]}; array_ns::map(nsp, [&](size_t i) { - auto p=new_sp->data()+offsets[i]; + auto p=new_sp[0].data()+offsets[i]; for (auto j=0; j(nsp); +// groupBarrier(); +// //asg_v(speciations, nsp, roundArray(mut_scale * density)); +// array_ns::map(nsp,[&](size_t i){ +// auto r=roundArray(mut_scale * density)[i]; +// speciations[i]=r; +// if (density[i]==0 && speciations[i]!=0) printf("aa %u %g %d\n",i,(mut_scale * density)[i],r); +// assert(density[i]>0 || r==0); +// }); +// groupBarrier(); +// array_ns::map(nsp, [&](size_t i){ +// if (density[i](nsp+1); +// sycl::joint_exclusive_scan(syclGroup(),speciations,speciations+nsp,offsets,sycl::plus()); +// array_ns::map(nsp, [&](size_t i){ +// if (density[i]=0)); +// +// groupBarrier(); +// GroupLocal new_sp(offsets[nsp], density.allocator()); +// +// //gen_index +// array_ns::map(nsp, [&](size_t i) { +// auto p=new_sp->data()+offsets[i]; +// for (auto j=0; j speciations=roundArray(mut_scale * density); @@ -661,6 +714,7 @@ unsigned SpatialModel::migrate() array capped_migration = merge(mm>cap,cap,mm); hostForAll([&,this](EcolabCell& c, size_t i) { + assert(all(c.density>=0)); /* loop over neighbours */ for (auto& n: c) { diff --git a/models/ecolab_model.h b/models/ecolab_model.h index 029f333..cf54d8b 100644 --- a/models/ecolab_model.h +++ b/models/ecolab_model.h @@ -83,7 +83,8 @@ class EcolabPoint void generate(unsigned niter, const ModelData&); void condense(const array& mask, size_t mask_true); using UnsignedArray=array>; - template UnsignedArray mutate(const E&); + using LocalArray=array>; + template 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); From 113482d902b467f06cd040a7468ad543320f8f14 Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Thu, 9 Jul 2026 18:25:36 +1000 Subject: [PATCH 05/28] Working on 12x12 on GPU. Need to check non-SYCL build. --- include/DeviceAllocator.h | 62 ++++++++++++----------- include/arrays.h | 11 ++--- include/ecolab.h | 60 +++++++++++++++++++---- include/sycl.h | 82 ------------------------------- models/ecolab_model.cc | 100 +++----------------------------------- models/spatial_ecolab.py | 3 -- 6 files changed, 95 insertions(+), 223 deletions(-) diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h index b8ee1fa..e6a6b3e 100644 --- a/include/DeviceAllocator.h +++ b/include/DeviceAllocator.h @@ -18,6 +18,11 @@ namespace ecolab /// memory allocated to each order, total memory allocated on device=poolSize*(maxOrder-minOrder) constexpr unsigned poolSize=32*1024*1024; + void setFatalError() { + *sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup())=~0U; + sycl::group_barrier(syclGroup()); + } + // Lock free circular buffer queue for SYCL template class Queue @@ -47,7 +52,7 @@ namespace ecolab template <> class DeviceAllocator { public: void init(sycl::queue& q) {} - void* allocate(size_t sz) {return nullptr;} + void* allocate(size_t sz) {setFatalError(); return nullptr;} void deallocate(void* p, size_t) {sycl::ext::oneapi::experimental::printf("%x leaked on device\n",p);} }; @@ -122,32 +127,23 @@ namespace ecolab template struct rebind {using other=HostSharedAllocator;}; }; + constexpr static unsigned LocalAllocatorSize=30*1024; // 32KiB = half typical local storage #ifdef __SYCL_DEVICE_ONLY__ - template - class LocalAllocatorImpl - { - static_assert((size&(size-1))==0); // assert power of 2 - unsigned next=0; - char buffer[size-1]; - public: - void* allocate(size_t n) { - unsigned offs=next; - sycl::group_barrier(syclGroup()); - if (syclGroup().leader()) next+=n; - return buffer+offs; - } - }; - + /** + 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 LocalAllocator { - constexpr static unsigned size=32*1024; - using BufferType=char[sizeof(LocalAllocatorImpl)]; - using LocalBufferType=decltype(sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup())); - // pointer is same for all instantiations in this group - LocalBufferType buffer=sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup()); - LocalAllocatorImpl& impl() const - {return reinterpret_cast&>(**buffer);} +// private: +// using LocalIndex=decltype(sycl::ext::oneapi::group_local_memory(syclGroup(),0)); +// //LocalIndex next=sycl::ext::oneapi::group_local_memory(syclGroup(),0); +// using BufferType=char[size]; +// using LocalBufferType=decltype(sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup())); +// // pointer is same for all instantiations in this group +// //LocalBufferType buffer=sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup()); public: using value_type=T; using pointer=T*; @@ -155,13 +151,23 @@ namespace ecolab using difference_type=std::ptrdiff_t; using propagate_on_container_move_assignment=std::true_type; - LocalAllocator() { - if (syclGroup().leader()) - new (*buffer) LocalAllocatorImpl(); + // no need for destructor, as Impl has nothing to tear down + __attribute__((always_inline)) T* allocate(size_t n) { + unsigned& next=*sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup()); + unsigned offs=next; + if (next==~0U) return nullptr; + if (offs+n*sizeof(T)>LocalAllocatorSize) + { + setFatalError(); + return nullptr; + } + //assert(offs+n*sizeof(T)<=LocalAllocatorSize); + sycl::group_barrier(syclGroup()); + if (syclGroup().leader()) next+=n*sizeof(T); sycl::group_barrier(syclGroup()); + return reinterpret_cast + ((*sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup()))+offs); } - // no need for destructor, as Impl has nothing to tear down - T* allocate(size_t n) {return reinterpret_cast(impl().allocate(n*sizeof(T)));} void deallocate(T*,size_t) {} // cleaned up when group exits template struct rebind {using other=LocalAllocator;}; }; diff --git a/include/arrays.h b/include/arrays.h index 7285c9c..86665d1 100644 --- a/include/arrays.h +++ b/include/arrays.h @@ -1477,7 +1477,10 @@ 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 *alloc(std::size_t n) { if (n==0) return nullptr; @@ -1495,9 +1498,6 @@ namespace ecolab #endif 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 @@ -1615,9 +1615,6 @@ 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()++; } diff --git a/include/ecolab.h b/include/ecolab.h index 6632aa1..d2600f9 100644 --- a/include/ecolab.h +++ b/include/ecolab.h @@ -240,17 +240,57 @@ namespace ecolab template 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(); - 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 (idxsize()) { - auto& cell=*(*this)[idx]->template as(); - f(cell,idx); - } + auto dev=syclQ().get_device(); + // 1. Max threads per single work-group + size_t max_wg_size = dev.get_info(); + + // 2. Hardware SIMD width (Sub-group size, usually 8, 16, or 32 on Intel) + std::vector sg_sizes = dev.get_info(); + 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(); + + // 4. Maximum number of compute units (Execution units / DSS count) + uint32_t max_compute_units = dev.get_info(); + + size_t workGroupSize=native_sg_size; +// 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< fatalError(false); + for (size_t cellStart=0; cellStartsize() && !*fatalError; 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) { + auto next=sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup()); + if (syclGroup().leader()) {*next = 0;} // reset local memory allocation + sycl::group_barrier(syclGroup()); + //sycl::ext::oneapi::group_local_memory_for_overwrite(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(); + f(cell,idx); + } + // flag fatal error to throw afterwards. + if (syclGroup().leader() && *next==~0U) *fatalError=true; + }); }); - }).wait_and_throw(); + syclQ().wait_and_throw(); + if (*fatalError) + throw std::runtime_error("Local Allocator Exhausted"); #else hostForAll(f); #endif diff --git a/include/sycl.h b/include/sycl.h index 28607ea..79ba281 100644 --- a/include/sycl.h +++ b/include/sycl.h @@ -58,80 +58,6 @@ namespace ecolab template struct rebind {using other=SyclQAllocator;}; }; -// template -// class GroupLocal -// { -//#ifdef __SYCL_DEVICE_ONLY__ -// using BufferType=char[sizeof(T)]; -// using LocalBufferType=decltype(sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup())); -// LocalBufferType buffer=sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup()); -// public: -// template -// GroupLocal(Args&&... args) { -// if (syclGroup().leader()) -// new (*buffer) T(std::forward(args)...); -// sycl::group_barrier(syclGroup()); -// } -// ~GroupLocal() { -// sycl::group_barrier(syclGroup()); -// if (syclGroup().leader()) -// (**this).~T(); -// } -// T& operator*() {return reinterpret_cast(**buffer);} -//#else -// T buffer; -// public: -// template -// GroupLocal(Args&&... args): buffer(std::forward(args)...) {} -// T& operator*() {return buffer;} -//#endif -// T* operator->() {return &**this;} -// }; -// -//#ifdef __SYCL_DEVICE_ONLY__ -// template -// T* getGroupBuffer() -// { -// if constexpr (N<=M) -// return *sycl::ext::oneapi::group_local_memory_for_overwrite(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 -// T* groupBuffer(size_t n) { -// // work out power of two >= n -// switch (sizeof(size_t)*8-sycl::clz(n-1)) -// { -// case 0: return getGroupBuffer(); -// case 1: return getGroupBuffer(); -// case 2: return getGroupBuffer(); -// case 3: return getGroupBuffer(); -// case 4: return getGroupBuffer(); -// case 5: return getGroupBuffer(); -// case 6: return getGroupBuffer(); -// case 7: return getGroupBuffer(); -// case 8: return getGroupBuffer(); -// case 9: return getGroupBuffer(); -// case 10: return getGroupBuffer(); -// case 11: return getGroupBuffer(); -// case 12: return getGroupBuffer(); -// case 13: return getGroupBuffer(); -// case 14: return getGroupBuffer(); -// case 15: return getGroupBuffer(); -// case 16: return getGroupBuffer(); -// default: assert(false && "Group buffer limit exceeded"); return nullptr; -// } -// } -//#endif - // random numbers template struct SyclRandomEngine @@ -186,11 +112,6 @@ namespace ecolab DeviceType(A... args): model(new M(std::forward(args)...)) {} #endif 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() { #ifdef SYCL_LANGUAGE_VERSION model->~M(); @@ -199,9 +120,6 @@ namespace ecolab delete model; #endif } -//#if defined(__GNUC__) && !defined(__clang__) -//#pragma GCC diagnostic pop -//#endif M& operator*() {return *model;} const M& operator*() const {return *model;} M* operator->() {return model;} diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc index aac4e8b..a7c6fe0 100644 --- a/models/ecolab_model.cc +++ b/models/ecolab_model.cc @@ -321,10 +321,10 @@ EcolabPoint::LocalArray EcolabPoint::mutate(const E& mut_scale) { /* calculate the number of mutants each species produces */ #ifdef __SYCL_DEVICE_ONLY__ - array> speciations=roundArray(mut_scale * density); + LocalArray speciations=roundArray(mut_scale * density); auto nsp=density.size(); // auto new_sp = gen_index(speciations); - array> offsets(density.size()+1); + LocalArray offsets(nsp+1); sycl::joint_exclusive_scan(syclGroup(),speciations.data(),speciations.data()+nsp, offsets.data(),sycl::plus()); groupBarrier(); @@ -332,63 +332,17 @@ EcolabPoint::LocalArray EcolabPoint::mutate(const E& mut_scale) offsets[nsp]=offsets[nsp-1]+speciations[nsp-1]; if (offsets[nsp]==0) return {}; - + density-=speciations; - vector> new_sp{offsets[nsp]}; - array_ns::map(nsp, [&](size_t i) { - auto p=new_sp[0].data()+offsets[i]; + LocalArray new_sp(offsets[nsp]); + array_ns::map(nsp, [offsets=offsets.data(),new_sp=new_sp.data()](size_t i) { + auto p=new_sp+offsets[i]; for (auto j=0; j(nsp); -// groupBarrier(); -// //asg_v(speciations, nsp, roundArray(mut_scale * density)); -// array_ns::map(nsp,[&](size_t i){ -// auto r=roundArray(mut_scale * density)[i]; -// speciations[i]=r; -// if (density[i]==0 && speciations[i]!=0) printf("aa %u %g %d\n",i,(mut_scale * density)[i],r); -// assert(density[i]>0 || r==0); -// }); -// groupBarrier(); -// array_ns::map(nsp, [&](size_t i){ -// if (density[i](nsp+1); -// sycl::joint_exclusive_scan(syclGroup(),speciations,speciations+nsp,offsets,sycl::plus()); -// array_ns::map(nsp, [&](size_t i){ -// if (density[i]=0)); -// -// groupBarrier(); -// GroupLocal new_sp(offsets[nsp], density.allocator()); -// -// //gen_index -// array_ns::map(nsp, [&](size_t i) { -// auto p=new_sp->data()+offsets[i]; -// for (auto j=0; j speciations=roundArray(mut_scale * density); auto new_sp = gen_index(speciations); @@ -789,19 +743,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(syclQ(),sycl::usm::alloc::shared)); -// create.allocator(falloc); -// repro_rate.allocator(falloc); -// mutation.allocator(falloc); -// migration.allocator(falloc); -// interaction.setAllocators -// (graphcode::Allocator(syclQ(),sycl::usm::alloc::shared),falloc); -// } -//#endif if (!species.size()) { species=pcoord(nsp); @@ -812,26 +753,6 @@ void ModelData::makeConsistent(size_t nsp) if (!migration.size()) migration.resize(species.size(),0); } -//void SpatialModel::setDensitiesShared() -//{ -//#ifdef SYCL_LANGUAGE_VERSION -// groupedForAll([=,this](EcolabCell& c) { -// c.memAlloc=sharedMemAlloc; -// c.density.allocator(c.allocator()); -// }); -//#endif -//} -// -//void SpatialModel::setDensitiesDevice() -//{ -//#ifdef SYCL_LANGUAGE_VERSION -// groupedForAll([=,this](EcolabCell& c) { -// c.memAlloc=deviceMemAlloc; -// c.density.allocator(c.allocator()); -// }); -//#endif -//} - void SpatialModel::makeConsistent() { // all cells must have same number of species. Pack out with zero density if necessary @@ -840,13 +761,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()); -// }); ModelData::makeConsistent(nsp); syncThreads(); } diff --git a/models/spatial_ecolab.py b/models/spatial_ecolab.py index d1de78d..fd05586 100644 --- a/models/spatial_ecolab.py +++ b/models/spatial_ecolab.py @@ -61,7 +61,6 @@ def randomList(num, min, max): extinctions=0 migrations=0 def stepImpl(): - #ecolab.setDensitiesDevice() ecolab.generate(100) ecolab.mutate() @@ -75,8 +74,6 @@ def stepImpl(): migrations+=ecolab.migrate() extinctions+=ecolab.condense() #ecolab.syncThreads() - #print(ecolab.nsp()()) - #ecolab.setDensitiesShared() #ecolab.gather() print(ecolab.nsp()()) From e02ad4ef3f82ac0bc4d86cab44b5874a57265fc8 Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Fri, 10 Jul 2026 13:30:21 +1000 Subject: [PATCH 06/28] bugfix: LocalAllocator working now - 12x12. --- include/DeviceAllocator.h | 51 +++++++++++++++++++++++++-------------- include/arrays.h | 3 ++- include/ecolab.h | 13 ++++++---- models/ecolab_model.cc | 6 +---- models/spatial_ecolab.py | 4 +-- 5 files changed, 46 insertions(+), 31 deletions(-) diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h index e6a6b3e..34fb97f 100644 --- a/include/DeviceAllocator.h +++ b/include/DeviceAllocator.h @@ -18,8 +18,14 @@ namespace ecolab /// memory allocated to each order, total memory allocated on device=poolSize*(maxOrder-minOrder) constexpr unsigned poolSize=32*1024*1024; + struct FatalErrorFlag + { + bool flag; + }; + void setFatalError() { - *sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup())=~0U; + printf("setting fatal error\n"); + sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup())->flag=true; sycl::group_barrier(syclGroup()); } @@ -52,7 +58,10 @@ namespace ecolab template <> class DeviceAllocator { public: void init(sycl::queue& q) {} - void* allocate(size_t sz) {setFatalError(); return nullptr;} + void* allocate(size_t sz) { + setFatalError(); + return nullptr; + } void deallocate(void* p, size_t) {sycl::ext::oneapi::experimental::printf("%x leaked on device\n",p);} }; @@ -128,7 +137,22 @@ namespace ecolab }; 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 avoid multiply defined symbols + */ + inline __attribute__((noinline)) LocalAllocatorBuffer& localAllocatorBuffer() + {return *sycl::ext::oneapi::group_local_memory_for_overwrite(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 @@ -137,13 +161,6 @@ namespace ecolab template class LocalAllocator { -// private: -// using LocalIndex=decltype(sycl::ext::oneapi::group_local_memory(syclGroup(),0)); -// //LocalIndex next=sycl::ext::oneapi::group_local_memory(syclGroup(),0); -// using BufferType=char[size]; -// using LocalBufferType=decltype(sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup())); -// // pointer is same for all instantiations in this group -// //LocalBufferType buffer=sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup()); public: using value_type=T; using pointer=T*; @@ -152,21 +169,19 @@ namespace ecolab using propagate_on_container_move_assignment=std::true_type; // no need for destructor, as Impl has nothing to tear down - __attribute__((always_inline)) T* allocate(size_t n) { - unsigned& next=*sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup()); - unsigned offs=next; - if (next==~0U) return nullptr; + T* allocate(size_t n) { + auto& b=localAllocatorBuffer(); + unsigned offs=b.next; if (offs+n*sizeof(T)>LocalAllocatorSize) { - setFatalError(); + setFatalError(); return nullptr; } - //assert(offs+n*sizeof(T)<=LocalAllocatorSize); sycl::group_barrier(syclGroup()); - if (syclGroup().leader()) next+=n*sizeof(T); + if (syclGroup().leader()) b.next+=n*sizeof(T); sycl::group_barrier(syclGroup()); - return reinterpret_cast - ((*sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup()))+offs); + char* alloc=b.buffer+offs; + return reinterpret_cast(alloc); } void deallocate(T*,size_t) {} // cleaned up when group exits template struct rebind {using other=LocalAllocator;}; diff --git a/include/arrays.h b/include/arrays.h index 86665d1..65ea4f9 100644 --- a/include/arrays.h +++ b/include/arrays.h @@ -1494,7 +1494,8 @@ namespace ecolab if (!p) { #ifdef __SYCL_DEVICE_ONLY__ - printf("failed to allocate %d bytes in array\n",sizeof(T)*n); + printf("failed to allocate %d bytes in array on group %u, thread %u\n", + sizeof(T)*n,syclGroup().get_group_linear_id(),syclGroup().get_local_linear_id()); #endif return nullptr; // SYCL allocator returns nullptr if not initialised } diff --git a/include/ecolab.h b/include/ecolab.h index d2600f9..cdd3d0b 100644 --- a/include/ecolab.h +++ b/include/ecolab.h @@ -267,16 +267,19 @@ namespace ecolab size_t num_work_groups = max_compute_units * wg_per_compute_unit; num_work_groups=std::min(num_work_groups,this->size()); - std::cout< fatalError(false); for (size_t cellStart=0; cellStartsize() && !*fatalError; 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) { - auto next=sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup()); - if (syclGroup().leader()) {*next = 0;} // reset local memory allocation - sycl::group_barrier(syclGroup()); + //auto next=sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup()); + auto fatalErrorFlag=sycl::ext::oneapi::group_local_memory(syclGroup(),false); //sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup()); + 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); @@ -285,7 +288,7 @@ namespace ecolab f(cell,idx); } // flag fatal error to throw afterwards. - if (syclGroup().leader() && *next==~0U) *fatalError=true; + if (syclGroup().leader() && fatalErrorFlag->flag) *fatalError=true; }); }); syclQ().wait_and_throw(); diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc index a7c6fe0..c0df504 100644 --- a/models/ecolab_model.cc +++ b/models/ecolab_model.cc @@ -124,12 +124,7 @@ void setArray(array>& x, const array(density.size()); -// array interactionResult(density.size()); -//#else array> interactionResult(density.size()); - //#endif for (unsigned step=0; step Date: Fri, 10 Jul 2026 14:21:24 +1000 Subject: [PATCH 07/28] bugfix: Fix non-sycl build --- include/DeviceAllocator.h | 2 -- include/sycl.h | 4 ++++ 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h index 34fb97f..253b67f 100644 --- a/include/DeviceAllocator.h +++ b/include/DeviceAllocator.h @@ -186,8 +186,6 @@ namespace ecolab void deallocate(T*,size_t) {} // cleaned up when group exits template struct rebind {using other=LocalAllocator;}; }; -#else - template using LocalAllocator=std::allocator; #endif } diff --git a/include/sycl.h b/include/sycl.h index 79ba281..394e17a 100644 --- a/include/sycl.h +++ b/include/sycl.h @@ -127,6 +127,10 @@ namespace ecolab operator bool() const {return true;} // always defined }; +#ifndef __SYCL_DEVICE_ONLY__ + template using LocalAllocator=std::allocator; +#endif + } From a2e4571a9d6786817912d4971d2fd651de248f74 Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Fri, 10 Jul 2026 15:19:35 +1000 Subject: [PATCH 08/28] chore: deadcode removal bugfix: fatalErrorFlag (SYCL) --- graphcode | 2 +- include/DeviceAllocator.h | 10 ++--- include/ecolab.h | 91 +-------------------------------------- models/ecolab_model.cc | 39 ++--------------- models/ecolab_model.h | 4 -- 5 files changed, 11 insertions(+), 135 deletions(-) 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 index 253b67f..eb7d868 100644 --- a/include/DeviceAllocator.h +++ b/include/DeviceAllocator.h @@ -23,10 +23,8 @@ namespace ecolab bool flag; }; - void setFatalError() { - printf("setting fatal error\n"); - sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup())->flag=true; - sycl::group_barrier(syclGroup()); + inline __attribute__((noinline)) bool& fatalErrorFlag() { + return sycl::ext::oneapi::group_local_memory(syclGroup(),false)->flag; } // Lock free circular buffer queue for SYCL @@ -59,7 +57,7 @@ namespace ecolab public: void init(sycl::queue& q) {} void* allocate(size_t sz) { - setFatalError(); + fatalErrorFlag()=true; return nullptr; } void deallocate(void* p, size_t) {sycl::ext::oneapi::experimental::printf("%x leaked on device\n",p);} @@ -174,7 +172,7 @@ namespace ecolab unsigned offs=b.next; if (offs+n*sizeof(T)>LocalAllocatorSize) { - setFatalError(); + fatalErrorFlag()=true; return nullptr; } sycl::group_barrier(syclGroup()); diff --git a/include/ecolab.h b/include/ecolab.h index cdd3d0b..3d22f28 100644 --- a/include/ecolab.h +++ b/include/ecolab.h @@ -102,93 +102,9 @@ namespace ecolab { 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 CellAllocator - { - public: - MemAllocator* memAlloc=nullptr; - template friend class Allocator; - CellAllocator()=default; - CellAllocator(MemAllocator* memAlloc): memAlloc(memAlloc) {} - template CellAllocator(const CellAllocator& x): memAlloc(x.memAlloc) {} - T* allocate(size_t sz) { - //#ifdef __SYCL_DEVICE_ONLY__ - if (memAlloc) { - auto r=reinterpret_cast(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(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 CellAllocator allocator() const { - return CellAllocator(memAlloc); - } -#else //!SYCL - template using CellAllocator=std::allocator; - template CellAllocator allocator() const {return CellAllocator();} -#endif - }; - - template - void printAllocator(const char* prefix, const T&) {} - - template - void printAllocator(const char* prefix, const CellBase::CellAllocator& 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(1, syclQ()); - new(sharedMemAlloc) MemAllocator; - sharedMemAlloc->init(syclQ()); -// deviceMemAlloc=sycl::malloc_shared(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 struct EcolabGraph: - public Exclude, public graphcode::Graph + template struct EcolabGraph: public graphcode::Graph { #ifdef SYCL_LANGUAGE_VERSION EcolabGraph(): graphcode::Graph @@ -273,9 +189,6 @@ namespace ecolab for (size_t cellStart=0; cellStartsize() && !*fatalError; 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) { - //auto next=sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup()); - auto fatalErrorFlag=sycl::ext::oneapi::group_local_memory(syclGroup(),false); - //sycl::ext::oneapi::group_local_memory_for_overwrite(syclGroup()); if (syclGroup().leader()) { localAllocatorBuffer().next = 0; } // reset local memory allocation @@ -288,7 +201,7 @@ namespace ecolab f(cell,idx); } // flag fatal error to throw afterwards. - if (syclGroup().leader() && fatalErrorFlag->flag) *fatalError=true; + *fatalError=fatalErrorFlag(); }); }); syclQ().wait_and_throw(); diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc index c0df504..2202343 100644 --- a/models/ecolab_model.cc +++ b/models/ecolab_model.cc @@ -87,41 +87,10 @@ template RoundArray EcolabPoint::roundArray(const E& expr) {return RoundArray(*this,expr);} -template <> void setArray(array>& x, const array& y) -{x=y;} -template <> array getArray(const array>& x) {return x;} - -#ifdef SYCL_LANGUAGE_VERSION -template <> -array getArray(const array>& x) -{ - DeviceType size; - DeviceType xData; - syclQ().single_task([size=&*size,xData=&*xData,x=&x](){ - *size=x->size(); - *xData=x->data(); - }).wait(); - array r(*size); - syclQ().copy(*xData,r.data(),*size); - return r; -} - -template <> -void setArray(array>& x, const array& y) -{ - auto size=y.size(); - DeviceType xData; - syclQ().single_task([size,x=&x,xData=&*xData](){ - //array> tmp(size,this->template allocator()); - //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 <> void setArray(array>& x, const array& y) +//{x=y;} +//template <> array getArray(const array>& x) {return x;} +// void EcolabPoint::generate(unsigned niter, const ModelData& model) { array> interactionResult(density.size()); diff --git a/models/ecolab_model.h b/models/ecolab_model.h index cf54d8b..1f6d627 100644 --- a/models/ecolab_model.h +++ b/models/ecolab_model.h @@ -118,7 +118,6 @@ struct PanmicticModel: public ModelData, public EcolabPoint, public ecolab::Mode struct EcolabCell: public EcolabPoint, public graphcode::Object { unsigned id=0; // stash the graphcode node id here so it is accessible from within the cell, rather than its pointer wrapper -// array> delta; }; class SpatialModel: public ModelData, public EcolabGraph, @@ -135,9 +134,6 @@ class SpatialModel: public ModelData, public EcolabGraph, EcolabCell& cell(size_t x, size_t y) { return *objects[makeId(x,y)]; } - /// on GPUs, set the -// void setDensitiesShared(); -// void setDensitiesDevice(); array nsp() const; void makeConsistent(); void seed(unsigned x) {forAll([=](EcolabCell& cell,size_t){cell.rand.seed(x);});} From 3d8753a6eea58b286bce507d7545e9392a92e420 Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Tue, 14 Jul 2026 11:59:08 +1000 Subject: [PATCH 09/28] bugfix: Increase working group size to 256. arrays assignment from expression should allocate on leader, the copy in parallel --- include/Makefile.config | 17 +++++++++++++++++ include/arrays.h | 15 +++++---------- include/ecolab.h | 2 +- models/ecolab_model.cc | 14 ++++++++------ models/spatial_ecolab.py | 2 +- 5 files changed, 32 insertions(+), 18 deletions(-) create mode 100644 include/Makefile.config diff --git a/include/Makefile.config b/include/Makefile.config new file mode 100644 index 0000000..6f8edc1 --- /dev/null +++ b/include/Makefile.config @@ -0,0 +1,17 @@ +TCL=1 +TK=1 +ZLIB=1 +READLINE=1 +XDR=1 +UNURAN=1 +GNUSL=1 +CAIRO=1 +BDB=1 +PANGO=1 +GDBM_COMPAT=1 +MPI= +PARALLEL= +OPENMP= +GCC= +NOGUI= +AQUA= diff --git a/include/arrays.h b/include/arrays.h index 65ea4f9..90216a8 100644 --- a/include/arrays.h +++ b/include/arrays.h @@ -1569,16 +1569,11 @@ namespace ecolab template void asgV(const A& alloc, size_t size, const E& x) { -//#ifdef __SYCL_DEVICE_ONLY__ -// GroupLocal tmp(size,alloc); -// array_ns::asg_v(tmp->dt->dt,size,x); -// groupBarrier(); -// if (syclGroup().leader()) swap(*tmp); -//#else - array tmp(size,alloc); - asg_v(tmp.data(),size,x); - swap(tmp); - //#endif +#ifdef __SYCL_DEVICE_ONLY__ + if (syclGroup().leader()) +#endif + resize(size); + asg_v(data(),size,x); } void copy() //any nonconst method needs to call this diff --git a/include/ecolab.h b/include/ecolab.h index 3d22f28..92a29e5 100644 --- a/include/ecolab.h +++ b/include/ecolab.h @@ -170,7 +170,7 @@ namespace ecolab // 4. Maximum number of compute units (Execution units / DSS count) uint32_t max_compute_units = dev.get_info(); - size_t workGroupSize=native_sg_size; + size_t workGroupSize=256; // if (workGroupSize > max_wg_size) // workGroupSize = max_wg_size; // Fallback for limited devices // diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc index 2202343..e91cdd8 100644 --- a/models/ecolab_model.cc +++ b/models/ecolab_model.cc @@ -94,10 +94,12 @@ RoundArray EcolabPoint::roundArray(const E& expr) void EcolabPoint::generate(unsigned niter, const ModelData& model) { array> interactionResult(density.size()); + array> lDensity(density); + for (unsigned step=0; step #endif (interactionResult[model.interaction.row[i]]) += - model.interaction.val[i]*density[model.interaction.col[i]]; + model.interaction.val[i]*lDensity[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])); + array_ns::map(lDensity.size(), [&](size_t i){ + lDensity[i]=ROUND(lDensity[i] + lDensity[i] * (model.repro_rate[i] + interactionResult[i])); }); - groupBarrier(); } + density=lDensity; assert(all(density>=0)); // // sequential/non-GPU version // for (unsigned i=0; i Date: Tue, 14 Jul 2026 12:32:26 +1000 Subject: [PATCH 10/28] bugfix: Added a groupBarrier to array::asgV, but left commented out. --- include/arrays.h | 1 + include/sycl.h | 6 +++++- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/include/arrays.h b/include/arrays.h index 90216a8..0272f25 100644 --- a/include/arrays.h +++ b/include/arrays.h @@ -1573,6 +1573,7 @@ namespace ecolab if (syclGroup().leader()) #endif resize(size); + //groupBarrier(); asg_v(data(),size,x); } diff --git a/include/sycl.h b/include/sycl.h index 394e17a..d901735 100644 --- a/include/sycl.h +++ b/include/sycl.h @@ -33,7 +33,11 @@ namespace ecolab sycl::queue& syclQ(); void* reallocSycl(void*,size_t); - inline void groupBarrier() {sycl::group_barrier(syclGroup());} + inline void groupBarrier() { +#ifdef __SYCL_DEVICE_ONLY__ + sycl::group_barrier(syclGroup()); +#endif + } inline void syncThreads() {syclQ().wait_and_throw();} From 74a953ebd065b9c85f2942dd9ab2072fac9703bb Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Tue, 14 Jul 2026 12:46:30 +1000 Subject: [PATCH 11/28] bugfix: Remove conditional expression around a groupBarrier. --- models/ecolab_model.cc | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc index e91cdd8..5f8f7f5 100644 --- a/models/ecolab_model.cc +++ b/models/ecolab_model.cc @@ -218,11 +218,7 @@ void SpatialModel::mutate() groupedForAll([newSp=newSp.data(),mut_scale=&*mut_scale,this](EcolabCell& c,size_t i) { assert(all(c.density>=0)); - auto tmp{c.mutate(*mut_scale)}; -#ifdef __SYCL_DEVICE_ONLY__ - if (syclGroup().leader() && tmp.size()) -#endif - newSp[i]=tmp; + newSp[i]=c.mutate(*mut_scale); }); hostForAll([](EcolabCell& c,size_t) { From 0b876c86efefd0dc1ee19faed7ea429efae74ab4 Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Wed, 15 Jul 2026 12:24:13 +1000 Subject: [PATCH 12/28] bugfix: more attempts to fix crashing problem. Reference counter updated solely ion one thread of work group, and all allocations called with same allocation on all threads, and return the same pointer to all threads. Still crashing after about 3E5 timesteps of spatial_ecolab.py. --- include/DeviceAllocator.h | 11 ++++-- include/arrays.h | 76 +++++++++++++++++++-------------------- include/sycl.h | 20 +++++++---- models/ecolab_model.cc | 22 ++++++------ models/spatial_ecolab.py | 14 ++++---- 5 files changed, 78 insertions(+), 65 deletions(-) diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h index eb7d868..13b6015 100644 --- a/include/DeviceAllocator.h +++ b/include/DeviceAllocator.h @@ -77,10 +77,15 @@ namespace ecolab }); nextAllocator.init(q); } + // all members of group get the same pointer void* allocate(size_t size) { if (size==0) return nullptr; if (size<=pageSize) { - auto offs=queue.dequeue(); + unsigned offs; + if (groupLeader()) offs=queue.dequeue(); +#ifdef __SYCL_DEVICE_ONLY__ + offs=sycl::group_broadcast(syclGroup(),offs); +#endif if (offs!=~0U) return memory+offs; } @@ -89,7 +94,9 @@ namespace ecolab void deallocate(void* p, size_t size) { if (!p) return; if (p>=memory && p(p)-memory)>>order); + groupBarrier(); + if (groupLeader()) + queue.enqueue((reinterpret_cast(p)-memory)>>order); return; } nextAllocator.deallocate(p,size); diff --git a/include/arrays.h b/include/arrays.h index 0272f25..0d9b758 100644 --- a/include/arrays.h +++ b/include/arrays.h @@ -1486,7 +1486,6 @@ namespace ecolab if (n==0) return nullptr; T *p; array_data *r; - //p = (char*)std::malloc((n-array_data::debug_display) * sizeof(T) + sizeof(array_data) + 16); // over allocate to allow for alignment and metadata auto allocation=n + (sizeof(array_data) + 16)/sizeof(T)+1-array_data::debug_display; p = m_allocator.allocate(allocation); @@ -1494,8 +1493,9 @@ namespace ecolab if (!p) { #ifdef __SYCL_DEVICE_ONLY__ - printf("failed to allocate %d bytes in array on group %u, thread %u\n", - sizeof(T)*n,syclGroup().get_group_linear_id(),syclGroup().get_local_linear_id()); + if (groupLeader()) + printf("failed to allocate %d bytes in array on group %u, thread %u\n", + sizeof(T)*n,syclGroup().get_group_linear_id(),syclGroup().get_local_linear_id()); #endif return nullptr; // SYCL allocator returns nullptr if not initialised } @@ -1519,47 +1519,49 @@ namespace ecolab void free(array_data *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 - ; -#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>::value || -// sycl::get_pointer_type(dt,syclQ().get_context())==sycl::usm::alloc::device) return; -//#endif if (ref()==1) { free(dt); + return; } - else - ref()--; + decrRef(); } } @@ -1567,13 +1569,9 @@ namespace ecolab // implements copying data from x to this in a GPU friendly way template - void asgV(const A& alloc, size_t size, const E& x) + void asgV(size_t size, const E& x) { -#ifdef __SYCL_DEVICE_ONLY__ - if (syclGroup().leader()) -#endif - resize(size); - //groupBarrier(); + resize(size); asg_v(data(),size,x); } @@ -1582,7 +1580,8 @@ namespace ecolab if (dt && ref()>1) { array_data* 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); @@ -1613,7 +1612,8 @@ namespace ecolab array(const array& x): m_allocator(x.m_allocator) { dt=x.dt; - if (dt) ref()++; + incrRef(); + } template @@ -1629,7 +1629,7 @@ 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()); + asgV(size(), data()); return m_allocator; } @@ -1673,7 +1673,7 @@ namespace ecolab release(); m_allocator=x.m_allocator; dt=x.dt; - if (dt) ref()++; + incrRef(); return *this; } @@ -1681,7 +1681,7 @@ namespace ecolab enable_if, 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 typename diff --git a/include/sycl.h b/include/sycl.h index d901735..de164d4 100644 --- a/include/sycl.h +++ b/include/sycl.h @@ -33,12 +33,6 @@ namespace ecolab sycl::queue& syclQ(); void* reallocSycl(void*,size_t); - inline void groupBarrier() { -#ifdef __SYCL_DEVICE_ONLY__ - sycl::group_barrier(syclGroup()); -#endif - } - inline void syncThreads() {syclQ().wait_and_throw();} /// SyclType is a pointer type, that when allocated via new is allocated from USM @@ -134,6 +128,20 @@ namespace ecolab #ifndef __SYCL_DEVICE_ONLY__ template using LocalAllocator=std::allocator; #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; + } + } diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc index 5f8f7f5..432b25d 100644 --- a/models/ecolab_model.cc +++ b/models/ecolab_model.cc @@ -170,7 +170,7 @@ unsigned PanmicticModel::condense() unsigned SpatialModel::condense() { - array total_density(species.size()); + array total_density(species.size(),0); for (auto& i: *this) total_density+=i->as()->density; // TODO #ifdef MPI_SUPPORT array recv(total_density.size()); @@ -240,16 +240,10 @@ void SpatialModel::mutate() // deallocate on device groupedForAll([newSp=newSp.data()](EcolabCell& c,size_t i) { - #ifdef __SYCL_DEVICE_ONLY__ - if (syclGroup().leader()) -#endif - { - newSp[i].clear(); - assert(newSp[i].refCnt()==0); - } + newSp[i].clear(); + assert(newSp[i].refCnt()==0); }); - if (new_sp.size()==0) return; #ifdef MPI_SUPPORT MPIbuf b; b<> (*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; + // 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,size_t) { + 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()); }); } @@ -290,7 +287,7 @@ EcolabPoint::LocalArray EcolabPoint::mutate(const E& mut_scale) sycl::joint_exclusive_scan(syclGroup(),speciations.data(),speciations.data()+nsp, offsets.data(),sycl::plus()); groupBarrier(); - if (syclGroup().leader()) + if (groupLeader()) offsets[nsp]=offsets[nsp-1]+speciations[nsp-1]; groupBarrier(); @@ -305,6 +302,7 @@ EcolabPoint::LocalArray EcolabPoint::mutate(const E& mut_scale) p[j]=i; }); + groupBarrier(); return new_sp; #else array speciations=roundArray(mut_scale * density); diff --git a/models/spatial_ecolab.py b/models/spatial_ecolab.py index 20ba90c..09f0317 100644 --- a/models/spatial_ecolab.py +++ b/models/spatial_ecolab.py @@ -28,10 +28,10 @@ def randomList(num, min, max): ecolab.species(range(nsp)) -numX=12 -numY=12 -#numX=1 -#numY=1 +numX=8 +numY=8 +#numX=2 +#numY=2 ecolab.setGrid(numX,numY) ecolab.partitionObjects() @@ -71,7 +71,7 @@ def stepImpl(): ecolab.migration([x/mut_factor for x in ecolab.migration()]) global extinctions, migrations - migrations+=ecolab.migrate() + #migrations+=ecolab.migrate() extinctions+=ecolab.condense() #ecolab.syncThreads() #ecolab.gather() @@ -90,13 +90,13 @@ def step(): migrations=0 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}') 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('No. species by cell',ecolab.tstep(),ecolab.nsp()()) plot('Extinctions',ecolab.tstep(),extinctions) plot('Migration',ecolab.tstep(),migrations) # for i in range(numX): From a14c382ec0957203a6b39cfbf61d55e7ab816982 Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Wed, 15 Jul 2026 18:21:13 +1000 Subject: [PATCH 13/28] bugfix: Change binop/unop and RVIndex expressions to copies rather than references. --- include/arrays.h | 10 +++++----- include/non-sycl.h | 1 - 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/include/arrays.h b/include/arrays.h index 0d9b758..0bfddb2 100644 --- a/include/arrays.h +++ b/include/arrays.h @@ -1137,8 +1137,8 @@ namespace ecolab template class RVindex { - const E& expr; - const I& idx; + E expr; + I idx; void operator=(const RVindex&); public: typedef typename E::value_type value_type; @@ -1246,7 +1246,7 @@ namespace ecolab class unop { public: - const E& e; + E e; Op op; unop(const E& expr): e(expr) {} @@ -1270,8 +1270,8 @@ namespace ecolab class binop { public: - const E1& e1; - const E2& e2; + E1 e1; + E2 e2; Op op; binop(const E1& ex1, const E2& ex2): e1(ex1), e2(ex2) {conformance_check(e1,e2);} diff --git a/include/non-sycl.h b/include/non-sycl.h index 0a0eb89..3a8a74c 100644 --- a/include/non-sycl.h +++ b/include/non-sycl.h @@ -11,7 +11,6 @@ namespace ecolab { template struct SyclQAllocator: public std::allocator {}; - inline void groupBarrier() {} inline void syncThreads() {} /// Non SYCL version of SyclType From 7db366417bd88a736a202f25575c7e4150e4e116 Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Wed, 15 Jul 2026 18:40:13 +1000 Subject: [PATCH 14/28] bugfix - reverted back tom how things were before the last lot of brain farts. --- include/arrays.h | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/include/arrays.h b/include/arrays.h index 0bfddb2..0f2fd12 100644 --- a/include/arrays.h +++ b/include/arrays.h @@ -1137,8 +1137,8 @@ namespace ecolab template class RVindex { - E expr; - I idx; + const E& expr; + const I& idx; void operator=(const RVindex&); public: typedef typename E::value_type value_type; @@ -1246,7 +1246,7 @@ namespace ecolab class unop { public: - E e; + const E& e; Op op; unop(const E& expr): e(expr) {} @@ -1270,8 +1270,8 @@ namespace ecolab class binop { public: - E1 e1; - E2 e2; + const E1& e1; + const E2& e2; Op op; binop(const E1& ex1, const E2& ex2): e1(ex1), e2(ex2) {conformance_check(e1,e2);} @@ -1571,8 +1571,10 @@ namespace ecolab template void asgV(size_t size, const E& x) { - resize(size); - asg_v(data(),size,x); + // copy into temporary data, as E may contain references to this + array tmp(size); + asg_v(tmp.data(),size,x); + swap(tmp); } void copy() //any nonconst method needs to call this From 03f614b410fa16cb57719ae2fb40979b7633c45f Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Wed, 15 Jul 2026 19:00:13 +1000 Subject: [PATCH 15/28] bugfix: For code running on device, we need a slightly different allocation policy. --- include/arrays.h | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/include/arrays.h b/include/arrays.h index 0f2fd12..ffbdad9 100644 --- a/include/arrays.h +++ b/include/arrays.h @@ -1572,9 +1572,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> tmp(size); + asg_v(tmp.data(),size,x); + resize(size); + asg_v(data(),size,tmp); +#else array tmp(size); asg_v(tmp.data(),size,x); swap(tmp); +#endif } void copy() //any nonconst method needs to call this From 454ceabc8004eec9d79266e774820657ac0466b7 Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Wed, 15 Jul 2026 19:36:39 +1000 Subject: [PATCH 16/28] chore: optimise array(expr) --- include/arrays.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/include/arrays.h b/include/arrays.h index ffbdad9..5b7f16a 100644 --- a/include/arrays.h +++ b/include/arrays.h @@ -1630,7 +1630,8 @@ namespace ecolab typename enable_if< is_expression, void*>::T dummy=0): m_allocator(alloc) { set_size(e.size()); - operator=(e); + asg_v(data(),e.size(),e); + //operator=(e); } ~array() {release();} From b8bc64afce4b6608088b11c342b784fd9670adbb Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Fri, 17 Jul 2026 16:55:29 +1000 Subject: [PATCH 17/28] bugfix: deallocating pointers in DeviceAllocator was incorrectly shifting the offsets --- include/DeviceAllocator.h | 10 ++++++++-- include/arrays.h | 5 +---- models/spatial_ecolab.py | 8 ++++---- 3 files changed, 13 insertions(+), 10 deletions(-) diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h index 13b6015..67e4b52 100644 --- a/include/DeviceAllocator.h +++ b/include/DeviceAllocator.h @@ -58,6 +58,11 @@ namespace ecolab void init(sycl::queue& q) {} void* allocate(size_t sz) { fatalErrorFlag()=true; + if (groupLeader()) + sycl::ext::oneapi::experimental::printf("failed to allocate %ul bytes\n",sz); +#ifndef __SYCL_DEVICE_ONLY__ + throw std::bad_alloc(); +#endif return nullptr; } void deallocate(void* p, size_t) {sycl::ext::oneapi::experimental::printf("%x leaked on device\n",p);} @@ -73,7 +78,8 @@ namespace ecolab public: void init(sycl::queue& q) { q.parallel_for(numPages, [this](auto i) { - queue.enqueue(i<=memory && p(p)-memory)>>order); + queue.enqueue(reinterpret_cast(p)-memory); return; } nextAllocator.deallocate(p,size); diff --git a/include/arrays.h b/include/arrays.h index 5b7f16a..fa2187c 100644 --- a/include/arrays.h +++ b/include/arrays.h @@ -1492,11 +1492,8 @@ namespace ecolab if (!p) { -#ifdef __SYCL_DEVICE_ONLY__ if (groupLeader()) - printf("failed to allocate %d bytes in array on group %u, thread %u\n", - sizeof(T)*n,syclGroup().get_group_linear_id(),syclGroup().get_local_linear_id()); -#endif + printf("failed to allocate %d bytes in array\n",sizeof(T)*n); return nullptr; // SYCL allocator returns nullptr if not initialised } diff --git a/models/spatial_ecolab.py b/models/spatial_ecolab.py index 09f0317..365428d 100644 --- a/models/spatial_ecolab.py +++ b/models/spatial_ecolab.py @@ -28,8 +28,8 @@ 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) @@ -71,7 +71,7 @@ def stepImpl(): ecolab.migration([x/mut_factor for x in ecolab.migration()]) global extinctions, migrations - #migrations+=ecolab.migrate() + migrations+=ecolab.migrate() extinctions+=ecolab.condense() #ecolab.syncThreads() #ecolab.gather() @@ -96,7 +96,7 @@ def step(): 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('No. species by cell',ecolab.tstep(),ecolab.nsp()()) plot('Extinctions',ecolab.tstep(),extinctions) plot('Migration',ecolab.tstep(),migrations) # for i in range(numX): From b9f37fc490287b2dfb91cc13549bbed783e62f31 Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Sat, 18 Jul 2026 09:42:24 +1000 Subject: [PATCH 18/28] refactor: Optimise generate step by precomputing row offsets --- models/ecolab_model.cc | 49 +++++++++++++++++++++++++++++------------- models/ecolab_model.h | 10 ++++++--- 2 files changed, 41 insertions(+), 18 deletions(-) diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc index 432b25d..678b014 100644 --- a/models/ecolab_model.cc +++ b/models/ecolab_model.cc @@ -93,27 +93,36 @@ RoundArray EcolabPoint::roundArray(const E& expr) // void EcolabPoint::generate(unsigned niter, const ModelData& model) { - array> interactionResult(density.size()); + //array> interactionResult(density.size()); array> lDensity(density); - + for (unsigned step=0; step -#endif - (interactionResult[model.interaction.row[i]]) += - model.interaction.val[i]*lDensity[model.interaction.col[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]]; +// for (unsigned k=0; k +//#endif +// (interactionResult[model.interaction.row[i]]) += +// model.interaction.val[i]*lDensity[model.interaction.col[i]]; +// }); +// +// groupBarrier(); - groupBarrier(); - array_ns::map(lDensity.size(), [&](size_t i){ - lDensity[i]=ROUND(lDensity[i] + lDensity[i] * (model.repro_rate[i] + interactionResult[i])); - }); +// array_ns::map(lDensity.size(), [&](size_t i){ +// lDensity[i]=ROUND(lDensity[i] + lDensity[i] * (model.repro_rate[i] + interactionResult[i])); +// }); } density=lDensity; assert(all(density>=0)); @@ -266,6 +275,8 @@ void SpatialModel::mutate() ModelData::mutate(new_sp); #endif + 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,this](EcolabCell& c,size_t) { @@ -712,6 +723,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 ModelData::computeODiagIdx() +{ + oDiagIdx.clear(); oDiagIdx.resize(species.size()); + for (unsigned i=0; i> species; array> create, repro_rate, mutation, migration; sparse_mat interaction; + using UnsignedArray=array>; + std::vector> oDiagIdx; sparse_mat_graph foodweb; unsigned long long tstep=0, last_mut_tstep=0, last_mig_tstep=0; //mutation parameters @@ -55,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& mask, size_t mask_true); void mutate(const array&); @@ -72,18 +75,19 @@ template array getArray(const array&); class EcolabPoint { public: - Float salt; /* random no. used for migration */ #ifdef SYCL_LANGUAGE_VERSION template using Allocator=GlobalDeviceAllocator; #else template using Allocator=std::allocator; #endif + using UnsignedArray=array>; + using LocalArray=array>; + + Float salt; /* random no. used for migration */ array> density; void generate(unsigned niter, const ModelData&); void condense(const array& mask, size_t mask_true); - using UnsignedArray=array>; - using LocalArray=array>; template 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 From 11fa6a8a4726e3e09789bf6071115dacb37407a9 Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Sat, 18 Jul 2026 10:21:20 +1000 Subject: [PATCH 19/28] chore: dead code removal --- models/ecolab_model.cc | 20 -------------------- 1 file changed, 20 deletions(-) diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc index 678b014..a2b87cb 100644 --- a/models/ecolab_model.cc +++ b/models/ecolab_model.cc @@ -93,7 +93,6 @@ RoundArray EcolabPoint::roundArray(const E& expr) // void EcolabPoint::generate(unsigned niter, const ModelData& model) { - //array> interactionResult(density.size()); array> lDensity(density); for (unsigned step=0; step -//#endif -// (interactionResult[model.interaction.row[i]]) += -// model.interaction.val[i]*lDensity[model.interaction.col[i]]; -// }); -// -// groupBarrier(); - -// array_ns::map(lDensity.size(), [&](size_t i){ -// lDensity[i]=ROUND(lDensity[i] + lDensity[i] * (model.repro_rate[i] + interactionResult[i])); -// }); } density=lDensity; assert(all(density>=0)); From c1cd4deffa59ff589e9132652e71154efd473d5c Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Sat, 18 Jul 2026 10:45:11 +1000 Subject: [PATCH 20/28] bugfix: missing oDiagIdx calculation in condense --- include/Makefile.config | 17 ----------------- models/ecolab_model.cc | 9 +-------- 2 files changed, 1 insertion(+), 25 deletions(-) delete mode 100644 include/Makefile.config diff --git a/include/Makefile.config b/include/Makefile.config deleted file mode 100644 index 6f8edc1..0000000 --- a/include/Makefile.config +++ /dev/null @@ -1,17 +0,0 @@ -TCL=1 -TK=1 -ZLIB=1 -READLINE=1 -XDR=1 -UNURAN=1 -GNUSL=1 -CAIRO=1 -BDB=1 -PANGO=1 -GDBM_COMPAT=1 -MPI= -PARALLEL= -OPENMP= -GCC= -NOGUI= -AQUA= diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc index a2b87cb..8d60170 100644 --- a/models/ecolab_model.cc +++ b/models/ecolab_model.cc @@ -144,6 +144,7 @@ void ModelData::condense(const array& 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; } @@ -201,19 +202,11 @@ void SpatialModel::mutate() vector> newSp(size()); - hostForAll([](EcolabCell& c,size_t) { - assert(all(c.density>=0)); - }); - 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); }); - hostForAll([](EcolabCell& c,size_t) { - assert(all(c.density>=0)); - }); - array new_sp; DeviceType>> cell_ids; syncThreads(); From cbdab62f4a447eb972800abecae1943a109a5447 Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Sat, 18 Jul 2026 17:11:20 +1000 Subject: [PATCH 21/28] bugfix: printf format string needs a %zu for size_t --- include/arrays.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/arrays.h b/include/arrays.h index fa2187c..d13d58a 100644 --- a/include/arrays.h +++ b/include/arrays.h @@ -1493,7 +1493,7 @@ namespace ecolab if (!p) { if (groupLeader()) - printf("failed to allocate %d bytes in array\n",sizeof(T)*n); + printf("failed to allocate %zu bytes in array\n",sizeof(T)*n); return nullptr; // SYCL allocator returns nullptr if not initialised } From 8342fad021ec34a5e9c167751ce3b948d9bdeec8 Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Sat, 18 Jul 2026 18:23:57 +1000 Subject: [PATCH 22/28] bugfix: update pred-prey test to current libs --- test/00/predator-prey.sh | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/test/00/predator-prey.sh b/test/00/predator-prey.sh index 999a8b9..9ca09d1 100644 --- a/test/00/predator-prey.sh +++ b/test/00/predator-prey.sh @@ -54,16 +54,16 @@ traceOption=0 case $traceOption in 0) cat >out1.dat < Date: Sun, 19 Jul 2026 09:49:05 +1000 Subject: [PATCH 23/28] Fix Queue enqueue/dequeue with bounded MPMC per-slot sequencing --- include/DeviceAllocator.h | 73 ++++++++++++++++++++++++++++++++++----- 1 file changed, 65 insertions(+), 8 deletions(-) diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h index 67e4b52..59a8f93 100644 --- a/include/DeviceAllocator.h +++ b/include/DeviceAllocator.h @@ -27,27 +27,84 @@ namespace ecolab return sycl::ext::oneapi::group_local_memory(syclGroup(),false)->flag; } - // Lock free circular buffer queue for SYCL + // Bounded MPMC circular buffer queue for SYCL using per-slot sequence numbers. + // dequeue() returns ~0U when queue appears empty (non-blocking empty signal). template class Queue { static_assert((size&(size-1))==0,"size must be power of two"); constexpr static unsigned mask=size-1; - unsigned queue[size]; + + struct Slot + { + unsigned seq; + unsigned value; + }; + + Slot slots[size]; unsigned head=0, tail=0; + using Atomic=sycl::atomic_ref; + using AtomicAcquire=sycl::atomic_ref; + using AtomicRelease=sycl::atomic_ref; public: + Queue() + { + for (unsigned i=0; i(syclGroup());} From c9eabc084c596aa48e448783fd3c619c74298435 Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Sun, 19 Jul 2026 10:11:24 +1000 Subject: [PATCH 24/28] Fix groupedForAll fatal-error tracking with atomic fetch_or and post-wait check --- include/ecolab.h | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/include/ecolab.h b/include/ecolab.h index 92a29e5..a395200 100644 --- a/include/ecolab.h +++ b/include/ecolab.h @@ -185,8 +185,8 @@ namespace ecolab num_work_groups=std::min(num_work_groups,this->size()); //std::cout< fatalError(false); - for (size_t cellStart=0; cellStartsize() && !*fatalError; cellStart+=num_work_groups) + DeviceType fatalError(0); + for (size_t cellStart=0; cellStartsize(); 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()) { @@ -201,7 +201,8 @@ namespace ecolab f(cell,idx); } // flag fatal error to throw afterwards. - *fatalError=fatalErrorFlag(); + if (fatalErrorFlag()) + sycl::atomic_ref(*fatalError).fetch_or(1); }); }); syclQ().wait_and_throw(); From 5e1b8720dababce93dc08fb931cf6a25c59490e2 Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Sun, 19 Jul 2026 13:26:45 +1000 Subject: [PATCH 25/28] chore: address code review comments. --- include/DeviceAllocator.h | 13 ++++++++----- include/arrays.h | 4 +++- include/ecolab.h | 2 +- include/sycl.h | 8 ++++++-- models/ecolab_model.cc | 6 +++--- models/ecolab_model.h | 2 +- 6 files changed, 22 insertions(+), 13 deletions(-) diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h index 59a8f93..f86be2f 100644 --- a/include/DeviceAllocator.h +++ b/include/DeviceAllocator.h @@ -45,8 +45,8 @@ namespace ecolab unsigned head=0, tail=0; using Atomic=sycl::atomic_ref; - using AtomicAcquire=sycl::atomic_ref; - using AtomicRelease=sycl::atomic_ref; + using AtomicAcquire=sycl::atomic_ref; + using AtomicRelease=sycl::atomic_ref; public: Queue() @@ -116,13 +116,13 @@ namespace ecolab void* allocate(size_t sz) { fatalErrorFlag()=true; if (groupLeader()) - sycl::ext::oneapi::experimental::printf("failed to allocate %ul bytes\n",sz); + sycl::ext::oneapi::experimental::printf("failed to allocate %zu bytes\n",sz); #ifndef __SYCL_DEVICE_ONLY__ throw std::bad_alloc(); #endif return nullptr; } - void deallocate(void* p, size_t) {sycl::ext::oneapi::experimental::printf("%x leaked on device\n",p);} + void deallocate(void* p, size_t) {sycl::ext::oneapi::experimental::printf("%p leaked on device\n",p);} }; template class DeviceAllocator @@ -170,7 +170,10 @@ namespace ecolab inline DeviceAllocator<>& deviceAllocator() { static DeviceType> deviceAllocator; - static int initialised=(deviceAllocator->init(syclQ()),0); + static int initialised=( + deviceAllocator->init(syclQ()), + syclQ().wait_and_throw(), + 0); return *deviceAllocator; } diff --git a/include/arrays.h b/include/arrays.h index d13d58a..8ad4bc3 100644 --- a/include/arrays.h +++ b/include/arrays.h @@ -1636,7 +1636,9 @@ namespace ecolab const Allocator& allocator() const {return m_allocator;} const Allocator& allocator(const Allocator& alloc) { if (alloc==m_allocator) return m_allocator; - asgV(size(), data()); + array tmp(size(),alloc); + tmp.asgV(size(), data()); + swap(tmp); return m_allocator; } diff --git a/include/ecolab.h b/include/ecolab.h index a395200..07e4b3a 100644 --- a/include/ecolab.h +++ b/include/ecolab.h @@ -170,7 +170,7 @@ namespace ecolab // 4. Maximum number of compute units (Execution units / DSS count) uint32_t max_compute_units = dev.get_info(); - size_t workGroupSize=256; + size_t workGroupSize=native_sg_size;//256; // if (workGroupSize > max_wg_size) // workGroupSize = max_wg_size; // Fallback for limited devices // diff --git a/include/sycl.h b/include/sycl.h index de164d4..654958e 100644 --- a/include/sycl.h +++ b/include/sycl.h @@ -21,6 +21,10 @@ #error "EcoLab requires OneAPI compiler for some experimental functions" #endif +#ifdef _OPENMP +#include +#endif + namespace ecolab { using sycl::ext::oneapi::experimental::printf; @@ -53,7 +57,7 @@ namespace ecolab struct SyclQAllocator: public graphcode::Allocator { SyclQAllocator(): graphcode::Allocator(syclQ(), UA) {} - template struct rebind {using other=SyclQAllocator;}; + template struct rebind {using other=SyclQAllocator;}; }; // random numbers @@ -79,7 +83,7 @@ namespace ecolab #ifdef __SYCL_DEVICE_ONLY__ return rngs[syclItem().get_global_linear_id() % rngs.size()](); #elif defined(_OPENMP) - return rngs[omp_get_thread_num()%numRngs](); + return rngs[omp_get_thread_num()%rngs.size()](); #else return rngs[0](); #endif diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc index 8d60170..ab7f2b1 100644 --- a/models/ecolab_model.cc +++ b/models/ecolab_model.cc @@ -263,6 +263,7 @@ template 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(); @@ -281,9 +282,8 @@ EcolabPoint::LocalArray EcolabPoint::mutate(const E& mut_scale) LocalArray new_sp(offsets[nsp]); array_ns::map(nsp, [offsets=offsets.data(),new_sp=new_sp.data()](size_t i) { - auto p=new_sp+offsets[i]; - for (auto j=0; j, } array nsp() const; void makeConsistent(); - void seed(unsigned x) {forAll([=](EcolabCell& cell,size_t){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 From b088debc60ac5ab605c7fe849d0e9f393ab0d67b Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Sun, 19 Jul 2026 16:09:24 +1000 Subject: [PATCH 26/28] bugfix: fixed refactor of DeviceAllocator --- include/DeviceAllocator.h | 49 +++++++++++++++------------------------ 1 file changed, 19 insertions(+), 30 deletions(-) diff --git a/include/DeviceAllocator.h b/include/DeviceAllocator.h index f86be2f..b4402a5 100644 --- a/include/DeviceAllocator.h +++ b/include/DeviceAllocator.h @@ -42,17 +42,17 @@ namespace ecolab }; Slot slots[size]; - unsigned head=0, tail=0; + unsigned head=size, tail=0; - using Atomic=sycl::atomic_ref; - using AtomicAcquire=sycl::atomic_ref; - using AtomicRelease=sycl::atomic_ref; + using Atomic=sycl::atomic_ref; public: Queue() { - for (unsigned i=0; i class DeviceAllocator { public: - void init(sycl::queue& q) {} void* allocate(size_t sz) { - fatalErrorFlag()=true; if (groupLeader()) sycl::ext::oneapi::experimental::printf("failed to allocate %zu bytes\n",sz); -#ifndef __SYCL_DEVICE_ONLY__ +#ifdef __SYCL_DEVICE_ONLY__ + fatalErrorFlag()=true; +#else throw std::bad_alloc(); #endif return nullptr; @@ -133,13 +133,6 @@ namespace ecolab char memory[poolSize]; DeviceAllocator nextAllocator; // next size up allocator public: - void init(sycl::queue& q) { - q.parallel_for(numPages, [this](auto i) { - for (unsigned j=i; j=memory && p(p)-memory); + queue.enqueue((reinterpret_cast(p)-memory)>>order); return; } nextAllocator.deallocate(p,size); @@ -170,10 +163,6 @@ namespace ecolab inline DeviceAllocator<>& deviceAllocator() { static DeviceType> deviceAllocator; - static int initialised=( - deviceAllocator->init(syclQ()), - syclQ().wait_and_throw(), - 0); return *deviceAllocator; } From ae4e3448ae0444c9919e2ebcfd43cee99a339173 Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Sun, 19 Jul 2026 16:42:44 +1000 Subject: [PATCH 27/28] bugfix: synchronous update was corrupted by prefious refactor. --- models/ecolab_model.cc | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/models/ecolab_model.cc b/models/ecolab_model.cc index ab7f2b1..c90e4c0 100644 --- a/models/ecolab_model.cc +++ b/models/ecolab_model.cc @@ -93,7 +93,7 @@ RoundArray EcolabPoint::roundArray(const E& expr) // void EcolabPoint::generate(unsigned niter, const ModelData& model) { - array> lDensity(density); + array> lDensity(density), tmp(density.size()); for (unsigned step=0; step=0)); From 2d45497bcd931fbc01318b5425707a53e7a1100e Mon Sep 17 00:00:00 2001 From: Russell Standish Date: Sun, 19 Jul 2026 16:57:52 +1000 Subject: [PATCH 28/28] Revert "bugfix: update pred-prey test to current libs" This reverts commit 8342fad021ec34a5e9c167751ce3b948d9bdeec8. --- test/00/predator-prey.sh | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/test/00/predator-prey.sh b/test/00/predator-prey.sh index 9ca09d1..999a8b9 100644 --- a/test/00/predator-prey.sh +++ b/test/00/predator-prey.sh @@ -54,16 +54,16 @@ traceOption=0 case $traceOption in 0) cat >out1.dat <