X-Git-Url: http://review.tizen.org/git/?a=blobdiff_plain;f=isl_tab_pip.c;h=41319feac3577c3dc46c788d775314f0dcc82e3a;hb=c266d821693f7c6b9ff24e73a4b5f661108a019d;hp=4bc8782475a5eef6aa803a23ebc0a346b1a0ff88;hpb=efeec2bdbfb2869aa8d509ed8cd494a3ee2b87ec;p=platform%2Fupstream%2Fisl.git diff --git a/isl_tab_pip.c b/isl_tab_pip.c index 4bc8782..41319fe 100644 --- a/isl_tab_pip.c +++ b/isl_tab_pip.c @@ -1,6 +1,16 @@ +/* + * Copyright 2008-2009 Katholieke Universiteit Leuven + * + * 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 + */ + #include "isl_map_private.h" #include "isl_seq.h" #include "isl_tab.h" +#include "isl_sample.h" /* * The implementation of parametric integer linear programming in this file @@ -31,13 +41,88 @@ * Taking as initial smaple value x' = 0 corresponds to x = -M, * which is always smaller than any possible value of x. * - * We use the big parameter trick both in the main tableau and - * the context tableau, each of course having its own big parameter. + * The big parameter trick is used in the main tableau and + * also in the context tableau if isl_context_lex is used. + * In this case, each tableaus has its own big parameter. * Before doing any real work, we check if all the parameters * happen to be non-negative. If so, we drop the column corresponding * to M from the initial context tableau. + * If isl_context_gbr is used, then the big parameter trick is only + * used in the main tableau. */ +struct isl_context; +struct isl_context_op { + /* detect nonnegative parameters in context and mark them in tab */ + struct isl_tab *(*detect_nonnegative_parameters)( + struct isl_context *context, struct isl_tab *tab); + /* return temporary reference to basic set representation of context */ + struct isl_basic_set *(*peek_basic_set)(struct isl_context *context); + /* return temporary reference to tableau representation of context */ + struct isl_tab *(*peek_tab)(struct isl_context *context); + /* add equality; check is 1 if eq may not be valid; + * update is 1 if we may want to call ineq_sign on context later. + */ + void (*add_eq)(struct isl_context *context, isl_int *eq, + int check, int update); + /* add inequality; check is 1 if ineq may not be valid; + * update is 1 if we may want to call ineq_sign on context later. + */ + void (*add_ineq)(struct isl_context *context, isl_int *ineq, + int check, int update); + /* check sign of ineq based on previous information. + * strict is 1 if saturation should be treated as a positive sign. + */ + enum isl_tab_row_sign (*ineq_sign)(struct isl_context *context, + isl_int *ineq, int strict); + /* check if inequality maintains feasibility */ + int (*test_ineq)(struct isl_context *context, isl_int *ineq); + /* return index of a div that corresponds to "div" */ + int (*get_div)(struct isl_context *context, struct isl_tab *tab, + struct isl_vec *div); + /* add div "div" to context and return non-negativity */ + int (*add_div)(struct isl_context *context, struct isl_vec *div); + int (*detect_equalities)(struct isl_context *context, + struct isl_tab *tab); + /* return row index of "best" split */ + int (*best_split)(struct isl_context *context, struct isl_tab *tab); + /* check if context has already been determined to be empty */ + int (*is_empty)(struct isl_context *context); + /* check if context is still usable */ + int (*is_ok)(struct isl_context *context); + /* save a copy/snapshot of context */ + void *(*save)(struct isl_context *context); + /* restore saved context */ + void (*restore)(struct isl_context *context, void *); + /* invalidate context */ + void (*invalidate)(struct isl_context *context); + /* free context */ + void (*free)(struct isl_context *context); +}; + +struct isl_context { + struct isl_context_op *op; +}; + +struct isl_context_lex { + struct isl_context context; + struct isl_tab *tab; +}; + +struct isl_partial_sol { + int level; + struct isl_basic_set *dom; + struct isl_mat *M; + + struct isl_partial_sol *next; +}; + +struct isl_sol; +struct isl_sol_callback { + struct isl_tab_callback callback; + struct isl_sol *sol; +}; + /* isl_sol is an interface for constructing a solution to * a parametric integer linear programming problem. * Every time the algorithm reaches a state where a solution @@ -55,55 +140,225 @@ * the solution. */ struct isl_sol { - struct isl_tab *context_tab; - struct isl_sol *(*add)(struct isl_sol *sol, struct isl_tab *tab); + int error; + int rational; + int level; + int max; + int n_out; + struct isl_context *context; + struct isl_partial_sol *partial; + void (*add)(struct isl_sol *sol, + struct isl_basic_set *dom, struct isl_mat *M); + void (*add_empty)(struct isl_sol *sol, struct isl_basic_set *bset); void (*free)(struct isl_sol *sol); + struct isl_sol_callback dec_level; }; static void sol_free(struct isl_sol *sol) { + struct isl_partial_sol *partial, *next; if (!sol) return; + for (partial = sol->partial; partial; partial = next) { + next = partial->next; + isl_basic_set_free(partial->dom); + isl_mat_free(partial->M); + free(partial); + } sol->free(sol); } -struct isl_sol_map { - struct isl_sol sol; - struct isl_map *map; - struct isl_set *empty; - int max; -}; - -static void sol_map_free(struct isl_sol_map *sol_map) +/* Push a partial solution represented by a domain and mapping M + * onto the stack of partial solutions. + */ +static void sol_push_sol(struct isl_sol *sol, + struct isl_basic_set *dom, struct isl_mat *M) { - isl_tab_free(sol_map->sol.context_tab); - isl_map_free(sol_map->map); - isl_set_free(sol_map->empty); - free(sol_map); + struct isl_partial_sol *partial; + + if (sol->error || !dom) + goto error; + + partial = isl_alloc_type(dom->ctx, struct isl_partial_sol); + if (!partial) + goto error; + + partial->level = sol->level; + partial->dom = dom; + partial->M = M; + partial->next = sol->partial; + + sol->partial = partial; + + return; +error: + isl_basic_set_free(dom); + sol->error = 1; } -static void sol_map_free_wrap(struct isl_sol *sol) +/* Pop one partial solution from the partial solution stack and + * pass it on to sol->add or sol->add_empty. + */ +static void sol_pop_one(struct isl_sol *sol) { - sol_map_free((struct isl_sol_map *)sol); + struct isl_partial_sol *partial; + + partial = sol->partial; + sol->partial = partial->next; + + if (partial->M) + sol->add(sol, partial->dom, partial->M); + else + sol->add_empty(sol, partial->dom); + free(partial); } -static struct isl_sol_map *add_empty(struct isl_sol_map *sol) +/* Return a fresh copy of the domain represented by the context tableau. + */ +static struct isl_basic_set *sol_domain(struct isl_sol *sol) { struct isl_basic_set *bset; - if (!sol->empty) - return sol; - sol->empty = isl_set_grow(sol->empty, 1); - bset = isl_basic_set_copy(sol->sol.context_tab->bset); - bset = isl_basic_set_simplify(bset); - bset = isl_basic_set_finalize(bset); - sol->empty = isl_set_add(sol->empty, bset); - if (!sol->empty) - goto error; - return sol; -error: - sol_map_free(sol); - return NULL; + if (sol->error) + return NULL; + + bset = isl_basic_set_dup(sol->context->op->peek_basic_set(sol->context)); + bset = isl_basic_set_update_from_tab(bset, + sol->context->op->peek_tab(sol->context)); + + return bset; +} + +/* Check whether two partial solutions have the same mapping, where n_div + * is the number of divs that the two partial solutions have in common. + */ +static int same_solution(struct isl_partial_sol *s1, struct isl_partial_sol *s2, + unsigned n_div) +{ + int i; + unsigned dim; + + if (!s1->M != !s2->M) + return 0; + if (!s1->M) + return 1; + + dim = isl_basic_set_total_dim(s1->dom) - s1->dom->n_div; + + for (i = 0; i < s1->M->n_row; ++i) { + if (isl_seq_first_non_zero(s1->M->row[i]+1+dim+n_div, + s1->M->n_col-1-dim-n_div) != -1) + return 0; + if (isl_seq_first_non_zero(s2->M->row[i]+1+dim+n_div, + s2->M->n_col-1-dim-n_div) != -1) + return 0; + if (!isl_seq_eq(s1->M->row[i], s2->M->row[i], 1+dim+n_div)) + return 0; + } + return 1; +} + +/* Pop all solutions from the partial solution stack that were pushed onto + * the stack at levels that are deeper than the current level. + * If the two topmost elements on the stack have the same level + * and represent the same solution, then their domains are combined. + * This combined domain is the same as the current context domain + * as sol_pop is called each time we move back to a higher level. + */ +static void sol_pop(struct isl_sol *sol) +{ + struct isl_partial_sol *partial; + unsigned n_div; + + if (sol->error) + return; + + if (sol->level == 0) { + for (partial = sol->partial; partial; partial = sol->partial) + sol_pop_one(sol); + return; + } + + partial = sol->partial; + if (!partial) + return; + + if (partial->level <= sol->level) + return; + + if (partial->next && partial->next->level == partial->level) { + n_div = isl_basic_set_dim( + sol->context->op->peek_basic_set(sol->context), + isl_dim_div); + + if (!same_solution(partial, partial->next, n_div)) { + sol_pop_one(sol); + sol_pop_one(sol); + } else { + struct isl_basic_set *bset; + + bset = sol_domain(sol); + + isl_basic_set_free(partial->next->dom); + partial->next->dom = bset; + partial->next->level = sol->level; + + sol->partial = partial->next; + isl_basic_set_free(partial->dom); + isl_mat_free(partial->M); + free(partial); + } + } else + sol_pop_one(sol); +} + +static void sol_dec_level(struct isl_sol *sol) +{ + if (sol->error) + return; + + sol->level--; + + sol_pop(sol); +} + +static int sol_dec_level_wrap(struct isl_tab_callback *cb) +{ + struct isl_sol_callback *callback = (struct isl_sol_callback *)cb; + + sol_dec_level(callback->sol); + + return callback->sol->error ? -1 : 0; +} + +/* Move down to next level and push callback onto context tableau + * to decrease the level again when it gets rolled back across + * the current state. That is, dec_level will be called with + * the context tableau in the same state as it is when inc_level + * is called. + */ +static void sol_inc_level(struct isl_sol *sol) +{ + struct isl_tab *tab; + + if (sol->error) + return; + + sol->level++; + tab = sol->context->op->peek_tab(sol->context); + if (isl_tab_push_callback(tab, &sol->dec_level.callback) < 0) + sol->error = 1; +} + +static void scale_rows(struct isl_mat *mat, isl_int m, int n_row) +{ + int i; + + if (isl_int_is_one(m)) + return; + + for (i = 0; i < n_row; ++i) + isl_seq_scale(mat->row[i], mat->row[i], m, mat->n_col); } /* Add the solution identified by the tableau and the context tableau. @@ -119,23 +374,23 @@ error: * dimensions in the input map * tab->n_div is equal to the number of divs in the context * - * If there is no solution, then the basic set corresponding to the - * context tableau is added to the set "empty". + * If there is no solution, then call add_empty with a basic set + * that corresponds to the context tableau. (If add_empty is NULL, + * then do nothing). * - * Otherwise, a basic map is constructed with the same parameters - * and divs as the context, the dimensions of the context as input - * dimensions and a number of output dimensions that is equal to - * the number of output dimensions in the input map. + * If there is a solution, then first construct a matrix that maps + * all dimensions of the context to the output variables, i.e., + * the output dimensions in the input map. * The divs in the input map (if any) that do not correspond to any * div in the context do not appear in the solution. * The algorithm will make sure that they have an integer value, * but these values themselves are of no interest. + * We have to be careful not to drop or rearrange any divs in the + * context because that would change the meaning of the matrix. * - * The constraints and divs of the context are simply copied - * fron context_tab->bset. * To extract the value of the output variables, it should be noted - * that we always use a big parameter M and so the variable stored - * in the tableau is not an output variable x itself, but + * that we always use a big parameter M in the main tableau and so + * the variable stored in this tableau is not an output variable x itself, but * x' = M + x (in case of minimization) * or * x' = M - x (in case of maximization) @@ -146,178 +401,265 @@ error: * are bounded, so this cannot occur. * Similarly, when x' appears in a row, then the coefficient of M in that * row is necessarily 1. - * If the row represents + * If the row in the tableau represents * d x' = c + d M + e(y) - * then, in case of minimization, an equality - * c + e(y) - d x' = 0 - * is added, and in case of maximization, - * c + e(y) + d x' = 0 + * then, in case of minimization, the corresponding row in the matrix + * will be + * a c + a e(y) + * with a d = m, the (updated) common denominator of the matrix. + * In case of maximization, the row will be + * -a c - a e(y) */ -static struct isl_sol_map *sol_map_add(struct isl_sol_map *sol, - struct isl_tab *tab) +static void sol_add(struct isl_sol *sol, struct isl_tab *tab) +{ + struct isl_basic_set *bset = NULL; + struct isl_mat *mat = NULL; + unsigned off; + int row, i; + isl_int m; + + if (sol->error || !tab) + goto error; + + if (tab->empty && !sol->add_empty) + return; + + bset = sol_domain(sol); + + if (tab->empty) { + sol_push_sol(sol, bset, NULL); + return; + } + + off = 2 + tab->M; + + mat = isl_mat_alloc(tab->mat->ctx, 1 + sol->n_out, + 1 + tab->n_param + tab->n_div); + if (!mat) + goto error; + + isl_int_init(m); + + isl_seq_clr(mat->row[0] + 1, mat->n_col - 1); + isl_int_set_si(mat->row[0][0], 1); + for (row = 0; row < sol->n_out; ++row) { + int i = tab->n_param + row; + int r, j; + + 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); + 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); + 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); + isl_int_divexact(m, mat->row[0][0], tab->mat->row[r][0]); + isl_int_mul(mat->row[1 + row][0], m, tab->mat->row[r][1]); + for (j = 0; j < tab->n_param; ++j) { + int col; + if (tab->var[j].is_row) + continue; + col = tab->var[j].index; + isl_int_mul(mat->row[1 + row][1 + j], m, + tab->mat->row[r][off + col]); + } + for (j = 0; j < tab->n_div; ++j) { + int col; + if (tab->var[tab->n_var - tab->n_div+j].is_row) + continue; + col = tab->var[tab->n_var - tab->n_div+j].index; + isl_int_mul(mat->row[1 + row][1 + tab->n_param + j], m, + tab->mat->row[r][off + col]); + } + if (sol->max) + isl_seq_neg(mat->row[1 + row], mat->row[1 + row], + mat->n_col); + } + + isl_int_clear(m); + + sol_push_sol(sol, bset, mat); + return; +error2: + isl_int_clear(m); +error: + isl_basic_set_free(bset); + isl_mat_free(mat); + sol_free(sol); +} + +struct isl_sol_map { + struct isl_sol sol; + struct isl_map *map; + struct isl_set *empty; +}; + +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); + isl_set_free(sol_map->empty); + free(sol_map); +} + +static void sol_map_free_wrap(struct isl_sol *sol) +{ + sol_map_free((struct isl_sol_map *)sol); +} + +/* 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_map_add_empty(struct isl_sol_map *sol, + struct 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, isl_basic_set_copy(bset)); + if (!sol->empty) + goto error; + isl_basic_set_free(bset); + return; +error: + isl_basic_set_free(bset); + sol->sol.error = 1; +} + +static void sol_map_add_empty_wrap(struct isl_sol *sol, + struct isl_basic_set *bset) +{ + 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 + * and divs as the context, the dimensions of the context as input + * dimensions and a number of output dimensions that is equal to + * the number of output dimensions in the input map. + * + * The constraints and divs of the context are simply copied + * from "dom". For each row + * x = c + e(y) + * an equality + * c + e(y) - d x = 0 + * is added, with d the common denominator of M. + */ +static void sol_map_add(struct isl_sol_map *sol, + struct isl_basic_set *dom, struct isl_mat *M) { int i; struct isl_basic_map *bmap = NULL; - struct isl_tab *context_tab; + isl_basic_set *context_bset; unsigned n_eq; unsigned n_ineq; unsigned nparam; unsigned total; unsigned n_div; unsigned n_out; - unsigned off; - if (!sol || !tab) + if (sol->sol.error || !dom || !M) goto error; - if (tab->empty) - return add_empty(sol); - - context_tab = sol->sol.context_tab; - off = 2 + tab->M; - n_out = isl_map_dim(sol->map, isl_dim_out); - n_eq = context_tab->bset->n_eq + n_out; - n_ineq = context_tab->bset->n_ineq; - nparam = tab->n_param; + n_out = sol->sol.n_out; + n_eq = dom->n_eq + n_out; + n_ineq = dom->n_ineq; + 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), - tab->n_div, n_eq, 2 * tab->n_div + n_ineq); + n_div, n_eq, 2 * n_div + n_ineq); if (!bmap) goto error; - n_div = tab->n_div; - if (tab->rational) + if (sol->sol.rational) ISL_F_SET(bmap, ISL_BASIC_MAP_RATIONAL); - for (i = 0; i < context_tab->bset->n_div; ++i) { + for (i = 0; i < dom->n_div; ++i) { int k = isl_basic_map_alloc_div(bmap); if (k < 0) goto error; - isl_seq_cpy(bmap->div[k], - context_tab->bset->div[i], 1 + 1 + nparam); + isl_seq_cpy(bmap->div[k], dom->div[i], 1 + 1 + nparam); isl_seq_clr(bmap->div[k] + 1 + 1 + nparam, total - nparam); isl_seq_cpy(bmap->div[k] + 1 + 1 + total, - context_tab->bset->div[i] + 1 + 1 + nparam, i); + dom->div[i] + 1 + 1 + nparam, i); } - for (i = 0; i < context_tab->bset->n_eq; ++i) { + for (i = 0; i < dom->n_eq; ++i) { int k = isl_basic_map_alloc_equality(bmap); if (k < 0) goto error; - isl_seq_cpy(bmap->eq[k], context_tab->bset->eq[i], 1 + nparam); + isl_seq_cpy(bmap->eq[k], dom->eq[i], 1 + nparam); isl_seq_clr(bmap->eq[k] + 1 + nparam, total - nparam); isl_seq_cpy(bmap->eq[k] + 1 + total, - context_tab->bset->eq[i] + 1 + nparam, n_div); + dom->eq[i] + 1 + nparam, n_div); } - for (i = 0; i < context_tab->bset->n_ineq; ++i) { + for (i = 0; i < dom->n_ineq; ++i) { int k = isl_basic_map_alloc_inequality(bmap); if (k < 0) goto error; - isl_seq_cpy(bmap->ineq[k], - context_tab->bset->ineq[i], 1 + nparam); + isl_seq_cpy(bmap->ineq[k], dom->ineq[i], 1 + nparam); isl_seq_clr(bmap->ineq[k] + 1 + nparam, total - nparam); isl_seq_cpy(bmap->ineq[k] + 1 + total, - context_tab->bset->ineq[i] + 1 + nparam, n_div); + dom->ineq[i] + 1 + nparam, n_div); } - for (i = tab->n_param; i < total; ++i) { + for (i = 0; i < M->n_row - 1; ++i) { int k = isl_basic_map_alloc_equality(bmap); if (k < 0) goto error; - isl_seq_clr(bmap->eq[k] + 1, isl_basic_map_total_dim(bmap)); - if (!tab->var[i].is_row) { - /* no unbounded */ - isl_assert(bmap->ctx, !tab->M, goto error); - isl_int_set_si(bmap->eq[k][0], 0); - if (sol->max) - isl_int_set_si(bmap->eq[k][1 + i], 1); - else - isl_int_set_si(bmap->eq[k][1 + i], -1); - } else { - int row, j; - row = tab->var[i].index; - /* no unbounded */ - if (tab->M) - isl_assert(bmap->ctx, - isl_int_eq(tab->mat->row[row][2], - tab->mat->row[row][0]), - goto error); - isl_int_set(bmap->eq[k][0], tab->mat->row[row][1]); - for (j = 0; j < tab->n_param; ++j) { - int col; - if (tab->var[j].is_row) - continue; - col = tab->var[j].index; - isl_int_set(bmap->eq[k][1 + j], - tab->mat->row[row][off + col]); - } - for (j = 0; j < tab->n_div; ++j) { - int col; - if (tab->var[tab->n_var - tab->n_div+j].is_row) - continue; - col = tab->var[tab->n_var - tab->n_div+j].index; - isl_int_set(bmap->eq[k][1 + total + j], - tab->mat->row[row][off + col]); - } - if (sol->max) - isl_int_set(bmap->eq[k][1 + i], - tab->mat->row[row][0]); - else - isl_int_neg(bmap->eq[k][1 + i], - tab->mat->row[row][0]); - } + isl_seq_cpy(bmap->eq[k], M->row[1 + i], 1 + nparam); + isl_seq_clr(bmap->eq[k] + 1 + nparam, n_out); + isl_int_neg(bmap->eq[k][1 + nparam + i], M->row[0][0]); + isl_seq_cpy(bmap->eq[k] + 1 + nparam + n_out, + M->row[1 + i] + 1 + nparam, n_div); } - bmap = isl_basic_map_gauss(bmap, NULL); - bmap = isl_basic_map_normalize_constraints(bmap); + bmap = isl_basic_map_simplify(bmap); bmap = isl_basic_map_finalize(bmap); sol->map = isl_map_grow(sol->map, 1); - sol->map = isl_map_add(sol->map, bmap); + sol->map = isl_map_add_basic_map(sol->map, bmap); if (!sol->map) goto error; - return sol; + isl_basic_set_free(dom); + isl_mat_free(M); + return; error: + isl_basic_set_free(dom); + isl_mat_free(M); isl_basic_map_free(bmap); - sol_free(&sol->sol); - return NULL; + sol->sol.error = 1; } -static struct isl_sol *sol_map_add_wrap(struct isl_sol *sol, - struct isl_tab *tab) +static void sol_map_add_wrap(struct isl_sol *sol, + struct isl_basic_set *dom, struct isl_mat *M) { - return (struct isl_sol *)sol_map_add((struct isl_sol_map *)sol, tab); -} - - -static struct isl_basic_set *isl_basic_set_add_ineq(struct isl_basic_set *bset, - isl_int *ineq) -{ - int k; - - bset = isl_basic_set_extend_constraints(bset, 0, 1); - if (!bset) - return NULL; - k = isl_basic_set_alloc_inequality(bset); - if (k < 0) - goto error; - isl_seq_cpy(bset->ineq[k], ineq, 1 + isl_basic_set_total_dim(bset)); - return bset; -error: - isl_basic_set_free(bset); - return NULL; -} - -static struct isl_basic_set *isl_basic_set_add_eq(struct isl_basic_set *bset, - isl_int *eq) -{ - int k; - - bset = isl_basic_set_extend_constraints(bset, 1, 0); - if (!bset) - return NULL; - k = isl_basic_set_alloc_equality(bset); - if (k < 0) - goto error; - isl_seq_cpy(bset->eq[k], eq, 1 + isl_basic_set_total_dim(bset)); - return bset; -error: - isl_basic_set_free(bset); - return NULL; + sol_map_add((struct isl_sol_map *)sol, dom, M); } @@ -476,6 +818,9 @@ static struct isl_vec *ineq_for_div(struct isl_basic_set *bset, unsigned div) unsigned div_pos; struct isl_vec *ineq; + if (!bset) + return NULL; + total = isl_basic_set_total_dim(bset); div_pos = 1 + total - bset->n_div + div; @@ -503,19 +848,25 @@ static struct isl_vec *ineq_for_div(struct isl_basic_set *bset, unsigned div) */ static struct isl_tab *set_row_cst_to_div(struct isl_tab *tab, int row, int div) { - int col; - unsigned off = 2 + tab->M; - isl_seq_fdiv_q(tab->mat->row[row] + 1, tab->mat->row[row] + 1, tab->mat->row[row][0], 1 + tab->M + tab->n_col); isl_int_set_si(tab->mat->row[row][0], 1); - isl_assert(tab->mat->ctx, - !tab->var[tab->n_var - tab->n_div + div].is_row, goto error); + if (tab->var[tab->n_var - tab->n_div + div].is_row) { + int drow = tab->var[tab->n_var - tab->n_div + div].index; + + isl_assert(tab->mat->ctx, + isl_int_is_one(tab->mat->row[drow][0]), goto error); + isl_seq_combine(tab->mat->row[row] + 1, + tab->mat->ctx->one, tab->mat->row[row] + 1, + tab->mat->ctx->one, tab->mat->row[drow] + 1, + 1 + tab->M + tab->n_col); + } else { + int dcol = tab->var[tab->n_var - tab->n_div + div].index; - col = tab->var[tab->n_var - tab->n_div + div].index; - isl_int_set_si(tab->mat->row[row][off + col], 1); + isl_int_set_si(tab->mat->row[row][2 + tab->M + dcol], 1); + } return tab; error: @@ -735,8 +1086,11 @@ static int first_neg(struct isl_tab *tab) for (row = tab->n_redundant; row < tab->n_row; ++row) { if (!isl_tab_var_from_row(tab, row)->is_nonneg) continue; - if (isl_int_is_neg(tab->mat->row[row][2])) - return row; + if (!isl_int_is_neg(tab->mat->row[row][2])) + continue; + if (tab->row_sign) + tab->row_sign[row] = isl_tab_row_neg; + return row; } for (row = tab->n_redundant; row < tab->n_row; ++row) { if (!isl_tab_var_from_row(tab, row)->is_nonneg) @@ -760,6 +1114,7 @@ static int first_neg(struct isl_tab *tab) * smallest increment in the sample point. If there is no such column * then the tableau is infeasible. */ +static struct isl_tab *restore_lexmin(struct isl_tab *tab) WARN_UNUSED; static struct isl_tab *restore_lexmin(struct isl_tab *tab) { int row, col; @@ -770,11 +1125,15 @@ static struct isl_tab *restore_lexmin(struct isl_tab *tab) return tab; while ((row = first_neg(tab)) != -1) { col = lexmin_pivot_col(tab, row); - if (col >= tab->n_col) - return isl_tab_mark_empty(tab); + if (col >= tab->n_col) { + if (isl_tab_mark_empty(tab) < 0) + goto error; + return tab; + } if (col < 0) goto error; - isl_tab_pivot(tab, row, col); + if (isl_tab_pivot(tab, row, col) < 0) + goto error; } return tab; error: @@ -848,16 +1207,20 @@ static struct isl_tab *add_lexmin_valid_eq(struct isl_tab *tab, isl_int *eq) i = last_var_col_or_int_par_col(tab, r); if (i < 0) { tab->con[r].is_nonneg = 1; - isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]); + if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0) + goto error; isl_seq_neg(eq, eq, 1 + tab->n_var); r = isl_tab_add_row(tab, eq); if (r < 0) goto error; tab->con[r].is_nonneg = 1; - isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]); + if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0) + goto error; } else { - isl_tab_pivot(tab, r, i); - isl_tab_kill_col(tab, i); + if (isl_tab_pivot(tab, r, i) < 0) + goto error; + if (isl_tab_kill_col(tab, i) < 0) + goto error; tab->n_eq++; tab = restore_lexmin(tab); @@ -888,30 +1251,33 @@ static int is_constant(struct isl_tab *tab, int row) * In the end we try to use one of the two constraints to eliminate * a column. */ +static struct isl_tab *add_lexmin_eq(struct isl_tab *tab, isl_int *eq) WARN_UNUSED; static struct isl_tab *add_lexmin_eq(struct isl_tab *tab, isl_int *eq) { int r1, r2; int row; + struct isl_tab_undo *snap; if (!tab) return NULL; - if (tab->bset) { - tab->bset = isl_basic_set_add_eq(tab->bset, eq); - isl_tab_push(tab, isl_tab_undo_bset_eq); - if (!tab->bset) - goto error; - } + snap = isl_tab_snap(tab); r1 = isl_tab_add_row(tab, eq); if (r1 < 0) goto error; tab->con[r1].is_nonneg = 1; - isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r1]); + if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r1]) < 0) + goto error; row = tab->con[r1].index; if (is_constant(tab, row)) { if (!isl_int_is_zero(tab->mat->row[row][1]) || - (tab->M && !isl_int_is_zero(tab->mat->row[row][2]))) - return isl_tab_mark_empty(tab); + (tab->M && !isl_int_is_zero(tab->mat->row[row][2]))) { + if (isl_tab_mark_empty(tab) < 0) + goto error; + return tab; + } + if (isl_tab_rollback(tab, snap) < 0) + goto error; return tab; } @@ -925,28 +1291,46 @@ static struct isl_tab *add_lexmin_eq(struct isl_tab *tab, isl_int *eq) if (r2 < 0) goto error; tab->con[r2].is_nonneg = 1; - isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r2]); + if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r2]) < 0) + goto error; tab = restore_lexmin(tab); if (!tab || tab->empty) return tab; - if (!tab->con[r1].is_row) - isl_tab_kill_col(tab, tab->con[r1].index); - else if (!tab->con[r2].is_row) - isl_tab_kill_col(tab, tab->con[r2].index); - else if (isl_int_is_zero(tab->mat->row[tab->con[r1].index][1])) { + if (!tab->con[r1].is_row) { + if (isl_tab_kill_col(tab, tab->con[r1].index) < 0) + goto error; + } else if (!tab->con[r2].is_row) { + if (isl_tab_kill_col(tab, tab->con[r2].index) < 0) + goto error; + } else if (isl_int_is_zero(tab->mat->row[tab->con[r1].index][1])) { unsigned off = 2 + tab->M; int i; int row = tab->con[r1].index; i = isl_seq_first_non_zero(tab->mat->row[row]+off+tab->n_dead, tab->n_col - tab->n_dead); if (i != -1) { - isl_tab_pivot(tab, row, tab->n_dead + i); - isl_tab_kill_col(tab, tab->n_dead + i); + if (isl_tab_pivot(tab, row, tab->n_dead + i) < 0) + goto error; + if (isl_tab_kill_col(tab, tab->n_dead + i) < 0) + goto error; } } + if (tab->bmap) { + tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq); + if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0) + goto error; + isl_seq_neg(eq, eq, 1 + tab->n_var); + tab->bmap = isl_basic_map_add_ineq(tab->bmap, eq); + isl_seq_neg(eq, eq, 1 + tab->n_var); + if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0) + goto error; + if (!tab->bmap) + goto error; + } + return tab; error: isl_tab_free(tab); @@ -962,26 +1346,30 @@ static struct isl_tab *add_lexmin_ineq(struct isl_tab *tab, isl_int *ineq) if (!tab) return NULL; - if (tab->bset) { - tab->bset = isl_basic_set_add_ineq(tab->bset, ineq); - isl_tab_push(tab, isl_tab_undo_bset_ineq); - if (!tab->bset) + if (tab->bmap) { + tab->bmap = isl_basic_map_add_ineq(tab->bmap, ineq); + if (isl_tab_push(tab, isl_tab_undo_bmap_ineq) < 0) + goto error; + if (!tab->bmap) goto error; } r = isl_tab_add_row(tab, ineq); if (r < 0) goto error; tab->con[r].is_nonneg = 1; - isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]); + if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0) + goto error; if (isl_tab_row_is_redundant(tab, tab->con[r].index)) { - isl_tab_mark_redundant(tab, tab->con[r].index); + if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0) + goto error; return tab; } tab = restore_lexmin(tab); if (tab && !tab->empty && tab->con[r].is_row && isl_tab_row_is_redundant(tab, tab->con[r].index)) - isl_tab_mark_redundant(tab, tab->con[r].index); + if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0) + goto error; return tab; error: isl_tab_free(tab); @@ -1023,7 +1411,7 @@ static int integer_variable(struct isl_tab *tab, int row) int i; unsigned off = 2 + tab->M; - for (i = 0; i < tab->n_col; ++i) { + for (i = tab->n_dead; i < tab->n_col; ++i) { if (tab->col_var[i] >= 0 && (tab->col_var[i] < tab->n_param || tab->col_var[i] >= tab->n_var - tab->n_div)) @@ -1047,8 +1435,9 @@ static int integer_constant(struct isl_tab *tab, int row) #define I_PAR 1 << 1 #define I_VAR 1 << 2 -/* Check for first (non-parameter) variable that is non-integer and - * therefore requires a cut. +/* Check for next (non-parameter) variable after "var" (first if var == -1) + * that is non-integer and therefore requires a cut and return + * the index of the variable. * For parametric tableaus, there are three parts in a row, * the constant, the coefficients of the parameters and the rest. * For each part, we check whether the coefficients in that part @@ -1057,16 +1446,16 @@ static int integer_constant(struct isl_tab *tab, int row) * current sample value is integral and no cut is required * (irrespective of whether the variable part is integral). */ -static int first_non_integer(struct isl_tab *tab, int *f) +static int next_non_integer_var(struct isl_tab *tab, int var, int *f) { - int i; + var = var < 0 ? tab->n_param : var + 1; - for (i = tab->n_param; i < tab->n_var - tab->n_div; ++i) { + for (; var < tab->n_var - tab->n_div; ++var) { int flags = 0; int row; - if (!tab->var[i].is_row) + if (!tab->var[var].is_row) continue; - row = tab->var[i].index; + row = tab->var[var].index; if (integer_constant(tab, row)) ISL_FL_SET(flags, I_CST); if (integer_parameter(tab, row)) @@ -1076,11 +1465,28 @@ static int first_non_integer(struct isl_tab *tab, int *f) if (integer_variable(tab, row)) ISL_FL_SET(flags, I_VAR); *f = flags; - return row; + return var; } return -1; } +/* Check for first (non-parameter) variable that is non-integer and + * therefore requires a cut and return the corresponding row. + * For parametric tableaus, there are three parts in a row, + * the constant, the coefficients of the parameters and the rest. + * For each part, we check whether the coefficients in that part + * are all integral and if so, set the corresponding flag in *f. + * If the constant and the parameter part are integral, then the + * current sample value is integral and no cut is required + * (irrespective of whether the variable part is integral). + */ +static int first_non_integer_row(struct isl_tab *tab, int *f) +{ + int var = next_non_integer_var(tab, -1, f); + + return var < 0 ? -1 : tab->var[var].index; +} + /* Add a (non-parametric) cut to cut away the non-integral sample * value of the given row. * @@ -1128,7 +1534,8 @@ static int add_cut(struct isl_tab *tab, int row) tab->mat->row[row][off + i], tab->mat->row[row][0]); tab->con[r].is_nonneg = 1; - isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]); + if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0) + return -1; if (tab->row_sign) tab->row_sign[tab->con[r].index] = isl_tab_row_neg; @@ -1139,15 +1546,17 @@ static int add_cut(struct isl_tab *tab, int row) * sample point is obtained or until the tableau is determined * to be integer infeasible. * As long as there is any non-integer value in the sample point, - * we add an appropriate cut, if possible and resolve the violated - * cut constraint using restore_lexmin. + * we add appropriate cuts, if possible, for each of these + * non-integer values and then resolve the violated + * cut constraints using restore_lexmin. * If one of the corresponding rows is equal to an integral * combination of variables/constraints plus a non-integral constant, - * then there is no way to obtain an integer point an we return + * then there is no way to obtain an integer point and we return * a tableau that is marked empty. */ static struct isl_tab *cut_to_integer_lexmin(struct isl_tab *tab) { + int var; int row; int flags; @@ -1156,12 +1565,18 @@ static struct isl_tab *cut_to_integer_lexmin(struct isl_tab *tab) if (tab->empty) return tab; - while ((row = first_non_integer(tab, &flags)) != -1) { - if (ISL_FL_ISSET(flags, I_VAR)) - return isl_tab_mark_empty(tab); - row = add_cut(tab, row); - if (row < 0) - goto error; + while ((var = next_non_integer_var(tab, -1, &flags)) != -1) { + do { + if (ISL_FL_ISSET(flags, I_VAR)) { + if (isl_tab_mark_empty(tab) < 0) + goto error; + return tab; + } + row = tab->var[var].index; + row = add_cut(tab, row); + if (row < 0) + goto error; + } while ((var = next_non_integer_var(tab, var, &flags)) != -1); tab = restore_lexmin(tab); if (!tab || tab->empty) break; @@ -1172,16 +1587,6 @@ error: return NULL; } -static struct isl_tab *drop_sample(struct isl_tab *tab, int s) -{ - if (s != tab->n_outside) - isl_mat_swap_rows(tab->samples, tab->n_outside, s); - tab->n_outside++; - isl_tab_push(tab, isl_tab_undo_drop_sample); - - return tab; -} - /* Check whether all the currently active samples also satisfy the inequality * "ineq" (treated as an equality if eq is set). * Remove those samples that do not. @@ -1194,7 +1599,7 @@ static struct isl_tab *check_samples(struct isl_tab *tab, isl_int *ineq, int eq) if (!tab) return NULL; - isl_assert(tab->mat->ctx, tab->bset, goto error); + isl_assert(tab->mat->ctx, tab->bmap, goto error); isl_assert(tab->mat->ctx, tab->samples, goto error); isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, goto error); @@ -1206,7 +1611,7 @@ static struct isl_tab *check_samples(struct isl_tab *tab, isl_int *ineq, int eq) sgn = isl_int_sgn(v); if (eq ? (sgn == 0) : (sgn >= 0)) continue; - tab = drop_sample(tab, i); + tab = isl_tab_drop_sample(tab, i); if (!tab) break; } @@ -1242,74 +1647,58 @@ static int sample_is_finite(struct isl_tab *tab) } /* Check if the context tableau of sol has any integer points. - * Returns -1 if an error occurred. + * Leave tab in empty state if no integer point can be found. * If an integer point can be found and if moreover it is finite, * then it is added to the list of sample values. * * This function is only called when none of the currently active sample * values satisfies the most recently added constraint. */ -static int context_is_feasible(struct isl_sol *sol) +static struct isl_tab *check_integer_feasible(struct isl_tab *tab) { struct isl_tab_undo *snap; - struct isl_tab *tab; int feasible; - if (!sol || !sol->context_tab) - return -1; + if (!tab) + return NULL; - snap = isl_tab_snap(sol->context_tab); - isl_tab_push_basis(sol->context_tab); + snap = isl_tab_snap(tab); + if (isl_tab_push_basis(tab) < 0) + goto error; - sol->context_tab = cut_to_integer_lexmin(sol->context_tab); - if (!sol->context_tab) + tab = cut_to_integer_lexmin(tab); + if (!tab) goto error; - tab = sol->context_tab; if (!tab->empty && sample_is_finite(tab)) { struct isl_vec *sample; - tab->samples = isl_mat_extend(tab->samples, - tab->n_sample + 1, tab->samples->n_col); - if (!tab->samples) - goto error; - sample = isl_tab_get_sample_value(tab); - if (!sample) - goto error; - isl_seq_cpy(tab->samples->row[tab->n_sample], - sample->el, sample->size); - isl_vec_free(sample); - tab->n_sample++; + + tab = isl_tab_add_sample(tab, sample); } - feasible = !sol->context_tab->empty; - if (isl_tab_rollback(sol->context_tab, snap) < 0) + if (!tab->empty && isl_tab_rollback(tab, snap) < 0) goto error; - return feasible; + return tab; error: - isl_tab_free(sol->context_tab); - sol->context_tab = NULL; - return -1; + isl_tab_free(tab); + return NULL; } -/* First check if any of the currently active sample values satisfies +/* Check if any of the currently active sample values satisfies * the inequality "ineq" (an equality if eq is set). - * If not, continue with check_integer_feasible. */ -static int context_valid_sample_or_feasible(struct isl_sol *sol, - isl_int *ineq, int eq) +static int tab_has_valid_sample(struct isl_tab *tab, isl_int *ineq, int eq) { int i; isl_int v; - struct isl_tab *tab; - if (!sol || !sol->context_tab) + if (!tab) return -1; - tab = sol->context_tab; - isl_assert(tab->mat->ctx, tab->bset, return -1); + isl_assert(tab->mat->ctx, tab->bmap, return -1); isl_assert(tab->mat->ctx, tab->samples, return -1); isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, return -1); @@ -1324,51 +1713,39 @@ static int context_valid_sample_or_feasible(struct isl_sol *sol, } isl_int_clear(v); - if (i < tab->n_sample) - return 1; - - return context_is_feasible(sol); + return i < tab->n_sample; } -/* For a div d = floor(f/m), add the constraints - * - * f - m d >= 0 - * -(f-(m-1)) + m d >= 0 - * - * Note that the second constraint is the negation of - * - * f - m d >= m +/* Add a div specifed by "div" to the tableau "tab" and return + * 1 if the div is obviously non-negative. */ -static struct isl_tab *add_div_constraints(struct isl_tab *tab, unsigned div) +static int context_tab_add_div(struct isl_tab *tab, struct isl_vec *div, + int (*add_ineq)(void *user, isl_int *), void *user) { - unsigned total; - unsigned div_pos; - struct isl_vec *ineq; - - if (!tab) - return NULL; - - total = isl_basic_set_total_dim(tab->bset); - div_pos = 1 + total - tab->bset->n_div + div; - - ineq = ineq_for_div(tab->bset, div); - if (!ineq) - goto error; - - tab = add_lexmin_ineq(tab, ineq->el); + int i; + int r; + struct isl_mat *samples; + int nonneg; - isl_seq_neg(ineq->el, tab->bset->div[div] + 1, 1 + total); - isl_int_set(ineq->el[div_pos], tab->bset->div[div][0]); - isl_int_add(ineq->el[0], ineq->el[0], ineq->el[div_pos]); - isl_int_sub_ui(ineq->el[0], ineq->el[0], 1); - tab = add_lexmin_ineq(tab, ineq->el); + r = isl_tab_add_div(tab, div, add_ineq, user); + if (r < 0) + return -1; + nonneg = tab->var[r].is_nonneg; + tab->var[r].frozen = 1; - isl_vec_free(ineq); + samples = isl_mat_extend(tab->samples, + tab->n_sample, 1 + tab->n_var); + tab->samples = samples; + if (!samples) + return -1; + for (i = tab->n_outside; i < samples->n_row; ++i) { + isl_seq_inner_product(div->el + 1, samples->row[i], + div->size - 1, &samples->row[i][samples->n_col - 1]); + isl_int_fdiv_q(samples->row[i][samples->n_col - 1], + samples->row[i][samples->n_col - 1], div->el[0]); + } - return tab; -error: - isl_tab_free(tab); - return NULL; + return nonneg; } /* Add a div specified by "div" to both the main tableau and @@ -1377,43 +1754,16 @@ error: * need to express the meaning of the div. * Return the index of the div or -1 if anything went wrong. */ -static int add_div(struct isl_tab *tab, struct isl_tab **context_tab, +static int add_div(struct isl_tab *tab, struct isl_context *context, struct isl_vec *div) { - int i; int r; - int k; - struct isl_mat *samples; - - if (isl_tab_extend_vars(*context_tab, 1) < 0) - goto error; - r = isl_tab_allocate_var(*context_tab); - if (r < 0) - goto error; - (*context_tab)->var[r].is_nonneg = 1; - (*context_tab)->var[r].frozen = 1; + int nonneg; - samples = isl_mat_extend((*context_tab)->samples, - (*context_tab)->n_sample, 1 + (*context_tab)->n_var); - (*context_tab)->samples = samples; - if (!samples) + if ((nonneg = context->op->add_div(context, div)) < 0) goto error; - for (i = (*context_tab)->n_outside; i < samples->n_row; ++i) { - isl_seq_inner_product(div->el + 1, samples->row[i], - div->size - 1, &samples->row[i][samples->n_col - 1]); - isl_int_fdiv_q(samples->row[i][samples->n_col - 1], - samples->row[i][samples->n_col - 1], div->el[0]); - } - (*context_tab)->bset = isl_basic_set_extend_dim((*context_tab)->bset, - isl_basic_set_get_dim((*context_tab)->bset), 1, 0, 2); - k = isl_basic_set_alloc_div((*context_tab)->bset); - if (k < 0) - goto error; - isl_seq_cpy((*context_tab)->bset->div[k], div->el, div->size); - isl_tab_push((*context_tab), isl_tab_undo_bset_div); - *context_tab = add_div_constraints(*context_tab, k); - if (!*context_tab) + if (!context->op->is_ok(context)) goto error; if (isl_tab_extend_vars(tab, 1) < 0) @@ -1421,27 +1771,26 @@ static int add_div(struct isl_tab *tab, struct isl_tab **context_tab, r = isl_tab_allocate_var(tab); if (r < 0) goto error; - if (!(*context_tab)->M) + if (nonneg) tab->var[r].is_nonneg = 1; tab->var[r].frozen = 1; tab->n_div++; return tab->n_div - 1; error: - isl_tab_free(*context_tab); - *context_tab = NULL; + context->op->invalidate(context); return -1; } static int find_div(struct isl_tab *tab, isl_int *div, isl_int denom) { int i; - unsigned total = isl_basic_set_total_dim(tab->bset); + unsigned total = isl_basic_map_total_dim(tab->bmap); - for (i = 0; i < tab->bset->n_div; ++i) { - if (isl_int_ne(tab->bset->div[i][0], denom)) + for (i = 0; i < tab->bmap->n_div; ++i) { + if (isl_int_ne(tab->bmap->div[i][0], denom)) continue; - if (!isl_seq_eq(tab->bset->div[i] + 1, div, total)) + if (!isl_seq_eq(tab->bmap->div[i] + 1, div, 1 + total)) continue; return i; } @@ -1451,16 +1800,20 @@ static int find_div(struct isl_tab *tab, isl_int *div, isl_int denom) /* Return the index of a div that corresponds to "div". * We first check if we already have such a div and if not, we create one. */ -static int get_div(struct isl_tab *tab, struct isl_tab **context_tab, +static int get_div(struct isl_tab *tab, struct isl_context *context, struct isl_vec *div) { int d; + struct isl_tab *context_tab = context->op->peek_tab(context); + + if (!context_tab) + return -1; - d = find_div(*context_tab, div->el + 1, div->el[0]); + d = find_div(context_tab, div->el + 1, div->el[0]); if (d != -1) return d; - return add_div(tab, context_tab, div); + return add_div(tab, context, div); } /* Add a parametric cut to cut away the non-integral sample value @@ -1487,7 +1840,7 @@ static int get_div(struct isl_tab *tab, struct isl_tab **context_tab, * Return the row of the cut or -1. */ static int add_parametric_cut(struct isl_tab *tab, int row, - struct isl_tab **context_tab) + struct isl_context *context) { struct isl_vec *div; int d; @@ -1495,21 +1848,20 @@ static int add_parametric_cut(struct isl_tab *tab, int row, int r; isl_int *r_row; int col; + int n; unsigned off = 2 + tab->M; - if (!*context_tab) - goto error; - - if (isl_tab_extend_cons(*context_tab, 3) < 0) - goto error; + if (!context) + return -1; div = get_row_parameter_div(tab, row); if (!div) return -1; - d = get_div(tab, context_tab, div); + n = tab->n_div; + d = context->op->get_div(context, tab, div); if (d < 0) - goto error; + return -1; if (isl_tab_extend_cons(tab, 1) < 0) return -1; @@ -1568,17 +1920,19 @@ static int add_parametric_cut(struct isl_tab *tab, int row, } tab->con[r].is_nonneg = 1; - isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]); + if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r]) < 0) + return -1; if (tab->row_sign) tab->row_sign[tab->con[r].index] = isl_tab_row_neg; isl_vec_free(div); - return tab->con[r].index; -error: - isl_tab_free(*context_tab); - *context_tab = NULL; - return -1; + row = tab->con[r].index; + + if (d >= n && context->op->detect_equalities(context, tab) < 0) + return -1; + + return row; } /* Construct a tableau for bmap that can be used for computing @@ -1613,8 +1967,11 @@ static struct isl_tab *tab_for_lexmin(struct isl_basic_map *bmap, if (!tab->row_sign) goto error; } - if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY)) - return isl_tab_mark_empty(tab); + if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY)) { + if (isl_tab_mark_empty(tab) < 0) + goto error; + return tab; + } for (i = tab->n_param; i < tab->n_var - tab->n_div; ++i) { tab->var[i].is_nonneg = 1; @@ -1652,28 +2009,1199 @@ error: return NULL; } -static struct isl_tab *context_tab_for_lexmin(struct isl_basic_set *bset) +/* Given a main tableau where more than one row requires a split, + * determine and return the "best" row to split on. + * + * Given two rows in the main tableau, if the inequality corresponding + * to the first row is redundant with respect to that of the second row + * in the current tableau, then it is better to split on the second row, + * since in the positive part, both row will be positive. + * (In the negative part a pivot will have to be performed and just about + * anything can happen to the sign of the other row.) + * + * As a simple heuristic, we therefore select the row that makes the most + * of the other rows redundant. + * + * Perhaps it would also be useful to look at the number of constraints + * that conflict with any given constraint. + */ +static int best_split(struct isl_tab *tab, struct isl_tab *context_tab) { - struct isl_tab *tab; + struct isl_tab_undo *snap; + int split; + int row; + int best = -1; + int best_r; - bset = isl_basic_set_cow(bset); - if (!bset) - return NULL; - tab = tab_for_lexmin((struct isl_basic_map *)bset, NULL, 1, 0); - if (!tab) - goto error; - tab->bset = bset; - tab->n_sample = 0; - tab->n_outside = 0; - tab->samples = isl_mat_alloc(bset->ctx, 1, 1 + tab->n_var); - if (!tab->samples) - goto error; - return tab; + if (isl_tab_extend_cons(context_tab, 2) < 0) + return -1; + + snap = isl_tab_snap(context_tab); + + for (split = tab->n_redundant; split < tab->n_row; ++split) { + struct isl_tab_undo *snap2; + struct isl_vec *ineq = NULL; + int r = 0; + int ok; + + if (!isl_tab_var_from_row(tab, split)->is_nonneg) + continue; + if (tab->row_sign[split] != isl_tab_row_any) + continue; + + ineq = get_row_parameter_ineq(tab, split); + if (!ineq) + return -1; + ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0; + isl_vec_free(ineq); + if (!ok) + return -1; + + snap2 = isl_tab_snap(context_tab); + + for (row = tab->n_redundant; row < tab->n_row; ++row) { + struct isl_tab_var *var; + + if (row == split) + continue; + if (!isl_tab_var_from_row(tab, row)->is_nonneg) + continue; + if (tab->row_sign[row] != isl_tab_row_any) + continue; + + ineq = get_row_parameter_ineq(tab, row); + if (!ineq) + return -1; + ok = isl_tab_add_ineq(context_tab, ineq->el) >= 0; + isl_vec_free(ineq); + if (!ok) + return -1; + var = &context_tab->con[context_tab->n_con - 1]; + if (!context_tab->empty && + !isl_tab_min_at_most_neg_one(context_tab, var)) + r++; + if (isl_tab_rollback(context_tab, snap2) < 0) + return -1; + } + if (best == -1 || r > best_r) { + best = split; + best_r = r; + } + if (isl_tab_rollback(context_tab, snap) < 0) + return -1; + } + + return best; +} + +static struct isl_basic_set *context_lex_peek_basic_set( + struct isl_context *context) +{ + struct isl_context_lex *clex = (struct isl_context_lex *)context; + if (!clex->tab) + return NULL; + return isl_tab_peek_bset(clex->tab); +} + +static struct isl_tab *context_lex_peek_tab(struct isl_context *context) +{ + struct isl_context_lex *clex = (struct isl_context_lex *)context; + return clex->tab; +} + +static void context_lex_extend(struct isl_context *context, int n) +{ + struct isl_context_lex *clex = (struct isl_context_lex *)context; + if (!clex->tab) + return; + if (isl_tab_extend_cons(clex->tab, n) >= 0) + return; + isl_tab_free(clex->tab); + clex->tab = NULL; +} + +static void context_lex_add_eq(struct isl_context *context, isl_int *eq, + int check, int update) +{ + struct isl_context_lex *clex = (struct isl_context_lex *)context; + if (isl_tab_extend_cons(clex->tab, 2) < 0) + goto error; + clex->tab = add_lexmin_eq(clex->tab, eq); + if (check) { + int v = tab_has_valid_sample(clex->tab, eq, 1); + if (v < 0) + goto error; + if (!v) + clex->tab = check_integer_feasible(clex->tab); + } + if (update) + clex->tab = check_samples(clex->tab, eq, 1); + return; +error: + isl_tab_free(clex->tab); + clex->tab = NULL; +} + +static void context_lex_add_ineq(struct isl_context *context, isl_int *ineq, + int check, int update) +{ + struct isl_context_lex *clex = (struct isl_context_lex *)context; + if (isl_tab_extend_cons(clex->tab, 1) < 0) + goto error; + clex->tab = add_lexmin_ineq(clex->tab, ineq); + if (check) { + int v = tab_has_valid_sample(clex->tab, ineq, 0); + if (v < 0) + goto error; + if (!v) + clex->tab = check_integer_feasible(clex->tab); + } + if (update) + clex->tab = check_samples(clex->tab, ineq, 0); + return; +error: + isl_tab_free(clex->tab); + clex->tab = NULL; +} + +static int context_lex_add_ineq_wrap(void *user, isl_int *ineq) +{ + struct isl_context *context = (struct isl_context *)user; + context_lex_add_ineq(context, ineq, 0, 0); + return context->op->is_ok(context) ? 0 : -1; +} + +/* Check which signs can be obtained by "ineq" on all the currently + * active sample values. See row_sign for more information. + */ +static enum isl_tab_row_sign tab_ineq_sign(struct isl_tab *tab, isl_int *ineq, + int strict) +{ + int i; + int sgn; + isl_int tmp; + enum isl_tab_row_sign res = isl_tab_row_unknown; + + isl_assert(tab->mat->ctx, tab->samples, return isl_tab_row_unknown); + isl_assert(tab->mat->ctx, tab->samples->n_col == 1 + tab->n_var, + return isl_tab_row_unknown); + + isl_int_init(tmp); + for (i = tab->n_outside; i < tab->n_sample; ++i) { + isl_seq_inner_product(tab->samples->row[i], ineq, + 1 + tab->n_var, &tmp); + sgn = isl_int_sgn(tmp); + if (sgn > 0 || (sgn == 0 && strict)) { + if (res == isl_tab_row_unknown) + res = isl_tab_row_pos; + if (res == isl_tab_row_neg) + res = isl_tab_row_any; + } + if (sgn < 0) { + if (res == isl_tab_row_unknown) + res = isl_tab_row_neg; + if (res == isl_tab_row_pos) + res = isl_tab_row_any; + } + if (res == isl_tab_row_any) + break; + } + isl_int_clear(tmp); + + return res; +} + +static enum isl_tab_row_sign context_lex_ineq_sign(struct isl_context *context, + isl_int *ineq, int strict) +{ + struct isl_context_lex *clex = (struct isl_context_lex *)context; + return tab_ineq_sign(clex->tab, ineq, strict); +} + +/* Check whether "ineq" can be added to the tableau without rendering + * it infeasible. + */ +static int context_lex_test_ineq(struct isl_context *context, isl_int *ineq) +{ + struct isl_context_lex *clex = (struct isl_context_lex *)context; + struct isl_tab_undo *snap; + int feasible; + + if (!clex->tab) + return -1; + + if (isl_tab_extend_cons(clex->tab, 1) < 0) + return -1; + + snap = isl_tab_snap(clex->tab); + if (isl_tab_push_basis(clex->tab) < 0) + return -1; + clex->tab = add_lexmin_ineq(clex->tab, ineq); + clex->tab = check_integer_feasible(clex->tab); + if (!clex->tab) + return -1; + feasible = !clex->tab->empty; + if (isl_tab_rollback(clex->tab, snap) < 0) + return -1; + + return feasible; +} + +static int context_lex_get_div(struct isl_context *context, struct isl_tab *tab, + struct isl_vec *div) +{ + return get_div(tab, context, div); +} + +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, + context_lex_add_ineq_wrap, context); +} + +static int context_lex_detect_equalities(struct isl_context *context, + struct isl_tab *tab) +{ + return 0; +} + +static int context_lex_best_split(struct isl_context *context, + struct isl_tab *tab) +{ + struct isl_context_lex *clex = (struct isl_context_lex *)context; + struct isl_tab_undo *snap; + int r; + + snap = isl_tab_snap(clex->tab); + if (isl_tab_push_basis(clex->tab) < 0) + return -1; + r = best_split(tab, clex->tab); + + if (isl_tab_rollback(clex->tab, snap) < 0) + return -1; + + return r; +} + +static int context_lex_is_empty(struct isl_context *context) +{ + struct isl_context_lex *clex = (struct isl_context_lex *)context; + if (!clex->tab) + return -1; + return clex->tab->empty; +} + +static void *context_lex_save(struct isl_context *context) +{ + struct isl_context_lex *clex = (struct isl_context_lex *)context; + struct isl_tab_undo *snap; + + snap = isl_tab_snap(clex->tab); + if (isl_tab_push_basis(clex->tab) < 0) + return NULL; + if (isl_tab_save_samples(clex->tab) < 0) + return NULL; + + return snap; +} + +static void context_lex_restore(struct isl_context *context, void *save) +{ + struct isl_context_lex *clex = (struct isl_context_lex *)context; + if (isl_tab_rollback(clex->tab, (struct isl_tab_undo *)save) < 0) { + isl_tab_free(clex->tab); + clex->tab = NULL; + } +} + +static int context_lex_is_ok(struct isl_context *context) +{ + struct isl_context_lex *clex = (struct isl_context_lex *)context; + return !!clex->tab; +} + +/* For each variable in the context tableau, check if the variable can + * only attain non-negative values. If so, mark the parameter as non-negative + * in the main tableau. This allows for a more direct identification of some + * cases of violated constraints. + */ +static struct isl_tab *tab_detect_nonnegative_parameters(struct isl_tab *tab, + struct isl_tab *context_tab) +{ + int i; + struct isl_tab_undo *snap; + struct isl_vec *ineq = NULL; + struct isl_tab_var *var; + int n; + + if (context_tab->n_var == 0) + return tab; + + ineq = isl_vec_alloc(tab->mat->ctx, 1 + context_tab->n_var); + if (!ineq) + goto error; + + if (isl_tab_extend_cons(context_tab, 1) < 0) + goto error; + + snap = isl_tab_snap(context_tab); + + n = 0; + isl_seq_clr(ineq->el, ineq->size); + for (i = 0; i < context_tab->n_var; ++i) { + isl_int_set_si(ineq->el[1 + i], 1); + if (isl_tab_add_ineq(context_tab, ineq->el) < 0) + goto error; + var = &context_tab->con[context_tab->n_con - 1]; + if (!context_tab->empty && + !isl_tab_min_at_most_neg_one(context_tab, var)) { + int j = i; + if (i >= tab->n_param) + j = i - tab->n_param + tab->n_var - tab->n_div; + tab->var[j].is_nonneg = 1; + n++; + } + isl_int_set_si(ineq->el[1 + i], 0); + if (isl_tab_rollback(context_tab, snap) < 0) + goto error; + } + + if (context_tab->M && n == context_tab->n_var) { + context_tab->mat = isl_mat_drop_cols(context_tab->mat, 2, 1); + context_tab->M = 0; + } + + isl_vec_free(ineq); + return tab; +error: + isl_vec_free(ineq); + isl_tab_free(tab); + return NULL; +} + +static struct isl_tab *context_lex_detect_nonnegative_parameters( + struct isl_context *context, struct isl_tab *tab) +{ + struct isl_context_lex *clex = (struct isl_context_lex *)context; + struct isl_tab_undo *snap; + + snap = isl_tab_snap(clex->tab); + if (isl_tab_push_basis(clex->tab) < 0) + goto error; + + tab = tab_detect_nonnegative_parameters(tab, clex->tab); + + if (isl_tab_rollback(clex->tab, snap) < 0) + goto error; + + return tab; +error: + isl_tab_free(tab); + return NULL; +} + +static void context_lex_invalidate(struct isl_context *context) +{ + struct isl_context_lex *clex = (struct isl_context_lex *)context; + isl_tab_free(clex->tab); + clex->tab = NULL; +} + +static void context_lex_free(struct isl_context *context) +{ + struct isl_context_lex *clex = (struct isl_context_lex *)context; + isl_tab_free(clex->tab); + free(clex); +} + +struct isl_context_op isl_context_lex_op = { + context_lex_detect_nonnegative_parameters, + context_lex_peek_basic_set, + context_lex_peek_tab, + context_lex_add_eq, + context_lex_add_ineq, + context_lex_ineq_sign, + context_lex_test_ineq, + context_lex_get_div, + context_lex_add_div, + context_lex_detect_equalities, + context_lex_best_split, + context_lex_is_empty, + context_lex_is_ok, + context_lex_save, + context_lex_restore, + context_lex_invalidate, + context_lex_free, +}; + +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); + if (!tab) + goto error; + if (isl_tab_track_bset(tab, bset) < 0) + goto error; + tab = isl_tab_init_samples(tab); + return tab; +error: + isl_basic_set_free(bset); + return NULL; +} + +static struct isl_context *isl_context_lex_alloc(struct isl_basic_set *dom) +{ + struct isl_context_lex *clex; + + if (!dom) + return NULL; + + clex = isl_alloc_type(dom->ctx, struct isl_context_lex); + if (!clex) + return NULL; + + clex->context.op = &isl_context_lex_op; + + clex->tab = context_tab_for_lexmin(isl_basic_set_copy(dom)); + clex->tab = restore_lexmin(clex->tab); + clex->tab = check_integer_feasible(clex->tab); + if (!clex->tab) + goto error; + + return &clex->context; +error: + clex->context.op->free(&clex->context); + return NULL; +} + +struct isl_context_gbr { + struct isl_context context; + struct isl_tab *tab; + struct isl_tab *shifted; + struct isl_tab *cone; +}; + +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; + return tab_detect_nonnegative_parameters(tab, cgbr->tab); +} + +static struct isl_basic_set *context_gbr_peek_basic_set( + struct isl_context *context) +{ + struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context; + if (!cgbr->tab) + return NULL; + return isl_tab_peek_bset(cgbr->tab); +} + +static struct isl_tab *context_gbr_peek_tab(struct isl_context *context) +{ + struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context; + return cgbr->tab; +} + +/* Initialize the "shifted" tableau of the context, which + * contains the constraints of the original tableau shifted + * by the sum of all negative coefficients. This ensures + * that any rational point in the shifted tableau can + * be rounded up to yield an integer point in the original tableau. + */ +static void gbr_init_shifted(struct isl_context_gbr *cgbr) +{ + int i, j; + struct isl_vec *cst; + struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab); + unsigned dim = isl_basic_set_total_dim(bset); + + cst = isl_vec_alloc(cgbr->tab->mat->ctx, bset->n_ineq); + if (!cst) + return; + + for (i = 0; i < bset->n_ineq; ++i) { + isl_int_set(cst->el[i], bset->ineq[i][0]); + for (j = 0; j < dim; ++j) { + if (!isl_int_is_neg(bset->ineq[i][1 + j])) + continue; + isl_int_add(bset->ineq[i][0], bset->ineq[i][0], + bset->ineq[i][1 + j]); + } + } + + cgbr->shifted = isl_tab_from_basic_set(bset); + + for (i = 0; i < bset->n_ineq; ++i) + isl_int_set(bset->ineq[i][0], cst->el[i]); + + isl_vec_free(cst); +} + +/* Check if the shifted tableau is non-empty, and if so + * use the sample point to construct an integer point + * of the context tableau. + */ +static struct isl_vec *gbr_get_shifted_sample(struct isl_context_gbr *cgbr) +{ + struct isl_vec *sample; + + if (!cgbr->shifted) + gbr_init_shifted(cgbr); + if (!cgbr->shifted) + return NULL; + if (cgbr->shifted->empty) + return isl_vec_alloc(cgbr->tab->mat->ctx, 0); + + sample = isl_tab_get_sample_value(cgbr->shifted); + sample = isl_vec_ceil(sample); + + return sample; +} + +static struct isl_basic_set *drop_constant_terms(struct isl_basic_set *bset) +{ + int i; + + if (!bset) + return NULL; + + for (i = 0; i < bset->n_eq; ++i) + isl_int_set_si(bset->eq[i][0], 0); + + for (i = 0; i < bset->n_ineq; ++i) + isl_int_set_si(bset->ineq[i][0], 0); + + return bset; +} + +static int use_shifted(struct isl_context_gbr *cgbr) +{ + return cgbr->tab->bmap->n_eq == 0 && cgbr->tab->bmap->n_div == 0; +} + +static struct isl_vec *gbr_get_sample(struct isl_context_gbr *cgbr) +{ + struct isl_basic_set *bset; + struct isl_basic_set *cone; + + if (isl_tab_sample_is_integer(cgbr->tab)) + return isl_tab_get_sample_value(cgbr->tab); + + if (use_shifted(cgbr)) { + struct isl_vec *sample; + + sample = gbr_get_shifted_sample(cgbr); + if (!sample || sample->size > 0) + return sample; + + isl_vec_free(sample); + } + + if (!cgbr->cone) { + bset = isl_tab_peek_bset(cgbr->tab); + 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) + return NULL; + } + if (isl_tab_detect_implicit_equalities(cgbr->cone) < 0) + return NULL; + + if (cgbr->cone->n_dead == cgbr->cone->n_col) { + struct isl_vec *sample; + struct isl_tab_undo *snap; + + if (cgbr->tab->basis) { + if (cgbr->tab->basis->n_col != 1 + cgbr->tab->n_var) { + isl_mat_free(cgbr->tab->basis); + cgbr->tab->basis = NULL; + } + cgbr->tab->n_zero = 0; + cgbr->tab->n_unbounded = 0; + } + + snap = isl_tab_snap(cgbr->tab); + + sample = isl_tab_sample(cgbr->tab); + + if (isl_tab_rollback(cgbr->tab, snap) < 0) { + isl_vec_free(sample); + return NULL; + } + + return sample; + } + + cone = isl_basic_set_dup(isl_tab_peek_bset(cgbr->cone)); + cone = drop_constant_terms(cone); + cone = isl_basic_set_update_from_tab(cone, cgbr->cone); + cone = isl_basic_set_underlying_set(cone); + cone = isl_basic_set_gauss(cone, NULL); + + bset = isl_basic_set_dup(isl_tab_peek_bset(cgbr->tab)); + bset = isl_basic_set_update_from_tab(bset, cgbr->tab); + bset = isl_basic_set_underlying_set(bset); + bset = isl_basic_set_gauss(bset, NULL); + + return isl_basic_set_sample_with_cone(bset, cone); +} + +static void check_gbr_integer_feasible(struct isl_context_gbr *cgbr) +{ + struct isl_vec *sample; + + if (!cgbr->tab) + return; + + if (cgbr->tab->empty) + return; + + sample = gbr_get_sample(cgbr); + if (!sample) + goto error; + + if (sample->size == 0) { + isl_vec_free(sample); + if (isl_tab_mark_empty(cgbr->tab) < 0) + goto error; + return; + } + + cgbr->tab = isl_tab_add_sample(cgbr->tab, sample); + + return; +error: + isl_tab_free(cgbr->tab); + cgbr->tab = NULL; +} + +static struct isl_tab *add_gbr_eq(struct isl_tab *tab, isl_int *eq) +{ + int r; + + if (!tab) + return NULL; + + if (isl_tab_extend_cons(tab, 2) < 0) + goto error; + + if (isl_tab_add_eq(tab, eq) < 0) + goto error; + + return tab; +error: + isl_tab_free(tab); + return NULL; +} + +static void context_gbr_add_eq(struct isl_context *context, isl_int *eq, + int check, int update) +{ + struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context; + + cgbr->tab = add_gbr_eq(cgbr->tab, eq); + + if (cgbr->cone && cgbr->cone->n_col != cgbr->cone->n_dead) { + if (isl_tab_extend_cons(cgbr->cone, 2) < 0) + goto error; + if (isl_tab_add_eq(cgbr->cone, eq) < 0) + goto error; + } + + if (check) { + int v = tab_has_valid_sample(cgbr->tab, eq, 1); + if (v < 0) + goto error; + if (!v) + check_gbr_integer_feasible(cgbr); + } + if (update) + cgbr->tab = check_samples(cgbr->tab, eq, 1); + return; +error: + isl_tab_free(cgbr->tab); + cgbr->tab = NULL; +} + +static void add_gbr_ineq(struct isl_context_gbr *cgbr, isl_int *ineq) +{ + if (!cgbr->tab) + return; + + if (isl_tab_extend_cons(cgbr->tab, 1) < 0) + goto error; + + if (isl_tab_add_ineq(cgbr->tab, ineq) < 0) + goto error; + + if (cgbr->shifted && !cgbr->shifted->empty && use_shifted(cgbr)) { + int i; + unsigned dim; + dim = isl_basic_map_total_dim(cgbr->tab->bmap); + + if (isl_tab_extend_cons(cgbr->shifted, 1) < 0) + goto error; + + for (i = 0; i < dim; ++i) { + if (!isl_int_is_neg(ineq[1 + i])) + continue; + isl_int_add(ineq[0], ineq[0], ineq[1 + i]); + } + + if (isl_tab_add_ineq(cgbr->shifted, ineq) < 0) + goto error; + + for (i = 0; i < dim; ++i) { + if (!isl_int_is_neg(ineq[1 + i])) + continue; + isl_int_sub(ineq[0], ineq[0], ineq[1 + i]); + } + } + + if (cgbr->cone && cgbr->cone->n_col != cgbr->cone->n_dead) { + if (isl_tab_extend_cons(cgbr->cone, 1) < 0) + goto error; + if (isl_tab_add_ineq(cgbr->cone, ineq) < 0) + goto error; + } + + return; +error: + isl_tab_free(cgbr->tab); + cgbr->tab = NULL; +} + +static void context_gbr_add_ineq(struct isl_context *context, isl_int *ineq, + int check, int update) +{ + struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context; + + add_gbr_ineq(cgbr, ineq); + if (!cgbr->tab) + return; + + if (check) { + int v = tab_has_valid_sample(cgbr->tab, ineq, 0); + if (v < 0) + goto error; + if (!v) + check_gbr_integer_feasible(cgbr); + } + if (update) + cgbr->tab = check_samples(cgbr->tab, ineq, 0); + return; +error: + isl_tab_free(cgbr->tab); + cgbr->tab = NULL; +} + +static int context_gbr_add_ineq_wrap(void *user, isl_int *ineq) +{ + struct isl_context *context = (struct isl_context *)user; + context_gbr_add_ineq(context, ineq, 0, 0); + return context->op->is_ok(context) ? 0 : -1; +} + +static enum isl_tab_row_sign context_gbr_ineq_sign(struct isl_context *context, + isl_int *ineq, int strict) +{ + struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context; + return tab_ineq_sign(cgbr->tab, ineq, strict); +} + +/* Check whether "ineq" can be added to the tableau without rendering + * it infeasible. + */ +static int context_gbr_test_ineq(struct isl_context *context, isl_int *ineq) +{ + struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context; + struct isl_tab_undo *snap; + struct isl_tab_undo *shifted_snap = NULL; + struct isl_tab_undo *cone_snap = NULL; + int feasible; + + if (!cgbr->tab) + return -1; + + if (isl_tab_extend_cons(cgbr->tab, 1) < 0) + return -1; + + snap = isl_tab_snap(cgbr->tab); + if (cgbr->shifted) + shifted_snap = isl_tab_snap(cgbr->shifted); + if (cgbr->cone) + cone_snap = isl_tab_snap(cgbr->cone); + add_gbr_ineq(cgbr, ineq); + check_gbr_integer_feasible(cgbr); + if (!cgbr->tab) + return -1; + feasible = !cgbr->tab->empty; + if (isl_tab_rollback(cgbr->tab, snap) < 0) + return -1; + if (shifted_snap) { + if (isl_tab_rollback(cgbr->shifted, shifted_snap)) + return -1; + } else if (cgbr->shifted) { + isl_tab_free(cgbr->shifted); + cgbr->shifted = NULL; + } + if (cone_snap) { + if (isl_tab_rollback(cgbr->cone, cone_snap)) + return -1; + } else if (cgbr->cone) { + isl_tab_free(cgbr->cone); + cgbr->cone = NULL; + } + + return feasible; +} + +/* Return the column of the last of the variables associated to + * a column that has a non-zero coefficient. + * This function is called in a context where only coefficients + * of parameters or divs can be non-zero. + */ +static int last_non_zero_var_col(struct isl_tab *tab, isl_int *p) +{ + int i; + int col; + unsigned dim = tab->n_var - tab->n_param - tab->n_div; + + if (tab->n_var == 0) + return -1; + + for (i = tab->n_var - 1; i >= 0; --i) { + if (i >= tab->n_param && i < tab->n_var - tab->n_div) + continue; + if (tab->var[i].is_row) + continue; + col = tab->var[i].index; + if (!isl_int_is_zero(p[col])) + return col; + } + + return -1; +} + +/* Look through all the recently added equalities in the context + * to see if we can propagate any of them to the main tableau. + * + * The newly added equalities in the context are encoded as pairs + * of inequalities starting at inequality "first". + * + * We tentatively add each of these equalities to the main tableau + * and if this happens to result in a row with a final coefficient + * that is one or negative one, we use it to kill a column + * in the main tableau. Otherwise, we discard the tentatively + * added row. + */ +static void propagate_equalities(struct isl_context_gbr *cgbr, + struct isl_tab *tab, unsigned first) +{ + int i; + struct isl_vec *eq = NULL; + + eq = isl_vec_alloc(tab->mat->ctx, 1 + tab->n_var); + if (!eq) + goto error; + + if (isl_tab_extend_cons(tab, (cgbr->tab->bmap->n_ineq - first)/2) < 0) + goto error; + + isl_seq_clr(eq->el + 1 + tab->n_param, + tab->n_var - tab->n_param - tab->n_div); + for (i = first; i < cgbr->tab->bmap->n_ineq; i += 2) { + int j; + int r; + struct isl_tab_undo *snap; + snap = isl_tab_snap(tab); + + isl_seq_cpy(eq->el, cgbr->tab->bmap->ineq[i], 1 + tab->n_param); + isl_seq_cpy(eq->el + 1 + tab->n_var - tab->n_div, + cgbr->tab->bmap->ineq[i] + 1 + tab->n_param, + tab->n_div); + + r = isl_tab_add_row(tab, eq->el); + if (r < 0) + goto error; + r = tab->con[r].index; + j = last_non_zero_var_col(tab, tab->mat->row[r] + 2 + tab->M); + if (j < 0 || j < tab->n_dead || + !isl_int_is_one(tab->mat->row[r][0]) || + (!isl_int_is_one(tab->mat->row[r][2 + tab->M + j]) && + !isl_int_is_negone(tab->mat->row[r][2 + tab->M + j]))) { + if (isl_tab_rollback(tab, snap) < 0) + goto error; + continue; + } + if (isl_tab_pivot(tab, r, j) < 0) + goto error; + if (isl_tab_kill_col(tab, j) < 0) + goto error; + + tab = restore_lexmin(tab); + } + + isl_vec_free(eq); + + return; +error: + isl_vec_free(eq); + isl_tab_free(cgbr->tab); + cgbr->tab = NULL; +} + +static int context_gbr_detect_equalities(struct isl_context *context, + struct isl_tab *tab) +{ + struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context; + struct isl_ctx *ctx; + int i; + enum isl_lp_result res; + unsigned n_ineq; + + ctx = cgbr->tab->mat->ctx; + + if (!cgbr->cone) { + struct isl_basic_set *bset = isl_tab_peek_bset(cgbr->tab); + 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) + goto error; + } + if (isl_tab_detect_implicit_equalities(cgbr->cone) < 0) + goto error; + + 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) + propagate_equalities(cgbr, tab, n_ineq); + + return 0; +error: + isl_tab_free(cgbr->tab); + cgbr->tab = NULL; + return -1; +} + +static int context_gbr_get_div(struct isl_context *context, struct isl_tab *tab, + struct isl_vec *div) +{ + return get_div(tab, context, div); +} + +static int context_gbr_add_div(struct isl_context *context, struct isl_vec *div) +{ + struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context; + if (cgbr->cone) { + int k; + + if (isl_tab_extend_cons(cgbr->cone, 3) < 0) + return -1; + if (isl_tab_extend_vars(cgbr->cone, 1) < 0) + return -1; + 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); + k = isl_basic_map_alloc_div(cgbr->cone->bmap); + if (k < 0) + return -1; + isl_seq_cpy(cgbr->cone->bmap->div[k], div->el, div->size); + if (isl_tab_push(cgbr->cone, isl_tab_undo_bmap_div) < 0) + return -1; + } + return context_tab_add_div(cgbr->tab, div, + context_gbr_add_ineq_wrap, context); +} + +static int context_gbr_best_split(struct isl_context *context, + struct isl_tab *tab) +{ + struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context; + struct isl_tab_undo *snap; + int r; + + snap = isl_tab_snap(cgbr->tab); + r = best_split(tab, cgbr->tab); + + if (isl_tab_rollback(cgbr->tab, snap) < 0) + return -1; + + return r; +} + +static int context_gbr_is_empty(struct isl_context *context) +{ + struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context; + if (!cgbr->tab) + return -1; + return cgbr->tab->empty; +} + +struct isl_gbr_tab_undo { + struct isl_tab_undo *tab_snap; + struct isl_tab_undo *shifted_snap; + struct isl_tab_undo *cone_snap; +}; + +static void *context_gbr_save(struct isl_context *context) +{ + struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context; + struct isl_gbr_tab_undo *snap; + + snap = isl_alloc_type(cgbr->tab->mat->ctx, struct isl_gbr_tab_undo); + if (!snap) + return NULL; + + snap->tab_snap = isl_tab_snap(cgbr->tab); + if (isl_tab_save_samples(cgbr->tab) < 0) + goto error; + + if (cgbr->shifted) + snap->shifted_snap = isl_tab_snap(cgbr->shifted); + else + snap->shifted_snap = NULL; + + if (cgbr->cone) + snap->cone_snap = isl_tab_snap(cgbr->cone); + else + snap->cone_snap = NULL; + + return snap; +error: + free(snap); + return NULL; +} + +static void context_gbr_restore(struct isl_context *context, void *save) +{ + struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context; + struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save; + if (!snap) + goto error; + if (isl_tab_rollback(cgbr->tab, snap->tab_snap) < 0) { + isl_tab_free(cgbr->tab); + cgbr->tab = NULL; + } + + if (snap->shifted_snap) { + if (isl_tab_rollback(cgbr->shifted, snap->shifted_snap) < 0) + goto error; + } else if (cgbr->shifted) { + isl_tab_free(cgbr->shifted); + cgbr->shifted = NULL; + } + + if (snap->cone_snap) { + if (isl_tab_rollback(cgbr->cone, snap->cone_snap) < 0) + goto error; + } else if (cgbr->cone) { + isl_tab_free(cgbr->cone); + cgbr->cone = NULL; + } + + free(snap); + + return; +error: + free(snap); + isl_tab_free(cgbr->tab); + cgbr->tab = NULL; +} + +static int context_gbr_is_ok(struct isl_context *context) +{ + struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context; + return !!cgbr->tab; +} + +static void context_gbr_invalidate(struct isl_context *context) +{ + struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context; + isl_tab_free(cgbr->tab); + cgbr->tab = NULL; +} + +static void context_gbr_free(struct isl_context *context) +{ + struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context; + isl_tab_free(cgbr->tab); + isl_tab_free(cgbr->shifted); + isl_tab_free(cgbr->cone); + free(cgbr); +} + +struct isl_context_op isl_context_gbr_op = { + context_gbr_detect_nonnegative_parameters, + context_gbr_peek_basic_set, + context_gbr_peek_tab, + context_gbr_add_eq, + context_gbr_add_ineq, + context_gbr_ineq_sign, + context_gbr_test_ineq, + context_gbr_get_div, + context_gbr_add_div, + context_gbr_detect_equalities, + context_gbr_best_split, + context_gbr_is_empty, + context_gbr_is_ok, + context_gbr_save, + context_gbr_restore, + context_gbr_invalidate, + context_gbr_free, +}; + +static struct isl_context *isl_context_gbr_alloc(struct isl_basic_set *dom) +{ + struct isl_context_gbr *cgbr; + + if (!dom) + return NULL; + + cgbr = isl_calloc_type(dom->ctx, struct isl_context_gbr); + if (!cgbr) + return NULL; + + cgbr->context.op = &isl_context_gbr_op; + + cgbr->shifted = NULL; + cgbr->cone = NULL; + cgbr->tab = isl_tab_from_basic_set(dom); + 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; error: - isl_basic_set_free(bset); + cgbr->context.op->free(&cgbr->context); return NULL; } +static struct isl_context *isl_context_alloc(struct isl_basic_set *dom) +{ + if (!dom) + return NULL; + + if (dom->ctx->opt->context == ISL_CONTEXT_LEXMIN) + return isl_context_lex_alloc(dom); + else + return isl_context_gbr_alloc(dom); +} + /* Construct an isl_sol_map 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. @@ -1684,27 +3212,30 @@ error: 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_tab *context_tab; - int f; + struct isl_sol_map *sol_map = NULL; - sol_map = isl_calloc_type(bset->ctx, struct isl_sol_map); + if (!bmap) + goto error; + + sol_map = isl_calloc_type(bmap->ctx, struct isl_sol_map); if (!sol_map) goto error; - sol_map->max = max; + sol_map->sol.rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL); + sol_map->sol.dec_level.callback.run = &sol_dec_level_wrap; + sol_map->sol.dec_level.sol = &sol_map->sol; + sol_map->sol.max = max; + sol_map->sol.n_out = isl_basic_map_dim(bmap, isl_dim_out); 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, ISL_MAP_DISJOINT); if (!sol_map->map) goto error; - context_tab = context_tab_for_lexmin(isl_basic_set_copy(dom)); - context_tab = restore_lexmin(context_tab); - sol_map->sol.context_tab = context_tab; - f = context_is_feasible(&sol_map->sol); - if (f < 0) + sol_map->sol.context = isl_context_alloc(dom); + if (!sol_map->sol.context) goto error; if (track_empty) { @@ -1722,70 +3253,6 @@ error: return NULL; } -/* For each variable in the context tableau, check if the variable can - * only attain non-negative values. If so, mark the parameter as non-negative - * in the main tableau. This allows for a more direct identification of some - * cases of violated constraints. - */ -static struct isl_tab *tab_detect_nonnegative_parameters(struct isl_tab *tab, - struct isl_tab *context_tab) -{ - int i; - struct isl_tab_undo *snap, *snap2; - struct isl_vec *ineq = NULL; - struct isl_tab_var *var; - int n; - - if (context_tab->n_var == 0) - return tab; - - ineq = isl_vec_alloc(tab->mat->ctx, 1 + context_tab->n_var); - if (!ineq) - goto error; - - if (isl_tab_extend_cons(context_tab, 1) < 0) - goto error; - - snap = isl_tab_snap(context_tab); - isl_tab_push_basis(context_tab); - - snap2 = isl_tab_snap(context_tab); - - n = 0; - isl_seq_clr(ineq->el, ineq->size); - for (i = 0; i < context_tab->n_var; ++i) { - isl_int_set_si(ineq->el[1 + i], 1); - context_tab = isl_tab_add_ineq(context_tab, ineq->el); - var = &context_tab->con[context_tab->n_con - 1]; - if (!context_tab->empty && - !isl_tab_min_at_most_neg_one(context_tab, var)) { - int j = i; - if (i >= tab->n_param) - j = i - tab->n_param + tab->n_var - tab->n_div; - tab->var[j].is_nonneg = 1; - n++; - } - isl_int_set_si(ineq->el[1 + i], 0); - if (isl_tab_rollback(context_tab, snap2) < 0) - goto error; - } - - if (isl_tab_rollback(context_tab, snap) < 0) - goto error; - - if (n == context_tab->n_var) { - context_tab->mat = isl_mat_drop_cols(context_tab->mat, 2, 1); - context_tab->M = 0; - } - - isl_vec_free(ineq); - return tab; -error: - isl_vec_free(ineq); - isl_tab_free(tab); - return NULL; -} - /* Check whether all coefficients of (non-parameter) variables * are non-positive, meaning that no pivots can be performed on the row. */ @@ -1886,18 +3353,14 @@ static int is_strict(struct isl_vec *vec) * >=0 ? Y N * any neg */ -static int row_sign(struct isl_tab *tab, struct isl_sol *sol, int row) +static enum isl_tab_row_sign row_sign(struct isl_tab *tab, + struct isl_sol *sol, int row) { - int i; - struct isl_tab_undo *snap = NULL; struct isl_vec *ineq = NULL; - int res = isl_tab_row_unknown; + enum isl_tab_row_sign res = isl_tab_row_unknown; int critical; int strict; - int sgn; int row2; - isl_int tmp; - struct isl_tab *context_tab = sol->context_tab; if (tab->row_sign[row] != isl_tab_row_unknown) return tab->row_sign[row]; @@ -1912,44 +3375,14 @@ static int row_sign(struct isl_tab *tab, struct isl_sol *sol, int row) critical = is_critical(tab, row); - isl_assert(tab->mat->ctx, context_tab->samples, goto error); - isl_assert(tab->mat->ctx, context_tab->samples->n_col == 1 + context_tab->n_var, goto error); - ineq = get_row_parameter_ineq(tab, row); if (!ineq) goto error; strict = is_strict(ineq); - isl_int_init(tmp); - for (i = context_tab->n_outside; i < context_tab->n_sample; ++i) { - isl_seq_inner_product(context_tab->samples->row[i], ineq->el, - ineq->size, &tmp); - sgn = isl_int_sgn(tmp); - if (sgn > 0 || (sgn == 0 && (critical || strict))) { - if (res == isl_tab_row_unknown) - res = isl_tab_row_pos; - if (res == isl_tab_row_neg) - res = isl_tab_row_any; - } - if (sgn < 0) { - if (res == isl_tab_row_unknown) - res = isl_tab_row_neg; - if (res == isl_tab_row_pos) - res = isl_tab_row_any; - } - if (res == isl_tab_row_any) - break; - } - isl_int_clear(tmp); - - if (res != isl_tab_row_any) { - if (isl_tab_extend_cons(context_tab, 1) < 0) - goto error; - - snap = isl_tab_snap(context_tab); - isl_tab_push_basis(context_tab); - } + res = sol->context->op->ineq_sign(sol->context, ineq->el, + critical || strict); if (res == isl_tab_row_unknown || res == isl_tab_row_pos) { /* test for negative values */ @@ -1957,20 +3390,14 @@ static int row_sign(struct isl_tab *tab, struct isl_sol *sol, int row) isl_seq_neg(ineq->el, ineq->el, ineq->size); isl_int_sub_ui(ineq->el[0], ineq->el[0], 1); - isl_tab_push_basis(context_tab); - sol->context_tab = add_lexmin_ineq(sol->context_tab, ineq->el); - feasible = context_is_feasible(sol); + feasible = sol->context->op->test_ineq(sol->context, ineq->el); if (feasible < 0) goto error; - context_tab = sol->context_tab; if (!feasible) res = isl_tab_row_pos; else res = (res == isl_tab_row_unknown) ? isl_tab_row_neg : isl_tab_row_any; - if (isl_tab_rollback(context_tab, snap) < 0) - goto error; - if (res == isl_tab_row_neg) { isl_seq_neg(ineq->el, ineq->el, ineq->size); isl_int_sub_ui(ineq->el[0], ineq->el[0], 1); @@ -1983,26 +3410,21 @@ static int row_sign(struct isl_tab *tab, struct isl_sol *sol, int row) if (!critical && !strict) isl_int_sub_ui(ineq->el[0], ineq->el[0], 1); - isl_tab_push_basis(context_tab); - sol->context_tab = add_lexmin_ineq(sol->context_tab, ineq->el); - feasible = context_is_feasible(sol); + feasible = sol->context->op->test_ineq(sol->context, ineq->el); if (feasible < 0) goto error; - context_tab = sol->context_tab; if (feasible) res = isl_tab_row_any; - if (isl_tab_rollback(context_tab, snap) < 0) - goto error; } isl_vec_free(ineq); return res; error: isl_vec_free(ineq); - return 0; + return isl_tab_row_unknown; } -static struct isl_sol *find_solutions(struct isl_sol *sol, struct isl_tab *tab); +static void find_solutions(struct isl_sol *sol, struct isl_tab *tab); /* Find solutions for values of the parameters that satisfy the given * inequality. @@ -2018,152 +3440,58 @@ static struct isl_sol *find_solutions(struct isl_sol *sol, struct isl_tab *tab); * and that we need to do this before saving the current basis * such that the basis has been restore before we restore the row signs. */ -static struct isl_sol *find_in_pos(struct isl_sol *sol, - struct isl_tab *tab, isl_int *ineq) +static void find_in_pos(struct isl_sol *sol, struct isl_tab *tab, isl_int *ineq) { - struct isl_tab_undo *snap; + void *saved; - snap = isl_tab_snap(sol->context_tab); - isl_tab_push_basis(sol->context_tab); - if (isl_tab_extend_cons(sol->context_tab, 1) < 0) + if (!sol->context) goto error; + saved = sol->context->op->save(sol->context); tab = isl_tab_dup(tab); if (!tab) goto error; - sol->context_tab = add_lexmin_ineq(sol->context_tab, ineq); - sol->context_tab = check_samples(sol->context_tab, ineq, 0); + sol->context->op->add_ineq(sol->context, ineq, 0, 1); - sol = find_solutions(sol, tab); + find_solutions(sol, tab); - isl_tab_rollback(sol->context_tab, snap); - return sol; + sol->context->op->restore(sol->context, saved); + return; error: - isl_tab_rollback(sol->context_tab, snap); - sol_free(sol); - return NULL; + sol->error = 1; } /* Record the absence of solutions for those values of the parameters * that do not satisfy the given inequality with equality. */ -static struct isl_sol *no_sol_in_strict(struct isl_sol *sol, +static void no_sol_in_strict(struct isl_sol *sol, struct isl_tab *tab, struct isl_vec *ineq) { int empty; - int f; - struct isl_tab_undo *snap; - snap = isl_tab_snap(sol->context_tab); - isl_tab_push_basis(sol->context_tab); - if (isl_tab_extend_cons(sol->context_tab, 1) < 0) + void *saved; + + if (!sol->context) goto error; + saved = sol->context->op->save(sol->context); isl_int_sub_ui(ineq->el[0], ineq->el[0], 1); - sol->context_tab = add_lexmin_ineq(sol->context_tab, ineq->el); - f = context_valid_sample_or_feasible(sol, ineq->el, 0); - if (f < 0) + sol->context->op->add_ineq(sol->context, ineq->el, 1, 0); + if (!sol->context) goto error; empty = tab->empty; tab->empty = 1; - sol = sol->add(sol, tab); + sol_add(sol, tab); tab->empty = empty; isl_int_add_ui(ineq->el[0], ineq->el[0], 1); - if (isl_tab_rollback(sol->context_tab, snap) < 0) - goto error; - return sol; + sol->context->op->restore(sol->context, saved); + return; error: - sol_free(sol); - return NULL; -} - -/* Given a main tableau where more than one row requires a split, - * determine and return the "best" row to split on. - * - * Given two rows in the main tableau, if the inequality corresponding - * to the first row is redundant with respect to that of the second row - * in the current tableau, then it is better to split on the second row, - * since in the positive part, both row will be positive. - * (In the negative part a pivot will have to be performed and just about - * anything can happen to the sign of the other row.) - * - * As a simple heuristic, we therefore select the row that makes the most - * of the other rows redundant. - * - * Perhaps it would also be useful to look at the number of constraints - * that conflict with any given constraint. - */ -static int best_split(struct isl_tab *tab, struct isl_tab *context_tab) -{ - struct isl_tab_undo *snap, *snap2; - int split; - int row; - int best = -1; - int best_r; - - if (isl_tab_extend_cons(context_tab, 2) < 0) - return -1; - - snap = isl_tab_snap(context_tab); - isl_tab_push_basis(context_tab); - snap2 = isl_tab_snap(context_tab); - - for (split = tab->n_redundant; split < tab->n_row; ++split) { - struct isl_tab_undo *snap3; - struct isl_vec *ineq = NULL; - int r = 0; - - if (!isl_tab_var_from_row(tab, split)->is_nonneg) - continue; - if (tab->row_sign[split] != isl_tab_row_any) - continue; - - ineq = get_row_parameter_ineq(tab, split); - if (!ineq) - return -1; - context_tab = isl_tab_add_ineq(context_tab, ineq->el); - isl_vec_free(ineq); - - snap3 = isl_tab_snap(context_tab); - - for (row = tab->n_redundant; row < tab->n_row; ++row) { - struct isl_tab_var *var; - - if (row == split) - continue; - if (!isl_tab_var_from_row(tab, row)->is_nonneg) - continue; - if (tab->row_sign[row] != isl_tab_row_any) - continue; - - ineq = get_row_parameter_ineq(tab, row); - if (!ineq) - return -1; - context_tab = isl_tab_add_ineq(context_tab, ineq->el); - isl_vec_free(ineq); - var = &context_tab->con[context_tab->n_con - 1]; - if (!context_tab->empty && - !isl_tab_min_at_most_neg_one(context_tab, var)) - r++; - if (isl_tab_rollback(context_tab, snap3) < 0) - return -1; - } - if (best == -1 || r > best_r) { - best = split; - best_r = r; - } - if (isl_tab_rollback(context_tab, snap2) < 0) - return -1; - } - - if (isl_tab_rollback(context_tab, snap) < 0) - return -1; - - return best; + sol->error = 1; } /* Compute the lexicographic minimum of the set represented by the main @@ -2260,24 +3588,24 @@ static int best_split(struct isl_tab *tab, struct isl_tab *context_tab) * In the part of the context where this inequality does not hold, the * main tableau is marked as being empty. */ -static struct isl_sol *find_solutions(struct isl_sol *sol, struct isl_tab *tab) +static void find_solutions(struct isl_sol *sol, struct isl_tab *tab) { - struct isl_tab **context_tab; + struct isl_context *context; - if (!tab || !sol) + if (!tab || sol->error) goto error; - context_tab = &sol->context_tab; + context = sol->context; if (tab->empty) goto done; - if ((*context_tab)->empty) + if (context->op->is_empty(context)) goto done; for (; tab && !tab->empty; tab = restore_lexmin(tab)) { int flags; int row; - int sgn; + enum isl_tab_row_sign sgn; int split = -1; int n_split = 0; @@ -2300,7 +3628,7 @@ static struct isl_sol *find_solutions(struct isl_sol *sol, struct isl_tab *tab) if (split != -1) { struct isl_vec *ineq; if (n_split != 1) - split = best_split(tab, *context_tab); + split = context->op->best_split(context, tab); if (split < 0) goto error; ineq = get_row_parameter_ineq(tab, split); @@ -2314,26 +3642,27 @@ static struct isl_sol *find_solutions(struct isl_sol *sol, struct isl_tab *tab) tab->row_sign[row] = isl_tab_row_unknown; } tab->row_sign[split] = isl_tab_row_pos; - sol = find_in_pos(sol, tab, ineq->el); + sol_inc_level(sol); + find_in_pos(sol, tab, ineq->el); tab->row_sign[split] = isl_tab_row_neg; row = split; isl_seq_neg(ineq->el, ineq->el, ineq->size); isl_int_sub_ui(ineq->el[0], ineq->el[0], 1); - *context_tab = add_lexmin_ineq(*context_tab, ineq->el); - *context_tab = check_samples(*context_tab, ineq->el, 0); + context->op->add_ineq(context, ineq->el, 0, 1); isl_vec_free(ineq); - if (!sol) + if (sol->error) goto error; continue; } if (tab->rational) break; - row = first_non_integer(tab, &flags); + row = first_non_integer_row(tab, &flags); if (row < 0) break; if (ISL_FL_ISSET(flags, I_PAR)) { if (ISL_FL_ISSET(flags, I_VAR)) { - tab = isl_tab_mark_empty(tab); + if (isl_tab_mark_empty(tab) < 0) + goto error; break; } row = add_cut(tab, row); @@ -2341,37 +3670,36 @@ static struct isl_sol *find_solutions(struct isl_sol *sol, struct isl_tab *tab) struct isl_vec *div; struct isl_vec *ineq; int d; - if (isl_tab_extend_cons(*context_tab, 3) < 0) - goto error; div = get_row_split_div(tab, row); if (!div) goto error; - d = get_div(tab, context_tab, div); + d = context->op->get_div(context, tab, div); isl_vec_free(div); if (d < 0) goto error; - ineq = ineq_for_div((*context_tab)->bset, d); - sol = no_sol_in_strict(sol, tab, ineq); + ineq = ineq_for_div(context->op->peek_basic_set(context), d); + sol_inc_level(sol); + no_sol_in_strict(sol, tab, ineq); isl_seq_neg(ineq->el, ineq->el, ineq->size); - *context_tab = add_lexmin_ineq(*context_tab, ineq->el); - *context_tab = check_samples(*context_tab, ineq->el, 0); + context->op->add_ineq(context, ineq->el, 1, 1); isl_vec_free(ineq); - if (!sol) + if (sol->error || !context->op->is_ok(context)) goto error; tab = set_row_cst_to_div(tab, row, d); + if (context->op->is_empty(context)) + break; } else - row = add_parametric_cut(tab, row, context_tab); + row = add_parametric_cut(tab, row, context); if (row < 0) goto error; } done: - sol = sol->add(sol, tab); + sol_add(sol, tab); isl_tab_free(tab); - return sol; + return; error: isl_tab_free(tab); sol_free(sol); - return NULL; } /* Compute the lexicographic minimum of the set represented by the main @@ -2385,11 +3713,12 @@ error: * In parts of the context where the added equality does not hold, * the main tableau is marked as being empty. */ -static struct isl_sol *find_solutions_main(struct isl_sol *sol, - struct isl_tab *tab) +static void find_solutions_main(struct isl_sol *sol, struct isl_tab *tab) { int row; + sol->level = 0; + for (row = tab->n_redundant; row < tab->n_row; ++row) { int p; struct isl_vec *eq; @@ -2405,47 +3734,47 @@ static struct isl_sol *find_solutions_main(struct isl_sol *sol, p = tab->row_var[row] + tab->n_param - (tab->n_var - tab->n_div); - if (isl_tab_extend_cons(sol->context_tab, 2) < 0) - goto error; - eq = isl_vec_alloc(tab->mat->ctx, 1+tab->n_param+tab->n_div); 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); - sol = no_sol_in_strict(sol, tab, eq); + sol_inc_level(sol); + no_sol_in_strict(sol, tab, eq); isl_seq_neg(eq->el, eq->el, eq->size); - sol = no_sol_in_strict(sol, tab, eq); + sol_inc_level(sol); + no_sol_in_strict(sol, tab, eq); isl_seq_neg(eq->el, eq->el, eq->size); - sol->context_tab = add_lexmin_eq(sol->context_tab, eq->el); - context_valid_sample_or_feasible(sol, eq->el, 1); - sol->context_tab = check_samples(sol->context_tab, eq->el, 1); + sol->context->op->add_eq(sol->context, eq->el, 1, 1); isl_vec_free(eq); - isl_tab_mark_redundant(tab, row); - - if (!sol->context_tab) + if (isl_tab_mark_redundant(tab, row) < 0) goto error; - if (sol->context_tab->empty) + + if (sol->context->op->is_empty(sol->context)) break; row = tab->n_redundant - 1; } - return find_solutions(sol, tab); + find_solutions(sol, tab); + + sol->level = 0; + sol_pop(sol); + + return; error: isl_tab_free(tab); sol_free(sol); - return NULL; } -static struct isl_sol_map *sol_map_find_solutions(struct isl_sol_map *sol_map, +static void sol_map_find_solutions(struct isl_sol_map *sol_map, struct isl_tab *tab) { - return (struct isl_sol_map *)find_solutions_main(&sol_map->sol, tab); + find_solutions_main(&sol_map->sol, tab); } /* Check if integer division "div" of "dom" also occurs in "bmap". @@ -2546,6 +3875,8 @@ struct isl_map *isl_tab_basic_map_partial_lexopt( 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; @@ -2555,7 +3886,10 @@ struct isl_map *isl_tab_basic_map_partial_lexopt( isl_assert(bmap->ctx, isl_basic_map_compatible_domain(bmap, dom), goto error); - bmap = isl_basic_map_detect_equalities(bmap); + 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); @@ -2565,28 +3899,29 @@ struct isl_map *isl_tab_basic_map_partial_lexopt( if (!sol_map) goto error; - if (isl_basic_set_fast_is_empty(sol_map->sol.context_tab->bset)) + context = sol_map->sol.context; + if (isl_basic_set_fast_is_empty(context->op->peek_basic_set(context))) /* nothing */; else if (isl_basic_map_fast_is_empty(bmap)) - sol_map = add_empty(sol_map); + sol_map_add_empty_if_needed(sol_map, + isl_basic_set_copy(context->op->peek_basic_set(context))); else { tab = tab_for_lexmin(bmap, - sol_map->sol.context_tab->bset, 1, max); - tab = tab_detect_nonnegative_parameters(tab, - sol_map->sol.context_tab); - sol_map = sol_map_find_solutions(sol_map, tab); - if (!sol_map) - goto error; + context->op->peek_basic_set(context), 1, max); + tab = context->op->detect_nonnegative_parameters(context, tab); + sol_map_find_solutions(sol_map, tab); } + if (sol_map->sol.error) + goto error; result = isl_map_copy(sol_map->map); if (empty) *empty = isl_set_copy(sol_map->empty); - sol_map_free(sol_map); + sol_free(&sol_map->sol); isl_basic_map_free(bmap); return result; error: - sol_map_free(sol_map); + sol_free(&sol_map->sol); isl_basic_map_free(bmap); return NULL; } @@ -2596,12 +3931,12 @@ struct isl_sol_for { int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_mat *map, void *user); void *user; - int max; }; static void sol_for_free(struct isl_sol_for *sol_for) { - isl_tab_free(sol_for->sol.context_tab); + if (sol_for->sol.context) + sol_for->sol.context->op->free(sol_for->sol.context); free(sol_for); } @@ -2612,7 +3947,7 @@ static void sol_for_free_wrap(struct isl_sol *sol) /* Add the solution identified by the tableau and the context tableau. * - * See documentation of sol_map_add for more details. + * See documentation of sol_add for more details. * * Instead of constructing a basic map, this function calls a user * defined function with the current context as a basic set and @@ -2624,89 +3959,31 @@ static void sol_for_free_wrap(struct isl_sol *sol) * may refer to the divs, the basic set is not simplified. * (Simplification may reorder or remove divs.) */ -static struct isl_sol_for *sol_for_add(struct isl_sol_for *sol, - struct isl_tab *tab) +static void sol_for_add(struct isl_sol_for *sol, + struct isl_basic_set *dom, struct isl_mat *M) { - struct isl_tab *context_tab; - struct isl_basic_set *bset; - struct isl_mat *mat = NULL; - unsigned n_out; - unsigned off; - int row, i; - - if (!sol || !tab) - goto error; - - if (tab->empty) - return sol; - - off = 2 + tab->M; - context_tab = sol->sol.context_tab; - - n_out = tab->n_var - tab->n_param - tab->n_div; - mat = isl_mat_alloc(tab->mat->ctx, 1 + n_out, 1 + tab->n_param + tab->n_div); - if (!mat) + if (sol->sol.error || !dom || !M) goto error; - isl_seq_clr(mat->row[0] + 1, mat->n_col - 1); - isl_int_set_si(mat->row[0][0], 1); - for (row = 0; row < n_out; ++row) { - int i = tab->n_param + row; - int r, j; - - isl_seq_clr(mat->row[1 + row], mat->n_col); - if (!tab->var[i].is_row) - 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 error); - isl_int_set(mat->row[1 + row][0], tab->mat->row[r][1]); - for (j = 0; j < tab->n_param; ++j) { - int col; - if (tab->var[j].is_row) - continue; - col = tab->var[j].index; - isl_int_set(mat->row[1 + row][1 + j], - tab->mat->row[r][off + col]); - } - for (j = 0; j < tab->n_div; ++j) { - int col; - if (tab->var[tab->n_var - tab->n_div+j].is_row) - continue; - col = tab->var[tab->n_var - tab->n_div+j].index; - isl_int_set(mat->row[1 + row][1 + tab->n_param + j], - tab->mat->row[r][off + col]); - } - if (!isl_int_is_one(tab->mat->row[r][0])) - isl_seq_scale_down(mat->row[1 + row], mat->row[1 + row], - tab->mat->row[r][0], mat->n_col); - if (sol->max) - isl_seq_neg(mat->row[1 + row], mat->row[1 + row], - mat->n_col); - } - - bset = isl_basic_set_dup(context_tab->bset); - bset = isl_basic_set_finalize(bset); + dom = isl_basic_set_simplify(dom); + dom = isl_basic_set_finalize(dom); - if (sol->fn(bset, isl_mat_copy(mat), sol->user) < 0) + if (sol->fn(isl_basic_set_copy(dom), isl_mat_copy(M), sol->user) < 0) goto error; - isl_mat_free(mat); - return sol; + isl_basic_set_free(dom); + isl_mat_free(M); + return; error: - isl_mat_free(mat); - sol_free(&sol->sol); - return NULL; + isl_basic_set_free(dom); + isl_mat_free(M); + sol->sol.error = 1; } -static struct isl_sol *sol_for_add_wrap(struct isl_sol *sol, - struct isl_tab *tab) +static void sol_for_add_wrap(struct isl_sol *sol, + struct isl_basic_set *dom, struct isl_mat *M) { - return (struct isl_sol *)sol_for_add((struct isl_sol_for *)sol, tab); + 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, @@ -2717,27 +3994,27 @@ static struct isl_sol_for *sol_for_init(struct isl_basic_map *bmap, int max, struct isl_sol_for *sol_for = NULL; struct isl_dim *dom_dim; struct isl_basic_set *dom = NULL; - struct isl_tab *context_tab; - int f; - sol_for = isl_calloc_type(bset->ctx, struct isl_sol_for); + 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 = 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->max = max; + 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; - context_tab = context_tab_for_lexmin(isl_basic_set_copy(dom)); - context_tab = restore_lexmin(context_tab); - sol_for->sol.context_tab = context_tab; - f = context_is_feasible(&sol_for->sol); - if (f < 0) + sol_for->sol.context = isl_context_alloc(dom); + if (!sol_for->sol.context) goto error; isl_basic_set_free(dom); @@ -2748,10 +4025,10 @@ error: return NULL; } -static struct isl_sol_for *sol_for_find_solutions(struct isl_sol_for *sol_for, +static void sol_for_find_solutions(struct isl_sol_for *sol_for, struct isl_tab *tab) { - return (struct isl_sol_for *)find_solutions_main(&sol_for->sol, tab); + find_solutions_main(&sol_for->sol, tab); } int isl_basic_map_foreach_lexopt(__isl_keep isl_basic_map *bmap, int max, @@ -2772,20 +4049,20 @@ int isl_basic_map_foreach_lexopt(__isl_keep isl_basic_map *bmap, int max, /* nothing */; else { struct isl_tab *tab; + struct isl_context *context = sol_for->sol.context; tab = tab_for_lexmin(bmap, - sol_for->sol.context_tab->bset, 1, max); - tab = tab_detect_nonnegative_parameters(tab, - sol_for->sol.context_tab); - sol_for = sol_for_find_solutions(sol_for, tab); - if (!sol_for) + 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_for_free(sol_for); + sol_free(&sol_for->sol); isl_basic_map_free(bmap); return 0; error: - sol_for_free(sol_for); + sol_free(&sol_for->sol); isl_basic_map_free(bmap); return -1; }