Accelerate

[Manuel M T Chakravarty]

Hi Ben,

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

The Trac server was only temporarily down for a few hours.

> 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 didn't respond earlier as I wanted to implement it first and see whether everything works out.  I'm still implementing it, but I'll summarise the state of the play below.

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

Yes, that is basically what I'm implementing right now.  Here are the types:

> -- |Boundary condition specification for stencil operations.
> --
> data Boundary a = Clamp               -- ^clamp coordinates to the extent of the array
>                 | Mirror              -- ^mirror coordinates beyond the array extent
>                 | Wrap                -- ^wrap coordinates around on each dimension
>                 | Constant a          -- ^use a constant value for outlying coordinates class Stencil dim a stencil
> 
> class Stencil dim a stencil
> 
> -- |Map a stencil over an array.  In contrast to 'map', the domain of a stencil function is an
> --  entire /neighbourhood/ of each array element.  Neighbourhoods are rectangular sub-arrays 
> --  centred around a focal point — i.e., the array position that is determined by the stencil.
> --  For those array positions where the neighbourhood extends past the boundaries of the source
> --  array, a boundary condition determines the contents of the out-of-bounds neighbourhood
> --  positions.
> --
> stencil :: (Ix dim, Elem a, Elem b, Stencil dim a stencil)
>         => (stencil -> Exp b)          -- ^stencil function
>         -> Boundary a                  -- ^boundary condition
>         -> Acc (Array dim a)           -- ^source array
>         -> Acc (Array dim b)           -- ^destination array

As you discussed with Trevor, it is better to pull the backpermute out of the stencil operation to ensure that we can (a) implement blocked stencil computations and (b) use shared memory on GPUs.

Boundary conditions are just an enumeration as you suggested.

Concerning the focal point of a stencil and unused positions in the stencil neighbourhood, I came to the conclusion that there doesn't seem to be an immediate benefit to be clever about either.  Instead, I think, we can require the following:

* The focal point is always in the centre of the stencil (and we only use tuples with an odd number of arguments, so that the centre is uniquely defined).
* We don't identify unused positions in the stencil in the type at all.

These restrictions make everything considerably easier.  I think that this doesn't really restrict stencil application in any serious form.  If we need a 2x2 stencil with the focal point in the lower right corner, we just define the function as a 3x3 stencil and ignore the positions to the right and bottom of the centre; ie,

  stencil ( (a, b, _)
          , (c, d, _)   -- 'd' is the focal point
          , (_, _, _)
          ) = ....

Personally, I think the underscores introduce less clutter than a special notation to identify the focal point, and it is easier to implement in the frontend, too.

Let me know if you are having trouble with any of your examples given these restrictions.

> 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

Yes, that is what I had in mind, too, but I haven't actually looked at it.  I want to implement the simple version first.

The simpler stencil types allow us to predefine common stencil types, whose definition I append, together with the examples I had earlier (but using the new types).

Cheers,
Manuel

-=-

> -- DIM0
> --
> type Stencil0 a = ()
> 
> -- DIM1
> type Stencil1 a = Exp a
> type Stencil3 a = (Exp a, Exp a, Exp a)
> type Stencil5 a = (Exp a, Exp a, Exp a, Exp a, Exp a)
> 
> -- DIM2
> type Stencil3x3 a = (Stencil3 a, Stencil3 a, Stencil3 a)
> type Stencil5x3 a = (Stencil5 a, Stencil5 a, Stencil5 a)
> type Stencil3x5 a = (Stencil3 a, Stencil3 a, Stencil3 a, Stencil3 a, Stencil3 a)
> type Stencil5x5 a = (Stencil5 a, Stencil5 a, Stencil5 a, Stencil5 a, Stencil5 a)
> 
> -- DIM3
> type Stencil3x3x3 a = (Stencil3x3 a, Stencil3x3 a, Stencil3x3 a)
> type Stencil5x3x3 a = (Stencil5x3 a, Stencil5x3 a, Stencil5x3 a)
> type Stencil3x5x3 a = (Stencil3x5 a, Stencil3x5 a, Stencil3x5 a)
> type Stencil3x3x5 a = (Stencil3x3 a, Stencil3x3 a, Stencil3x3 a, Stencil3x3 a, Stencil3x3 a)
> type Stencil5x5x3 a = (Stencil5x5 a, Stencil5x5 a, Stencil5x5 a)
> type Stencil5x3x5 a = (Stencil5x3 a, Stencil5x3 a, Stencil5x3 a, Stencil5x3 a, Stencil5x3 a)
> type Stencil3x5x5 a = (Stencil3x5 a, Stencil3x5 a, Stencil3x5 a, Stencil3x5 a, Stencil3x5 a)
> type Stencil5x5x5 a = (Stencil5x5 a, Stencil5x5 a, Stencil5x5 a, Stencil5x5 a, Stencil5x5 a)

> stencil1D :: (Elem a, IsFloating a) 
>           => Stencil3 a -> Exp a
> stencil1D (x, y, z) = (x + z - 2 * y) / 2
> 
> stencil2D5 :: (Elem a, IsFloating a) 
>            => Stencil3x3 a -> Exp a
> stencil2D5 ( (_, t, _)
>            , (l, m, r)
>            , (_, b, _)
>            ) 
>            = (t + l + r + b - 4 * m) / 4
> 
> stencil2D :: (Elem a, IsFloating a) 
>           => Stencil3x3 a -> Exp a
> stencil2D ( (t1, t2, t3)
>           , (l , 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)
>           => Stencil3x3x3 a -> Exp a
> stencil3D (front, back, _) =      -- 'b4' is the focal point
>   let ((f1, f2, _),
>        (f3, f4, _),
>        _          ) = front
>       ((b1, b2, _),
>        (b3, b4, _),
>        _          ) = back
>   in
>   f1 + f2 + f3 + f4 + b1 + b2 + b3 + b4