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",

examples/performance/00_overview.ipynb

@@ -508,7 +508,7 @@
       "}\n",
       "STOP(section0,timers)\n",
       "\n",
-      "float r1 = 1.0F/h_y;\n",
+      "const float r1 = 1.0F/h_y;\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",
@@ -572,7 +572,13 @@
       "+            u[t1][x + 4][y + 4][z + 4] = (f[x + 1][y + 1][z + 1]*f[x + 1][y + 1][z + 1])*((-6.66666667e-1F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 1][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 2][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 4][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 5][z + 4]) + (-8.33333333e-2F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 4][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 5][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 7][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 8][z + 4]) + (8.33333333e-2F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 1][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 3][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 4][z + 4]) + (6.66666667e-1F*r0)*(8.33333333e-2F*r0*u[t0][x + 4][y + 3][z + 4] - 6.66666667e-1F*r0*u[t0][x + 4][y + 4][z + 4] + 6.66666667e-1F*r0*u[t0][x + 4][y + 6][z + 4] - 8.33333333e-2F*r0*u[t0][x + 4][y + 7][z + 4]))*sinf(f[x + 1][y + 1][z + 1]);\n",
       "           }\n",
       "         }\n",
-      "       }\n",
+      "       }\n"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
       "\n"
      ]
     }
@@ -652,7 +658,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 12,
    "metadata": {},
    "outputs": [
     {
@@ -712,7 +718,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 14,
+   "execution_count": 13,
    "metadata": {},
    "outputs": [
     {
@@ -753,7 +759,14 @@
       "    }\n",
       "  }\n",
       "  STOP(section0,timers)\n",
-      "}\n"
+      "}"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "\n"
      ]
     }
    ],
@@ -772,7 +785,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": 14,
    "metadata": {},
    "outputs": [
     {
@@ -863,7 +876,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 15,
    "metadata": {},
    "outputs": [
     {
@@ -919,7 +932,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 17,
+   "execution_count": 16,
    "metadata": {},
    "outputs": [
     {
@@ -976,7 +989,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 18,
+   "execution_count": 17,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -994,7 +1007,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 19,
+   "execution_count": 18,
    "metadata": {},
    "outputs": [
     {
@@ -1044,7 +1057,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 20,
+   "execution_count": 19,
    "metadata": {},
    "outputs": [
     {
@@ -1112,7 +1125,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 21,
+   "execution_count": 20,
    "metadata": {},
    "outputs": [
     {
@@ -1192,7 +1205,7 @@
       "  }\n",
       "  STOP(section0,timers)\n",
       "\n",
-      "  float r1 = 1.0F/h_y;\n",
+      "  const float r1 = 1.0F/h_y;\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",
@@ -1279,7 +1292,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 22,
+   "execution_count": 21,
    "metadata": {},
    "outputs": [
     {
@@ -1304,7 +1317,7 @@
       "}\n",
       "STOP(section0,timers)\n",
       "\n",
-      "float r1 = 1.0F/h_y;\n",
+      "const float r1 = 1.0F/h_y;\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",
@@ -1369,7 +1382,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 23,
+   "execution_count": 22,
    "metadata": {},
    "outputs": [
     {
@@ -1404,7 +1417,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 24,
+   "execution_count": 23,
    "metadata": {},
    "outputs": [
     {
@@ -1483,7 +1496,7 @@
       "  }\n",
       "  STOP(section0,timers)\n",
       "\n",
-      "  float r1 = 1.0F/h_y;\n",
+      "  const float r1 = 1.0F/h_y;\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",
@@ -1546,7 +1559,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 25,
+   "execution_count": 24,
    "metadata": {},
    "outputs": [
     {
@@ -1624,8 +1637,8 @@
       "  }\n",
       "  STOP(section0,timers)\n",
       "\n",
-      "  float r1 = 1.0F/h_x;\n",
-      "  float r2 = 1.0F/h_y;\n",
+      "  const float r1 = 1.0F/h_x;\n",
+      "  const float r2 = 1.0F/h_y;\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",
@@ -1718,7 +1731,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.12"
+   "version": "3.13.11"
   }
  },
  "nbformat": 4,