From 88750b3f84bd7f76b168550186d820d7d478124c Mon Sep 17 00:00:00 2001 From: mloubout Date: Mon, 26 Oct 2026 14:21:12 -0400 Subject: [PATCH 1/2] compiler: Key the injective storage exemption on a zero distance An injective Dimension was exempted from `is_storage_related` outright. That is too strong: injectivity says distinct iterations map to distinct indices, so only the *same* iteration can reach a given location -- it says nothing about two different offsets off that Dimension, which do collide. `is_storage_related` also gates fusion and lifting, so the blanket exemption split `f[ci]` and `f[ci+1]` into separate loops, and the second overwrote what the first produced instead of interleaving per iteration: Eq(g, i), Eq(f, g), Eq(f[ci+1], g, implicit_dims=cend) # f == [0, 0, 4, 8, 12, 16], one slot off throughout Require a zero distance for the exemption, which is exactly the subsampled `Inc` it was added for, and leave every other dependence as before. --- devito/ir/support/basic.py | 20 ++++++++++--------- tests/test_subsampling.py | 41 ++++++++++++++++++++++++++++++++++++-- 2 files changed, 50 insertions(+), 11 deletions(-) diff --git a/devito/ir/support/basic.py b/devito/ir/support/basic.py index 326562b5c9..0a0d092a55 100644 --- a/devito/ir/support/basic.py +++ b/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 diff --git a/tests/test_subsampling.py b/tests/test_subsampling.py index 07ff1a5a1f..033237c193 100644 --- a/tests/test_subsampling.py +++ b/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), From b45f9fcca19b8ff83f41309729dfad2c3cbafd82 Mon Sep 17 00:00:00 2001 From: mloubout Date: Mon, 26 Oct 2026 16:32:21 -0400 Subject: [PATCH 2/2] compiler: Don't re-lower an already-lowered DummyEq `DummyEq(lhs, rhs)` wrapped its arguments in a DSL `Eq` and pushed them back through `LoweredEq`'s lowering path. But its callers are compiler passes holding already-lowered IR, so the DSL-to-IR rewrites ran a second time. For most LHSs that is a no-op, which is why it went unnoticed. A ConditionalDimension LHS, though, becomes `IntDiv(parent, factor)` -- an expression, not an assignment target. A backend that emits an index assignment then produces IntDiv(time, ctf) = ct_m + (int)idx[0] and `LoweredEq.dtype` raises `AttributeError: 'IntDiv' object has no attribute 'dtype'`. The SYCL parallelizer does exactly this for its parallel_for indices, where every other backend prints the loop variable by name and so never looks at the lowered form. Construct the ClusterizedEq directly instead, as the class docstring already describes. `LoweredEq`'s own two-argument form bypasses lowering for the same reason. --- devito/ir/equations/equation.py | 11 +++++++++-- tests/test_ir.py | 21 ++++++++++++++++++++- 2 files changed, 29 insertions(+), 3 deletions(-) diff --git a/devito/ir/equations/equation.py b/devito/ir/equations/equation.py index 3139cb3354..f84e24c58c 100644 --- a/devito/ir/equations/equation.py +++ b/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) diff --git a/tests/test_ir.py b/tests/test_ir.py index 9a7b2e9fea..f4fd0d2635 100644 --- a/tests/test_ir.py +++ b/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]') ])