Re: [pooma-dev] missing reductions

[Stephen A. Smith]

The comment refers to the error term, and we actually do
have reductions now.  (By reductions, we mean getting
the sum of elements in an array, the maximum element
in an array, etc.)  The loop could now be replaced with:

    error_m(I,J) = -b_m(I,J) + e0 * x0_m(I,J) -
                e1 * (x0_m(I-1,J) + x0_m(I+1,J) +
                             x0_m(I,J+1) + x0_m(I,J-1)) -
                e2 * (x0_m(I-1,J-1) + x0_m(I+1,J-1) +
                             x0_m(I-1,J+1) + x0_m(I+1,J+1));

     err = sum(error_m*error_m);

provided the array error_m is added as a member.

We can implement the stencil object you suggest as
well, and this code actually is intended to illustrate
the advantage of doing so.  The code in the class
JacobiInP2UnOpt should be rewritten as above,
while the code in the class JacobiInP2Opt should be
written as:

    WeightedAverage wa(e0, e1, e2);
    Stencil<WeightedAverage> average(wa);

    error_m(I,J) = -b_m(I,J) + average(x0_m);

     err = sum(error_m*error_m);

if you provide a WeightedAverage object like the
NinePoint object already there.

I agree that it would be good to update the
examples and benchmarks to use the best
features of POOMA.  It might worth asking,
if you're looking at the solver benchmarks,
if all those solver variations are really useful.

Eventually, POOMA should provide more predefined
stencil objects.   Since they're easy to write, and
Blanca likes to write their own stencils anyway,
the contract should probably focus on the functionality
we provide through the stencil engines and the evaluator.
For example, one feature that has yet to be implemented
is multi-argument stencils for fields and arrays, so you
could write code for non-constant weights:

x = weightedAverage(w0, w1, y);

to implement:

x(i) = w1(i) * y(i -1) + w0(i) * y(i) + w1(i) * y(i + 1);

or similar things.  This may not be high-priority, however.
I think John Hall told me this was not needed for their
current code.  In the long term, they will need something
like this functionality under the hood to represent differential
operators on non-uniform meshes, where the coefficients are
computed from arrays instead of constants.

    Stephen


Allan Stokes wrote:

> In benchmarks/Solvers/Jacobi/JacobiInP2.h there is this comment:
>
>   // Since we don't have reductions yet, do this by hand.
>
> The hand-wrought code looks like a simple nine-point stencil reduction with
> multiplicative co-efficients for each cell.  It doesn't look like something
> that Pooma couldn't implement.
>
> Is there a technical reason this operation is not yet provided?
>
> Allan