tosquash: Docs and autotune, notebooks

devito/core/cpu.py

@@ -317,7 +317,7 @@ def _make_iet_passes_mapper(cls, **kwargs):
         return {
             'denormals': avoid_denormals,
             'optcomms': optimize_halospots,
-            'blocking': partial(finalize_loop_bounds, sregistry=sregistry),
+            'blocking': partial(finalize_loop_bounds),
             'parallel': parizer.make_parallel,
             'openmp': parizer.make_parallel,
             'mpi': partial(mpiize, mode=options['mpi']),

devito/passes/iet/misc.py

@@ -88,10 +88,10 @@ def finalize_loop_bounds(iet, **kwargs):
 
         outer, inner = split(iterations, lambda i: not i.dim.parent.is_Incr)
 
-        # Get symbolic_max out of each outer dimension
+        # Get root symbolic_max out of each outer dimension
         roots_max = {i.dim.root: i.symbolic_max for i in outer}
 
-        # A dictionary to map maximum of parent dimensions
+        # A dictionary to map maximum of processed parent dimensions
         # useful for hierarchical blocking
         proc_parents_max = {}
 
@@ -100,50 +100,56 @@ def finalize_loop_bounds(iet, **kwargs):
             assert i.direction is Forward
             if (i.dim.parent in proc_parents_max and
                i.symbolic_size == i.dim.parent.step):
-                # Special case: A parent dimension may have already been processed
-                # as a member of 'inner' iterations. In this case we can use parent's
-                # Max. Usually encountered in hierarchical blocking (BLOCKLEVELS > 1)
-                it_max = proc_parents_max[i.dim.parent]
+                # In hierarchical blocking (BLOCKLEVELS > 1) a parent dimension may
+                # have already been processed as an 'inner' iteration. Since in
+                # hierarchical blocking, block sizes in lower levels always perfectly
+                # divide block sizes of upper levels we can use parent's iteration
+                # maximum.
+                iter_max = proc_parents_max[i.dim.parent]
             else:
+                # Candidate 1: symbolic_max of current iteration
                 # Most of the cases pass though this code:
-                # Candidates for iteration's max are:
+                # Candidates for iteration's maximum are:
                 # Candidate 1: symbolic_max of current iteration
                 # e.g. `i.symbolic_max = x0_blk0 + x0_blk0_size`
                 symbolic_max = i.symbolic_max
 
-                # Candidate 2: The domain max. Usualy it is the max of parent/root
-                # dimension.
-                # e.g. `x_M`
-                # This may not always be true as the symbolic_size of an Iteration may
-                # exceed the size of a parent's block size (e.g. after CIRE passes)
+                # Candidate 2: maximum of the root iteration
+                # Candidate 2a: Usualy it is the max of parent/root
+                # dimension, e.g. `x_M`
+                root_max = roots_max[i.dim.root]
+
+                # Candidate 2b:
+                # the symbolic_size of an Iteration may exceed the size of a parent's
+                # block size (e.g. after CIRE passes)
                 # e.g.
                 # `i.dim.parent.step = x0_blk1_size`
                 # `i.symbolic_size = x0_blk1_size + 4`
-
                 # For this case, proper margin should be allowed
-                # in order not to drop iterations. So Candidate 2 is the maximum of the
-                # root's max and the current iteration's required max
-                # and instead of `x_M` we may need `x_M + 1` or `x_M + 2`
-                root_max = roots_max[i.dim.root]
-
+                # in order not to omit loop iterations.
                 lb_margin = i.symbolic_min - i.dim.symbolic_min
                 size_margin = i.symbolic_size - i.dim.parent.step
-                upper_margin = i.dim.parent.symbolic_max + size_margin + lb_margin
+                ub_margin = i.dim.parent.symbolic_max + size_margin + lb_margin
+
+                # So Candidate 2 is the maximum of the
+                # root's max (2a) and the current iteration's required max (2b)
+                # and instead of `x_M` we may need `x_M + 1` or `x_M + 2`
 
-                # Domain max candidate
+                # So candidate 2 is
                 # e.g. `domain_max = Max(x_M + 1, x_M)``
-                domain_max = Max(upper_margin, root_max)
+                domain_max = Max(ub_margin, root_max)
 
-                # Finally the iteration's maximum is the minimum of the candidates
-                # e.g. `it_max = Min(x0_blk0 + x0_blk0_size, domain_max)``
-                it_max = Min(symbolic_max, domain_max)
+                # Finally the iteration's maximum is the minimum of the
+                # candidates 1 and 2
+                # e.g. `iter_max = Min(x0_blk0 + x0_blk0_size, domain_max)``
+                iter_max = Min(symbolic_max, domain_max)
 
             # Store the selected maximum of this iteration's dimension for
             # possible reference in case of children iterations
             # Usually encountered in subdims and hierarchical blocking
-            proc_parents_max[i.dim] = it_max
+            proc_parents_max[i.dim] = iter_max
 
-            mapper[i] = i._rebuild(limits=(i.symbolic_min, it_max, i.step))
+            mapper[i] = i._rebuild(limits=(i.symbolic_min, iter_max, i.step))
 
     iet = Transformer(mapper, nested=True).visit(iet)
 

examples/cfd/01_convection.ipynb

@@ -460,7 +460,7 @@
       "    {\n",
       "      for (int y = y_m; y <= y_M; y += 1)\n",
       "      {\n",
-      "        u[t1][x + 1][y + 1] = dt*(1.0F*u[t0][x + 1][y]/h_y - 1.0F*u[t0][x + 1][y + 1]/h_y + 1.0F*u[t0][x][y + 1]/h_x - 1.0F*u[t0][x + 1][y + 1]/h_x + u[t0][x + 1][y + 1]/dt);\n",
+      "        u[t1][x + 1][y + 1] = dt*(-(-u[t0][x][y + 1]/h_x + u[t0][x + 1][y + 1]/h_x) - (-u[t0][x + 1][y]/h_y + u[t0][x + 1][y + 1]/h_y) + u[t0][x + 1][y + 1]/dt);\n",
       "      }\n",
       "    }\n",
       "    STOP_TIMER(section0,timers)\n",

examples/cfd/01_convection_revisited.ipynb

@@ -435,7 +435,7 @@
       "    {\n",
       "      for (int i0y = i0y_ltkn + y_m; i0y <= -i0y_rtkn + y_M; i0y += 1)\n",
       "      {\n",
-      "        u[t1][i0x + 1][i0y + 1] = dt*(1.0F*u[t0][i0x + 1][i0y]/h_y - 1.0F*u[t0][i0x + 1][i0y + 1]/h_y + 1.0F*u[t0][i0x][i0y + 1]/h_x - 1.0F*u[t0][i0x + 1][i0y + 1]/h_x + u[t0][i0x + 1][i0y + 1]/dt);\n",
+      "        u[t1][i0x + 1][i0y + 1] = dt*(-(-u[t0][i0x][i0y + 1]/h_x + u[t0][i0x + 1][i0y + 1]/h_x) - (-u[t0][i0x + 1][i0y]/h_y + u[t0][i0x + 1][i0y + 1]/h_y) + u[t0][i0x + 1][i0y + 1]/dt);\n",
       "      }\n",
       "    }\n",
       "    STOP_TIMER(section0,timers)\n",