diff --git a/devito/finite_differences/finite_difference.py b/devito/finite_differences/finite_difference.py
index cd428641735..44859b09038 100644
--- a/devito/finite_differences/finite_difference.py
+++ b/devito/finite_differences/finite_difference.py
@@ -239,7 +239,10 @@ def make_derivative(expr, dim, fd_order, deriv_order, side, matvec, x0, symbolic
             # For homogeneity, always generate e.g. `x + i0` rather than `x - i0`
             # for transpose and `x + i0` for direct
             indices = indices.transpose()
-            weights = weights._subs(indices.free_dim, -indices.free_dim)
+
+            # Do the same for the Weights, though this is more than just a
+            # transposition, we also must switch to the transposed StencilDimension
+            weights = weights._subs(weights.dimension, -indices.free_dim)
 
         # Inject the StencilDimension
         # E.g. `x + i*h_x` into `f(x)` s.t. `f(x + i*h_x)`
diff --git a/devito/finite_differences/tools.py b/devito/finite_differences/tools.py
index 5ec3629403d..e1f209c4d6f 100644
--- a/devito/finite_differences/tools.py
+++ b/devito/finite_differences/tools.py
@@ -203,10 +203,7 @@ def transpose(self):
         """
         indices = tuple(reversed(self))
 
-        free_dim = StencilDimension(self.free_dim.name,
-                                    -self.free_dim._max,
-                                    -self.free_dim._min,
-                                    backward=True)
+        free_dim = self.free_dim.transpose()
 
         try:
             expr = self.expr._subs(self.free_dim, -free_dim)
diff --git a/devito/ir/clusters/algorithms.py b/devito/ir/clusters/algorithms.py
index 6a9cd2e3c70..296196d0096 100644
--- a/devito/ir/clusters/algorithms.py
+++ b/devito/ir/clusters/algorithms.py
@@ -64,7 +64,7 @@ def impose_total_ordering(clusters):
         try:
             relations = {tuple(sorted(c.ispace.itdims, key=key))}
         except ValueError:
-            # See issue #X
+            # See issue #2204
             relations = c.ispace.relations
         ispace = c.ispace.reorder(relations=relations, mode='total')
 
diff --git a/devito/ir/clusters/cluster.py b/devito/ir/clusters/cluster.py
index 6d5808caa0c..f9f0f3be91a 100644
--- a/devito/ir/clusters/cluster.py
+++ b/devito/ir/clusters/cluster.py
@@ -484,12 +484,12 @@ def tailor_properties(properties, ispace):
     Create a new Properties object off `properties` that retains all and only
     the iteration dimensions in `ispace`.
     """
-    for d in ispace.itdims:
-        properties = properties.add(d)
-
     for i in properties:
         for d in as_tuple(i):
             if d not in ispace.itdims:
                 properties = properties.drop(d)
 
+    for d in ispace.itdims:
+        properties = properties.add(d)
+
     return properties
diff --git a/devito/passes/clusters/aliases.py b/devito/passes/clusters/aliases.py
index 5db4a787c26..bfa10147488 100644
--- a/devito/passes/clusters/aliases.py
+++ b/devito/passes/clusters/aliases.py
@@ -996,7 +996,8 @@ def pick_best(variants):
 
         # Magic sauce
         # The coefficients were obtained empirically studying the behaviour
-        # of different variants in several kernels
+        # of different variants in several kernels and platforms
+        # Intuitively, it's like trading 70 operations for 1 temporary
         ws_curve = lambda x: (-1 / 70)*x
 
         if delta_ws <= ws_curve(delta_flops):
diff --git a/devito/types/dimension.py b/devito/types/dimension.py
index bc60998fe82..15f53325fc8 100644
--- a/devito/types/dimension.py
+++ b/devito/types/dimension.py
@@ -1376,9 +1376,9 @@ class StencilDimension(BasicDimension):
     is_Stencil = True
 
     __rargs__ = BasicDimension.__rargs__ + ('_min', '_max')
-    __rkwargs__ = BasicDimension.__rkwargs__ + ('backward',)
+    __rkwargs__ = BasicDimension.__rkwargs__ + ('step',)
 
-    def __init_finalize__(self, name, _min, _max, spacing=None, backward=False,
+    def __init_finalize__(self, name, _min, _max, spacing=None, step=1,
                           **kwargs):
         self._spacing = sympy.sympify(spacing) or sympy.S.One
 
@@ -1388,20 +1388,25 @@ def __init_finalize__(self, name, _min, _max, spacing=None, backward=False,
             raise ValueError("Expected integer `max` (got %s)" % _max)
         if not is_integer(self._spacing):
             raise ValueError("Expected integer `spacing` (got %s)" % self._spacing)
+        if not is_integer(step):
+            raise ValueError("Expected integer `step` (got %s)" % step)
 
         self._min = _min
         self._max = _max
+        self._step = step
 
         self._size = _max - _min + 1
 
-        self._backward = backward
-
         if self._size < 1:
             raise ValueError("Expected size greater than 0 (got %s)" % self._size)
 
-    @cached_property
+    @property
+    def step(self):
+        return self._step
+
+    @property
     def backward(self):
-        return self._backward
+        return self.step < 0
 
     @cached_property
     def symbolic_size(self):
@@ -1419,6 +1424,9 @@ def symbolic_max(self):
     def range(self):
         return range(self._min, self._max + 1)
 
+    def transpose(self):
+        return StencilDimension(self.name, -self._max, -self._min, step=-1)
+
     @property
     def _arg_names(self):
         return ()