compiler: prevent hosted per-thread arrays are dereferenced within partree at read

Author: mloubout

devito/finite_differences/differentiable.py

@@ -259,6 +259,9 @@ def __rfloordiv__(self, other):
         from .elementary import floor
         return floor(other / self)
 
+    def safe_inv(self, ref):
+        return SafeInv(self, ref or self)
+
     def __mod__(self, other):
         return Mod(self, other)
 

devito/passes/iet/parpragma.py

@@ -317,6 +317,9 @@ def _make_parregion(self, partree, parrays):
             i = n.write
             if not (i.is_Array or i.is_TempFunction):
                 continue
+            elif partree.dim in i.dimensions:
+                # Non-local Array (full iteration space): no need to vector-expand
+                continue
             elif i in parrays:
                 pi = parrays[i]
             else:

devito/symbolics/inspection.py

@@ -4,6 +4,7 @@
 from sympy import (Function, Indexed, Integer, Mul, Number,
                    Pow, S, Symbol, Tuple)
 from sympy.core.numbers import ImaginaryUnit
+from sympy.core.function import Application
 
 from devito.finite_differences import Derivative
 from devito.finite_differences.differentiable import IndexDerivative
@@ -116,7 +117,7 @@ def estimate_cost(exprs, estimate=False):
 estimate_values = {
     'elementary': 100,
     'pow': 50,
-    'SafeInv': 10,
+    'SafeInv': 50,
     'div': 5,
     'Abs': 5,
     'floor': 1,
@@ -211,6 +212,7 @@ def _(expr, estimate, seen):
 
 
 @_estimate_cost.register(Function)
+@_estimate_cost.register(Application)
 def _(expr, estimate, seen):
     if q_routine(expr):
         flops, _ = zip(*[_estimate_cost(a, estimate, seen) for a in expr.args])

devito/types/basic.py

@@ -989,7 +989,7 @@ def c0(self):
     def _eval_deriv(self):
         return self
 
-    @cached_property
+    @property
     def _grid_map(self):
         """
         Mapper of off-grid interpolation points indices for each dimension.
@@ -1049,14 +1049,13 @@ def _evaluate(self, **kwargs):
             return self
 
         io = self.interp_order
+        retval = self.subs({i.subs(subs): self.indices_ref[d]
+                            for d, i in mapper.items()})
         if self._avg_mode == 'harmonic':
-            retval = 1 / self
-        else:
-            retval = self
+            retval = retval.safe_inv(retval)
 
         # Apply interpolation from inner most dim
         for d, i in mapper.items():
-            retval = retval._subs(i.subs(subs), self.indices_ref[d])
             retval = retval.diff(d, deriv_order=0, fd_order=io, x0={d: i})
 
         # Evaluate. Since we used `self.function` it will be on the grid when
@@ -1066,8 +1065,7 @@ def _evaluate(self, **kwargs):
 
         # If harmonic averaging, invert at the end
         if self._avg_mode == 'harmonic':
-            from devito.finite_differences.differentiable import SafeInv
-            retval = SafeInv(retval, self.function.subs(subs))
+            retval = retval.safe_inv(self.function.subs(subs))
 
         return retval