         std::map<ptrdiff_t,int>::const_reverse_iterator rit = strides.rbegin();
         for (int i = 0; i < N; i++, rit++)
           order.index[i] = rit->second;
-        if (!space.dense()) 
+        
+        if (fold)
         {
-          // The index space is not dense so launch a kernel for each rectangle
-          Realm::IndexSpaceIterator<N,T> iterator(space);
-          while (true)
+          const Realm::AffineAccessor<typename REDOP::RHS,N,T> 
+            dst_accessor(dst, field_id, space.bounds);
+          if (!space.dense()) 
           {
-            const size_t volume = iterator.rect.volume();
+            // The index space is not dense so launch a kernel for each rect
+            Realm::IndexSpaceIterator<N,T> iterator(space);
+            while (true)
+            {
+              const size_t volume = iterator.rect.volume();
+              const size_t blocks = 
+                (volume + THREADS_PER_BLOCK - 1) / THREADS_PER_BLOCK;
+              Point<N,T> extents;
+              size_t pitch = 1;
+              for (std::map<ptrdiff_t,int>::const_iterator it = 
+                    strides.begin(); it != strides.end(); it++)
+              {
+                extents[it->second] = pitch;
+                pitch *= ((iterator.rect.hi[it->second] - 
+                            iterator.rect.lo[it->second]) + 1);
+              }
+              if (exclusive)
+                fold_kernel<REDOP,N,T,true><<<blocks,THREADS_PER_BLOCK>>>(
+                    dst_accessor, src_accessor, order, 
+                    iterator.rect.lo, extents, volume); 
+              else
+                fold_kernel<REDOP,N,T,false><<<blocks,THREADS_PER_BLOCK>>>(
+                    dst_accessor, src_accessor, order, 
+                    iterator.rect.lo, extents, volume);
+              if (!iterator.step())
+                break;
+            }
+          }
+          else
+          {
+            // Space is dense so we just need a single kernel launch here
+            const size_t volume = space.bounds.volume();
             const size_t blocks = 
               (volume + THREADS_PER_BLOCK - 1) / THREADS_PER_BLOCK;
             Point<N,T> extents;

        offset = offset % extents[index];
      }
      REDOP::template apply<EXCLUSIVE>(dst[point], src[point]);
    }

    template<typename REDOP>
    struct ReductionRunner {

    template<typename REDOP, int N, typename T, bool EXCLUSIVE>
    __global__ void
    __launch_bounds__(THREADS_PER_BLOCK,MIN_BLOCKS_PER_SM)
    fold_kernel(const Realm::AffineAccessor<typename REDOP::RHS,N,T> dst,
                const Realm::AffineAccessor<typename REDOP::RHS,N,T> src,
                const DimOrder<N> order,
                const Realm::Point<N,T> lo,
                const Realm::Point<N,T> extents,
                const size_t max_offset)
    {
      size_t offset = blockIdx.x * blockDim.x + threadIdx.x;
      if (offset >= max_offset)
        return;
      Point<N,T> point = lo;
      for (int i = 0; i < N; i++)
      {
        const int index = order.index[i];
        point[index] += (offset / extents[index]);
        offset = offset % extents[index];
      }
      REDOP::template fold<EXCLUSIVE>(dst[point], src[point]);
    }

    template<typename REDOP, int N, typename T, bool EXCLUSIVE>
    __global__ void
    __launch_bounds__(THREADS_PER_BLOCK,MIN_BLOCKS_PER_SM)
    apply_kernel(const Realm::AffineAccessor<typename REDOP::LHS,N,T> dst,
                 const Realm::AffineAccessor<typename REDOP::RHS,N,T> src,
                 const DimOrder<N> order,
                 const Realm::Point<N,T> lo,
                 const Realm::Point<N,T> extents,
                 const size_t max_offset)
    {
      size_t offset = blockIdx.x * blockDim.x + threadIdx.x;
      if (offset >= max_offset)
