compiler: Patch race conditions due to storage-related dependencies

devito/ir/support/basic.py

@@ -100,8 +100,6 @@ def aindices(self):
                 retval.append(dims.pop())
             elif isinstance(i, Dimension):
                 retval.append(i)
-            elif q_constant(i):
-                retval.append(fi)
             else:
                 retval.append(None)
         return DimensionTuple(*retval, getters=self.findices)
@@ -262,10 +260,14 @@ def is_regular(self):
         # space Dimensions
         positions = []
         for d in self.aindices:
-            for n, i in enumerate(self.intervals):
-                if i.dim._defines & d._defines:
-                    positions.append(n)
-                    break
+            try:
+                for n, i in enumerate(self.intervals):
+                    if i.dim._defines & d._defines:
+                        positions.append(n)
+                        break
+            except AttributeError:
+                # `d is None` due to e.g. constant access
+                continue
         return positions == sorted(positions)
 
     def __lt__(self, other):
@@ -548,6 +550,15 @@ def is_cross(self):
     def is_local(self):
         return self.function.is_Symbol
 
+    @memoized_meth
+    def is_const(self, dim):
+        """
+        True if a constant depedence, that is no Dimensions involved, False otherwise.
+        """
+        return (self.source.aindices[dim] is None and
+                self.sink.aindices[dim] is None and
+                self.distance_mapper[dim] == 0)
+
     @memoized_meth
     def is_carried(self, dim=None):
         """Return True if definitely a dimension-carried dependence, False otherwise."""
@@ -623,9 +634,10 @@ def is_storage_related(self, dims=None):
         cause the access of the same memory location, False otherwise.
         """
         for d in self.findices:
-            if (d._defines & set(as_tuple(dims)) and
-                    any(i.is_NonlinearDerived for i in d._defines)):
-                return True
+            if d._defines & set(as_tuple(dims)):
+                if any(i.is_NonlinearDerived for i in d._defines) or \
+                   self.is_const(d):
+                    return True
         return False
 
 

devito/mpi/halo_scheme.py

@@ -373,7 +373,7 @@ def classify(exprs, ispace):
                         v[(d, LEFT)] = STENCIL
                         v[(d, RIGHT)] = STENCIL
                 else:
-                    v[(d, i.aindices[d])] = NONE
+                    v[(d, i[d])] = NONE
 
             # Does `i` actually require a halo exchange?
             if not any(hl is STENCIL for hl in v.values()):
@@ -426,7 +426,12 @@ def classify(exprs, ispace):
                 func = Max
             candidates = [i for i in aindices if not is_integer(i)]
             candidates = {(i.origin if d.is_Stepping else i) - d: i for i in candidates}
-            loc_indices[d] = candidates[func(*candidates.keys())]
+            try:
+                loc_indices[d] = candidates[func(*candidates.keys())]
+            except KeyError:
+                # E.g., `aindices = [0, 1, d+1]` -- it doesn't really matter
+                # what we put here, so we place 0 as it's the old behaviour
+                loc_indices[d] = 0
 
         mapper[f] = HaloSchemeEntry(frozendict(loc_indices), frozenset(halos))
 

devito/passes/clusters/asynchrony.py

@@ -91,7 +91,7 @@ def callback(self, clusters, prefix):
                         else:
                             # Functions over non-stepping Dimensions need no lock
                             continue
-                    except KeyError:
+                    except (AttributeError, KeyError):
                         # Would degenerate to a scalar, but we rather use a lock
                         # of size 1 for simplicity
                         ld = CustomDimension(name='ld', symbolic_size=1)

devito/passes/clusters/buffering.py

@@ -5,7 +5,8 @@
 import numpy as np
 
 from devito.ir import (Cluster, Forward, GuardBound, Interval, IntervalGroup,
-                       IterationSpace, PARALLEL, Queue, Vector, lower_exprs, vmax, vmin)
+                       IterationSpace, PARALLEL, Queue, SEQUENTIAL, Vector,
+                       lower_exprs, normalize_properties, vmax, vmin)
 from devito.exceptions import InvalidOperator
 from devito.logger import warning
 from devito.symbolics import retrieve_function_carriers, uxreplace
@@ -207,7 +208,16 @@ def callback(self, clusters, prefix, cache=None):
                 expr = lower_exprs(uxreplace(Eq(lhs, rhs), b.subdims_mapper))
                 ispace = b.written
 
-                processed.append(c.rebuild(exprs=expr, ispace=ispace))
+                # Buffering creates a storage-related dependence along the
+                # contracted dimensions
+                properties = dict(c.properties)
+                for d in b.contraction_mapper:
+                    d = ispace[d].dim  # E.g., `time_sub -> time`
+                    properties[d] = normalize_properties(properties[d], {SEQUENTIAL})
+
+                processed.append(
+                    c.rebuild(exprs=expr, ispace=ispace, properties=properties)
+                )
 
             # Substitute buffered Functions with the newly created buffers
             exprs = [uxreplace(e, subs) for e in c.exprs]
@@ -233,7 +243,16 @@ def callback(self, clusters, prefix, cache=None):
                 expr = lower_exprs(uxreplace(Eq(lhs, rhs), b.subdims_mapper))
                 ispace = b.written
 
-                processed.append(c.rebuild(exprs=expr, ispace=ispace))
+                # Buffering creates a storage-related dependence along the
+                # contracted dimensions
+                properties = dict(c.properties)
+                for d in b.contraction_mapper:
+                    d = ispace[d].dim  # E.g., `time_sub -> time`
+                    properties[d] = normalize_properties(properties[d], {SEQUENTIAL})
+
+                processed.append(
+                    c.rebuild(exprs=expr, ispace=ispace, properties=properties)
+                )
 
         return processed
 

tests/test_buffering.py

@@ -603,6 +603,25 @@ def test_multi_access():
     assert np.all(w.data == w1.data)
 
 
+def test_issue_1901():
+    grid = Grid(shape=(2, 2))
+    time = grid.time_dim
+    x, y = grid.dimensions
+
+    usave = TimeFunction(name='usave', grid=grid, save=10)
+    v = TimeFunction(name='v', grid=grid)
+
+    eq = [Eq(v[time, x, y], usave)]
+
+    op = Operator(eq, opt='buffering')
+
+    trees = retrieve_iteration_tree(op)
+    assert len(trees) == 2
+    assert trees[1].root.dim is time
+    assert not trees[1].root.is_Parallel
+    assert trees[1].root.is_Sequential  # Obv
+
+
 def test_everything():
     nt = 50
     grid = Grid(shape=(6, 6))

tests/test_dle.py

@@ -667,7 +667,7 @@ def test_incs_no_atomic(self):
     @pytest.mark.parametrize('exprs,simd_level,expected', [
         (['Eq(y.symbolic_max, g[0, x], implicit_dims=(t, x))',
          'Inc(h1[0, 0], 1, implicit_dims=(t, x, y))'],
-         2, [6, 0, 0]),
+         None, [6, 0, 0]),
         (['Eq(y.symbolic_max, g[0, x], implicit_dims=(t, x))',
          'Eq(h1[0, y], y, implicit_dims=(t, x, y))'],  # 1695
          2, [0, 1, 2]),

tests/test_ir.py

@@ -707,7 +707,7 @@ def test_multiple_eqs(self, exprs, expected, ti0, ti1, ti3, fa):
         assert len(expected) == 0
 
 
-class TestAnalysis(object):
+class TestParallelismAnalysis(object):
 
     @pytest.mark.parametrize('exprs,atomic,parallel', [
         (['Inc(u[gp[p, 0]+rx, gp[p, 1]+ry], cx*cy*src)'],
@@ -726,7 +726,10 @@ class TestAnalysis(object):
         (['Eq(v.forward, v[t+1, x, y-1]+v[t, x, y]+v[t, x, y-1])'],
          [], ['x']),
         (['Eq(v.forward, v+1)', 'Inc(u, v)'],
-         [], ['x', 'y'])
+         [], ['x', 'y']),
+        # Test for issue #1902
+        (['Eq(u[0, y], v)', 'Eq(w, u[0, y])'],
+         [], ['y'])
     ])
     def test_iteration_parallelism_2d(self, exprs, atomic, parallel):
         """Tests detection of PARALLEL_* properties."""
@@ -742,6 +745,7 @@ def test_iteration_parallelism_2d(self, exprs, atomic, parallel):
 
         u = Function(name='u', grid=grid)  # noqa
         v = TimeFunction(name='v', grid=grid, save=None)  # noqa
+        w = TimeFunction(name='w', grid=grid, save=None)  # noqa
 
         cx = Function(name='coeff_x', dimensions=(p, rx), shape=(1, 2))  # noqa
         cy = Function(name='coeff_y', dimensions=(p, ry), shape=(1, 2))  # noqa