diff --git a/devito/ir/clusters/algorithms.py b/devito/ir/clusters/algorithms.py
index 979d13ee8f..066cd029b3 100644
--- a/devito/ir/clusters/algorithms.py
+++ b/devito/ir/clusters/algorithms.py
@@ -156,9 +156,10 @@ def callback(self, clusters, prefix, backlog=None, known_break=None):
 
         # Schedule Clusters over different IterationSpaces if this increases
         # parallelism
-        for i in range(1, len(clusters)):
+        for i in reversed(range(len(clusters) - 1)):
             if self._break_for_parallelism(scope, candidates, i):
-                return self.callback(clusters[:i], prefix, clusters[i:] + backlog,
+                return self.callback(clusters[:i+1], prefix,
+                                     clusters[i+1:] + backlog,
                                      candidates | known_break)
 
         # Compute iteration direction
diff --git a/tests/test_fission.py b/tests/test_fission.py
index 9f7c74b6a6..bf3651ede1 100644
--- a/tests/test_fission.py
+++ b/tests/test_fission.py
@@ -2,7 +2,7 @@
 
 from conftest import assert_structure
 from devito import (Eq, Inc, Grid, Function, TimeFunction, SubDimension, SubDomain,
-                    Operator, solve)
+                    Operator, VectorTimeFunction, Buffer, curl, solve)
 
 
 def test_issue_1725():
@@ -125,3 +125,21 @@ def test_issue_1921():
     op1.apply(time_m=1, time_M=5, g=g1)
 
     assert np.all(g.data == g1.data)
+
+
+def test_fission_largest_cluster():
+    grid = Grid(shape=(10, 10, 10))
+    x, y, z = grid.dimensions
+
+    A = VectorTimeFunction(name='A', grid=grid, save=Buffer(1), space_order=2,
+                           time_order=1, staggered=(x, y, z))
+    B = VectorTimeFunction(name='B', grid=grid, save=Buffer(1), space_order=2,
+                           time_order=1, staggered=((y, z), (x, z), (x, y)))
+    f = Function(name='f', grid=grid, space_order=2)
+
+    eqs = [Eq(B.forward, solve(Eq(curl(A), -B.dt), B.forward)),
+           Eq(A.forward, solve(Eq(curl(B.forward), f*A.dt), A.forward))]
+
+    op = Operator(eqs, opt='fission')
+
+    assert_structure(op, ['t,x,y,z', 't,x,y,z'], 't,x,y,z,x,y,z')