compiler: Shed reduction-only halos when lowering injections

After lifting lower_sparse_ops out of mpiize, an injection nest is
turned into a Call before optimize_halospots runs, so
_drop_reduction_halospots can no longer detect that the wrapping
HaloSpot's entry for the injected field is reduction-only. The stale
entry was left in place, and on save=True (no modulo buffering) the
hoist pass propagated the loop iteration variable out of the time
loop, producing an undeclared 'time' reference. Drop those entries at
lowering time so the resulting HaloSpot only carries entries with a
genuine read at the IET level.

Author: mloubout

devito/ir/cgen/printer.py

@@ -212,7 +212,12 @@ def _print_Pow(self, expr):
         suffix = self.func_literal(expr)
         base = self._print(expr.base)
         if equal_valued(expr.exp, -1):
-            return self._print_Float(Float(1.0)) + '/' + \
+            # Pick the literal precision from this Pow's dtype rather than
+            # the printer default, so e.g. ``Pow(DOUBLE(h), -1)`` emits
+            # ``1.0/(double)h`` not ``1.0F/(double)h``. This branch only
+            # fires when the surrounding Mul printer chose not to group
+            # the Pow into a denominator (e.g. lone reciprocal).
+            return f'1.0{self.prec_literal(expr)}/' + \
                 self.parenthesize(expr.base, PREC)
         elif equal_valued(expr.exp, 0.5):
             return f'{self.ns(expr)}sqrt{suffix}({base})'

devito/operations/interpolators.py

@@ -14,6 +14,7 @@
 from devito.finite_differences.elementary import floor
 from devito.logger import warning
 from devito.symbolics import INT, retrieve_function_carriers, retrieve_functions
+from devito.symbolics.extended_dtypes import DOUBLE
 from devito.tools import as_tuple, filter_ordered, memoized_meth
 from devito.types import (
     Eq, Inc, IncrInterpolation, Injection, Interpolation, SubFunction, Symbol
@@ -257,9 +258,44 @@ def _augment_implicit_dims(self, implicit_dims, extras=None):
     def _coeff_temps(self, implicit_dims, shifts=None):
         return []
 
+    @memoized_meth
+    def _raw_pos_symbols(self, shifts=None):
+        """
+        Per-Dimension Symbol holding the unrounded grid-relative position
+        ``(coord - origin - shift)/h``. Both the integer position
+        (``floor(...)``) and the linear-interp fractional part
+        (``... - floor(...)``) reuse this Symbol so the divide-and-shift
+        expression is emitted only once per sparse point.
+        """
+        dtype = self.sfunction.coordinates.dtype
+        symbols = []
+        for d, s in zip(self.grid.dimensions,
+                        shifts or (0,) * len(self.grid.dimensions),
+                        strict=True):
+            suffix = '_s1' if s != 0 else ''
+            symbols.append(Symbol(name=f'rpos{d}{suffix}', dtype=dtype))
+        return DimensionTuple(*symbols, getters=self.grid.dimensions)
+
     def _positions(self, implicit_dims, shifts=None):
-        return [Eq(v, INT(floor(k)), implicit_dims=implicit_dims)
-                for k, v in self.sfunction._position_map(shifts=shifts).items()]
+        # The ``(coord - origin)/h`` subtract is the only step that can lose
+        # precision to catastrophic cancellation when ``coord`` and ``origin``
+        # are large and close to each other (e.g. an origin-shifted survey).
+        # Promote ``origin`` and ``h`` to float64 so the subtract and divide
+        # happen in double precision in C (one cast operand promotes the
+        # whole expression); the result narrows to the field dtype on store
+        # to ``rpos*`` so downstream ``floor`` / fractional math stays in
+        # the field dtype.
+        rposs = self._raw_pos_symbols(shifts=shifts)
+        subs = {o: DOUBLE(o) for o in self.grid.origin_symbols}
+        subs.update({d.spacing: DOUBLE(d.spacing) for d in self._gdims})
+        return [Eq(rposs[d], k.xreplace(subs), implicit_dims=implicit_dims)
+                for d, k in zip(self._gdims,
+                                self.sfunction._position_map(shifts=shifts),
+                                strict=True)] + \
+               [Eq(v, INT(floor(rposs[d])), implicit_dims=implicit_dims)
+                for d, v in zip(self._gdims,
+                                self.sfunction._position_map(shifts=shifts).values(),
+                                strict=True)]
 
     def sparse_temps(self, rhs, implicit_dims, field=None):
         """
@@ -458,13 +494,14 @@ def _point_symbols(self, shifts=None):
         return DimensionTuple(*symbols, getters=self.grid.dimensions)
 
     def _coeff_temps(self, implicit_dims, shifts=None):
-        # Positions
-        pmap = self.sfunction._position_map(shifts=shifts)
+        # The fractional part of the unrounded position; reuse the
+        # ``rpos*`` Symbols emitted by ``_positions`` rather than the full
+        # ``(c - o)/h`` expression so the divide is computed only once.
+        rposs = self._raw_pos_symbols(shifts=shifts)
         psyms = self._point_symbols(shifts)
-        poseq = [Eq(psyms[d], pos - floor(pos),
-                    implicit_dims=implicit_dims)
-                 for (d, pos) in zip(self._gdims, pmap.keys(), strict=True)]
-        return poseq
+        return [Eq(psyms[d], rposs[d] - floor(rposs[d]),
+                   implicit_dims=implicit_dims)
+                for d in self._gdims]
 
 
 class PrecomputedInterpolator(WeightedInterpolator):

devito/passes/iet/parpragma.py

@@ -130,6 +130,16 @@ def _make_simd(self, iet):
                 if any(i.is_Indexed for i in reductions):
                     continue
 
+            # A scalar reduction in the loop body (e.g. ``acc += ...`` with
+            # ``acc`` a Symbol, as emitted by the interpolation accumulator
+            # pattern) is a cross-iteration dependence the SIMD pragma alone
+            # can't express -- it would need a ``reduction(+:acc)`` clause,
+            # which we don't emit here. Without it, ``_Complex`` accumulators
+            # are miscompiled by some gcc releases.
+            exprs = FindNodes(Expression).visit(candidate)
+            if any(e.is_reduction and not e.output.is_Indexed for e in exprs):
+                continue
+
             # Add SIMD pragma
             simd = self._make_simd_pragma(candidate)
             pragmas = candidate.pragmas + simd

devito/passes/iet/sparse.py

@@ -69,13 +69,16 @@ def lower_sparse_ops(iet, sregistry=None, **kwargs):
             continue
         groups.setdefault(nest, []).append(expr)
 
-    # If a sparse-op nest sits inside a HaloSpot whose halo scheme is
-    # void (the reduction-only halo got dropped by
-    # ``_drop_reduction_halospots``), replace the HaloSpot rather than
-    # just the nest so we don't leave behind an empty HaloSpot — the
-    # MPI overlap machinery would otherwise try to wrap our Call with
-    # its own dynamic-args plumbing.
-    parents = {nest: _enclosing_void_halospot(iet, nest) for nest in groups}
+    # ``lower_sparse_ops`` runs before ``optimize_halospots``, so the
+    # halo-exchange optimiser hasn't yet had a chance to drop the
+    # reduction-only halo entries that the IR scheduler put around an
+    # injection nest (e.g. an entry for ``u`` at ``loc_indices={time:
+    # time+1}`` wrapping ``u[time+1] += ...``). Once the nest becomes a
+    # Call those expressions are no longer visible to
+    # ``_drop_reduction_halospots``, so we shed those entries here -- and
+    # if that empties the HaloSpot, replace it whole so the MPI overlap
+    # machinery doesn't wrap our Call with stale dynamic-args plumbing.
+    parents = {nest: _enclosing_halospot(iet, nest) for nest in groups}
 
     mapper = {}
     efuncs = []
@@ -87,7 +90,26 @@ def lower_sparse_ops(iet, sregistry=None, **kwargs):
         efunc = make_callable(sregistry.make_name(prefix=prefix), new_nest)
         efuncs.append(efunc)
 
-        mapper[parents[nest] or nest] = Call(efunc.name, list(efunc.parameters))
+        call = Call(efunc.name, list(efunc.parameters))
+        parent = parents[nest]
+        if parent is None:
+            mapper[nest] = call
+            continue
+
+        # Drop fields that the (now-opaque) Call only writes/increments,
+        # since the wrapping HaloSpot's purpose was to ensure read-side
+        # coherency for them and the read no longer exists at the IET
+        # level. Interpolation reads its target field, so its entries
+        # stay.
+        reduced = {e.expr.lhs.function for e in exprs
+                   if isinstance(e.expr, InjectionMixin)}
+        hs = parent.halo_scheme.drop(reduced) if reduced else parent.halo_scheme
+        if hs.is_void:
+            mapper[parent] = call
+        elif hs is parent.halo_scheme:
+            mapper[nest] = call
+        else:
+            mapper[parent] = parent._rebuild(halo_scheme=hs, body=call)
 
     if not mapper:
         return iet, {}
@@ -107,14 +129,12 @@ def _find_outer_iteration(iet, expr):
     return None
 
 
-def _enclosing_void_halospot(iet, nest):
+def _enclosing_halospot(iet, nest):
     """
-    Return the HaloSpot directly wrapping ``nest`` if it carries an
-    empty (void) HaloScheme, otherwise None. Such HaloSpots are leftover
-    after ``_drop_reduction_halospots`` cleared all entries.
+    Return the HaloSpot directly wrapping ``nest``, if any.
     """
     for hs in FindNodes(HaloSpot).visit(iet):
-        if hs.is_void and nest in FindNodes(Iteration).visit(hs):
+        if nest in FindNodes(Iteration).visit(hs):
             return hs
     return None
 

devito/symbolics/extended_sympy.py

@@ -380,6 +380,8 @@ class UnaryOp(Expr, Pickable, BasicWrapperMixin):
 
     _op = ''
 
+    is_commutative = True
+
     __rargs__ = ('base',)
 
     def __new__(cls, base, **kwargs):

devito/symbolics/inspection.py

@@ -321,6 +321,14 @@ def sympy_dtype(expr, base=None, default=None, smin=None):
         with suppress(AttributeError):
             dtypes.add(i.dtype)
 
+    # A ``Cast`` overrides the dtype of the symbol(s) it wraps -- e.g.
+    # ``DOUBLE(o_x)`` is observably double-typed even though ``o_x``
+    # itself is float-typed. Without this, the C printer would pick
+    # float-precision literals (``1.0F``) when emitting a ``Pow(DOUBLE(h),
+    # -1)`` inside a wider expression.
+    for c in expr.atoms(Cast):
+        dtypes.add(c.dtype)
+
     if not dtypes or not np.issubdtype(base, np.complexfloating):
         dtypes.update({base} - {None})
 

devito/types/dense.py

@@ -1633,8 +1633,11 @@ def _arg_values(self, **kwargs):
             return self._arg_defaults(alias=self)
 
     def _arg_apply(self, *args, **kwargs):
+        # Parent-owned SubFunction data is computed once and read by the
+        # Operator; the parent gathers its own data via its own parameter
+        # entry, so there's nothing for the SubFunction to do here.
         if self._parent is not None:
-            return self._parent._arg_apply(*args, **kwargs)
+            return
         return super()._arg_apply(*args, **kwargs)
 
     @property