compiler: prevent undefined Temp through global init

Author: mloubout

devito/passes/clusters/aliases.py

@@ -12,6 +12,7 @@
     Interval, IntervalGroup, IterationSpace, LabeledVector, Queue, Vector, extrema,
     maximum, minimum, normalize_properties, relax_properties, unbounded, vmax, vmin
 )
+from devito.ir.support import null_ispace
 from devito.passes.clusters.cse import _cse
 from devito.symbolics import (
     Uxmapper, estimate_cost, retrieve_functions, reuse_if_untouched, search, sympy_dtype,
@@ -131,10 +132,11 @@ def _aliases_from_clusters(self, cgroup, exclude, meta):
             # AliasList -> Schedule
             schedule = lower_aliases(aliases, meta, self.opt_maxpar)
 
-            variants.append(Variant(schedule, exprs))
+            if schedule:
+                variants.append(Variant(schedule, exprs))
 
         if not variants:
-            return []
+            return [], []
 
         # [Schedule]_m -> Schedule (s.t. best memory/flops trade-off)
         schedule, exprs = self._select(variants)
@@ -144,9 +146,9 @@ def _aliases_from_clusters(self, cgroup, exclude, meta):
             schedule = optimize_schedule_rotations(schedule, self.sregistry)
 
         # Schedule -> [Clusters]_k
-        processed, subs = lower_schedule(schedule, meta, self.sregistry,
-                                         self.opt_ftemps, self.opt_min_dtype,
-                                         self.opt_minmem)
+        processed, subs, inits = lower_schedule(schedule, meta, self.sregistry,
+                                                self.opt_ftemps, self.opt_min_dtype,
+                                                self.opt_minmem)
 
         # [Clusters]_k -> [Clusters]_k (optimization)
         if self.opt_multisubdomain:
@@ -166,7 +168,7 @@ def _aliases_from_clusters(self, cgroup, exclude, meta):
 
         assert len(exprs) == 0
 
-        return processed
+        return processed, inits
 
     def process(self, clusters):
         raise NotImplementedError
@@ -302,12 +304,13 @@ def callback(self, clusters, prefix, xtracted=None):
 
         key = lambda c: self._lookup_key(c, d)
         processed = list(clusters)
+        inits = []
         for ak, group in as_mapper(clusters, key=key).items():
             g = [c for c in group if c.is_dense and c not in xtracted]
             if not g:
                 continue
 
-            made = self._aliases_from_clusters(ClusterGroup(g), exclude, ak)
+            made, cinits = self._aliases_from_clusters(ClusterGroup(g), exclude, ak)
 
             if made:
                 idx = processed.index(g[0])
@@ -317,7 +320,10 @@ def callback(self, clusters, prefix, xtracted=None):
 
                 xtracted.extend(made)
 
-        return processed
+            if cinits:
+                inits.extend(cinits)
+
+        return inits + processed
 
     def _lookup_key(self, c, d):
         ispace = c.ispace.reset()
@@ -390,7 +396,7 @@ def process(self, clusters):
             # TODO: to process third- and higher-order derivatives, we could
             # extend this by calling `_aliases_from_clusters` repeatedly until
             # `made` is empty. To be investigated
-            made = self._aliases_from_clusters(
+            made, _ = self._aliases_from_clusters(
                 ClusterGroup(c), exclude, self._lookup_key(c)
             )
 
@@ -860,6 +866,7 @@ def lower_schedule(schedule, meta, sregistry, opt_ftemps, opt_min_dtype,
     make = TempFunction if opt_ftemps else TempArray
 
     clusters = []
+    inits = []
     subs = {}
     for pivot, writeto, ispace, aliaseds, indicess in schedule:
         name = sregistry.make_name()
@@ -928,8 +935,11 @@ def lower_schedule(schedule, meta, sregistry, opt_ftemps, opt_min_dtype,
             assert writeto.size == 0
 
             dtype = sympy_dtype(pivot, base=meta.dtype, smin=opt_min_dtype)
-            obj = Temp(name=name, dtype=dtype)
+            const = not meta.guards
+            obj = Temp(name=name, dtype=dtype, is_const=const)
             expression = Eq(obj, uxreplace(pivot, subs))
+            if not const:
+                inits.append(Eq(obj, 0))
 
             callback = lambda idx: obj  # noqa: B023
 
@@ -959,7 +969,14 @@ def lower_schedule(schedule, meta, sregistry, opt_ftemps, opt_min_dtype,
         # Finally, build the alias Cluster
         clusters.append(Cluster(expression, ispace, meta.guards, properties))
 
-    return clusters, subs
+    if inits:
+        # To avoid undefined variables when an constant (Temp) alias is used
+        # within different guards/loop, we need to initialize it outside of the loops
+        # so that it's globally defined.
+        # See tests/test_dse.py::TestAliases::test_split_cond
+        inits = [Cluster(inits, null_ispace)]
+
+    return clusters, subs, inits
 
 
 def optimize_clusters_msds(clusters):

devito/passes/clusters/cse.py

@@ -40,6 +40,12 @@ def cse_dtype(exprdtype, cdtype):
     Return the dtype of a CSE temporary given the dtype of the expression to be
     captured and the cluster's dtype.
     """
+    if np.issubdtype(cdtype, np.floating) and np.issubdtype(exprdtype, np.integer):
+        # Integer expression and floating-point cluster: promote to the floating point
+        # np.promote_types upcast integers (e.g int32 -> Float64) so we
+        # need to ensure that the promoted type is not larger than the cluster's dtype
+        return cdtype
+
     if np.issubdtype(cdtype, np.complexfloating):
         return np.promote_types(exprdtype, cdtype(0).real.__class__).type
     else:
@@ -97,8 +103,9 @@ def cse(cluster, sregistry=None, options=None, **kwargs):
     if cluster.is_fence:
         return cluster
 
-    make_dtype = lambda e: cse_dtype(e.dtype, dtype)
-    make = lambda e: CTemp(name=sregistry.make_name(), dtype=make_dtype(e))
+    def make(e):
+        edtype = cse_dtype(e.dtype, dtype)
+        return CTemp(name=sregistry.make_name(), dtype=edtype)
 
     exprs = _cse(cluster, make, min_cost=min_cost, mode=mode)
 

devito/types/basic.py

@@ -1589,7 +1589,7 @@ def _rebuild(self, *args, **kwargs):
         comps = [f.func(*args, name=f.name.replace(self.name, newname), **kwargs)
                  for f in self.flat()]
         # Rebuild the matrix with the new components
-        return self._new(comps)
+        return self._new(*self.shape, comps)
 
     func = _rebuild
 

examples/mpi/overview.ipynb

@@ -460,7 +460,7 @@
        "  _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);\n",
        "  _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);\n",
        "\n",
-       "  float r0 = 1.0F/h_x;\n",
+       "  const float r0 = 1.0F/h_x;\n",
        "\n",
        "  for (int time = time_m, t0 = (time)%(2), t1 = (time + 1)%(2); time <= time_M; time += 1, t0 = (time)%(2), t1 = (time + 1)%(2))\n",
        "  {\n",

tests/test_dse.py

@@ -2685,6 +2685,52 @@ def test_space_and_time_invariant_together(self):
             'tx0_blk0y0_blk0xyzyz'
         )
 
+    def test_split_cond(self):
+        grid = Grid((11, 11))
+        time = grid.time_dim
+
+        u = TimeFunction(name='u', grid=grid, time_order=2, space_order=2)
+
+        ct = ConditionalDimension(name='ct', parent=time, factor=2)
+        ct2 = ConditionalDimension(name='ct2', parent=time, factor=4)
+
+        eq0 = Eq(u.forward, u + cos(time), implicit_dims=ct)
+        eq1 = Eq(u.forward, u.forward + 1, implicit_dims=ct2)
+        eq2 = Eq(u.forward, u.forward + cos(time), implicit_dims=ct)
+
+        op = Operator([eq0, eq1, eq2])
+        cond = FindNodes(Conditional).visit(op)
+        assert len(cond) == 3
+        assert str(cond[0].args['then_body'][0].exprs[0]) == 'r0 = cos(time);'
+        assert str(op.body.body[0].body[0].body[0]) == 'float r0 = 0;'
+
+    def test_split_cond_multi_alias(self):
+        grid = Grid((11, 11))
+        time = grid.time_dim
+
+        u = TimeFunction(name='u', grid=grid, time_order=2, space_order=2)
+
+        ct = ConditionalDimension(name='ct', parent=time, factor=2)
+        ct2 = ConditionalDimension(name='ct2', parent=time, factor=4)
+
+        eq0 = Eq(u.forward, u + cos(time) + sin(time), implicit_dims=ct)
+        eq1 = Eq(u.forward, u.forward + 1, implicit_dims=ct2)
+        eq2 = Eq(u.forward, u.forward + cos(time) - sin(time), implicit_dims=ct)
+
+        op = Operator([eq0, eq1, eq2])
+        cond = FindNodes(Conditional).visit(op)
+        assert len(cond) == 3
+        print(op)
+        assert str(cond[0].args['then_body'][0].exprs[0]) == 'float r3 = cos(time);'
+        assert str(cond[0].args['then_body'][0].exprs[1]) == 'float r4 = sin(time);'
+        assert str(cond[0].args['then_body'][0].exprs[2]) == 'r0 = r3 + r4;'
+        assert str(cond[0].args['then_body'][0].exprs[3]) == 'r1 = r3;'
+        assert str(cond[0].args['then_body'][0].exprs[4]) == 'r2 = r4;'
+
+        assert str(op.body.body[0].body[0].body[0]) == 'float r0 = 0;'
+        assert str(op.body.body[0].body[0].body[1]) == 'float r1 = 0;'
+        assert str(op.body.body[0].body[0].body[2]) == 'float r2 = 0;'
+
 
 class TestIsoAcoustic: