compiler: Key the injective storage exemption on a zero distance

devito/ir/equations/equation.py

@@ -393,8 +393,15 @@ def __new__(cls, *args, **kwargs):
             input_expr = args[0]
             assert isinstance(input_expr, Eq)
             obj = LoweredEq(input_expr)
+            return ClusterizedEq.__new__(cls, obj, ispace=obj.ispace)
         elif len(args) == 2:
-            obj = LoweredEq(Eq(*args, evaluate=False))
+            # `lhs` and `rhs` are already lowered -- a pass built them out of
+            # IR objects -- so go straight to a ClusterizedEq. Lowering them
+            # again would re-apply the DSL-to-IR rewrites: a ConditionalDimension
+            # LHS would become `IntDiv(parent, factor)`, which no backend can
+            # assign to
+            return ClusterizedEq.__new__(cls, *args, ispace=IterationSpace([]),
+                                         conditionals={}, implicit_dims=(),
+                                         operation=None)
         else:
             raise ValueError(f"Cannot construct DummyEq from args={str(args)}")
-        return ClusterizedEq.__new__(cls, obj, ispace=obj.ispace)

devito/ir/support/basic.py

@@ -869,15 +869,17 @@ def is_storage_related(self, dims=None):
         """
         for d in self.findices:
             if d._defines & set(as_tuple(dims)):  # noqa: SIM102
-                # An injective Dimension maps distinct iterations to distinct
-                # locations by construction -- e.g. a factor-based
-                # ConditionalDimension, whose guard admits only
-                # `parent % factor == 0` and whose `i -> i//factor` is
-                # one-to-one on those. Unlike, say, a SteppingDimension, where
-                # `time % 2` deliberately revisits the same slot
-                nonlinear = [i for i in d._defines
-                             if i.is_NonlinearDerived and not i.is_Injective]
-                if nonlinear or self.is_const(d):
+                if self.is_const(d):
+                    return True
+                for i in d._defines:
+                    if not i.is_NonlinearDerived:
+                        continue
+                    # An injective Dimension maps distinct iterations to
+                    # distinct locations, so only the *same* iteration can
+                    # reach the same location -- unlike a SteppingDimension,
+                    # where `time % 2` has distinct iterations revisit a slot
+                    if i.is_Injective and not self.distance_mapper.get(i):
+                        continue
                     return True
         return False
 

tests/test_ir.py

@@ -10,7 +10,7 @@
 from devito.finite_differences.differentiable import IndexSum, LocalSum
 from devito.ir.cgen import ccode
 from devito.ir.clusters import Cluster, ClusterGroup
-from devito.ir.equations import LoweredEq
+from devito.ir.equations import DummyEq, LoweredEq
 from devito.ir.equations.algorithms import dimension_sort
 from devito.ir.iet import FindNodes, Iteration
 from devito.ir.stree import stree_build
@@ -1575,6 +1575,25 @@ def test_iteration_parallelism_3d(self, exprs, atomic, parallel):
 
 class TestEquationAlgorithms:
 
+    def test_dummyeq_does_not_relower(self):
+        """
+        `DummyEq(lhs, rhs)` builds already-lowered IR, so it must not re-apply
+        the DSL-to-IR rewrites. Re-lowering turns a ConditionalDimension LHS
+        into `IntDiv(parent, factor)`, which is not assignable: a backend that
+        emits an index assignment -- as the SYCL one does for the parallel_for
+        indices -- then produces `IntDiv(x, ctf) = ...` and trips over
+        `LoweredEq.dtype`.
+        """
+        grid = Grid(shape=(16,))
+        x, = grid.dimensions
+        cd = ConditionalDimension(name='ct', parent=x, factor=4)
+
+        for d in (cd, x):
+            eq = DummyEq(d, 7)
+            assert eq.lhs is d
+            # A valid assignment target carries a dtype
+            assert eq.dtype is not None
+
     @pytest.mark.parametrize('expr,expected', [
         ('Eq(a[time, p], b[time, c[p, 0]+r, c[p, 1]] * f[p, r])', '[time, p, r, d]')
     ])

tests/test_subsampling.py

@@ -5,8 +5,8 @@
 from sympy import Or
 
 from devito import (
-    CondEq, ConditionalDimension, Constant, Eq, Function, Grid, Inc, Operator, SubGrid,
-    TimeFunction, switchconfig
+    CondEq, ConditionalDimension, Constant, Eq, Function, Grid, Inc, Operator,
+    SpaceDimension, SubGrid, TimeFunction, switchconfig
 )
 from devito.ir.iet import Conditional, Expression, FindNodes, Iteration
 from devito.symbolics import retrieve_indexed
@@ -218,6 +218,43 @@ class TestNumerics:
 
     """Results must match an independent reference."""
 
+    def test_or_relation_stays_fused(self):
+        """
+        Two `Or`-related guards writing overlapping slots of the same Function
+        must stay in one loop. Splitting them into separate loops reorders the
+        writes: `f[ci]` and `f[ci+1]` alias, so the second loop would overwrite
+        what the first produced instead of interleaving per iteration.
+
+        Injectivity says distinct iterations map to distinct indices, so only
+        the *same* iteration can reach the same location. It does not say two
+        different offsets off that Dimension cannot collide, which is why the
+        exemption is keyed on a zero distance.
+        """
+        size, factor = 18, 4
+        i = SpaceDimension(name='i')
+        grid = Grid(shape=(size,), dimensions=(i,))
+
+        ci = ConditionalDimension(name='ci', parent=i, factor=factor,
+                                  relation=Or)
+        cend = ConditionalDimension(name='cend', parent=i,
+                                    condition=CondEq(i, size - 1), relation=Or)
+
+        g = Function(name='g', grid=grid, shape=(size,), dimensions=(i,))
+        f = Function(name='f', grid=grid, shape=(size//factor + 2,),
+                     dimensions=(ci,))
+
+        op = Operator([Eq(g, i), Eq(f, g),
+                       Eq(f[ci + 1], g, implicit_dims=cend)])
+        op.apply()
+
+        # One loop, not one per guard
+        assert len(FindNodes(Iteration).visit(op)) == 1
+
+        assert np.array_equal(g.data, np.arange(size, dtype=g.dtype))
+        for j in range(size//factor + 1):
+            assert f.data[j] == g.data[j*factor]
+        assert f.data[-1] == g.data[-1]
+
     @pytest.mark.parametrize('shape,factor,ndim', [
         ((12,), 2, 1), ((12,), 3, 1), ((13,), 4, 1),
         ((12, 12), 2, 2), ((12, 12), 3, 2), ((11, 9), 4, 2),