How should domain trees be deduced/specified?

inducer commented 4 years ago

@isuruf showed this kernel, which can be reproduced approximately by running

 SUMPY_NO_CACHE=1 pycl -m pudb test_fmm.py 'test_sumpy_fmm(cl._csc, LaplaceKernel(2), VolumeTaylorLocalExpansion,VolumeTaylorMultipoleExpansion)'

in https://github.com/inducer/sumpy.

The if is clearly emitted prematurely: It should only cover statements in isrc, but here it seems to cover the accumulator initialization and the result store.

inducer commented 4 years ago

I think I've made some progress here. I've got a workaround and some design questions to make this better going forward.

The domains in this kernel are given as:

            "{[itgt_box]: 0 <= itgt_box < ntgt_boxes}",
            "{[isrc_box]: isrc_box_start <= isrc_box < isrc_box_end}",
            "{[itgt, isrc, idim]: \
                itgt_start <= itgt < itgt_end and \
                isrc_start <= isrc < isrc_end and \
                0 <= idim < dim}",

The above domains are also responsible for the placement of the conditional: As soon as the itgt loop is entered, that last domain gets to enforce its constraints. Since the accumulator is inside itgt, and since that domain (projected onto itgt) is empty if isrc_start == isrc_end, Loopy is in principle correct in emitting the conditional where it does.

So in this case, the kernel is incorrect and should be fixed.

The natural fix (just make separate domains for itgt and isrc):

[ntgt_boxes] -> { [itgt_box] : 0 <= itgt_box < ntgt_boxes }                                                             
  [isrc_box_end, isrc_box_start] -> { [isrc_box] : isrc_box_start <= isrc_box < isrc_box_end }                          
  [itgt_end, itgt_start] -> { [itgt] : itgt_start <= itgt < itgt_end }                                                  
  [isrc_end, isrc_start] -> { [isrc] : isrc_start <= isrc < isrc_end }                                                  
{ [idim, idim_0] : 0 <= idim <= 1 and 0 <= idim_0 <= 1 }

fails because Loopy gets the domain tree wrong and then complains about needing to branch the domain tree in order to satisfy the inames needed by various statements.

Helping it along by (kind of unnaturally) forcing the right nesting of domains produces the desired kernel without the conditional:

            "{[itgt_box]: 0 <= itgt_box < ntgt_boxes}",
            "{[isrc_box]: isrc_box_start <= isrc_box < isrc_box_end}",
            "[itgt_start,itgt_end,isrc_box] -> {[itgt]:  itgt_start <= itgt < itgt_end}",
            "[isrc_start, isrc_end, itgt] -> {[isrc]: isrc_start <= isrc < isrc_end}",
            "{[idim]: 0 <= idim < dim}"

This works, but it's stupid, which begs the question of what we should do instead. Some thoughts on that follow.

One could ask whether the itgt loop should be inside or outside the isrc_box loop. As far as I can tell, both orderings are justifiable. Bonus points for a solution that makes it so that the "math-level" input kernel doesn't (have to) prescribe this nesting.
Is the current data structure of a "tree" of domains appropriate? Why can't a domain have multiple parents?
The "parent-child" relationship doesn't seem totally flawed to me: A domain is a child of a parent if it has parameters that are inames defined by or data fetched in loops controlled by the parent.
Based on this definition, it seems that we can (as we currently do) determine parent-child relationships automatically. (though we need to do better, as illustrated above) But are there situations where this fails (i.e. where only user input can help us find a parent-child relationship)?
I think we can be a bit more liberal with "branching", however a free-for-all can't work. Consider the domains
```
A: {[i]: 0<=i<10}
B: {[j]: 0<=i<j and i > 5}
C: {[k]: 0<=i<j and i <= 5}
```
B and C are naturally both children of A. The A & B and A & C intersections both work nicely, however if someone were to try a A & B & C nesting (i.e. a statement inside i, j, and k), they would get no statement instances at all. A potential test for whether this is OK might be whether A & B projected onto the domain of just A equals A (and likewise for C), i.e. neither child "substantially modifies" A.

cc'ing the semantics squad @mattwala @jdsteve2 and @kaushikcfd for general Loopy interest.

mattwala commented 4 years ago

inducer / loopy

How should domain trees be deduced/specified? #127