compiler: Fix par-tile unfolding behaviour

devito/finite_differences/differentiable.py

@@ -1111,6 +1111,58 @@ def dimension(self):
 
     weights = Array.initvalue
 
+    @cached_property
+    def symmetry(self):
+        """
+        The sign relating mirrored coefficients: `1` for symmetric weights,
+        `-1` for antisymmetric weights, or `None` otherwise.
+
+        An odd antisymmetric sequence must have a zero center coefficient.
+        """
+        weights = self.weights
+        if not weights:
+            return None
+
+        middle = len(weights) // 2
+        pairs = tuple(zip(weights[:middle],
+                          reversed(weights[middle + len(weights) % 2:]),
+                          strict=True))
+
+        if all(left.equals(right) is True for left, right in pairs):
+            return sympy.S.One
+        if all(left.equals(-right) is True for left, right in pairs) and \
+           (len(weights) % 2 == 0 or weights[middle].equals(0) is True):
+            return sympy.S.NegativeOne
+
+        return None
+
+    @cached_property
+    def coefficient_span(self):
+        """
+        Paired weights and their offset within the original stencil.
+
+        Symmetric or antisymmetric weights retain their full span, including
+        symmetric zero padding. Otherwise, trim known zeros at both ends and
+        return rebuilt weights together with their starting offset. The
+        resulting weights may still have no symmetry.
+        """
+        if self.symmetry is not None:
+            return self, 0
+
+        weights = self.weights
+        first, last = 0, len(weights)
+        while first < last and weights[first].is_zero is True:
+            first += 1
+        while last > first and weights[last - 1].is_zero is True:
+            last -= 1
+
+        if first == 0 and last == len(weights):
+            return self, 0
+
+        d = self.dimension
+        dimension = d._rebuild(_min=d._min + first, _max=d._min + last - 1)
+        return self._rebuild(dimensions=dimension, initvalue=weights[first:last]), first
+
     def _xreplace(self, rule):
         if self in rule:
             return rule[self], True

devito/passes/clusters/blocking.py

@@ -541,7 +541,6 @@ def schedule(self, dims, clusters):
 
         else:
             umt = self.umt
-            unbound = self.unbound and umt.is_multi
 
         umt.iter()
 

devito/symbolics/unevaluation.py

@@ -1,6 +1,6 @@
 import sympy
 
-__all__ = ['Add', 'Mul', 'Pow']
+__all__ = ['Add', 'Mod', 'Mul', 'Pow']
 
 
 class UnevaluableMixin:
@@ -13,6 +13,10 @@ class Add(sympy.Add, UnevaluableMixin):
     __new__ = UnevaluableMixin.__new__
 
 
+class Mod(sympy.Mod, UnevaluableMixin):
+    __new__ = UnevaluableMixin.__new__
+
+
 class Mul(sympy.Mul, UnevaluableMixin):
     __new__ = UnevaluableMixin.__new__
 

tests/test_derivatives.py

@@ -941,6 +941,54 @@ def test_unexpand_space_interp_w_saved_timefunc(self):
         assert_structure(op, ['t,x,y,z', 't,x,y,z,i1', 't,x,y,z,i1,i0'])
 
 
+class TestWeights:
+
+    @pytest.mark.parametrize('coeffs,expected', [
+        ((1, 2, -2, -1), -1),
+        ((1, 2, 0, -2, -1), -1),
+        ((1, 2, 2, 1), 1),
+        ((1, 2, 3, 2, 1), 1),
+        ((1, 2, 3, -2, -1), None),
+        ((0, 0, 0), 1),
+        ((0,), 1),
+        ((3,), 1),
+        ((Symbol('a'), 0, -Symbol('a')), -1),
+        ((Symbol('a'), Symbol('b'), Symbol('a')), 1),
+        ((Symbol('a'), Symbol('b')), None),
+    ])
+    def test_symmetry(self, coeffs, expected):
+        i = StencilDimension('i', 0, len(coeffs) - 1)
+        weights = Weights(name='w', dimensions=i, initvalue=coeffs)
+
+        assert weights.symmetry == expected
+        assert weights._rebuild().symmetry == expected
+
+    @pytest.mark.parametrize('coeffs,expected,offset,sign', [
+        ((0, 1, 2, -2, -1), (1, 2, -2, -1), 1, -1),
+        ((1, 2, -2, -1, 0), (1, 2, -2, -1), 0, -1),
+        ((0, 0, 1, 2, 2, 1, 0), (1, 2, 2, 1), 2, 1),
+        ((0, 1, 2, -2, -1, 0), (0, 1, 2, -2, -1, 0), 0, -1),
+        ((0, 1, 2, 3, 0), (1, 2, 3), 1, None),
+        ((1, 2, 3), (1, 2, 3), 0, None),
+        ((0, 0, 0), (0, 0, 0), 0, 1),
+        ((0, Symbol('a'), -Symbol('a')), (Symbol('a'), -Symbol('a')), 1, -1),
+    ])
+    def test_coefficient_span(self, coeffs, expected, offset, sign):
+        i = StencilDimension('i', -3, len(coeffs) - 4)
+        weights = Weights(name='w', dimensions=i, initvalue=coeffs)
+        span, start = weights.coefficient_span
+
+        assert span.weights == expected
+        assert span.symmetry == sign
+        assert start == offset
+        assert span.dimension._min == i._min + offset
+        assert span.dimension._size == len(expected)
+        assert weights.coefficient_span[0] is span
+        assert weights._rebuild().coefficient_span == (span, start)
+        if coeffs == expected:
+            assert span is weights
+
+
 class TestTwoStageEvaluation:
 
     def test_exceptions(self):

tests/test_gpu_openacc.py

@@ -144,11 +144,10 @@ def test_multiple_tile_sizes(self, par_tile):
         assert trees[3][1].pragmas[0].ccode.value ==\
             f'acc parallel loop {sclause} present(src,src_gp,src_wx,src_wy,src_wz,u)'
 
-    def test_short_multi_tile_keeps_outer_dim_blocked(self):
+    def test_short_multi_tile_keeps_outer_dim_on_device(self):
         """
-        A multi `par-tile` entry shorter than the nest it lands on must not cost
-        the outermost Dimension its BlockDimension: on a device, dropping it
-        would leave `x` iterated outside the offloaded nest.
+        A short multi `par-tile` leaves the outermost Dimension unblocked.
+        Its full loop must stay inside the OpenACC region.
         """
         grid = Grid(shape=(8, 8, 8))
 
@@ -166,9 +165,21 @@ def test_short_multi_tile_keeps_outer_dim_blocked(self):
                           'advanced',
                           {'par-tile': par_tile, 'blocklevels': 1, 'blockinner': True}))
 
-        bns, _ = assert_blocking(op, {'x0_blk0', 'x1_blk0'})
+        bns, _ = assert_blocking(op, {'x0_blk0', 'y1_blk0'})
 
-        expected = ((4, 4, 32), (4, 4, 16))
+        trees = retrieve_iteration_tree(op)
+        assert len(trees) == 2
+        assert all(tree[0].dim is grid.time_dim for tree in trees)
+        assert trees[0][1].pragmas[0].ccode.value ==\
+            'acc parallel loop tile(32,4,4) present(u)'
+
+        x = grid.dimensions[0]
+        assert trees[1][1].dim is x
+        assert trees[1][1].limits == (x.symbolic_min, x.symbolic_max, 1)
+        assert trees[1][1].pragmas[0].ccode.value ==\
+            'acc parallel loop tile(16,4,4) present(u,v)'
+
+        expected = ((4, 4, 32), (4, 16))
         for root, v in zip(bns.values(), expected, strict=True):
             iters = FindNodes(Iteration).visit(root)
             iters = [i for i in iters if i.dim.is_Block and i.dim._depth == 1]

tests/test_pickle.py

@@ -22,6 +22,7 @@
     CallFromPointer, Cast, DefFunction, FieldFromPointer, IntDiv, ListInitializer, SizeOf,
     indexify, pow_to_mul
 )
+from devito.symbolics.unevaluation import Mod
 from devito.tools import EnrichedTuple
 from devito.types import (
     Array, ComponentAccess, CustomDimension, DefaultDimension, DeviceID, FIndexed,
@@ -720,6 +721,15 @@ def test_local_sum(self, pickle):
             tuple(map(str, reduction.conditionals.values()))
         assert rebuilt.free_symbols.isdisjoint(rebuilt.bound_symbols)
 
+    @pytest.mark.parametrize('args', [(0, 3), (1, 3), (3, 3), (Symbol('i'), 3)])
+    def test_unevaluated_mod(self, pickle, args):
+        expr = Mod(*args)
+
+        rebuilt = pickle.loads(pickle.dumps(expr))
+
+        assert rebuilt == expr
+        assert rebuilt.args == args
+
 
 class TestAdvanced: