BZDEV: Re: Thoughts about stencils

[Steve K. <kar...@la...>]

Something like this can definitely work, and will have all the important
benefits that Todd says.  There are a bunch of possible tweaks to the user
interface.

For example, I would like to be able to do the Laplacian operator as a
stencil object and use it in larger expressions like:

	P3 = 2 * P2 + c * Laplacian(P2) - P1;

where now this is one data parallel expression, but the Laplacian operator
is a stencil which stitches itself into the expression template.  This
would get a bunch of the benefits Todd mentions.  This means that the
stencil wouldn't have the equals inside it.  I consider that to be an
advantage, actually.

Composition would also be important with these.  I'd want to be able to
construct an Acoustic stencil using Laplacian and then say:

	P3 = Acoustic(c,P1,P2);

How would something like this actually get put together?  Laplacian and
Acoustic would have to return someting that the rest of the expression
template machinery understands.  In POOMA II that means it would return an
Array with some engine.

I think the way to construct this is built up in much the same way that
Todd shows: some little gadget that knows it is operating elementwise, and
another that takes that elementwise thing and applies it all over the place.

The elementwise gadget would look very much like Todd's:

template<class A>
A::Element_t
laplace(const A& a) 
{
  return a(-1,0,0)+a(1,0,0)+a(0,1,0)+a(0,-1,0)+a(0,0,1)+a(0,0,-1)-6*a(0,0,0);
}

You could the apply this using an apply_stencil gadget:

  P3 = 2 * P2 + c * apply_stencil(laplace,P2) - P1;

Then you could write a function like:

template<class T, class Engine>
Array<3,T,StencilEngine<3,T,Engine,laplace> >
laplace(const Array<3,T,Engine>& a)
{
  return apply_stencil(laplace,a);
}

Then you would say

  P3 = 2 * P2 + c * laplace(P2) - P1;

Partial ordering would get the right laplace called at the right place.

These are then composable.  You could build the elementwise acoustic
stencil by saying

template<class A, class B, class C>
A::Element_t
acoustic(const A& a, const B& b, const C& c) 
{
  return 2 * b + c * laplace(b) - a;
}

Then you could build a data parallel 'acoustic' operator the same way as
with Laplace.

Unfortunately I don't see how to compose these things at the data parallel
level with something like:

template<class T1, class E1, class T2, class E2, class C>
Array<3,T,????>
acoustic(const Array<3,T1,E1>& a, const Array<3,T2,E2>& b, const C& c) 
{
  return 2 * b + c * laplace(b) - a;
}

because of the gnarly return engine type.

For simple stencils like laplace it would certainly be possible to give it
some sort of proxy class which would go through the expression and return
the footprint of the stencil.  More complex stencils could have
conditionals inside which vary the footprint depending on run-time info.  I
think that means that the user would have to indicate in some way the
footprint.  There are a couple of ways that could be handled.

One is to explicitly say something inside the function:

template<class A>
A::Element_t
laplace(const A& a) 
{
  Footprint(a,Interval<1>(-1,1),Interval<1>(-1,1));
  return a(-1,0,0)+a(1,0,0)+a(0,1,0)+a(0,-1,0)+a(0,0,1)+a(0,0,-1)-6*a(0,0,0);
}

Then the expression itself wouldn't do anything to a, but the Footprint
statement would record it.

Another way would be to make the elementwise operator a functor instead of
a function.  Then other member functions could be used to interrogate it
for the footprint.  I tend to lean that direction, also because plugging
the functor into the StencilEngine is easier if it is a class than if it is
a function.

-Steve Karmesin

======================================================================
Pooma Team Leader
kar...@la...
phone: (505)665-6019
fax  : (505)667-4939
For PGP public key, send message with subject line: SEND PUBLIC KEY
======================================================================
--------------------- blitz-dev list --------------------------------
* To subscribe/unsubscribe: mail to maj...@oo..., with
"subscribe blitz-dev" or "unsubscribe blitz-dev" in the body of the message
* Blitz++ web page: http://oonumerics.org/blitz/