X-Git-Url: http://review.tizen.org/git/?a=blobdiff_plain;f=isl_tab_pip.c;h=1f2ab1b410becb2eca67593aa0d1815f82a888b5;hb=78b3cfae0a9cad8fbca88a9b5ce79337c2ba01c3;hp=b616a335e402c33bea1e259ae2ee1e416d2d741a;hpb=3dd551186a6f18654214338a83d141c80885ddf6;p=platform%2Fupstream%2Fisl.git diff --git a/isl_tab_pip.c b/isl_tab_pip.c index b616a33..1f2ab1b 100644 --- a/isl_tab_pip.c +++ b/isl_tab_pip.c @@ -1,16 +1,20 @@ /* * Copyright 2008-2009 Katholieke Universiteit Leuven + * Copyright 2010 INRIA Saclay * * Use of this software is governed by the GNU LGPLv2.1 license * * Written by Sven Verdoolaege, K.U.Leuven, Departement * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium + * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite, + * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France */ #include "isl_map_private.h" -#include "isl_seq.h" +#include #include "isl_tab.h" #include "isl_sample.h" +#include /* * The implementation of parametric integer linear programming in this file @@ -35,10 +39,10 @@ * then the initial sample value may be chosen equal to zero. * However, we will not make this assumption. Instead, we apply * the "big parameter" trick. Any variable x is then not directly - * used in the tableau, but instead it its represented by another + * used in the tableau, but instead it is represented by another * variable x' = M + x, where M is an arbitrarily large (positive) * value. x' is therefore always non-negative, whatever the value of x. - * Taking as initial smaple value x' = 0 corresponds to x = -M, + * Taking as initial sample value x' = 0 corresponds to x = -M, * which is always smaller than any possible value of x. * * The big parameter trick is used in the main tableau and @@ -448,17 +452,17 @@ static void sol_add(struct isl_sol *sol, struct isl_tab *tab) isl_seq_clr(mat->row[1 + row], mat->n_col); if (!tab->var[i].is_row) { - /* no unbounded */ - isl_assert(mat->ctx, !tab->M, goto error2); + if (tab->M) + isl_die(mat->ctx, isl_error_invalid, + "unbounded optimum", goto error2); continue; } r = tab->var[i].index; - /* no unbounded */ - if (tab->M) - isl_assert(mat->ctx, isl_int_eq(tab->mat->row[r][2], - tab->mat->row[r][0]), - goto error2); + if (tab->M && + isl_int_ne(tab->mat->row[r][2], tab->mat->row[r][0])) + isl_die(mat->ctx, isl_error_invalid, + "unbounded optimum", goto error2); isl_int_gcd(m, mat->row[0][0], tab->mat->row[r][0]); isl_int_divexact(m, tab->mat->row[r][0], m); scale_rows(mat, m, 1 + row); @@ -494,7 +498,7 @@ error2: error: isl_basic_set_free(bset); isl_mat_free(mat); - sol_free(sol); + sol->error = 1; } struct isl_sol_map { @@ -505,6 +509,8 @@ struct isl_sol_map { static void sol_map_free(struct isl_sol_map *sol_map) { + if (!sol_map) + return; if (sol_map->sol.context) sol_map->sol.context->op->free(sol_map->sol.context); isl_map_free(sol_map->map); @@ -1070,7 +1076,7 @@ error: } /* Return the first known violated constraint, i.e., a non-negative - * contraint that currently has an either obviously negative value + * constraint that currently has an either obviously negative value * or a previously determined to be negative value. * * If any constraint has a negative coefficient for the big parameter, @@ -1108,7 +1114,7 @@ static int first_neg(struct isl_tab *tab) /* Resolve all known or obviously violated constraints through pivoting. * In particular, as long as we can find any violated constraint, we - * look for a pivoting column that would result in the lexicographicallly + * look for a pivoting column that would result in the lexicographically * smallest increment in the sample point. If there is no such column * then the tableau is infeasible. */ @@ -1220,8 +1226,6 @@ static struct isl_tab *add_lexmin_valid_eq(struct isl_tab *tab, isl_int *eq) if (isl_tab_kill_col(tab, i) < 0) goto error; tab->n_eq++; - - tab = restore_lexmin(tab); } return tab; @@ -1714,7 +1718,7 @@ static int tab_has_valid_sample(struct isl_tab *tab, isl_int *ineq, int eq) return i < tab->n_sample; } -/* Add a div specifed by "div" to the tableau "tab" and return +/* Add a div specified by "div" to the tableau "tab" and return * 1 if the div is obviously non-negative. */ static int context_tab_add_div(struct isl_tab *tab, struct isl_vec *div, @@ -1988,6 +1992,8 @@ static struct isl_tab *tab_for_lexmin(struct isl_basic_map *bmap, if (!tab || tab->empty) return tab; } + if (bmap->n_eq) + tab = restore_lexmin(tab); for (i = 0; i < bmap->n_ineq; ++i) { if (max) isl_seq_neg(bmap->ineq[i] + 1 + tab->n_param, @@ -2251,11 +2257,25 @@ static int context_lex_get_div(struct isl_context *context, struct isl_tab *tab, return get_div(tab, context, div); } +/* Add a div specified by "div" to the context tableau and return + * 1 if the div is obviously non-negative. + * context_tab_add_div will always return 1, because all variables + * in a isl_context_lex tableau are non-negative. + * However, if we are using a big parameter in the context, then this only + * reflects the non-negativity of the variable used to _encode_ the + * div, i.e., div' = M + div, so we can't draw any conclusions. + */ static int context_lex_add_div(struct isl_context *context, struct isl_vec *div) { struct isl_context_lex *clex = (struct isl_context_lex *)context; - return context_tab_add_div(clex->tab, div, + int nonneg; + nonneg = context_tab_add_div(clex->tab, div, context_lex_add_ineq_wrap, context); + if (nonneg < 0) + return -1; + if (clex->tab->M) + return 0; + return nonneg; } static int context_lex_detect_equalities(struct isl_context *context, @@ -2276,7 +2296,7 @@ static int context_lex_best_split(struct isl_context *context, return -1; r = best_split(tab, clex->tab); - if (isl_tab_rollback(clex->tab, snap) < 0) + if (r >= 0 && isl_tab_rollback(clex->tab, snap) < 0) return -1; return r; @@ -2384,6 +2404,9 @@ static struct isl_tab *context_lex_detect_nonnegative_parameters( struct isl_context_lex *clex = (struct isl_context_lex *)context; struct isl_tab_undo *snap; + if (!tab) + return NULL; + snap = isl_tab_snap(clex->tab); if (isl_tab_push_basis(clex->tab) < 0) goto error; @@ -2488,6 +2511,8 @@ static struct isl_tab *context_gbr_detect_nonnegative_parameters( struct isl_context *context, struct isl_tab *tab) { struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context; + if (!tab) + return NULL; return tab_detect_nonnegative_parameters(tab, cgbr->tab); } @@ -3032,7 +3057,7 @@ static int context_gbr_best_split(struct isl_context *context, snap = isl_tab_snap(cgbr->tab); r = best_split(tab, cgbr->tab); - if (isl_tab_rollback(cgbr->tab, snap) < 0) + if (r >= 0 && isl_tab_rollback(cgbr->tab, snap) < 0) return -1; return r; @@ -3210,7 +3235,10 @@ static struct isl_context *isl_context_alloc(struct isl_basic_set *dom) static struct isl_sol_map *sol_map_init(struct isl_basic_map *bmap, struct isl_basic_set *dom, int track_empty, int max) { - struct isl_sol_map *sol_map; + struct isl_sol_map *sol_map = NULL; + + if (!bmap) + goto error; sol_map = isl_calloc_type(bmap->ctx, struct isl_sol_map); if (!sol_map) @@ -3451,7 +3479,8 @@ static void find_in_pos(struct isl_sol *sol, struct isl_tab *tab, isl_int *ineq) find_solutions(sol, tab); - sol->context->op->restore(sol->context, saved); + if (!sol->error) + sol->context->op->restore(sol->context, saved); return; error: sol->error = 1; @@ -3466,7 +3495,7 @@ static void no_sol_in_strict(struct isl_sol *sol, int empty; void *saved; - if (!sol->context) + if (!sol->context || sol->error) goto error; saved = sol->context->op->save(sol->context); @@ -3563,7 +3592,7 @@ error: * coefficient are integral, then there is nothing that can be done * and the tableau has no integral solution. * If, on the other hand, one or more of the other columns have rational - * coeffcients, but the parameter coefficients are all integral, then + * coefficients, but the parameter coefficients are all integral, then * we can perform a regular (non-parametric) cut. * Finally, if there is any parameter coefficient that is non-integral, * then we need to involve the context tableau. There are two cases here. @@ -3643,7 +3672,8 @@ static void find_solutions(struct isl_sol *sol, struct isl_tab *tab) row = split; isl_seq_neg(ineq->el, ineq->el, ineq->size); isl_int_sub_ui(ineq->el[0], ineq->el[0], 1); - context->op->add_ineq(context, ineq->el, 0, 1); + if (!sol->error) + context->op->add_ineq(context, ineq->el, 0, 1); isl_vec_free(ineq); if (sol->error) goto error; @@ -3673,6 +3703,8 @@ static void find_solutions(struct isl_sol *sol, struct isl_tab *tab) if (d < 0) goto error; ineq = ineq_for_div(context->op->peek_basic_set(context), d); + if (!ineq) + goto error; sol_inc_level(sol); no_sol_in_strict(sol, tab, ineq); isl_seq_neg(ineq->el, ineq->el, ineq->size); @@ -3694,7 +3726,7 @@ done: return; error: isl_tab_free(tab); - sol_free(sol); + sol->error = 1; } /* Compute the lexicographic minimum of the set represented by the main @@ -3712,6 +3744,9 @@ static void find_solutions_main(struct isl_sol *sol, struct isl_tab *tab) { int row; + if (!tab) + goto error; + sol->level = 0; for (row = tab->n_redundant; row < tab->n_row; ++row) { @@ -3730,6 +3765,8 @@ static void find_solutions_main(struct isl_sol *sol, struct isl_tab *tab) + tab->n_param - (tab->n_var - tab->n_div); eq = isl_vec_alloc(tab->mat->ctx, 1+tab->n_param+tab->n_div); + if (!eq) + goto error; get_row_parameter_line(tab, row, eq->el); isl_int_neg(eq->el[1 + p], tab->mat->row[row][0]); eq = isl_vec_normalize(eq); @@ -3763,7 +3800,7 @@ static void find_solutions_main(struct isl_sol *sol, struct isl_tab *tab) return; error: isl_tab_free(tab); - sol_free(sol); + sol->error = 1; } static void sol_map_find_solutions(struct isl_sol_map *sol_map, @@ -3848,43 +3885,21 @@ error: return NULL; } -/* Compute the lexicographic minimum (or maximum if "max" is set) - * of "bmap" over the domain "dom" and return the result as a map. - * 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. +/* Base case of isl_tab_basic_map_partial_lexopt, after removing + * some obvious symmetries. * - * 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 make sure the divs in the domain are properly order, + * We make sure the divs in the domain are properly ordered, * because they will be added one by one in the given order * during the construction of the solution map. */ -struct isl_map *isl_tab_basic_map_partial_lexopt( - struct isl_basic_map *bmap, struct isl_basic_set *dom, - struct isl_set **empty, int max) +static __isl_give isl_map *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_map *result = NULL; struct isl_tab *tab; - struct isl_map *result = NULL; struct isl_sol_map *sol_map = NULL; struct isl_context *context; - struct isl_basic_map *eq; - - if (empty) - *empty = NULL; - if (!bmap || !dom) - goto error; - - isl_assert(bmap->ctx, - isl_basic_map_compatible_domain(bmap, dom), goto error); - - eq = isl_basic_map_copy(bmap); - eq = isl_basic_map_intersect_domain(eq, isl_basic_set_copy(dom)); - eq = isl_basic_map_affine_hull(eq); - bmap = isl_basic_map_intersect(bmap, eq); if (dom->n_div) { dom = isl_basic_set_order_divs(dom); @@ -3921,6 +3936,512 @@ error: return NULL; } +/* Structure used during detection of parallel constraints. + * n_in: number of "input" variables: isl_dim_param + isl_dim_in + * n_out: number of "output" variables: isl_dim_out + isl_dim_div + * val: the coefficients of the output variables + */ +struct isl_constraint_equal_info { + isl_basic_map *bmap; + unsigned n_in; + unsigned n_out; + isl_int *val; +}; + +/* Check whether the coefficients of the output variables + * of the constraint in "entry" are equal to info->val. + */ +static int constraint_equal(const void *entry, const void *val) +{ + isl_int **row = (isl_int **)entry; + const struct isl_constraint_equal_info *info = val; + + return isl_seq_eq((*row) + 1 + info->n_in, info->val, info->n_out); +} + +/* Check whether "bmap" has a pair of constraints that have + * the same coefficients for the output variables. + * Note that the coefficients of the existentially quantified + * variables need to be zero since the existentially quantified + * of the result are usually not the same as those of the input. + * the isl_dim_out and isl_dim_div dimensions. + * If so, return 1 and return the row indices of the two constraints + * in *first and *second. + */ +static int parallel_constraints(__isl_keep isl_basic_map *bmap, + int *first, int *second) +{ + int i; + isl_ctx *ctx = isl_basic_map_get_ctx(bmap); + struct isl_hash_table *table = NULL; + struct isl_hash_table_entry *entry; + struct isl_constraint_equal_info info; + unsigned n_out; + unsigned n_div; + + ctx = isl_basic_map_get_ctx(bmap); + table = isl_hash_table_alloc(ctx, bmap->n_ineq); + if (!table) + goto error; + + info.n_in = isl_basic_map_dim(bmap, isl_dim_param) + + isl_basic_map_dim(bmap, isl_dim_in); + info.bmap = bmap; + n_out = isl_basic_map_dim(bmap, isl_dim_out); + n_div = isl_basic_map_dim(bmap, isl_dim_div); + info.n_out = n_out + n_div; + for (i = 0; i < bmap->n_ineq; ++i) { + uint32_t hash; + + info.val = bmap->ineq[i] + 1 + info.n_in; + if (isl_seq_first_non_zero(info.val, n_out) < 0) + continue; + if (isl_seq_first_non_zero(info.val + n_out, n_div) >= 0) + continue; + hash = isl_seq_get_hash(info.val, info.n_out); + entry = isl_hash_table_find(ctx, table, hash, + constraint_equal, &info, 1); + if (!entry) + goto error; + if (entry->data) + break; + entry->data = &bmap->ineq[i]; + } + + if (i < bmap->n_ineq) { + *first = ((isl_int **)entry->data) - bmap->ineq; + *second = i; + } + + isl_hash_table_free(ctx, table); + + return i < bmap->n_ineq; +error: + isl_hash_table_free(ctx, table); + return -1; +} + +/* 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 + * of the upper bounds is smaller than all the others and assign + * this upper bound to m. + * + * In particular, if there are n bounds b_i, then the result + * consists of n basic sets, each one of the form + * + * m = b_i + * 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, + __isl_take isl_mat *var) +{ + int i, j, k; + isl_basic_set *bset = NULL; + isl_ctx *ctx; + isl_set *set = NULL; + + if (!dim || !var) + goto error; + + ctx = isl_dim_get_ctx(dim); + set = isl_set_alloc_dim(isl_dim_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, + 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); + set = isl_set_add_basic_set(set, bset); + } + + isl_dim_free(dim); + isl_mat_free(var); + return set; +error: + isl_basic_set_free(bset); + isl_set_free(set); + isl_dim_free(dim); + isl_mat_free(var); + return NULL; +} + +/* Given that the last input variable of "bmap" represents the minimum + * of the bounds in "cst", check whether we need to split the domain + * based on which bound attains the minimum. + * + * A split is needed when the minimum appears in an integer division + * or in an equality. Otherwise, it is only needed if it appears in + * 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, + __isl_keep isl_mat *cst) +{ + int i, j; + unsigned total; + unsigned pos; + + pos = cst->n_col - 1; + total = isl_basic_map_dim(bmap, isl_dim_all); + + for (i = 0; i < bmap->n_div; ++i) + if (!isl_int_is_zero(bmap->div[i][2 + pos])) + return 1; + + for (i = 0; i < bmap->n_eq; ++i) + if (!isl_int_is_zero(bmap->eq[i][1 + pos])) + return 1; + + for (i = 0; i < bmap->n_ineq; ++i) { + if (isl_int_is_nonneg(bmap->ineq[i][1 + pos])) + continue; + if (!isl_int_is_negone(bmap->ineq[i][1 + pos])) + return 1; + if (isl_seq_first_non_zero(bmap->ineq[i] + 1 + pos + 1, + total - pos - 1) >= 0) + return 1; + + for (j = 0; j < cst->n_row; ++j) + if (isl_seq_eq(bmap->ineq[i], cst->row[j], cst->n_col)) + break; + if (j >= cst->n_row) + return 1; + } + + return 0; +} + +static int need_split_set(__isl_keep isl_basic_set *bset, + __isl_keep isl_mat *cst) +{ + return need_split_map((isl_basic_map *)bset, cst); +} + +/* Given a set of which the last set variable is the minimum + * of the bounds in "cst", split each basic set in the set + * in pieces where one of the bounds is (strictly) smaller than the others. + * This subdivision is given in "min_expr". + * The variable is subsequently projected out. + * + * We only do the split when it is needed. + * For example if the last input variable m = min(a,b) and the only + * constraints in the given basic set are lower bounds on m, + * i.e., l <= m = min(a,b), then we can simply project out m + * to obtain l <= a and l <= b, without having to split on whether + * m is equal to a or b. + */ +static __isl_give isl_set *split(__isl_take isl_set *empty, + __isl_take isl_set *min_expr, __isl_take isl_mat *cst) +{ + int n_in; + int i; + isl_dim *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); + 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)) + set = isl_set_intersect(set, isl_set_copy(min_expr)); + set = isl_set_remove_dims(set, isl_dim_set, n_in - 1, 1); + + res = isl_set_union_disjoint(res, set); + } + + isl_set_free(empty); + isl_set_free(min_expr); + isl_mat_free(cst); + return res; +error: + isl_set_free(empty); + isl_set_free(min_expr); + isl_mat_free(cst); + return NULL; +} + +/* Given a map of which the last input variable is the minimum + * of the bounds in "cst", split each basic set in the set + * in pieces where one of the bounds is (strictly) smaller than the others. + * This subdivision is given in "min_expr". + * The variable is subsequently projected out. + * + * The implementation is essentially the same as that of "split". + */ +static __isl_give isl_map *split_domain(__isl_take isl_map *opt, + __isl_take isl_set *min_expr, __isl_take isl_mat *cst) +{ + int n_in; + int i; + isl_dim *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); + 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)) + map = isl_map_intersect_domain(map, + isl_set_copy(min_expr)); + map = isl_map_remove_dims(map, isl_dim_in, n_in - 1, 1); + + res = isl_map_union_disjoint(res, map); + } + + isl_map_free(opt); + isl_set_free(min_expr); + isl_mat_free(cst); + return res; +error: + isl_map_free(opt); + isl_set_free(min_expr); + isl_mat_free(cst); + return NULL; +} + +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); + +/* 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 + * constraints by a single constraint with as "input" part the minimum + * of the input parts of the list of constraints. Then, recursively call + * basic_map_partial_lexopt (possibly finding more parallel constraints) + * and plug in the definition of the minimum in the result. + * + * More specifically, given a set of constraints + * + * a x + b_i(p) >= 0 + * + * Replace this set by a single constraint + * + * a x + u >= 0 + * + * with u a new parameter with constraints + * + * u <= b_i(p) + * + * Any solution to the new system is also a solution for the original system + * since + * + * a x >= -u >= -b_i(p) + * + * 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( + __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom, + __isl_give isl_set **empty, int max, int first, int second) +{ + int i, n, k; + int *list = NULL; + unsigned n_in, n_out, n_div; + isl_ctx *ctx; + isl_vec *var = NULL; + isl_mat *cst = NULL; + isl_map *opt; + isl_set *min_expr; + isl_dim *map_dim, *set_dim; + + map_dim = isl_basic_map_get_dim(bmap); + set_dim = empty ? isl_basic_set_get_dim(dom) : NULL; + + n_in = isl_basic_map_dim(bmap, isl_dim_param) + + isl_basic_map_dim(bmap, isl_dim_in); + n_out = isl_basic_map_dim(bmap, isl_dim_all) - n_in; + + ctx = isl_basic_map_get_ctx(bmap); + list = isl_alloc_array(ctx, int, bmap->n_ineq); + var = isl_vec_alloc(ctx, n_out); + if (!list || !var) + goto error; + + list[0] = first; + list[1] = second; + isl_seq_cpy(var->el, bmap->ineq[first] + 1 + n_in, n_out); + for (i = second + 1, n = 2; i < bmap->n_ineq; ++i) { + if (isl_seq_eq(var->el, bmap->ineq[i] + 1 + n_in, n_out)) + list[n++] = i; + } + + cst = isl_mat_alloc(ctx, n, 1 + n_in); + if (!cst) + goto error; + + for (i = 0; i < n; ++i) + isl_seq_cpy(cst->row[i], bmap->ineq[list[i]], 1 + n_in); + + bmap = isl_basic_map_cow(bmap); + if (!bmap) + goto error; + for (i = n - 1; i >= 0; --i) + if (isl_basic_map_drop_inequality(bmap, list[i]) < 0) + goto error; + + bmap = isl_basic_map_add(bmap, isl_dim_in, 1); + bmap = isl_basic_map_extend_constraints(bmap, 0, 1); + k = isl_basic_map_alloc_inequality(bmap); + if (k < 0) + goto error; + isl_seq_clr(bmap->ineq[k], 1 + n_in); + isl_int_set_si(bmap->ineq[k][1 + n_in], 1); + isl_seq_cpy(bmap->ineq[k] + 1 + n_in + 1, var->el, n_out); + 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_extend_constraints(dom, 0, n); + for (i = 0; i < n; ++i) { + k = isl_basic_set_alloc_inequality(dom); + if (k < 0) + goto error; + isl_seq_cpy(dom->ineq[k], cst->row[i], 1 + n_in); + isl_int_set_si(dom->ineq[k][1 + n_in], -1); + 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; +error: + isl_dim_free(map_dim); + isl_dim_free(set_dim); + isl_mat_free(cst); + isl_vec_free(var); + free(list); + isl_basic_set_free(dom); + isl_basic_map_free(bmap); + return NULL; +} + +/* Recursive part of isl_tab_basic_map_partial_lexopt, 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 and then call + * this function recursively to look for more parallel constraints. + */ +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) +{ + 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(bmap, dom, empty, max); + + return basic_map_partial_lexopt_symm(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 map. + * 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. + */ +struct isl_map *isl_tab_basic_map_partial_lexopt( + struct isl_basic_map *bmap, struct isl_basic_set *dom, + struct 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); + + 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(bmap, dom, empty, max); +error: + isl_basic_set_free(dom); + isl_basic_map_free(bmap); + return NULL; +} + struct isl_sol_for { struct isl_sol sol; int (*fn)(__isl_take isl_basic_set *dom, @@ -3946,7 +4467,7 @@ static void sol_for_free_wrap(struct isl_sol *sol) * * Instead of constructing a basic map, this function calls a user * defined function with the current context as a basic set and - * an affine matrix reprenting the relation between the input and output. + * an affine matrix representing the relation between the input and output. * The number of rows in this matrix is equal to one plus the number * of output variables. The number of columns is equal to one plus * the total dimension of the context, i.e., the number of parameters,