Accelerate

[Manuel M T Chakravarty]

Ben,

> It's starting to look overly complicated - I feel there's a better way. Perhaps the stencil description should be extracted out as something separate given that it specifies so much?

Yes, I'd also like to abstract stencils out as a separate type.  Regarding a stencil as an array and a reduction function also has some disadvantages.  We don't always need all points of a regular array for a stencil (eg, in a 5-point stencil, we only use 5 of 9 points) and the reduction operation is somewhat cluttered as it also gets the indices.

I have been wondering whether we want to make stencil more explicit and represent them as proper functions whose pattern characterises the stencil and whose right-hand side implements the reduction.  Stencil patterns would be represented as nested tuples, which means that they can be irregular.  Still disregarding boundary conditions, the stencil function would become much more compact:

stencil
   :: (Ix dim, Ix dim', Elem a, Elem b, StencilPattern dim a pat)
   => (pat -> Exp b)					       -- ^Stencil function
   -> Exp dim'                                                 -- ^Extent of the destination array
   -> Acc (Array dim a)                                        -- ^Source array
   -> Acc (Array dim' b)                                       -- ^Destination array

The type class StencilPattern would ensure that only sensible types are used as stencil patterns.  I append some example stencil functions at the end of this message.

What do you think?  Does this make sense?  It is obviously biased towards smaller stencils, but AFAIK that is what is mostly used in practice.  In return, I believe the code is far more intuitive for smaller stencil.  It is also straight forward to extend to multiple stencils (ie, stencil2, etc.).

BTW, currently the dimensions of the source array (dim) and the result array (dim') are completely unconnected.  Is there any constraint?  In the two examples in the ticket, dim' is always greater (more dimensions) or equal dim.  Is that always the case?

Manuel

PS: The stencil functions below do type check.  I haven't tried to implement the type class 'StencilPattern' yet.  I have an idea how I would approach it, but it might turn out to be tricky.  I would just like to know what you think about this approach before pushing further in this direction.

PPS: I think everybody except BenLi is on the mailing list now.

-=-

data F a = F a    -- marks the focal point

stencil1D :: (Elem a, IsFloating a) 
          => (Exp a, F (Exp a), Exp a) -> Exp a
stencil1D (x, F y, z) = (x + z - 2 * y) / 2

stencil2D5 :: (Elem a, IsFloating a) 
           => (Exp a, (Exp a, F (Exp a), Exp a), Exp a) -> Exp a
stencil2D5 (      t
           , (l, F m, r)
           ,      b    ) = (t + l + r + b - 4 * m) / 4

stencil2D :: (Elem a, IsFloating a) 
          => ((Exp a,   (Exp a), Exp a),
              (Exp a, F (Exp a), Exp a),
              (Exp a,   (Exp a), Exp a)) -> Exp a
stencil2D ( (t1,  t2, t3)
          , (l , F m, r )
          , (b1,  b2, b3)
          ) 
          = (t1/2 + t2 + t3/2 + l + r + b1/2 + b2 + b3/2 - 4 * m) / 4

stencil3D :: (Elem a, IsNum a)
          => (((Exp a, Exp a), (Exp a, Exp a)), ((Exp a, Exp a), (Exp a, F (Exp a)))) -> Exp a
stencil3D (front, back) =
  let ((f1,   f2),
       (f3,   f4)) = front
      ((b1,   b2),
       (b3, F b4)) = back
  in
  f1 + f2 + f3 + f4 + b1 + b2 + b3 + b4