Conditionals not dealt with correctly

[wence-]

Consider the following:

from firedrake import *
parameters["coffee"]["O2"] = True
mesh = UnitIntervalMesh(4)
V = FunctionSpace(mesh, "DG", 0)
u = Function(V)
du = TrialFunction(V)
v = TestFunction(V)
bhp = Constant(2)
u.dat.data[...] = range(10)

Jc = conditional(ge(u - bhp, 0.0), 1.0, 0.0)*du*v*dx
Jm = conditional(gt(u - bhp, 0.0), 1.0, 0.0)*du*v*dx
print assemble(Jc).M.values.diagonal()
print assemble(Jm).M.values.diagonal()

I expect to see:

[ 0.    0.    0.25  0.25]
[ 0.    0.    0.    0.25]

and do so with coffee off.
When on, I instead see:

[ 0.25  0.25  0.25  0.25]
[ 0.  0.  0.  0.]

FWIW, I suspect this is somehow a problem that some of the conditional stuff in ffc is built as strings, rather than ast nodes.

[wence-]

A little bit more digging (using https://bitbucket.org/mapdes/ffc/pull-request/78/create-correct-nodes-for-max_value-and/diff) suggests this is a problem with the loop merger step of the rewriter. After licm and expression factorisation, you do:

            lm = SSALoopMerger(self.header, ew.expr_graph)
            merged_loops = lm.merge()
            for merged, merged_in in merged_loops:
                [self.hoisted.update_loop(l, merged_in) for l in merged]
            lm.simplify()

However, the expr_graph only tracks dependencies in the loop being optimised. The resulting code ends up doing:

double F0 = 0;
double F1 = 0;

double C0 = some_expression(F0, F1);

for ( .... ) {
   F0 += ...;
   F1 += ...;
}

// quad loop here

i.e. the read-after-write dependency between C and F0/F1 is missing. If I apply this patch in coffee:

diff --git a/coffee/optimizer.py b/coffee/optimizer.py
index a744f7e..3dd87d4 100644
--- a/coffee/optimizer.py
+++ b/coffee/optimizer.py
@@ -58,7 +58,7 @@ class LoopOptimizer(object):
         # Track nonzero regions accessed in the loop nest
         self.nonzero_info = {}
         # Track data dependencies
-        self.expr_graph = ExpressionGraph(loop)
+        self.expr_graph = ExpressionGraph(self.header)
         # Track hoisted expressions
         self.hoisted = StmtTracker()

then I get the right answer. But is this the correct thing to do?

[FabioLuporini]

I can't reproduce the error because of this:

u.dat.data[...] = range(10)

did you mean:

u.dat._data[:] = range(4)

or presumably larger...

As for your patch: it is actually probably the right thing to do to start tracking from the header. Let's make the example run on my machine, and then I can tell definitely what we should do.

[wence-]

u.dat.data[...] = range(4) is correct (not 10)

[FabioLuporini]

This seems to be the right fix. Plus, this also fixes firedrakeproject/firedrake#525.

I'm running the test suite right now. If they all pass, I'll just add this commit to the PR I'm about to make.

[wence-]

Please don't. It is a standalone fix that should be merged separately

[FabioLuporini]

So the patch works, for me we can merge it even now. If we do so, let's close both this issue and firedrakeproject/firedrake#525

[FabioLuporini]

Shall we hit merge?

[kynan]

Go ahead I think.

[FabioLuporini]

PRs merged, closing