I think the GPU implementation mutate is now working, but it is so slow compared with running on the host.

Author: highperformancecoder

include/arrays.h

@@ -1469,7 +1469,7 @@ namespace ecolab
     template <class T, class A=std::allocator<T>>
     class array
     {
-      array_data<T> *dt;
+      array_data<T> *dt=nullptr;
       A m_allocator;
 
 
@@ -1492,6 +1492,10 @@ 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
         size_t d = (size_t)(reinterpret_cast<array_data<T>*>(p)->dt);
@@ -1557,15 +1561,37 @@ namespace ecolab
 
     protected:
 
+      // implements copying data from x to this in a GPU friendly way
+      template <class E>
+      void asgV(const A& alloc, size_t size, const E& x)
+      {
+#ifdef __SYCL_DEVICE_ONLY__
+        GroupLocal<array> 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
+      }
+      
       void copy() //any nonconst method needs to call this
       {           // to implement copy-on-write semantics
         if (dt && ref()>1)
           {
+#ifdef __SYCL_DEVICE_ONLY__
+            syclPrintf("b4 asgV in copy\n");
+            asgV(m_allocator, size(), dt->dt);
+            syclPrintf("after asgV in copy\n");
+#else
             array_data<T>* oldData=dt;
             bool freeMem = ref()-- == 0;
-            dt=alloc(size()); 
+            dt=alloc(size());
             memcpy(dt->dt,oldData->dt,size()*sizeof(T));
             if (freeMem) free(oldData);
+#endif
           }
       }
 
@@ -1574,7 +1600,7 @@ namespace ecolab
       typedef size_t size_type; 
       using Allocator=A;
 
-      array(const Allocator& alloc={}): m_allocator(alloc) {set_size(0);}
+      array(const Allocator& alloc={}): m_allocator(alloc) {}
       explicit array(size_t s, const Allocator& alloc=Allocator()): m_allocator(alloc)
       {
         set_size(s);
@@ -1605,11 +1631,12 @@ namespace ecolab
       const Allocator& allocator() const {return m_allocator;}
       const Allocator& allocator(const Allocator& alloc) {
         if (alloc==m_allocator) return m_allocator;
-        array tmp(size(),alloc);
-        asg_v(tmp.data(),size(),data());
-        swap(tmp);
+        asgV(alloc, size(), data());
         return m_allocator;
       }
+
+      /// current value of the reference counter
+      unsigned refCnt() const {return dt->cnt;}
 
       /// resize array to \a s elements
       void resize(size_t s) {
@@ -1645,33 +1672,18 @@ namespace ecolab
 
       array& operator=(const array& x) {
         if (x.dt==dt) return *this;
-        if (m_allocator==x.m_allocator) { // shared data optimisation
-          release();
-          dt=x.dt;
-          if (dt) ref()++;
-          return *this;
-        } 
-        array tmp(x.size(),m_allocator);
-        array_ns::asg_v(tmp.data(),tmp.size(),x);
-        swap(tmp);
+        release();
+        m_allocator=x.m_allocator;
+        dt=x.dt;
+        if (dt) ref()++;
         return *this;
       }
 
       template <class expr> typename
       enable_if<is_expression<expr>, array&>::T 
       operator=(const expr& x) {
         if ((void*)(&x)==(void*)(this)) return *this;
-        // since expression x may contain a reference to this, assign to a temporary
-#ifdef __SYCL_DEVICE_ONLY__
-        GroupLocal<array> tmp(x.size(),m_allocator);
-        array_ns::asg_v(tmp.ref().data(),x.size(),x);
-        groupBarrier();
-        if (syclGroup().leader()) swap(tmp.ref());
-#else
-        array tmp(x.size(),m_allocator);
-        array_ns::asg_v(tmp.data(),tmp.size(),x);
-        swap(tmp);
-#endif
+        asgV(m_allocator, x.size(), x);
         return *this;
       }
       template <class expr> typename

include/ecolab.h

@@ -205,57 +205,52 @@ namespace ecolab
 
   struct CellBase
   {
+    size_t m_idx=0; // stash the position within the local node vector here
+    size_t idx() const {return m_idx;}
 #ifdef SYCL_LANGUAGE_VERSION
-    Ouro::SyclDesc<1>* desc=nullptr;
     using MemAllocator=Ouro::MultiOuroPQ;
     MemAllocator* memAlloc=nullptr;
-    const sycl::stream* out=nullptr;
-    size_t idx() const {
-      if (desc) return desc->item.get_global_linear_id();
-      else return 0;
-    }
     template <class T> class CellAllocator
     {
     public:
-      Ouro::SyclDesc<1>*const* desc=nullptr;
-      MemAllocator*const* memAlloc=nullptr;
+      MemAllocator* memAlloc=nullptr;
       template <class U> friend class Allocator;
       CellAllocator()=default;
-      CellAllocator(Ouro::SyclDesc<1>* const& desc, MemAllocator*const& memAlloc):
-        desc(&desc), memAlloc(&memAlloc) {
-      }
-      template <class U> CellAllocator(const CellAllocator<U>& x):
-        desc(x.desc), memAlloc(x.memAlloc) {}
+      CellAllocator(MemAllocator* memAlloc): memAlloc(memAlloc) {}
+      template <class U> CellAllocator(const CellAllocator<U>& x): memAlloc(x.memAlloc) {}
       T* allocate(size_t sz) {
 #ifdef  __SYCL_DEVICE_ONLY__
-        if (memAlloc && *memAlloc && desc && *desc)  {
-          auto r=reinterpret_cast<T*>((*memAlloc)->malloc(**desc,sz*sizeof(T)));
+        if (memAlloc)  {
+          auto r=reinterpret_cast<T*>(memAlloc->malloc(Ouro::SyclDesc<>(syclItem(),{}),sz*sizeof(T)));
           if (!r) syclPrintf("Mem allocation failed\n");
           return r;
         }
+        syclPrintf("Missing allocator memAlloc=%x\n",memAlloc);
         return nullptr; // TODO raise an error??
 #else
         return sycl::malloc_shared<T>(sz,syclQ());
 #endif
       }
-      void deallocate(T* p,size_t) {
+      void deallocate(T* p,size_t n) {
+        if (!p) return;
 #ifdef  __SYCL_DEVICE_ONLY__
-        if (memAlloc && *memAlloc && desc && *desc)
-          (*memAlloc)->free(**desc,p);
+        if (memAlloc && reinterpret_cast<Ouro::memory_t*>(p)>=memAlloc->memory.d_data &&
+            reinterpret_cast<Ouro::memory_t*>(p)<memAlloc->memory.d_data_end)
+          memAlloc->free(Ouro::SyclDesc<>(syclItem(),{}),p);
+        else
+          syclPrintf("leaked %d bytes\n",n*sizeof(T));
 #else
         sycl::free(p,syclQ());
 #endif
       }
-      bool operator==(const CellAllocator& x) const {return desc==x.desc && memAlloc==x.memAlloc;}
+      bool operator==(const CellAllocator& x) const {return memAlloc==x.memAlloc;}
     };
     template <class T> CellAllocator<T> allocator() const {
-      return CellAllocator<T>(desc,memAlloc);
+      return CellAllocator<T>(memAlloc);
     }
-#else
+#else //!SYCL
     template <class T> using CellAllocator=std::allocator<T>;
     template <class T> CellAllocator<T> allocator() const {return CellAllocator<T>();}
-    size_t m_idx=0; // stash the position within the local node vector here
-    size_t idx() const {return m_idx;}
 #endif
   };
 
@@ -304,6 +299,21 @@ namespace ecolab
        graphcode::Allocator<graphcode::ObjectPtr<Cell>>(syclQ(),sycl::usm::alloc::shared)) {}
     ~EcolabGraph() {syncThreads();}
 #endif
+    /// apply a functional to all local cells of this processor using the host processor
+    /// @param f 
+    template <class F>
+    void hostForAll(F f) {
+      auto sz=this->size();
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+      for (size_t idx=0; idx<sz; ++idx) {
+        auto& cell=*(*this)[idx]->template as<Cell>();
+        cell.m_idx=idx;
+        f(cell);
+      }
+    }
+    
     /// apply a functional to all local cells of this processor in parallel
     /// @param f 
     template <class F>
@@ -313,34 +323,17 @@ namespace ecolab
       size_t range=this->size()/workGroupSize;
       if (range*workGroupSize < this->size()) ++range;
       syclQ().submit([&](auto& h) {
-#ifndef NDEBUG
-        sycl::stream out(1000000,1000,h);
-#endif
         h.parallel_for(sycl::nd_range<1>(range*workGroupSize, workGroupSize), [=,this](auto i) {
           auto idx=i.get_global_linear_id();
           if (idx<this->size()) {
             auto& cell=*(*this)[idx]->template as<Cell>();
-            Ouro::SyclDesc<> desc(i,{});
-            cell.desc=&desc;
-#ifndef NDEBUG
-            cell.out=&out;
-#endif
+            cell.m_idx=idx;
             f(cell);
-            cell.desc=nullptr;
-            cell.out=nullptr;
           }
         });
       });
 #else
-      auto sz=this->size();
-#ifdef _OPENMP
-#pragma omp parallel for
-#endif
-      for (size_t idx=0; idx<sz; ++idx) {
-        auto& cell=*(*this)[idx]->template as<Cell>();
-        cell.m_idx=idx;
-        f(cell);
-      }
+      hostForAll(f);
 #endif
     }
 
@@ -351,28 +344,20 @@ namespace ecolab
     template <class F>
     void groupedForAll(F f) {
 #ifdef SYCL_LANGUAGE_VERSION
+      // TODO - pass in workGroupSize as an optional parameter??
       static size_t workGroupSize=32;//syclQ().get_device().get_info<sycl::info::device::max_work_group_size>();
       syclQ().submit([&](auto& h) {
-#ifndef NDEBUG
-        sycl::stream out(1000000,1000,h);
-#endif
         h.parallel_for(sycl::nd_range<1>(this->size()*workGroupSize, workGroupSize), [=,this](auto i) {
           auto idx=i.get_group_linear_id();
           if (idx<this->size()) {
             auto& cell=*(*this)[idx]->template as<Cell>();
-            Ouro::SyclDesc<> desc(i,{});
-            cell.desc=&desc;
-#ifndef NDEBUG
-            cell.out=&out;
-#endif
+            cell.m_idx=idx;
             f(cell);
-            cell.desc=nullptr;
-            cell.out=nullptr;
           }
         });
       });
 #else
-      forAll(f);
+      hostForAll(f);
 #endif
     }
 
@@ -423,10 +408,17 @@ namespace ecolab
     ~GroupLocal() {
       sycl::group_barrier(syclGroup());
       if (syclGroup().leader())
-        ref().~T();
+        (**this).~T();
     }
-    T& ref() {return reinterpret_cast<T&>(**buffer);}
+    T& operator*() {return reinterpret_cast<T&>(**buffer);}
+#else
+    T buffer;
+  public:
+    template <class... Args>
+    GroupLocal(Args&&... args): buffer(std::forward<Args>(args)...) {}
+    T& operator*() {return buffer;}
 #endif
+    T* operator->() {return &**this;}
   };
 
 #ifdef __SYCL_DEVICE_ONLY__