+ int i;
+ isl_ctx *ctx;
+ isl_local_space *ls;
+ isl_aff *aff;
+ isl_aff_list *list;
+
+ if (sol->sol.error || !dom || !M)
+ goto error;
+
+ ctx = isl_basic_set_get_ctx(dom);
+ ls = isl_basic_set_get_local_space(dom);
+ list = isl_aff_list_alloc(ctx, M->n_row - 1);
+ 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);
+ }
+ list = isl_aff_list_add(list, aff);
+ }
+ isl_local_space_free(ls);
+
+ dom = isl_basic_set_finalize(dom);
+
+ if (sol->fn(isl_basic_set_copy(dom), list, sol->user) < 0)
+ goto error;
+
+ isl_basic_set_free(dom);
+ isl_mat_free(M);
+ return;
+error:
+ isl_basic_set_free(dom);
+ isl_mat_free(M);
+ sol->sol.error = 1;
+}
+
+static void sol_for_add_wrap(struct isl_sol *sol,
+ struct isl_basic_set *dom, struct isl_mat *M)
+{
+ sol_for_add((struct isl_sol_for *)sol, dom, M);
+}
+
+static struct isl_sol_for *sol_for_init(struct isl_basic_map *bmap, int max,
+ int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_aff_list *list,
+ void *user),
+ void *user)
+{
+ struct isl_sol_for *sol_for = NULL;
+ 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_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);
+ sol_for->sol.dec_level.callback.run = &sol_dec_level_wrap;
+ sol_for->sol.dec_level.sol = &sol_for->sol;
+ sol_for->fn = fn;
+ sol_for->user = user;
+ sol_for->sol.max = max;
+ sol_for->sol.n_out = isl_basic_map_dim(bmap, isl_dim_out);
+ sol_for->sol.add = &sol_for_add_wrap;
+ sol_for->sol.add_empty = NULL;
+ sol_for->sol.free = &sol_for_free_wrap;
+
+ sol_for->sol.context = isl_context_alloc(dom);
+ if (!sol_for->sol.context)
+ goto error;
+
+ isl_basic_set_free(dom);
+ return sol_for;
+error:
+ isl_basic_set_free(dom);
+ sol_for_free(sol_for);
+ return NULL;
+}
+
+static void sol_for_find_solutions(struct isl_sol_for *sol_for,
+ struct isl_tab *tab)
+{
+ find_solutions_main(&sol_for->sol, tab);
+}
+
+int isl_basic_map_foreach_lexopt(__isl_keep isl_basic_map *bmap, int max,
+ int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_aff_list *list,
+ void *user),
+ void *user)
+{
+ struct isl_sol_for *sol_for = NULL;
+
+ bmap = isl_basic_map_copy(bmap);
+ if (!bmap)
+ return -1;
+
+ bmap = isl_basic_map_detect_equalities(bmap);
+ sol_for = sol_for_init(bmap, max, fn, user);
+
+ if (isl_basic_map_plain_is_empty(bmap))
+ /* nothing */;
+ else {
+ struct isl_tab *tab;
+ struct isl_context *context = sol_for->sol.context;
+ tab = tab_for_lexmin(bmap,
+ context->op->peek_basic_set(context), 1, max);
+ tab = context->op->detect_nonnegative_parameters(context, tab);
+ sol_for_find_solutions(sol_for, tab);
+ if (sol_for->sol.error)
+ goto error;
+ }
+
+ sol_free(&sol_for->sol);
+ isl_basic_map_free(bmap);
+ return 0;
+error:
+ sol_free(&sol_for->sol);
+ isl_basic_map_free(bmap);
+ return -1;
+}
+
+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(bset, max, fn, user);
+}
+
+/* Check if the given sequence of len variables starting at pos
+ * represents a trivial (i.e., zero) solution.
+ * The variables are assumed to be non-negative and to come in pairs,
+ * with each pair representing a variable of unrestricted sign.
+ * The solution is trivial if each such pair in the sequence consists
+ * of two identical values, meaning that the variable being represented
+ * has value zero.
+ */
+static int region_is_trivial(struct isl_tab *tab, int pos, int len)
+{
+ int i;
+
+ if (len == 0)
+ return 0;
+
+ for (i = 0; i < len; i += 2) {
+ int neg_row;
+ int pos_row;
+
+ neg_row = tab->var[pos + i].is_row ?
+ tab->var[pos + i].index : -1;
+ pos_row = tab->var[pos + i + 1].is_row ?
+ tab->var[pos + i + 1].index : -1;
+
+ if ((neg_row < 0 ||
+ isl_int_is_zero(tab->mat->row[neg_row][1])) &&
+ (pos_row < 0 ||
+ isl_int_is_zero(tab->mat->row[pos_row][1])))
+ continue;
+
+ if (neg_row < 0 || pos_row < 0)
+ return 0;
+ if (isl_int_ne(tab->mat->row[neg_row][1],
+ tab->mat->row[pos_row][1]))
+ return 0;
+ }
+
+ return 1;
+}
+
+/* Return the index of the first trivial region or -1 if all regions
+ * are non-trivial.
+ */
+static int first_trivial_region(struct isl_tab *tab,
+ int n_region, struct isl_region *region)
+{
+ int i;
+
+ for (i = 0; i < n_region; ++i) {
+ if (region_is_trivial(tab, region[i].pos, region[i].len))
+ return i;
+ }
+
+ return -1;
+}
+
+/* Check if the solution is optimal, i.e., whether the first
+ * n_op entries are zero.
+ */
+static int is_optimal(__isl_keep isl_vec *sol, int n_op)
+{
+ int i;
+
+ for (i = 0; i < n_op; ++i)
+ if (!isl_int_is_zero(sol->el[1 + i]))
+ return 0;
+ return 1;
+}
+
+/* Add constraints to "tab" that ensure that any solution is significantly
+ * better that that represented by "sol". That is, find the first
+ * relevant (within first n_op) non-zero coefficient and force it (along
+ * with all previous coefficients) to be zero.
+ * If the solution is already optimal (all relevant coefficients are zero),
+ * then just mark the table as empty.
+ */
+static int force_better_solution(struct isl_tab *tab,
+ __isl_keep isl_vec *sol, int n_op)
+{
+ int i;
+ isl_ctx *ctx;
+ isl_vec *v = NULL;
+
+ if (!sol)
+ return -1;
+
+ for (i = 0; i < n_op; ++i)
+ if (!isl_int_is_zero(sol->el[1 + i]))
+ break;
+
+ if (i == n_op) {
+ if (isl_tab_mark_empty(tab) < 0)
+ return -1;
+ return 0;
+ }
+
+ ctx = isl_vec_get_ctx(sol);
+ v = isl_vec_alloc(ctx, 1 + tab->n_var);
+ if (!v)
+ return -1;
+
+ for (; i >= 0; --i) {
+ v = isl_vec_clr(v);
+ isl_int_set_si(v->el[1 + i], -1);
+ if (add_lexmin_eq(tab, v->el) < 0)
+ goto error;
+ }
+
+ isl_vec_free(v);
+ return 0;
+error:
+ isl_vec_free(v);
+ return -1;
+}
+
+struct isl_trivial {
+ int update;
+ int region;
+ int side;
+ struct isl_tab_undo *snap;
+};
+
+/* Return the lexicographically smallest non-trivial solution of the
+ * given ILP problem.
+ *
+ * All variables are assumed to be non-negative.
+ *
+ * n_op is the number of initial coordinates to optimize.
+ * That is, once a solution has been found, we will only continue looking
+ * for solution that result in significantly better values for those
+ * initial coordinates. That is, we only continue looking for solutions
+ * that increase the number of initial zeros in this sequence.
+ *
+ * A solution is non-trivial, if it is non-trivial on each of the
+ * specified regions. Each region represents a sequence of pairs
+ * of variables. A solution is non-trivial on such a region if
+ * at least one of these pairs consists of different values, i.e.,
+ * such that the non-negative variable represented by the pair is non-zero.
+ *
+ * Whenever a conflict is encountered, all constraints involved are
+ * reported to the caller through a call to "conflict".
+ *
+ * We perform a simple branch-and-bound backtracking search.
+ * Each level in the search represents initially trivial region that is forced
+ * to be non-trivial.
+ * At each level we consider n cases, where n is the length of the region.
+ * In terms of the n/2 variables of unrestricted signs being encoded by
+ * the region, we consider the cases
+ * x_0 >= 1
+ * x_0 <= -1
+ * x_0 = 0 and x_1 >= 1
+ * x_0 = 0 and x_1 <= -1
+ * x_0 = 0 and x_1 = 0 and x_2 >= 1
+ * x_0 = 0 and x_1 = 0 and x_2 <= -1
+ * ...
+ * The cases are considered in this order, assuming that each pair
+ * x_i_a x_i_b represents the value x_i_b - x_i_a.
+ * That is, x_0 >= 1 is enforced by adding the constraint
+ * x_0_b - x_0_a >= 1
+ */
+__isl_give isl_vec *isl_tab_basic_set_non_trivial_lexmin(
+ __isl_take isl_basic_set *bset, int n_op, int n_region,
+ struct isl_region *region,
+ int (*conflict)(int con, void *user), void *user)
+{
+ int i, j;
+ int r;
+ isl_ctx *ctx = isl_basic_set_get_ctx(bset);
+ isl_vec *v = NULL;
+ isl_vec *sol = isl_vec_alloc(ctx, 0);
+ struct isl_tab *tab;
+ struct isl_trivial *triv = NULL;
+ int level, init;
+
+ tab = tab_for_lexmin(bset, NULL, 0, 0);
+ if (!tab)
+ goto error;
+ tab->conflict = conflict;
+ tab->conflict_user = user;
+
+ v = isl_vec_alloc(ctx, 1 + tab->n_var);
+ triv = isl_calloc_array(ctx, struct isl_trivial, n_region);
+ if (!v || !triv)
+ goto error;
+
+ level = 0;
+ init = 1;
+
+ while (level >= 0) {
+ int side, base;
+
+ if (init) {
+ tab = cut_to_integer_lexmin(tab);
+ if (!tab)
+ goto error;
+ if (tab->empty)
+ goto backtrack;
+ r = first_trivial_region(tab, n_region, region);
+ if (r < 0) {
+ for (i = 0; i < level; ++i)
+ triv[i].update = 1;
+ isl_vec_free(sol);
+ sol = isl_tab_get_sample_value(tab);
+ if (!sol)
+ goto error;
+ if (is_optimal(sol, n_op))
+ break;
+ goto backtrack;
+ }
+ if (level >= n_region)
+ isl_die(ctx, isl_error_internal,
+ "nesting level too deep", goto error);
+ if (isl_tab_extend_cons(tab,
+ 2 * region[r].len + 2 * n_op) < 0)
+ goto error;
+ triv[level].region = r;
+ triv[level].side = 0;
+ }
+
+ r = triv[level].region;
+ side = triv[level].side;
+ base = 2 * (side/2);
+
+ if (side >= region[r].len) {
+backtrack:
+ level--;
+ init = 0;
+ if (level >= 0)
+ if (isl_tab_rollback(tab, triv[level].snap) < 0)
+ goto error;
+ continue;
+ }
+
+ if (triv[level].update) {
+ if (force_better_solution(tab, sol, n_op) < 0)
+ goto error;
+ triv[level].update = 0;
+ }
+
+ if (side == base && base >= 2) {
+ for (j = base - 2; j < base; ++j) {
+ v = isl_vec_clr(v);
+ isl_int_set_si(v->el[1 + region[r].pos + j], 1);
+ if (add_lexmin_eq(tab, v->el) < 0)
+ goto error;
+ }
+ }
+
+ triv[level].snap = isl_tab_snap(tab);
+ if (isl_tab_push_basis(tab) < 0)
+ goto error;
+
+ v = isl_vec_clr(v);
+ isl_int_set_si(v->el[0], -1);
+ isl_int_set_si(v->el[1 + region[r].pos + side], -1);
+ isl_int_set_si(v->el[1 + region[r].pos + (side ^ 1)], 1);
+ tab = add_lexmin_ineq(tab, v->el);
+
+ triv[level].side++;
+ level++;
+ init = 1;
+ }
+
+ free(triv);
+ isl_vec_free(v);
+ isl_tab_free(tab);
+ isl_basic_set_free(bset);
+
+ return sol;
+error:
+ free(triv);
+ isl_vec_free(v);
+ isl_tab_free(tab);
+ isl_basic_set_free(bset);
+ isl_vec_free(sol);
+ return NULL;
+}
+
+/* Return the lexicographically smallest rational point in "bset",
+ * 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_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(bset, NULL, 0, 0);
+ if (!tab)
+ goto error;
+ if (tab->empty)
+ sol = isl_vec_alloc(ctx, 0);
+ else
+ sol = isl_tab_get_sample_value(tab);
+ isl_tab_free(tab);
+ isl_basic_set_free(bset);
+ return sol;
+error:
+ isl_tab_free(tab);
+ 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;