compiler: Patch race conditions due to storage-related dependencies

devito/arch/compiler.py

@@ -481,7 +481,13 @@ def __init__(self, *args, **kwargs):
         self.cflags.remove('-std=c99')
         self.cflags.remove('-O3')
         self.cflags.remove('-Wall')
-        self.cflags += ['-std=c++11', '-acc:gpu', '-gpu=pinned', '-mp']
+
+        self.cflags += ['-std=c++11', '-mp']
+
+        platform = kwargs.pop('platform', configuration['platform'])
+        if platform is NVIDIAX:
+            self.cflags += ['-acc:gpu', '-gpu=pinned']
+
         if not configuration['safe-math']:
             self.cflags.append('-fast')
         # Default PGI compile for a target is GPU and single threaded host.

devito/core/cpu.py

@@ -143,14 +143,15 @@ class Cpu64NoopOperator(Cpu64OperatorMixin, CoreOperator):
     def _specialize_iet(cls, graph, **kwargs):
         options = kwargs['options']
         platform = kwargs['platform']
+        compiler = kwargs['compiler']
         sregistry = kwargs['sregistry']
 
         # Distributed-memory parallelism
         mpiize(graph, sregistry=sregistry, options=options)
 
         # Shared-memory parallelism
         if options['openmp']:
-            parizer = cls._Target.Parizer(sregistry, options, platform)
+            parizer = cls._Target.Parizer(sregistry, options, platform, compiler)
             parizer.make_parallel(graph)
             parizer.initialize(graph)
 
@@ -213,6 +214,7 @@ def _specialize_clusters(cls, clusters, **kwargs):
     def _specialize_iet(cls, graph, **kwargs):
         options = kwargs['options']
         platform = kwargs['platform']
+        compiler = kwargs['compiler']
         sregistry = kwargs['sregistry']
 
         # Flush denormal numbers
@@ -225,7 +227,7 @@ def _specialize_iet(cls, graph, **kwargs):
         relax_incr_dimensions(graph)
 
         # Parallelism
-        parizer = cls._Target.Parizer(sregistry, options, platform)
+        parizer = cls._Target.Parizer(sregistry, options, platform, compiler)
         parizer.make_simd(graph)
         parizer.make_parallel(graph)
         parizer.initialize(graph)
@@ -304,9 +306,10 @@ def callback(f):
     def _make_iet_passes_mapper(cls, **kwargs):
         options = kwargs['options']
         platform = kwargs['platform']
+        compiler = kwargs['compiler']
         sregistry = kwargs['sregistry']
 
-        parizer = cls._Target.Parizer(sregistry, options, platform)
+        parizer = cls._Target.Parizer(sregistry, options, platform, compiler)
 
         return {
             'denormals': avoid_denormals,

devito/core/gpu.py

@@ -130,13 +130,14 @@ class DeviceNoopOperator(DeviceOperatorMixin, CoreOperator):
     def _specialize_iet(cls, graph, **kwargs):
         options = kwargs['options']
         platform = kwargs['platform']
+        compiler = kwargs['compiler']
         sregistry = kwargs['sregistry']
 
         # Distributed-memory parallelism
         mpiize(graph, sregistry=sregistry, options=options)
 
         # GPU parallelism
-        parizer = cls._Target.Parizer(sregistry, options, platform)
+        parizer = cls._Target.Parizer(sregistry, options, platform, compiler)
         parizer.make_parallel(graph)
         parizer.initialize(graph)
 
@@ -202,6 +203,7 @@ def _specialize_clusters(cls, clusters, **kwargs):
     def _specialize_iet(cls, graph, **kwargs):
         options = kwargs['options']
         platform = kwargs['platform']
+        compiler = kwargs['compiler']
         sregistry = kwargs['sregistry']
 
         # Distributed-memory parallelism
@@ -211,7 +213,7 @@ def _specialize_iet(cls, graph, **kwargs):
         relax_incr_dimensions(graph)
 
         # GPU parallelism
-        parizer = cls._Target.Parizer(sregistry, options, platform)
+        parizer = cls._Target.Parizer(sregistry, options, platform, compiler)
         parizer.make_parallel(graph)
         parizer.initialize(graph)
 
@@ -282,9 +284,10 @@ def callback(f):
     def _make_iet_passes_mapper(cls, **kwargs):
         options = kwargs['options']
         platform = kwargs['platform']
+        compiler = kwargs['compiler']
         sregistry = kwargs['sregistry']
 
-        parizer = cls._Target.Parizer(sregistry, options, platform)
+        parizer = cls._Target.Parizer(sregistry, options, platform, compiler)
         orchestrator = cls._Target.Orchestrator(sregistry)
 
         return {

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
 

devito/passes/iet/langbase.py

@@ -185,7 +185,7 @@ class LangTransformer(ABC):
     The constructs of the target language. To be specialized by a subclass.
     """
 
-    def __init__(self, key, sregistry, platform):
+    def __init__(self, key, sregistry, platform, compiler):
         """
         Parameters
         ----------
@@ -195,13 +195,16 @@ def __init__(self, key, sregistry, platform):
             The symbol registry, to access the symbols appearing in an IET.
         platform : Platform
             The underlying platform.
+        compiler : Compiler
+            The underlying JIT compiler.
         """
         if key is not None:
             self.key = key
         else:
             self.key = lambda i: False
         self.sregistry = sregistry
         self.platform = platform
+        self.compiler = compiler
 
     @iet_pass
     def make_parallel(self, iet):

devito/passes/iet/languages/openmp.py

@@ -1,7 +1,10 @@
+from distutils import version
+
 import cgen as c
 from sympy import Not
 
 from devito.arch import AMDGPUX, NVIDIAX, INTELGPUX
+from devito.arch.compiler import GNUCompiler
 from devito.ir import (Call, Conditional, List, Prodder, ParallelIteration,
                        ParallelBlock, PointerCast, While, FindSymbols)
 from devito.passes.iet.definitions import DataManager, DeviceAwareDataManager
@@ -186,8 +189,18 @@ class SimdOmpizer(PragmaSimdTransformer):
 
 
 class Ompizer(PragmaShmTransformer):
+
     lang = OmpBB
 
+    @classmethod
+    def _support_array_reduction(cls, compiler):
+        # Not all backend compilers support array reduction!
+        # Here are the known unsupported ones:
+        if isinstance(compiler, GNUCompiler) and \
+           compiler.version < version.StrictVersion("6.0"):
+            return False
+        return True
+
 
 class DeviceOmpizer(PragmaDeviceAwareTransformer):
     lang = DeviceOmpBB

devito/passes/iet/parpragma.py

@@ -95,7 +95,7 @@ class PragmaShmTransformer(PragmaSimdTransformer):
     and shared-memory-parallel IETs.
     """
 
-    def __init__(self, sregistry, options, platform):
+    def __init__(self, sregistry, options, platform, compiler):
         """
         Parameters
         ----------
@@ -116,9 +116,11 @@ def __init__(self, sregistry, options, platform):
                is greater than this threshold.
         platform : Platform
             The underlying platform.
+        compiler : Compiler
+            The underlying JIT compiler.
         """
         key = lambda i: i.is_ParallelRelaxed and not i.is_Vectorized
-        super().__init__(key, sregistry, platform)
+        super().__init__(key, sregistry, platform, compiler)
 
         self.collapse_ncores = options['par-collapse-ncores']
         self.collapse_work = options['par-collapse-work']
@@ -209,14 +211,22 @@ def _select_candidates(self, candidates):
 
         return root, list(collapsable)
 
+    @classmethod
+    def _support_array_reduction(cls, compiler):
+        return True
+
     def _make_reductions(self, partree):
         if not any(i.is_ParallelAtomic for i in partree.collapsed):
             return partree
 
         exprs = [i for i in FindNodes(Expression).visit(partree) if i.is_Increment]
         reduction = [i.output for i in exprs]
-        if all(i.is_Affine for i in partree.collapsed) or \
-           all(not i.is_Indexed for i in reduction):
+
+        test0 = all(not i.is_Indexed for i in reduction)
+        test1 = (self._support_array_reduction(self.compiler) and
+                 all(i.is_Affine for i in partree.collapsed))
+
+        if test0 or test1:
             # Implement reduction
             mapper = {partree.root: partree.root._rebuild(reduction=reduction)}
         else:
@@ -399,8 +409,8 @@ class PragmaDeviceAwareTransformer(DeviceAwareMixin, PragmaShmTransformer):
     shared-memory-parallel, and device-parallel IETs.
     """
 
-    def __init__(self, sregistry, options, platform):
-        super().__init__(sregistry, options, platform)
+    def __init__(self, sregistry, options, platform, compiler):
+        super().__init__(sregistry, options, platform, compiler)
 
         self.gpu_fit = options['gpu-fit']
         self.par_tile = options['par-tile']

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))