Accelerate

[Ben Lever]

Hi Manuel,

Trevor mentioned there were some problems with the trac server - that may explain why my previous response wasn't going through to the mailing list. Hopefully this one goes through :)

Rami and I have been talking with Trevor this morning and he's been filling us in on his discussion with you about the stencil proposals. I'll summarise below what we discussed and agreed on and hopefully that matches your expectations:

Here's our latest take on the function type:

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

Comments:

 *   Output array extent is the same as input array. This means stencil behaves more like a map operation - here the mapping function instead takes a set of gathered values (as described by the stencil pattern). We realised that our thoughts of adding a backpermute-type argument to stencil were flawed as ti removed all the benefits of stencil knowing statically where to draw input values from for a given output value.
 *   Trevor clarified that the dim in StencilPattern simply means that the stencil pattern is of the same dimensionality of the input (and output) array(s), but obviously can be different extents. This is fine.
 *   The above type doesn't show support for boundary conditions. As mentioned below, we think that this should be part of specifying a "stencil" (as opposed to an argument of the stencil function) - i.e. it is specified along with the pattern, etc. After chatting with Trevor, we think it's reasonable that a boundary condition enumeration (e.g. clamp, mirror, constant, wrap) be independent of the coordinate being fetched.

For completeness, we can extended to multiple stencils like so:

stencil2 :: (Ix dim, Elem a, Elem b, Elem c, StencilPattern dim a pat, Stencil Pattern dim b pat')
         => (pat -> pat' -> Exp c)  -- ^Stencil function
         -> Acc (Array dim a)       -- ^First source array
         -> Acc (Array dim b)       -- ^Second source array
         -> Acc (Array dim c)       -- ^Destination array

stencil3 :: (Ix dim, Elem a, Elem b, Elem c, Elem d, StencilPattern dim a pat, Stencil Pattern dim b pat', StencilPattern dim c pat'')
         => (pat -> pat' -> pat'' -> Exp c)  -- ^Stencil function
         -> Acc (Array dim a)                -- ^First source array
         -> Acc (Array dim b)                -- ^Second source array
         -> Acc (Array dim c)                -- ^Third source array
         -> Acc (Array dim d)                -- ^Destination array


So, from our point of view, we're happy with the above. Let us know if you think there should be any further additions/removals.


Regards,
Ben.


On 03/08/2010, at 9:18 AM, Ben Lever wrote:

Hi Manuel,

I'm not sure whether or not you got this e-mail I sent yesterday. There seems to be quite a delay in sending to the mailing list and it actually sending out emails to subscribers - at least for Rami, Trevor and I.

Also, we've just come across a paper (attached) that is very related to the Accelerate work. Do you have much insight into their work and how their framework compares to Accelerate? We're planning on looking into it more detail but would be interested to hear your perspective.


Regards,
Ben.

Begin forwarded message:

From: Ben Lever <Ben.Lever at nicta.com.au<mailto:Ben.Lever at nicta.com.au>>
Date: 2 August 2010 4:08:46 PM AEST
To: Manuel M T Chakravarty <chak at cse.unsw.edu.au<mailto:chak at cse.unsw.edu.au>>
Cc: "accelerate at projects.haskell.org<mailto:accelerate at projects.haskell.org>" <accelerate at projects.haskell.org<mailto:accelerate at projects.haskell.org>>
Subject: Re: Accelerate

Hi Manuel,

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.

This version of stencil looks much simpler, which is great. Some comments:


 *   We need some way of specifying which input array element is the focal point for a given output array element. That is, a backpermute type function - e.g. (dim' -> dim). This is necessary when the output extent is different to the input.
 *   From what we can see, the output extent can only have a dimension that is equal or less than that of the input. We do have some examples that could stencil to go up a dimension, but we believe they're not particularly compelling ones. So, we think this constraint should be ok.
 *   In "StencilPattern dim a pat" above, does the use of "dim" mean that the dimensionality of the stencil pattern is the same as the source array? Or is this a typo and the stencil pattern can be of any dimension?
 *   Boundary conditions seem to only be relevant to the stencil-type operation, so could it be something that is specified as part of StencilPattern? Perhaps it would be enough to simply enumerate how the stencil should handle an out-of-bounds coordinate, for example, (dim -> BoundCond e) where BoundCond could have data constructions like Clamp, Mirror, Wrap and Constant e. This would be enough to cover the all the cases we know about without having an effect on performance.


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.).

I agree that the code is very intuitive for smaller stencils. For larger stencils (e.g. 36x36),  we should be able to use all of Haskell to generate the resulting Exp type.


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?

In our latest "test" version of stencil, we opted for disconnecting dim' from dim. also We don't believe there is any need to make dim' a function of dim. Even if there is, the relationship is static and can be computed outside the call to stencil. As mentioned above, we believe you can't go up a dimension with your stencil function - is there a way to have the types enforce this?


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.

We wish you luck in figuring out how to implement it :)

Will the implementation provide a way to calculate the rectangular-bounds of the neighbourhood specified by a stencil pattern in order to determine what should be loaded into shared-memory on a GPU, etc?


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