--- isl-0.18/isl_options_private.h	2016-11-24 12:32:48.000000000 +0100
+++ isl-0.18/isl_options_private.h	2017-11-03 22:54:57.000000000 +0100
@@ -16,6 +16,7 @@
 
 	#define			ISL_CLOSURE_ISL		0
 	#define			ISL_CLOSURE_BOX		1
+        #define                 ISL_CLOSURE_ITER        2
 	unsigned		closure;
 
 	int			bound;
--- isl-0.18/isl_transitive_closure.c	2016-12-15 12:27:30.000000000 +0100
+++ isl-0.18/isl_transitive_closure.c	2017-11-05 21:42:09.000000000 +0100
@@ -2572,6 +2572,89 @@
 	return NULL;
 }
 
+isl_stat union_closure(__isl_take isl_basic_map *bmap, void *user)
+{
+        isl_map** union_map = (isl_map**)user;
+        int exact = 0;
+    
+        isl_map* rPlus = isl_map_from_basic_map(isl_basic_map_copy(bmap));
+        rPlus->ctx->opt->closure = ISL_CLOSURE_ISL;
+        rPlus = isl_map_transitive_closure(rPlus, &exact);
+        rPlus = isl_map_compute_divs(rPlus);
+        rPlus = isl_map_coalesce(isl_map_remove_redundancies(rPlus));
+
+        if (*union_map)
+                *union_map = isl_map_union(*union_map, exact ? rPlus : 
+                        isl_map_from_basic_map(bmap));
+        else
+                *union_map = exact ? rPlus : isl_map_from_basic_map(bmap);
+
+        *union_map = isl_map_compute_divs(*union_map);
+        *union_map = isl_map_coalesce(isl_map_remove_redundancies(*union_map));
+    
+        return isl_stat_ok;
+}
+
+/* Compute the transitive closure of union map incrementally by
+ * composing a union of transitive closures of each basic map 
+ * from union map. For each iteration we first extract a list of 
+ * basic maps from the union map. Then we calculate transitive 
+ * closure for each map separately and compose themselves. 
+ * If the result is exact, then *exact is set to 1, 
+ * 0 for underapproximation. 
+*/
+static __isl_give isl_map *transitive_closure_iterate(__isl_take isl_map* map, 
+        unsigned max_iterations, 
+        unsigned max_disjunctions, 
+        int* exact) 
+{
+        isl_map* out = NULL;
+        unsigned n = 0;
+        
+        isl_map* identity = isl_map_identity(isl_map_get_space(map));
+
+        isl_map_foreach_basic_map(map, union_closure, (void*)&out);
+        out = isl_map_union(out, isl_map_copy(identity));
+        
+        
+        isl_map* compose = isl_map_apply_range(isl_map_copy(out), 
+                isl_map_copy(out));
+        compose = isl_map_compute_divs(compose);
+        compose = isl_map_coalesce(isl_map_remove_redundancies(compose));
+
+        *exact = isl_map_is_equal(compose, out);
+               
+        while(!(*exact) && (compose->n < max_disjunctions) && 
+                (n < max_iterations)) 
+        {
+                compose = isl_map_subtract(compose, isl_map_copy(identity));
+                compose = isl_map_compute_divs(compose);
+                compose = isl_map_remove_redundancies(compose);
+                compose = isl_map_coalesce(compose);
+                
+                out = (isl_map_free(out), NULL);
+                isl_map_foreach_basic_map(compose, union_closure, (void*)&out);
+                out = isl_map_union(out, isl_map_copy(identity));
+    
+                compose = isl_map_union(compose, isl_map_copy(identity));
+                compose = isl_map_apply_range(compose, isl_map_copy(out));
+                compose = isl_map_remove_redundancies(compose);
+                compose = isl_map_coalesce(compose);
+            
+                *exact = isl_map_is_equal(compose, out);
+        }
+        
+        out = isl_map_subtract(out, identity);
+        out = isl_map_compute_divs(out);
+        out = isl_map_coalesce(isl_map_remove_redundancies(out));
+        
+        isl_map_free(map);
+        isl_map_free(compose);
+        
+        return out;
+}
+
+
 /* Compute the transitive closure  of "map", or an overapproximation.
  * If the result is exact, then *exact is set to 1.
  * Simply use map_power to compute the powers of map, but tell
@@ -2589,6 +2672,9 @@
 
 	if (map->ctx->opt->closure == ISL_CLOSURE_BOX)
 		return transitive_closure_omega(map, exact);
+        
+	if (map->ctx->opt->closure == ISL_CLOSURE_ITER)
+		return transitive_closure_iterate(map, 10, 100, exact);       
 
 	map = isl_map_compute_divs(map);
 	map = isl_map_coalesce(map);