X-Git-Url: http://review.tizen.org/git/?a=blobdiff_plain;f=isl_tab_pip.c;h=2592a285c0e179164254439f4b2c0e0d84bd12ed;hb=63fb8a7f484648c3caa25351c8c94ac2395ec563;hp=073380e021dee8761ae5971582dcefd9d217aa2f;hpb=bbf85ab5afe150c925b01e126c438784b015230a;p=platform%2Fupstream%2Fisl.git

diff --git a/isl_tab_pip.c b/isl_tab_pip.c
index 073380e..2592a28 100644
--- a/isl_tab_pip.c
+++ b/isl_tab_pip.c
@@ -2,7 +2,7 @@
  * Copyright 2008-2009 Katholieke Universiteit Leuven
  * Copyright 2010      INRIA Saclay
  *
- * Use of this software is governed by the GNU LGPLv2.1 license
+ * Use of this software is governed by the MIT license
  *
  * Written by Sven Verdoolaege, K.U.Leuven, Departement
  * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
@@ -17,6 +17,7 @@
 #include "isl_sample.h"
 #include <isl_mat_private.h>
 #include <isl_aff_private.h>
+#include <isl_options_private.h>
 #include <isl_config.h>
 
 /*
@@ -101,6 +102,8 @@ struct isl_context_op {
 	void *(*save)(struct isl_context *context);
 	/* restore saved context */
 	void (*restore)(struct isl_context *context, void *);
+	/* discard saved context */
+	void (*discard)(void *);
 	/* invalidate context */
 	void (*invalidate)(struct isl_context *context);
 	/* free context */
@@ -116,6 +119,14 @@ struct isl_context_lex {
 	struct isl_tab *tab;
 };
 
+/* A stack (linked list) of solutions of subtrees of the search space.
+ *
+ * "M" describes the solution in terms of the dimensions of "dom".
+ * The number of columns of "M" is one more than the total number
+ * of dimensions of "dom".
+ *
+ * If "M" is NULL, then there is no solution on "dom".
+ */
 struct isl_partial_sol {
 	int level;
 	struct isl_basic_set *dom;
@@ -200,6 +211,7 @@ static void sol_push_sol(struct isl_sol *sol,
 	return;
 error:
 	isl_basic_set_free(dom);
+	isl_mat_free(M);
 	sol->error = 1;
 }
 
@@ -303,13 +315,26 @@ static void sol_pop(struct isl_sol *sol)
 			sol_pop_one(sol);
 		} else {
 			struct isl_basic_set *bset;
+			isl_mat *M;
+			unsigned n;
 
+			n = isl_basic_set_dim(partial->next->dom, isl_dim_div);
+			n -= n_div;
 			bset = sol_domain(sol);
-
 			isl_basic_set_free(partial->next->dom);
 			partial->next->dom = bset;
+			M = partial->next->M;
+			if (M) {
+				M = isl_mat_drop_cols(M, M->n_col - n, n);
+				partial->next->M = M;
+				if (!M)
+					goto error;
+			}
 			partial->next->level = sol->level;
 
+			if (!bset)
+				goto error;
+
 			sol->partial = partial->next;
 			isl_basic_set_free(partial->dom);
 			isl_mat_free(partial->M);
@@ -317,6 +342,9 @@ static void sol_pop(struct isl_sol *sol)
 		}
 	} else
 		sol_pop_one(sol);
+
+	if (0)
+error:		sol->error = 1;
 }
 
 static void sol_dec_level(struct isl_sol *sol)
@@ -430,6 +458,8 @@ static void sol_add(struct isl_sol *sol, struct isl_tab *tab)
 
 	if (tab->empty && !sol->add_empty)
 		return;
+	if (sol->context->op->is_empty(sol->context))
+		return;
 
 	bset = sol_domain(sol);
 
@@ -556,18 +586,6 @@ static void sol_map_add_empty_wrap(struct isl_sol *sol,
 	sol_map_add_empty((struct isl_sol_map *)sol, bset);
 }
 
-/* Add bset to sol's empty, but only if we are actually collecting
- * the empty set.
- */
-static void sol_map_add_empty_if_needed(struct isl_sol_map *sol,
-	struct isl_basic_set *bset)
-{
-	if (sol->empty)
-		sol_map_add_empty(sol, bset);
-	else
-		isl_basic_set_free(bset);
-}
-
 /* Given a basic map "dom" that represents the context and an affine
  * matrix "M" that maps the dimensions of the context to the
  * output variables, construct a basic map with the same parameters
@@ -603,7 +621,7 @@ static void sol_map_add(struct isl_sol_map *sol,
 	n_div = dom->n_div;
 	nparam = isl_basic_set_total_dim(dom) - n_div;
 	total = isl_map_dim(sol->map, isl_dim_all);
-	bmap = isl_basic_map_alloc_dim(isl_map_get_dim(sol->map),
+	bmap = isl_basic_map_alloc_space(isl_map_get_space(sol->map),
 					n_div, n_eq, 2 * n_div + n_ineq);
 	if (!bmap)
 		goto error;
@@ -650,10 +668,10 @@ static void sol_map_add(struct isl_sol_map *sol,
 	bmap = isl_basic_map_finalize(bmap);
 	sol->map = isl_map_grow(sol->map, 1);
 	sol->map = isl_map_add_basic_map(sol->map, bmap);
-	if (!sol->map)
-		goto error;
 	isl_basic_set_free(dom);
 	isl_mat_free(M);
+	if (!sol->map)
+		sol->sol.error = 1;
 	return;
 error:
 	isl_basic_set_free(dom);
@@ -754,6 +772,30 @@ static struct isl_vec *get_row_parameter_ineq(struct isl_tab *tab, int row)
 	return ineq;
 }
 
+/* Normalize a div expression of the form
+ *
+ *	[(g*f(x) + c)/(g * m)]
+ *
+ * with c the constant term and f(x) the remaining coefficients, to
+ *
+ *	[(f(x) + [c/g])/m]
+ */
+static void normalize_div(__isl_keep isl_vec *div)
+{
+	isl_ctx *ctx = isl_vec_get_ctx(div);
+	int len = div->size - 2;
+
+	isl_seq_gcd(div->el + 2, len, &ctx->normalize_gcd);
+	isl_int_gcd(ctx->normalize_gcd, ctx->normalize_gcd, div->el[0]);
+
+	if (isl_int_is_one(ctx->normalize_gcd))
+		return;
+
+	isl_int_divexact(div->el[0], div->el[0], ctx->normalize_gcd);
+	isl_int_fdiv_q(div->el[1], div->el[1], ctx->normalize_gcd);
+	isl_seq_scale_down(div->el + 2, div->el + 2, ctx->normalize_gcd, len);
+}
+
 /* Return a integer division for use in a parametric cut based on the given row.
  * In particular, let the parametric constant of the row be
  *
@@ -774,8 +816,8 @@ static struct isl_vec *get_row_parameter_div(struct isl_tab *tab, int row)
 
 	isl_int_set(div->el[0], tab->mat->row[row][0]);
 	get_row_parameter_line(tab, row, div->el + 1);
-	div = isl_vec_normalize(div);
 	isl_seq_neg(div->el + 1, div->el + 1, div->size - 1);
+	normalize_div(div);
 	isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
 
 	return div;
@@ -802,7 +844,7 @@ static struct isl_vec *get_row_split_div(struct isl_tab *tab, int row)
 
 	isl_int_set(div->el[0], tab->mat->row[row][0]);
 	get_row_parameter_line(tab, row, div->el + 1);
-	div = isl_vec_normalize(div);
+	normalize_div(div);
 	isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
 
 	return div;
@@ -840,7 +882,7 @@ static struct isl_vec *ineq_for_div(struct isl_basic_set *bset, unsigned div)
 }
 
 /* Given a row in the tableau and a div that was created
- * using get_row_split_div and that been constrained to equality, i.e.,
+ * using get_row_split_div and that has been constrained to equality, i.e.,
  *
  *		d = floor(\sum_i {a_i} y_i) = \sum_i {a_i} y_i
  *
@@ -871,7 +913,8 @@ static struct isl_tab *set_row_cst_to_div(struct isl_tab *tab, int row, int div)
 	} else {
 		int dcol = tab->var[tab->n_var - tab->n_div + div].index;
 
-		isl_int_set_si(tab->mat->row[row][2 + tab->M + dcol], 1);
+		isl_int_add_ui(tab->mat->row[row][2 + tab->M + dcol],
+				tab->mat->row[row][2 + tab->M + dcol], 1);
 	}
 
 	return tab;
@@ -1161,8 +1204,8 @@ static int report_conflicting_constraint(struct isl_tab *tab, int con)
 
 /* Given a conflicting row in the tableau, report all constraints
  * involved in the row to the caller.  That is, the row itself
- * (if represents a constraint) and all constraint columns with
- * non-zero (and therefore negative) coefficient.
+ * (if it represents a constraint) and all constraint columns with
+ * non-zero (and therefore negative) coefficients.
  */
 static int report_conflict(struct isl_tab *tab, int row)
 {
@@ -1614,6 +1657,9 @@ static int add_cut(struct isl_tab *tab, int row)
 	return tab->con[r].index;
 }
 
+#define CUT_ALL 1
+#define CUT_ONE 0
+
 /* Given a non-parametric tableau, add cuts until an integer
  * sample point is obtained or until the tableau is determined
  * to be integer infeasible.
@@ -1625,8 +1671,12 @@ static int add_cut(struct isl_tab *tab, int row)
  * combination of variables/constraints plus a non-integral constant,
  * then there is no way to obtain an integer point and we return
  * a tableau that is marked empty.
+ * The parameter cutting_strategy controls the strategy used when adding cuts
+ * to remove non-integer points. CUT_ALL adds all possible cuts
+ * before continuing the search. CUT_ONE adds only one cut at a time.
  */
-static struct isl_tab *cut_to_integer_lexmin(struct isl_tab *tab)
+static struct isl_tab *cut_to_integer_lexmin(struct isl_tab *tab,
+	int cutting_strategy)
 {
 	int var;
 	int row;
@@ -1648,6 +1698,8 @@ static struct isl_tab *cut_to_integer_lexmin(struct isl_tab *tab)
 			row = add_cut(tab, row);
 			if (row < 0)
 				goto error;
+			if (cutting_strategy == CUT_ONE)
+				break;
 		} while ((var = next_non_integer_var(tab, var, &flags)) != -1);
 		if (restore_lexmin(tab) < 0)
 			goto error;
@@ -1738,7 +1790,7 @@ static struct isl_tab *check_integer_feasible(struct isl_tab *tab)
 	if (isl_tab_push_basis(tab) < 0)
 		goto error;
 
-	tab = cut_to_integer_lexmin(tab);
+	tab = cut_to_integer_lexmin(tab, CUT_ALL);
 	if (!tab)
 		goto error;
 
@@ -1932,6 +1984,7 @@ static int add_parametric_cut(struct isl_tab *tab, int row,
 
 	n = tab->n_div;
 	d = context->op->get_div(context, tab, div);
+	isl_vec_free(div);
 	if (d < 0)
 		return -1;
 
@@ -1997,8 +2050,6 @@ static int add_parametric_cut(struct isl_tab *tab, int row,
 	if (tab->row_sign)
 		tab->row_sign[tab->con[r].index] = isl_tab_row_neg;
 
-	isl_vec_free(div);
-
 	row = tab->con[r].index;
 
 	if (d >= n && context->op->detect_equalities(context, tab) < 0)
@@ -2393,6 +2444,10 @@ static void context_lex_restore(struct isl_context *context, void *save)
 	}
 }
 
+static void context_lex_discard(void *save)
+{
+}
+
 static int context_lex_is_ok(struct isl_context *context)
 {
 	struct isl_context_lex *clex = (struct isl_context_lex *)context;
@@ -2512,6 +2567,7 @@ struct isl_context_op isl_context_lex_op = {
 	context_lex_is_ok,
 	context_lex_save,
 	context_lex_restore,
+	context_lex_discard,
 	context_lex_invalidate,
 	context_lex_free,
 };
@@ -2520,7 +2576,6 @@ static struct isl_tab *context_tab_for_lexmin(struct isl_basic_set *bset)
 {
 	struct isl_tab *tab;
 
-	bset = isl_basic_set_cow(bset);
 	if (!bset)
 		return NULL;
 	tab = tab_for_lexmin((struct isl_basic_map *)bset, NULL, 1, 0);
@@ -2561,6 +2616,14 @@ error:
 	return NULL;
 }
 
+/* Representation of the context when using generalized basis reduction.
+ *
+ * "shifted" contains the offsets of the unit hypercubes that lie inside the
+ * context.  Any rational point in "shifted" can therefore be rounded
+ * up to an integer point in the context.
+ * If the context is constrained by any equality, then "shifted" is not used
+ * as it would be empty.
+ */
 struct isl_context_gbr {
 	struct isl_context context;
 	struct isl_tab *tab;
@@ -2619,7 +2682,7 @@ static void gbr_init_shifted(struct isl_context_gbr *cgbr)
 		}
 	}
 
-	cgbr->shifted = isl_tab_from_basic_set(bset);
+	cgbr->shifted = isl_tab_from_basic_set(bset, 0);
 
 	for (i = 0; i < bset->n_ineq; ++i)
 		isl_int_set(bset->ineq[i][0], cst->el[i]);
@@ -2692,7 +2755,8 @@ static struct isl_vec *gbr_get_sample(struct isl_context_gbr *cgbr)
 		cgbr->cone = isl_tab_from_recession_cone(bset, 0);
 		if (!cgbr->cone)
 			return NULL;
-		if (isl_tab_track_bset(cgbr->cone, isl_basic_set_dup(bset)) < 0)
+		if (isl_tab_track_bset(cgbr->cone,
+					isl_basic_set_copy(bset)) < 0)
 			return NULL;
 	}
 	if (isl_tab_detect_implicit_equalities(cgbr->cone) < 0)
@@ -2783,6 +2847,15 @@ error:
 	return NULL;
 }
 
+/* Add the equality described by "eq" to the context.
+ * If "check" is set, then we check if the context is empty after
+ * adding the equality.
+ * If "update" is set, then we check if the samples are still valid.
+ *
+ * We do not explicitly add shifted copies of the equality to
+ * cgbr->shifted since they would conflict with each other.
+ * Instead, we directly mark cgbr->shifted empty.
+ */
 static void context_gbr_add_eq(struct isl_context *context, isl_int *eq,
 		int check, int update)
 {
@@ -2790,6 +2863,11 @@ static void context_gbr_add_eq(struct isl_context *context, isl_int *eq,
 
 	cgbr->tab = add_gbr_eq(cgbr->tab, eq);
 
+	if (cgbr->shifted && !cgbr->shifted->empty && use_shifted(cgbr)) {
+		if (isl_tab_mark_empty(cgbr->shifted) < 0)
+			goto error;
+	}
+
 	if (cgbr->cone && cgbr->cone->n_col != cgbr->cone->n_dead) {
 		if (isl_tab_extend_cons(cgbr->cone, 2) < 0)
 			goto error;
@@ -3054,7 +3132,8 @@ static int context_gbr_detect_equalities(struct isl_context *context,
 		cgbr->cone = isl_tab_from_recession_cone(bset, 0);
 		if (!cgbr->cone)
 			goto error;
-		if (isl_tab_track_bset(cgbr->cone, isl_basic_set_dup(bset)) < 0)
+		if (isl_tab_track_bset(cgbr->cone,
+					isl_basic_set_copy(bset)) < 0)
 			goto error;
 	}
 	if (isl_tab_detect_implicit_equalities(cgbr->cone) < 0)
@@ -3062,7 +3141,9 @@ static int context_gbr_detect_equalities(struct isl_context *context,
 
 	n_ineq = cgbr->tab->bmap->n_ineq;
 	cgbr->tab = isl_tab_detect_equalities(cgbr->tab, cgbr->cone);
-	if (cgbr->tab && cgbr->tab->bmap->n_ineq > n_ineq)
+	if (!cgbr->tab)
+		return -1;
+	if (cgbr->tab->bmap->n_ineq > n_ineq)
 		propagate_equalities(cgbr, tab, n_ineq);
 
 	return 0;
@@ -3091,8 +3172,8 @@ static int context_gbr_add_div(struct isl_context *context, struct isl_vec *div)
 		if (isl_tab_allocate_var(cgbr->cone) <0)
 			return -1;
 
-		cgbr->cone->bmap = isl_basic_map_extend_dim(cgbr->cone->bmap,
-			isl_basic_map_get_dim(cgbr->cone->bmap), 1, 0, 2);
+		cgbr->cone->bmap = isl_basic_map_extend_space(cgbr->cone->bmap,
+			isl_basic_map_get_space(cgbr->cone->bmap), 1, 0, 2);
 		k = isl_basic_map_alloc_div(cgbr->cone->bmap);
 		if (k < 0)
 			return -1;
@@ -3199,6 +3280,12 @@ error:
 	cgbr->tab = NULL;
 }
 
+static void context_gbr_discard(void *save)
+{
+	struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save;
+	free(snap);
+}
+
 static int context_gbr_is_ok(struct isl_context *context)
 {
 	struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
@@ -3237,6 +3324,7 @@ struct isl_context_op isl_context_gbr_op = {
 	context_gbr_is_ok,
 	context_gbr_save,
 	context_gbr_restore,
+	context_gbr_discard,
 	context_gbr_invalidate,
 	context_gbr_free,
 };
@@ -3256,13 +3344,10 @@ static struct isl_context *isl_context_gbr_alloc(struct isl_basic_set *dom)
 
 	cgbr->shifted = NULL;
 	cgbr->cone = NULL;
-	cgbr->tab = isl_tab_from_basic_set(dom);
+	cgbr->tab = isl_tab_from_basic_set(dom, 1);
 	cgbr->tab = isl_tab_init_samples(cgbr->tab);
 	if (!cgbr->tab)
 		goto error;
-	if (isl_tab_track_bset(cgbr->tab,
-				isl_basic_set_cow(isl_basic_set_copy(dom))) < 0)
-		goto error;
 	check_gbr_integer_feasible(cgbr);
 
 	return &cgbr->context;
@@ -3289,7 +3374,7 @@ static struct isl_context *isl_context_alloc(struct isl_basic_set *dom)
  * a minimization problem, which means that the variables in the
  * tableau have value "M - x" rather than "M + x".
  */
-static struct isl_sol_map *sol_map_init(struct isl_basic_map *bmap,
+static struct isl_sol *sol_map_init(struct isl_basic_map *bmap,
 	struct isl_basic_set *dom, int track_empty, int max)
 {
 	struct isl_sol_map *sol_map = NULL;
@@ -3309,7 +3394,7 @@ static struct isl_sol_map *sol_map_init(struct isl_basic_map *bmap,
 	sol_map->sol.add = &sol_map_add_wrap;
 	sol_map->sol.add_empty = track_empty ? &sol_map_add_empty_wrap : NULL;
 	sol_map->sol.free = &sol_map_free_wrap;
-	sol_map->map = isl_map_alloc_dim(isl_basic_map_get_dim(bmap), 1,
+	sol_map->map = isl_map_alloc_space(isl_basic_map_get_space(bmap), 1,
 					    ISL_MAP_DISJOINT);
 	if (!sol_map->map)
 		goto error;
@@ -3319,14 +3404,14 @@ static struct isl_sol_map *sol_map_init(struct isl_basic_map *bmap,
 		goto error;
 
 	if (track_empty) {
-		sol_map->empty = isl_set_alloc_dim(isl_basic_set_get_dim(dom),
+		sol_map->empty = isl_set_alloc_space(isl_basic_set_get_space(dom),
 							1, ISL_SET_DISJOINT);
 		if (!sol_map->empty)
 			goto error;
 	}
 
 	isl_basic_set_free(dom);
-	return sol_map;
+	return &sol_map->sol;
 error:
 	isl_basic_set_free(dom);
 	sol_map_free(sol_map);
@@ -3538,6 +3623,8 @@ static void find_in_pos(struct isl_sol *sol, struct isl_tab *tab, isl_int *ineq)
 
 	if (!sol->error)
 		sol->context->op->restore(sol->context, saved);
+	else
+		sol->context->op->discard(saved);
 	return;
 error:
 	sol->error = 1;
@@ -3789,6 +3876,24 @@ error:
 	sol->error = 1;
 }
 
+/* Does "sol" contain a pair of partial solutions that could potentially
+ * be merged?
+ *
+ * We currently only check that "sol" is not in an error state
+ * and that there are at least two partial solutions of which the final two
+ * are defined at the same level.
+ */
+static int sol_has_mergeable_solutions(struct isl_sol *sol)
+{
+	if (sol->error)
+		return 0;
+	if (!sol->partial)
+		return 0;
+	if (!sol->partial->next)
+		return 0;
+	return sol->partial->level == sol->partial->next->level;
+}
+
 /* Compute the lexicographic minimum of the set represented by the main
  * tableau "tab" within the context "sol->context_tab".
  *
@@ -3799,10 +3904,20 @@ error:
  * corresponding rows may not be marked as being non-negative.
  * In parts of the context where the added equality does not hold,
  * the main tableau is marked as being empty.
+ *
+ * Before we embark on the actual computation, we save a copy
+ * of the context.  When we return, we check if there are any
+ * partial solutions that can potentially be merged.  If so,
+ * we perform a rollback to the initial state of the context.
+ * The merging of partial solutions happens inside calls to
+ * sol_dec_level that are pushed onto the undo stack of the context.
+ * If there are no partial solutions that can potentially be merged
+ * then the rollback is skipped as it would just be wasted effort.
  */
 static void find_solutions_main(struct isl_sol *sol, struct isl_tab *tab)
 {
 	int row;
+	void *saved;
 
 	if (!tab)
 		goto error;
@@ -3852,8 +3967,15 @@ static void find_solutions_main(struct isl_sol *sol, struct isl_tab *tab)
 		row = tab->n_redundant - 1;
 	}
 
+	saved = sol->context->op->save(sol->context);
+
 	find_solutions(sol, tab);
 
+	if (sol_has_mergeable_solutions(sol))
+		sol->context->op->restore(sol->context, saved);
+	else
+		sol->context->op->discard(saved);
+
 	sol->level = 0;
 	sol_pop(sol);
 
@@ -3863,12 +3985,6 @@ error:
 	sol->error = 1;
 }
 
-static void sol_map_find_solutions(struct isl_sol_map *sol_map,
-	struct isl_tab *tab)
-{
-	find_solutions_main(&sol_map->sol, tab);
-}
-
 /* Check if integer division "div" of "dom" also occurs in "bmap".
  * If so, return its position within the divs.
  * If not, return -1.
@@ -3877,8 +3993,8 @@ static int find_context_div(struct isl_basic_map *bmap,
 	struct isl_basic_set *dom, unsigned div)
 {
 	int i;
-	unsigned b_dim = isl_dim_total(bmap->dim);
-	unsigned d_dim = isl_dim_total(dom->dim);
+	unsigned b_dim = isl_space_dim(bmap->dim, isl_dim_all);
+	unsigned d_dim = isl_space_dim(dom->dim, isl_dim_all);
 
 	if (isl_int_is_zero(dom->div[div][0]))
 		return -1;
@@ -3923,7 +4039,7 @@ static struct isl_basic_map *align_context_divs(struct isl_basic_map *bmap,
 			common++;
 	other = bmap->n_div - common;
 	if (dom->n_div - common > 0) {
-		bmap = isl_basic_map_extend_dim(bmap, isl_dim_copy(bmap->dim),
+		bmap = isl_basic_map_extend_space(bmap, isl_space_copy(bmap->dim),
 				dom->n_div - common, 0, 0);
 		if (!bmap)
 			return NULL;
@@ -3952,48 +4068,72 @@ error:
  * because they will be added one by one in the given order
  * during the construction of the solution map.
  */
-static __isl_give isl_map *basic_map_partial_lexopt_base(
+static struct isl_sol *basic_map_partial_lexopt_base(
 	__isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
-	__isl_give isl_set **empty, int max)
+	__isl_give isl_set **empty, int max,
+	struct isl_sol *(*init)(__isl_keep isl_basic_map *bmap,
+		    __isl_take isl_basic_set *dom, int track_empty, int max))
 {
-	isl_map *result = NULL;
 	struct isl_tab *tab;
-	struct isl_sol_map *sol_map = NULL;
+	struct isl_sol *sol = NULL;
 	struct isl_context *context;
 
 	if (dom->n_div) {
 		dom = isl_basic_set_order_divs(dom);
 		bmap = align_context_divs(bmap, dom);
 	}
-	sol_map = sol_map_init(bmap, dom, !!empty, max);
-	if (!sol_map)
+	sol = init(bmap, dom, !!empty, max);
+	if (!sol)
 		goto error;
 
-	context = sol_map->sol.context;
+	context = sol->context;
 	if (isl_basic_set_plain_is_empty(context->op->peek_basic_set(context)))
 		/* nothing */;
-	else if (isl_basic_map_plain_is_empty(bmap))
-		sol_map_add_empty_if_needed(sol_map,
+	else if (isl_basic_map_plain_is_empty(bmap)) {
+		if (sol->add_empty)
+			sol->add_empty(sol,
 		    isl_basic_set_copy(context->op->peek_basic_set(context)));
-	else {
+	} else {
 		tab = tab_for_lexmin(bmap,
 				    context->op->peek_basic_set(context), 1, max);
 		tab = context->op->detect_nonnegative_parameters(context, tab);
-		sol_map_find_solutions(sol_map, tab);
+		find_solutions_main(sol, tab);
 	}
-	if (sol_map->sol.error)
+	if (sol->error)
 		goto error;
 
+	isl_basic_map_free(bmap);
+	return sol;
+error:
+	sol_free(sol);
+	isl_basic_map_free(bmap);
+	return NULL;
+}
+
+/* Base case of isl_tab_basic_map_partial_lexopt, after removing
+ * some obvious symmetries.
+ *
+ * We call basic_map_partial_lexopt_base and extract the results.
+ */
+static __isl_give isl_map *basic_map_partial_lexopt_base_map(
+	__isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
+	__isl_give isl_set **empty, int max)
+{
+	isl_map *result = NULL;
+	struct isl_sol *sol;
+	struct isl_sol_map *sol_map;
+
+	sol = basic_map_partial_lexopt_base(bmap, dom, empty, max,
+					    &sol_map_init);
+	if (!sol)
+		return NULL;
+	sol_map = (struct isl_sol_map *) sol;
+
 	result = isl_map_copy(sol_map->map);
 	if (empty)
 		*empty = isl_set_copy(sol_map->empty);
 	sol_free(&sol_map->sol);
-	isl_basic_map_free(bmap);
 	return result;
-error:
-	sol_free(&sol_map->sol);
-	isl_basic_map_free(bmap);
-	return NULL;
 }
 
 /* Structure used during detection of parallel constraints.
@@ -4081,6 +4221,43 @@ error:
 	return -1;
 }
 
+/* Given a set of upper bounds in "var", add constraints to "bset"
+ * that make the i-th bound smallest.
+ *
+ * In particular, if there are n bounds b_i, then add the constraints
+ *
+ *	b_i <= b_j	for j > i
+ *	b_i <  b_j	for j < i
+ */
+static __isl_give isl_basic_set *select_minimum(__isl_take isl_basic_set *bset,
+	__isl_keep isl_mat *var, int i)
+{
+	isl_ctx *ctx;
+	int j, k;
+
+	ctx = isl_mat_get_ctx(var);
+
+	for (j = 0; j < var->n_row; ++j) {
+		if (j == i)
+			continue;
+		k = isl_basic_set_alloc_inequality(bset);
+		if (k < 0)
+			goto error;
+		isl_seq_combine(bset->ineq[k], ctx->one, var->row[j],
+				ctx->negone, var->row[i], var->n_col);
+		isl_int_set_si(bset->ineq[k][var->n_col], 0);
+		if (j < i)
+			isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
+	}
+
+	bset = isl_basic_set_finalize(bset);
+
+	return bset;
+error:
+	isl_basic_set_free(bset);
+	return NULL;
+}
+
 /* Given a set of upper bounds on the last "input" variable m,
  * construct a set that assigns the minimal upper bound to m, i.e.,
  * construct a set that divides the space into cells where one
@@ -4094,10 +4271,10 @@ error:
  *	b_i <= b_j	for j > i
  *	b_i <  b_j	for j < i
  */
-static __isl_give isl_set *set_minimum(__isl_take isl_dim *dim,
+static __isl_give isl_set *set_minimum(__isl_take isl_space *dim,
 	__isl_take isl_mat *var)
 {
-	int i, j, k;
+	int i, k;
 	isl_basic_set *bset = NULL;
 	isl_ctx *ctx;
 	isl_set *set = NULL;
@@ -4105,43 +4282,29 @@ static __isl_give isl_set *set_minimum(__isl_take isl_dim *dim,
 	if (!dim || !var)
 		goto error;
 
-	ctx = isl_dim_get_ctx(dim);
-	set = isl_set_alloc_dim(isl_dim_copy(dim),
+	ctx = isl_space_get_ctx(dim);
+	set = isl_set_alloc_space(isl_space_copy(dim),
 				var->n_row, ISL_SET_DISJOINT);
 
 	for (i = 0; i < var->n_row; ++i) {
-		bset = isl_basic_set_alloc_dim(isl_dim_copy(dim), 0,
+		bset = isl_basic_set_alloc_space(isl_space_copy(dim), 0,
 					       1, var->n_row - 1);
 		k = isl_basic_set_alloc_equality(bset);
 		if (k < 0)
 			goto error;
 		isl_seq_cpy(bset->eq[k], var->row[i], var->n_col);
 		isl_int_set_si(bset->eq[k][var->n_col], -1);
-		for (j = 0; j < var->n_row; ++j) {
-			if (j == i)
-				continue;
-			k = isl_basic_set_alloc_inequality(bset);
-			if (k < 0)
-				goto error;
-			isl_seq_combine(bset->ineq[k], ctx->one, var->row[j],
-					ctx->negone, var->row[i],
-					var->n_col);
-			isl_int_set_si(bset->ineq[k][var->n_col], 0);
-			if (j < i)
-				isl_int_sub_ui(bset->ineq[k][0],
-					       bset->ineq[k][0], 1);
-		}
-		bset = isl_basic_set_finalize(bset);
+		bset = select_minimum(bset, var, i);
 		set = isl_set_add_basic_set(set, bset);
 	}
 
-	isl_dim_free(dim);
+	isl_space_free(dim);
 	isl_mat_free(var);
 	return set;
 error:
 	isl_basic_set_free(bset);
 	isl_set_free(set);
-	isl_dim_free(dim);
+	isl_space_free(dim);
 	isl_mat_free(var);
 	return NULL;
 }
@@ -4155,7 +4318,7 @@ error:
  * an upper bound that is different from the upper bounds on which it
  * is defined.
  */
-static int need_split_map(__isl_keep isl_basic_map *bmap,
+static int need_split_basic_map(__isl_keep isl_basic_map *bmap,
 	__isl_keep isl_mat *cst)
 {
 	int i, j;
@@ -4192,10 +4355,33 @@ static int need_split_map(__isl_keep isl_basic_map *bmap,
 	return 0;
 }
 
-static int need_split_set(__isl_keep isl_basic_set *bset,
+/* Given that the last set variable of "bset" represents the minimum
+ * of the bounds in "cst", check whether we need to split the domain
+ * based on which bound attains the minimum.
+ *
+ * We simply call need_split_basic_map here.  This is safe because
+ * the position of the minimum is computed from "cst" and not
+ * from "bmap".
+ */
+static int need_split_basic_set(__isl_keep isl_basic_set *bset,
 	__isl_keep isl_mat *cst)
 {
-	return need_split_map((isl_basic_map *)bset, cst);
+	return need_split_basic_map((isl_basic_map *)bset, cst);
+}
+
+/* Given that the last set variable of "set" represents the minimum
+ * of the bounds in "cst", check whether we need to split the domain
+ * based on which bound attains the minimum.
+ */
+static int need_split_set(__isl_keep isl_set *set, __isl_keep isl_mat *cst)
+{
+	int i;
+
+	for (i = 0; i < set->n; ++i)
+		if (need_split_basic_set(set->p[i], cst))
+			return 1;
+
+	return 0;
 }
 
 /* Given a set of which the last set variable is the minimum
@@ -4216,22 +4402,22 @@ static __isl_give isl_set *split(__isl_take isl_set *empty,
 {
 	int n_in;
 	int i;
-	isl_dim *dim;
+	isl_space *dim;
 	isl_set *res;
 
 	if (!empty || !min_expr || !cst)
 		goto error;
 
 	n_in = isl_set_dim(empty, isl_dim_set);
-	dim = isl_set_get_dim(empty);
-	dim = isl_dim_drop(dim, isl_dim_set, n_in - 1, 1);
+	dim = isl_set_get_space(empty);
+	dim = isl_space_drop_dims(dim, isl_dim_set, n_in - 1, 1);
 	res = isl_set_empty(dim);
 
 	for (i = 0; i < empty->n; ++i) {
 		isl_set *set;
 
 		set = isl_set_from_basic_set(isl_basic_set_copy(empty->p[i]));
-		if (need_split_set(empty->p[i], cst))
+		if (need_split_basic_set(empty->p[i], cst))
 			set = isl_set_intersect(set, isl_set_copy(min_expr));
 		set = isl_set_remove_dims(set, isl_dim_set, n_in - 1, 1);
 
@@ -4262,22 +4448,22 @@ static __isl_give isl_map *split_domain(__isl_take isl_map *opt,
 {
 	int n_in;
 	int i;
-	isl_dim *dim;
+	isl_space *dim;
 	isl_map *res;
 
 	if (!opt || !min_expr || !cst)
 		goto error;
 
 	n_in = isl_map_dim(opt, isl_dim_in);
-	dim = isl_map_get_dim(opt);
-	dim = isl_dim_drop(dim, isl_dim_in, n_in - 1, 1);
+	dim = isl_map_get_space(opt);
+	dim = isl_space_drop_dims(dim, isl_dim_in, n_in - 1, 1);
 	res = isl_map_empty(dim);
 
 	for (i = 0; i < opt->n; ++i) {
 		isl_map *map;
 
 		map = isl_map_from_basic_map(isl_basic_map_copy(opt->p[i]));
-		if (need_split_map(opt->p[i], cst))
+		if (need_split_basic_map(opt->p[i], cst))
 			map = isl_map_intersect_domain(map,
 						       isl_set_copy(min_expr));
 		map = isl_map_remove_dims(map, isl_dim_in, n_in - 1, 1);
@@ -4300,6 +4486,47 @@ static __isl_give isl_map *basic_map_partial_lexopt(
 	__isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
 	__isl_give isl_set **empty, int max);
 
+union isl_lex_res {
+	void *p;
+	isl_map *map;
+	isl_pw_multi_aff *pma;
+};
+
+/* This function is called from basic_map_partial_lexopt_symm.
+ * The last variable of "bmap" and "dom" corresponds to the minimum
+ * of the bounds in "cst".  "map_space" is the space of the original
+ * input relation (of basic_map_partial_lexopt_symm) and "set_space"
+ * is the space of the original domain.
+ *
+ * We recursively call basic_map_partial_lexopt and then plug in
+ * the definition of the minimum in the result.
+ */
+static __isl_give union isl_lex_res basic_map_partial_lexopt_symm_map_core(
+	__isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
+	__isl_give isl_set **empty, int max, __isl_take isl_mat *cst,
+	__isl_take isl_space *map_space, __isl_take isl_space *set_space)
+{
+	isl_map *opt;
+	isl_set *min_expr;
+	union isl_lex_res res;
+
+	min_expr = set_minimum(isl_basic_set_get_space(dom), isl_mat_copy(cst));
+
+	opt = basic_map_partial_lexopt(bmap, dom, empty, max);
+
+	if (empty) {
+		*empty = split(*empty,
+			       isl_set_copy(min_expr), isl_mat_copy(cst));
+		*empty = isl_set_reset_space(*empty, set_space);
+	}
+
+	opt = split_domain(opt, min_expr, cst);
+	opt = isl_map_reset_space(opt, map_space);
+
+	res.map = opt;
+	return res;
+}
+
 /* Given a basic map with at least two parallel constraints (as found
  * by the function parallel_constraints), first look for more constraints
  * parallel to the two constraint and replace the found list of parallel
@@ -4328,9 +4555,15 @@ static __isl_give isl_map *basic_map_partial_lexopt(
  * Moreover, m = min_i(b_i(p)) satisfies the constraints on u and can
  * therefore be plugged into the solution.
  */
-static __isl_give isl_map *basic_map_partial_lexopt_symm(
+static union isl_lex_res basic_map_partial_lexopt_symm(
 	__isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
-	__isl_give isl_set **empty, int max, int first, int second)
+	__isl_give isl_set **empty, int max, int first, int second,
+	__isl_give union isl_lex_res (*core)(__isl_take isl_basic_map *bmap,
+					    __isl_take isl_basic_set *dom,
+					    __isl_give isl_set **empty,
+					    int max, __isl_take isl_mat *cst,
+					    __isl_take isl_space *map_space,
+					    __isl_take isl_space *set_space))
 {
 	int i, n, k;
 	int *list = NULL;
@@ -4338,12 +4571,11 @@ static __isl_give isl_map *basic_map_partial_lexopt_symm(
 	isl_ctx *ctx;
 	isl_vec *var = NULL;
 	isl_mat *cst = NULL;
-	isl_map *opt;
-	isl_set *min_expr;
-	isl_dim *map_dim, *set_dim;
+	isl_space *map_space, *set_space;
+	union isl_lex_res res;
 
-	map_dim = isl_basic_map_get_dim(bmap);
-	set_dim = empty ? isl_basic_set_get_dim(dom) : NULL;
+	map_space = isl_basic_map_get_space(bmap);
+	set_space = empty ? isl_basic_set_get_space(dom) : NULL;
 
 	n_in = isl_basic_map_dim(bmap, isl_dim_param) +
 	       isl_basic_map_dim(bmap, isl_dim_in);
@@ -4388,7 +4620,7 @@ static __isl_give isl_map *basic_map_partial_lexopt_symm(
 	bmap = isl_basic_map_finalize(bmap);
 
 	n_div = isl_basic_set_dim(dom, isl_dim_div);
-	dom = isl_basic_set_add(dom, isl_dim_set, 1);
+	dom = isl_basic_set_add_dims(dom, isl_dim_set, 1);
 	dom = isl_basic_set_extend_constraints(dom, 0, n);
 	for (i = 0; i < n; ++i) {
 		k = isl_basic_set_alloc_inequality(dom);
@@ -4399,32 +4631,28 @@ static __isl_give isl_map *basic_map_partial_lexopt_symm(
 		isl_seq_clr(dom->ineq[k] + 1 + n_in + 1, n_div);
 	}
 
-	min_expr = set_minimum(isl_basic_set_get_dim(dom), isl_mat_copy(cst));
-
 	isl_vec_free(var);
 	free(list);
 
-	opt = basic_map_partial_lexopt(bmap, dom, empty, max);
-
-	if (empty) {
-		*empty = split(*empty,
-			       isl_set_copy(min_expr), isl_mat_copy(cst));
-		*empty = isl_set_reset_dim(*empty, set_dim);
-	}
-
-	opt = split_domain(opt, min_expr, cst);
-	opt = isl_map_reset_dim(opt, map_dim);
-
-	return opt;
+	return core(bmap, dom, empty, max, cst, map_space, set_space);
 error:
-	isl_dim_free(map_dim);
-	isl_dim_free(set_dim);
+	isl_space_free(map_space);
+	isl_space_free(set_space);
 	isl_mat_free(cst);
 	isl_vec_free(var);
 	free(list);
 	isl_basic_set_free(dom);
 	isl_basic_map_free(bmap);
-	return NULL;
+	res.p = NULL;
+	return res;
+}
+
+static __isl_give isl_map *basic_map_partial_lexopt_symm_map(
+	__isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
+	__isl_give isl_set **empty, int max, int first, int second)
+{
+	return basic_map_partial_lexopt_symm(bmap, dom, empty, max,
+		    first, second, &basic_map_partial_lexopt_symm_map_core).map;
 }
 
 /* Recursive part of isl_tab_basic_map_partial_lexopt, after detecting
@@ -4451,10 +4679,10 @@ static __isl_give isl_map *basic_map_partial_lexopt(
 	if (par < 0)
 		goto error;
 	if (!par)
-		return basic_map_partial_lexopt_base(bmap, dom, empty, max);
+		return basic_map_partial_lexopt_base_map(bmap, dom, empty, max);
 	
-	return basic_map_partial_lexopt_symm(bmap, dom, empty, max,
-					     first, second);
+	return basic_map_partial_lexopt_symm_map(bmap, dom, empty, max,
+						 first, second);
 error:
 	isl_basic_set_free(dom);
 	isl_basic_map_free(bmap);
@@ -4553,9 +4781,10 @@ static void sol_for_add(struct isl_sol_for *sol,
 	for (i = 1; i < M->n_row; ++i) {
 		aff = isl_aff_alloc(isl_local_space_copy(ls));
 		if (aff) {
-			isl_int_set_si(aff->v->el[0], 1);
+			isl_int_set(aff->v->el[0], M->row[0][0]);
 			isl_seq_cpy(aff->v->el + 1, M->row[i], M->n_col);
 		}
+		aff = isl_aff_normalize(aff);
 		list = isl_aff_list_add(list, aff);
 	}
 	isl_local_space_free(ls);
@@ -4586,14 +4815,14 @@ static struct isl_sol_for *sol_for_init(struct isl_basic_map *bmap, int max,
 	void *user)
 {
 	struct isl_sol_for *sol_for = NULL;
-	struct isl_dim *dom_dim;
+	isl_space *dom_dim;
 	struct isl_basic_set *dom = NULL;
 
 	sol_for = isl_calloc_type(bmap->ctx, struct isl_sol_for);
 	if (!sol_for)
 		goto error;
 
-	dom_dim = isl_dim_domain(isl_dim_copy(bmap->dim));
+	dom_dim = isl_space_domain(isl_space_copy(bmap->dim));
 	dom = isl_basic_set_universe(dom_dim);
 
 	sol_for->sol.rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
@@ -4633,11 +4862,13 @@ int isl_basic_map_foreach_lexopt(__isl_keep isl_basic_map *bmap, int max,
 	struct isl_sol_for *sol_for = NULL;
 
 	bmap = isl_basic_map_copy(bmap);
+	bmap = isl_basic_map_detect_equalities(bmap);
 	if (!bmap)
 		return -1;
 
-	bmap = isl_basic_map_detect_equalities(bmap);
 	sol_for = sol_for_init(bmap, max, fn, user);
+	if (!sol_for)
+		goto error;
 
 	if (isl_basic_map_plain_is_empty(bmap))
 		/* nothing */;
@@ -4661,20 +4892,12 @@ error:
 	return -1;
 }
 
-int isl_basic_map_foreach_lexmin(__isl_keep isl_basic_map *bmap,
+int isl_basic_set_foreach_lexopt(__isl_keep isl_basic_set *bset, int max,
 	int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_aff_list *list,
 		  void *user),
 	void *user)
 {
-	return isl_basic_map_foreach_lexopt(bmap, 0, fn, user);
-}
-
-int isl_basic_map_foreach_lexmax(__isl_keep isl_basic_map *bmap,
-	int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_aff_list *list,
-		  void *user),
-	void *user)
-{
-	return isl_basic_map_foreach_lexopt(bmap, 1, fn, user);
+	return isl_basic_map_foreach_lexopt(bset, max, fn, user);
 }
 
 /* Check if the given sequence of len variables starting at pos
@@ -4844,14 +5067,20 @@ __isl_give isl_vec *isl_tab_basic_set_non_trivial_lexmin(
 {
 	int i, j;
 	int r;
-	isl_ctx *ctx = isl_basic_set_get_ctx(bset);
+	isl_ctx *ctx;
 	isl_vec *v = NULL;
-	isl_vec *sol = isl_vec_alloc(ctx, 0);
+	isl_vec *sol = NULL;
 	struct isl_tab *tab;
 	struct isl_trivial *triv = NULL;
 	int level, init;
 
-	tab = tab_for_lexmin(isl_basic_map_from_range(bset), NULL, 0, 0);
+	if (!bset)
+		return NULL;
+
+	ctx = isl_basic_set_get_ctx(bset);
+	sol = isl_vec_alloc(ctx, 0);
+
+	tab = tab_for_lexmin(bset, NULL, 0, 0);
 	if (!tab)
 		goto error;
 	tab->conflict = conflict;
@@ -4869,7 +5098,7 @@ __isl_give isl_vec *isl_tab_basic_set_non_trivial_lexmin(
 		int side, base;
 
 		if (init) {
-			tab = cut_to_integer_lexmin(tab);
+			tab = cut_to_integer_lexmin(tab, CUT_ONE);
 			if (!tab)
 				goto error;
 			if (tab->empty)
@@ -4959,14 +5188,17 @@ error:
  * assuming that all variables are non-negative.
  * If "bset" is empty, then return a zero-length vector.
  */
- __isl_give isl_vec *isl_tab_basic_set_non_neg_lexmin(
+__isl_give isl_vec *isl_tab_basic_set_non_neg_lexmin(
 	__isl_take isl_basic_set *bset)
 {
 	struct isl_tab *tab;
 	isl_ctx *ctx = isl_basic_set_get_ctx(bset);
 	isl_vec *sol;
 
-	tab = tab_for_lexmin(isl_basic_map_from_range(bset), NULL, 0, 0);
+	if (!bset)
+		return NULL;
+
+	tab = tab_for_lexmin(bset, NULL, 0, 0);
 	if (!tab)
 		goto error;
 	if (tab->empty)
@@ -4981,3 +5213,448 @@ error:
 	isl_basic_set_free(bset);
 	return NULL;
 }
+
+struct isl_sol_pma {
+	struct isl_sol	sol;
+	isl_pw_multi_aff *pma;
+	isl_set *empty;
+};
+
+static void sol_pma_free(struct isl_sol_pma *sol_pma)
+{
+	if (!sol_pma)
+		return;
+	if (sol_pma->sol.context)
+		sol_pma->sol.context->op->free(sol_pma->sol.context);
+	isl_pw_multi_aff_free(sol_pma->pma);
+	isl_set_free(sol_pma->empty);
+	free(sol_pma);
+}
+
+/* This function is called for parts of the context where there is
+ * no solution, with "bset" corresponding to the context tableau.
+ * Simply add the basic set to the set "empty".
+ */
+static void sol_pma_add_empty(struct isl_sol_pma *sol,
+	__isl_take isl_basic_set *bset)
+{
+	if (!bset)
+		goto error;
+	isl_assert(bset->ctx, sol->empty, goto error);
+
+	sol->empty = isl_set_grow(sol->empty, 1);
+	bset = isl_basic_set_simplify(bset);
+	bset = isl_basic_set_finalize(bset);
+	sol->empty = isl_set_add_basic_set(sol->empty, bset);
+	if (!sol->empty)
+		sol->sol.error = 1;
+	return;
+error:
+	isl_basic_set_free(bset);
+	sol->sol.error = 1;
+}
+
+/* Given a basic map "dom" that represents the context and an affine
+ * matrix "M" that maps the dimensions of the context to the
+ * output variables, construct an isl_pw_multi_aff with a single
+ * cell corresponding to "dom" and affine expressions copied from "M".
+ */
+static void sol_pma_add(struct isl_sol_pma *sol,
+	__isl_take isl_basic_set *dom, __isl_take isl_mat *M)
+{
+	int i;
+	isl_local_space *ls;
+	isl_aff *aff;
+	isl_multi_aff *maff;
+	isl_pw_multi_aff *pma;
+
+	maff = isl_multi_aff_alloc(isl_pw_multi_aff_get_space(sol->pma));
+	ls = isl_basic_set_get_local_space(dom);
+	for (i = 1; i < M->n_row; ++i) {
+		aff = isl_aff_alloc(isl_local_space_copy(ls));
+		if (aff) {
+			isl_int_set(aff->v->el[0], M->row[0][0]);
+			isl_seq_cpy(aff->v->el + 1, M->row[i], M->n_col);
+		}
+		aff = isl_aff_normalize(aff);
+		maff = isl_multi_aff_set_aff(maff, i - 1, aff);
+	}
+	isl_local_space_free(ls);
+	isl_mat_free(M);
+	dom = isl_basic_set_simplify(dom);
+	dom = isl_basic_set_finalize(dom);
+	pma = isl_pw_multi_aff_alloc(isl_set_from_basic_set(dom), maff);
+	sol->pma = isl_pw_multi_aff_add_disjoint(sol->pma, pma);
+	if (!sol->pma)
+		sol->sol.error = 1;
+}
+
+static void sol_pma_free_wrap(struct isl_sol *sol)
+{
+	sol_pma_free((struct isl_sol_pma *)sol);
+}
+
+static void sol_pma_add_empty_wrap(struct isl_sol *sol,
+	__isl_take isl_basic_set *bset)
+{
+	sol_pma_add_empty((struct isl_sol_pma *)sol, bset);
+}
+
+static void sol_pma_add_wrap(struct isl_sol *sol,
+	__isl_take isl_basic_set *dom, __isl_take isl_mat *M)
+{
+	sol_pma_add((struct isl_sol_pma *)sol, dom, M);
+}
+
+/* Construct an isl_sol_pma structure for accumulating the solution.
+ * If track_empty is set, then we also keep track of the parts
+ * of the context where there is no solution.
+ * If max is set, then we are solving a maximization, rather than
+ * a minimization problem, which means that the variables in the
+ * tableau have value "M - x" rather than "M + x".
+ */
+static struct isl_sol *sol_pma_init(__isl_keep isl_basic_map *bmap,
+	__isl_take isl_basic_set *dom, int track_empty, int max)
+{
+	struct isl_sol_pma *sol_pma = NULL;
+
+	if (!bmap)
+		goto error;
+
+	sol_pma = isl_calloc_type(bmap->ctx, struct isl_sol_pma);
+	if (!sol_pma)
+		goto error;
+
+	sol_pma->sol.rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
+	sol_pma->sol.dec_level.callback.run = &sol_dec_level_wrap;
+	sol_pma->sol.dec_level.sol = &sol_pma->sol;
+	sol_pma->sol.max = max;
+	sol_pma->sol.n_out = isl_basic_map_dim(bmap, isl_dim_out);
+	sol_pma->sol.add = &sol_pma_add_wrap;
+	sol_pma->sol.add_empty = track_empty ? &sol_pma_add_empty_wrap : NULL;
+	sol_pma->sol.free = &sol_pma_free_wrap;
+	sol_pma->pma = isl_pw_multi_aff_empty(isl_basic_map_get_space(bmap));
+	if (!sol_pma->pma)
+		goto error;
+
+	sol_pma->sol.context = isl_context_alloc(dom);
+	if (!sol_pma->sol.context)
+		goto error;
+
+	if (track_empty) {
+		sol_pma->empty = isl_set_alloc_space(isl_basic_set_get_space(dom),
+							1, ISL_SET_DISJOINT);
+		if (!sol_pma->empty)
+			goto error;
+	}
+
+	isl_basic_set_free(dom);
+	return &sol_pma->sol;
+error:
+	isl_basic_set_free(dom);
+	sol_pma_free(sol_pma);
+	return NULL;
+}
+
+/* Base case of isl_tab_basic_map_partial_lexopt, after removing
+ * some obvious symmetries.
+ *
+ * We call basic_map_partial_lexopt_base and extract the results.
+ */
+static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_base_pma(
+	__isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
+	__isl_give isl_set **empty, int max)
+{
+	isl_pw_multi_aff *result = NULL;
+	struct isl_sol *sol;
+	struct isl_sol_pma *sol_pma;
+
+	sol = basic_map_partial_lexopt_base(bmap, dom, empty, max,
+					    &sol_pma_init);
+	if (!sol)
+		return NULL;
+	sol_pma = (struct isl_sol_pma *) sol;
+
+	result = isl_pw_multi_aff_copy(sol_pma->pma);
+	if (empty)
+		*empty = isl_set_copy(sol_pma->empty);
+	sol_free(&sol_pma->sol);
+	return result;
+}
+
+/* Given that the last input variable of "maff" represents the minimum
+ * of some bounds, check whether we need to plug in the expression
+ * of the minimum.
+ *
+ * In particular, check if the last input variable appears in any
+ * of the expressions in "maff".
+ */
+static int need_substitution(__isl_keep isl_multi_aff *maff)
+{
+	int i;
+	unsigned pos;
+
+	pos = isl_multi_aff_dim(maff, isl_dim_in) - 1;
+
+	for (i = 0; i < maff->n; ++i)
+		if (isl_aff_involves_dims(maff->p[i], isl_dim_in, pos, 1))
+			return 1;
+
+	return 0;
+}
+
+/* Given a set of upper bounds on the last "input" variable m,
+ * construct a piecewise affine expression that selects
+ * the minimal upper bound to m, i.e.,
+ * divide the space into cells where one
+ * of the upper bounds is smaller than all the others and select
+ * this upper bound on that cell.
+ *
+ * In particular, if there are n bounds b_i, then the result
+ * consists of n cell, each one of the form
+ *
+ *	b_i <= b_j	for j > i
+ *	b_i <  b_j	for j < i
+ *
+ * The affine expression on this cell is
+ *
+ *	b_i
+ */
+static __isl_give isl_pw_aff *set_minimum_pa(__isl_take isl_space *space,
+	__isl_take isl_mat *var)
+{
+	int i;
+	isl_aff *aff = NULL;
+	isl_basic_set *bset = NULL;
+	isl_ctx *ctx;
+	isl_pw_aff *paff = NULL;
+	isl_space *pw_space;
+	isl_local_space *ls = NULL;
+
+	if (!space || !var)
+		goto error;
+
+	ctx = isl_space_get_ctx(space);
+	ls = isl_local_space_from_space(isl_space_copy(space));
+	pw_space = isl_space_copy(space);
+	pw_space = isl_space_from_domain(pw_space);
+	pw_space = isl_space_add_dims(pw_space, isl_dim_out, 1);
+	paff = isl_pw_aff_alloc_size(pw_space, var->n_row);
+
+	for (i = 0; i < var->n_row; ++i) {
+		isl_pw_aff *paff_i;
+
+		aff = isl_aff_alloc(isl_local_space_copy(ls));
+		bset = isl_basic_set_alloc_space(isl_space_copy(space), 0,
+					       0, var->n_row - 1);
+		if (!aff || !bset)
+			goto error;
+		isl_int_set_si(aff->v->el[0], 1);
+		isl_seq_cpy(aff->v->el + 1, var->row[i], var->n_col);
+		isl_int_set_si(aff->v->el[1 + var->n_col], 0);
+		bset = select_minimum(bset, var, i);
+		paff_i = isl_pw_aff_alloc(isl_set_from_basic_set(bset), aff);
+		paff = isl_pw_aff_add_disjoint(paff, paff_i);
+	}
+
+	isl_local_space_free(ls);
+	isl_space_free(space);
+	isl_mat_free(var);
+	return paff;
+error:
+	isl_aff_free(aff);
+	isl_basic_set_free(bset);
+	isl_pw_aff_free(paff);
+	isl_local_space_free(ls);
+	isl_space_free(space);
+	isl_mat_free(var);
+	return NULL;
+}
+
+/* Given a piecewise multi-affine expression of which the last input variable
+ * is the minimum of the bounds in "cst", plug in the value of the minimum.
+ * This minimum expression is given in "min_expr_pa".
+ * The set "min_expr" contains the same information, but in the form of a set.
+ * The variable is subsequently projected out.
+ *
+ * The implementation is similar to those of "split" and "split_domain".
+ * If the variable appears in a given expression, then minimum expression
+ * is plugged in.  Otherwise, if the variable appears in the constraints
+ * and a split is required, then the domain is split.  Otherwise, no split
+ * is performed.
+ */
+static __isl_give isl_pw_multi_aff *split_domain_pma(
+	__isl_take isl_pw_multi_aff *opt, __isl_take isl_pw_aff *min_expr_pa,
+	__isl_take isl_set *min_expr, __isl_take isl_mat *cst)
+{
+	int n_in;
+	int i;
+	isl_space *space;
+	isl_pw_multi_aff *res;
+
+	if (!opt || !min_expr || !cst)
+		goto error;
+
+	n_in = isl_pw_multi_aff_dim(opt, isl_dim_in);
+	space = isl_pw_multi_aff_get_space(opt);
+	space = isl_space_drop_dims(space, isl_dim_in, n_in - 1, 1);
+	res = isl_pw_multi_aff_empty(space);
+
+	for (i = 0; i < opt->n; ++i) {
+		isl_pw_multi_aff *pma;
+
+		pma = isl_pw_multi_aff_alloc(isl_set_copy(opt->p[i].set),
+					 isl_multi_aff_copy(opt->p[i].maff));
+		if (need_substitution(opt->p[i].maff))
+			pma = isl_pw_multi_aff_substitute(pma,
+					isl_dim_in, n_in - 1, min_expr_pa);
+		else if (need_split_set(opt->p[i].set, cst))
+			pma = isl_pw_multi_aff_intersect_domain(pma,
+						       isl_set_copy(min_expr));
+		pma = isl_pw_multi_aff_project_out(pma,
+						    isl_dim_in, n_in - 1, 1);
+
+		res = isl_pw_multi_aff_add_disjoint(res, pma);
+	}
+
+	isl_pw_multi_aff_free(opt);
+	isl_pw_aff_free(min_expr_pa);
+	isl_set_free(min_expr);
+	isl_mat_free(cst);
+	return res;
+error:
+	isl_pw_multi_aff_free(opt);
+	isl_pw_aff_free(min_expr_pa);
+	isl_set_free(min_expr);
+	isl_mat_free(cst);
+	return NULL;
+}
+
+static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_pma(
+	__isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
+	__isl_give isl_set **empty, int max);
+
+/* This function is called from basic_map_partial_lexopt_symm.
+ * The last variable of "bmap" and "dom" corresponds to the minimum
+ * of the bounds in "cst".  "map_space" is the space of the original
+ * input relation (of basic_map_partial_lexopt_symm) and "set_space"
+ * is the space of the original domain.
+ *
+ * We recursively call basic_map_partial_lexopt and then plug in
+ * the definition of the minimum in the result.
+ */
+static __isl_give union isl_lex_res basic_map_partial_lexopt_symm_pma_core(
+	__isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
+	__isl_give isl_set **empty, int max, __isl_take isl_mat *cst,
+	__isl_take isl_space *map_space, __isl_take isl_space *set_space)
+{
+	isl_pw_multi_aff *opt;
+	isl_pw_aff *min_expr_pa;
+	isl_set *min_expr;
+	union isl_lex_res res;
+
+	min_expr = set_minimum(isl_basic_set_get_space(dom), isl_mat_copy(cst));
+	min_expr_pa = set_minimum_pa(isl_basic_set_get_space(dom),
+					isl_mat_copy(cst));
+
+	opt = basic_map_partial_lexopt_pma(bmap, dom, empty, max);
+
+	if (empty) {
+		*empty = split(*empty,
+			       isl_set_copy(min_expr), isl_mat_copy(cst));
+		*empty = isl_set_reset_space(*empty, set_space);
+	}
+
+	opt = split_domain_pma(opt, min_expr_pa, min_expr, cst);
+	opt = isl_pw_multi_aff_reset_space(opt, map_space);
+
+	res.pma = opt;
+	return res;
+}
+
+static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_symm_pma(
+	__isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
+	__isl_give isl_set **empty, int max, int first, int second)
+{
+	return basic_map_partial_lexopt_symm(bmap, dom, empty, max,
+		    first, second, &basic_map_partial_lexopt_symm_pma_core).pma;
+}
+
+/* Recursive part of isl_basic_map_partial_lexopt_pw_multi_aff, after detecting
+ * equalities and removing redundant constraints.
+ *
+ * We first check if there are any parallel constraints (left).
+ * If not, we are in the base case.
+ * If there are parallel constraints, we replace them by a single
+ * constraint in basic_map_partial_lexopt_symm_pma and then call
+ * this function recursively to look for more parallel constraints.
+ */
+static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_pma(
+	__isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
+	__isl_give isl_set **empty, int max)
+{
+	int par = 0;
+	int first, second;
+
+	if (!bmap)
+		goto error;
+
+	if (bmap->ctx->opt->pip_symmetry)
+		par = parallel_constraints(bmap, &first, &second);
+	if (par < 0)
+		goto error;
+	if (!par)
+		return basic_map_partial_lexopt_base_pma(bmap, dom, empty, max);
+	
+	return basic_map_partial_lexopt_symm_pma(bmap, dom, empty, max,
+						 first, second);
+error:
+	isl_basic_set_free(dom);
+	isl_basic_map_free(bmap);
+	return NULL;
+}
+
+/* Compute the lexicographic minimum (or maximum if "max" is set)
+ * of "bmap" over the domain "dom" and return the result as a piecewise
+ * multi-affine expression.
+ * If "empty" is not NULL, then *empty is assigned a set that
+ * contains those parts of the domain where there is no solution.
+ * If "bmap" is marked as rational (ISL_BASIC_MAP_RATIONAL),
+ * then we compute the rational optimum.  Otherwise, we compute
+ * the integral optimum.
+ *
+ * We perform some preprocessing.  As the PILP solver does not
+ * handle implicit equalities very well, we first make sure all
+ * the equalities are explicitly available.
+ *
+ * We also add context constraints to the basic map and remove
+ * redundant constraints.  This is only needed because of the
+ * way we handle simple symmetries.  In particular, we currently look
+ * for symmetries on the constraints, before we set up the main tableau.
+ * It is then no good to look for symmetries on possibly redundant constraints.
+ */
+__isl_give isl_pw_multi_aff *isl_basic_map_partial_lexopt_pw_multi_aff(
+	__isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
+	__isl_give isl_set **empty, int max)
+{
+	if (empty)
+		*empty = NULL;
+	if (!bmap || !dom)
+		goto error;
+
+	isl_assert(bmap->ctx,
+	    isl_basic_map_compatible_domain(bmap, dom), goto error);
+
+	if (isl_basic_set_dim(dom, isl_dim_all) == 0)
+		return basic_map_partial_lexopt_pma(bmap, dom, empty, max);
+
+	bmap = isl_basic_map_intersect_domain(bmap, isl_basic_set_copy(dom));
+	bmap = isl_basic_map_detect_equalities(bmap);
+	bmap = isl_basic_map_remove_redundancies(bmap);
+
+	return basic_map_partial_lexopt_pma(bmap, dom, empty, max);
+error:
+	isl_basic_set_free(dom);
+	isl_basic_map_free(bmap);
+	return NULL;
+}