#include "isl_tab.h"
#include "isl_sample.h"
#include <isl_mat_private.h>
+#include <isl_aff_private.h>
+#include <isl_config.h>
/*
* The implementation of parametric integer linear programming in this file
struct isl_basic_set *bset = NULL;
struct isl_mat *mat = NULL;
unsigned off;
- int row, i;
+ int row;
isl_int m;
if (sol->error || !tab)
{
int i;
struct isl_basic_map *bmap = NULL;
- isl_basic_set *context_bset;
unsigned n_eq;
unsigned n_ineq;
unsigned nparam;
* is satisfied. This function is not called from the current code
* but is useful during debugging.
*/
+static void check_lexpos(struct isl_tab *tab) __attribute__ ((unused));
static void check_lexpos(struct isl_tab *tab)
{
unsigned off = 2 + tab->M;
}
}
+/* Report to the caller that the given constraint is part of an encountered
+ * conflict.
+ */
+static int report_conflicting_constraint(struct isl_tab *tab, int con)
+{
+ return tab->conflict(con, tab->conflict_user);
+}
+
+/* Given a conflicting row in the tableau, report all constraints
+ * involved in the row to the caller. That is, the row itself
+ * (if represents a constraint) and all constraint columns with
+ * non-zero (and therefore negative) coefficient.
+ */
+static int report_conflict(struct isl_tab *tab, int row)
+{
+ int j;
+ isl_int *tr;
+
+ if (!tab->conflict)
+ return 0;
+
+ if (tab->row_var[row] < 0 &&
+ report_conflicting_constraint(tab, ~tab->row_var[row]) < 0)
+ return -1;
+
+ tr = tab->mat->row[row] + 2 + tab->M;
+
+ for (j = tab->n_dead; j < tab->n_col; ++j) {
+ if (tab->col_var[j] >= 0 &&
+ (tab->col_var[j] < tab->n_param ||
+ tab->col_var[j] >= tab->n_var - tab->n_div))
+ continue;
+
+ if (!isl_int_is_neg(tr[j]))
+ continue;
+
+ if (tab->col_var[j] < 0 &&
+ report_conflicting_constraint(tab, ~tab->col_var[j]) < 0)
+ return -1;
+ }
+
+ return 0;
+}
+
/* Resolve all known or obviously violated constraints through pivoting.
* In particular, as long as we can find any violated constraint, we
* look for a pivoting column that would result in the lexicographically
* 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)
+static int restore_lexmin(struct isl_tab *tab) WARN_UNUSED;
+static int restore_lexmin(struct isl_tab *tab)
{
int row, col;
if (!tab)
- return NULL;
+ return -1;
if (tab->empty)
- return tab;
+ return 0;
while ((row = first_neg(tab)) != -1) {
col = lexmin_pivot_col(tab, row);
if (col >= tab->n_col) {
+ if (report_conflict(tab, row) < 0)
+ return -1;
if (isl_tab_mark_empty(tab) < 0)
- goto error;
- return tab;
+ return -1;
+ return 0;
}
if (col < 0)
- goto error;
+ return -1;
if (isl_tab_pivot(tab, row, col) < 0)
- goto error;
+ return -1;
}
- return tab;
-error:
- isl_tab_free(tab);
- return NULL;
+ return 0;
}
/* Given a row that represents an equality, look for an appropriate
* 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)
+static int add_lexmin_eq(struct isl_tab *tab, isl_int *eq) WARN_UNUSED;
+static int add_lexmin_eq(struct isl_tab *tab, isl_int *eq)
{
int r1, r2;
int row;
struct isl_tab_undo *snap;
if (!tab)
- return NULL;
+ return -1;
snap = isl_tab_snap(tab);
r1 = isl_tab_add_row(tab, eq);
if (r1 < 0)
- goto error;
+ return -1;
tab->con[r1].is_nonneg = 1;
if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r1]) < 0)
- goto error;
+ return -1;
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]))) {
if (isl_tab_mark_empty(tab) < 0)
- goto error;
- return tab;
+ return -1;
+ return 0;
}
if (isl_tab_rollback(tab, snap) < 0)
- goto error;
- return tab;
+ return -1;
+ return 0;
}
- tab = restore_lexmin(tab);
- if (!tab || tab->empty)
- return tab;
+ if (restore_lexmin(tab) < 0)
+ return -1;
+ if (tab->empty)
+ return 0;
isl_seq_neg(eq, eq, 1 + tab->n_var);
r2 = isl_tab_add_row(tab, eq);
if (r2 < 0)
- goto error;
+ return -1;
tab->con[r2].is_nonneg = 1;
if (isl_tab_push_var(tab, isl_tab_undo_nonneg, &tab->con[r2]) < 0)
- goto error;
+ return -1;
- tab = restore_lexmin(tab);
- if (!tab || tab->empty)
- return tab;
+ if (restore_lexmin(tab) < 0)
+ return -1;
+ if (tab->empty)
+ return 0;
if (!tab->con[r1].is_row) {
if (isl_tab_kill_col(tab, tab->con[r1].index) < 0)
- goto error;
+ return -1;
} else if (!tab->con[r2].is_row) {
if (isl_tab_kill_col(tab, tab->con[r2].index) < 0)
- goto error;
+ return -1;
}
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;
+ return -1;
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;
+ return -1;
if (!tab->bmap)
- goto error;
+ return -1;
}
- return tab;
-error:
- isl_tab_free(tab);
- return NULL;
+ return 0;
}
/* Add an inequality to the tableau, resolving violations using
return tab;
}
- tab = restore_lexmin(tab);
- if (tab && !tab->empty && tab->con[r].is_row &&
+ if (restore_lexmin(tab) < 0)
+ goto error;
+ if (!tab->empty && tab->con[r].is_row &&
isl_tab_row_is_redundant(tab, tab->con[r].index))
if (isl_tab_mark_redundant(tab, tab->con[r].index) < 0)
goto error;
if (row < 0)
goto error;
} while ((var = next_non_integer_var(tab, var, &flags)) != -1);
- tab = restore_lexmin(tab);
- if (!tab || tab->empty)
+ if (restore_lexmin(tab) < 0)
+ goto error;
+ if (tab->empty)
break;
}
return tab;
static struct isl_tab *check_integer_feasible(struct isl_tab *tab)
{
struct isl_tab_undo *snap;
- int feasible;
if (!tab)
return NULL;
if (!tab || tab->empty)
return tab;
}
- if (bmap->n_eq)
- tab = restore_lexmin(tab);
+ if (bmap->n_eq && restore_lexmin(tab) < 0)
+ goto error;
for (i = 0; i < bmap->n_ineq; ++i) {
if (max)
isl_seq_neg(bmap->ineq[i] + 1 + tab->n_param,
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 (add_lexmin_eq(clex->tab, eq) < 0)
+ goto error;
if (check) {
int v = tab_has_valid_sample(clex->tab, eq, 1);
if (v < 0)
clex->context.op = &isl_context_lex_op;
clex->tab = context_tab_for_lexmin(isl_basic_set_copy(dom));
- clex->tab = restore_lexmin(clex->tab);
+ if (restore_lexmin(clex->tab) < 0)
+ goto error;
clex->tab = check_integer_feasible(clex->tab);
if (!clex->tab)
goto error;
static struct isl_tab *add_gbr_eq(struct isl_tab *tab, isl_int *eq)
{
- int r;
-
if (!tab)
return NULL;
{
int i;
int col;
- unsigned dim = tab->n_var - tab->n_param - tab->n_div;
if (tab->n_var == 0)
return -1;
if (isl_tab_kill_col(tab, j) < 0)
goto error;
- tab = restore_lexmin(tab);
+ if (restore_lexmin(tab) < 0)
+ goto error;
}
isl_vec_free(eq);
{
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;
static void find_solutions(struct isl_sol *sol, struct isl_tab *tab)
{
struct isl_context *context;
+ int r;
if (!tab || sol->error)
goto error;
if (context->op->is_empty(context))
goto done;
- for (; tab && !tab->empty; tab = restore_lexmin(tab)) {
+ for (r = 0; r >= 0 && tab && !tab->empty; r = restore_lexmin(tab)) {
int flags;
int row;
enum isl_tab_row_sign sgn;
if (row < 0)
goto error;
}
+ if (r < 0)
+ goto error;
done:
sol_add(sol, tab);
isl_tab_free(tab);
goto error;
context = sol_map->sol.context;
- if (isl_basic_set_fast_is_empty(context->op->peek_basic_set(context)))
+ if (isl_basic_set_plain_is_empty(context->op->peek_basic_set(context)))
/* nothing */;
- else if (isl_basic_map_fast_is_empty(bmap))
+ else if (isl_basic_map_plain_is_empty(bmap))
sol_map_add_empty_if_needed(sol_map,
isl_basic_set_copy(context->op->peek_basic_set(context)));
else {
struct isl_sol_for {
struct isl_sol sol;
int (*fn)(__isl_take isl_basic_set *dom,
- __isl_take isl_mat *map, void *user);
+ __isl_take isl_aff_list *list, void *user);
void *user;
};
*
* Instead of constructing a basic map, this function calls a user
* defined function with the current context as a basic set and
- * an affine matrix representing the relation between the input and output.
- * The number of rows in this matrix is equal to one plus the number
- * of output variables. The number of columns is equal to one plus
- * the total dimension of the context, i.e., the number of parameters,
- * input variables and divs. Since some of the columns in the matrix
- * may refer to the divs, the basic set is not simplified.
- * (Simplification may reorder or remove divs.)
+ * a list of affine expressions representing the relation between
+ * the input and output. The space over which the affine expressions
+ * are defined is the same as that of the domain. The number of
+ * affine expressions in the list is equal to the number of output variables.
*/
static void sol_for_add(struct isl_sol_for *sol,
struct isl_basic_set *dom, struct isl_mat *M)
{
+ int i;
+ isl_ctx *ctx;
+ isl_local_space *ls;
+ isl_aff *aff;
+ isl_aff_list *list;
+
if (sol->sol.error || !dom || !M)
goto error;
- dom = isl_basic_set_simplify(dom);
+ ctx = isl_basic_set_get_ctx(dom);
+ ls = isl_basic_set_get_local_space(dom);
+ list = isl_aff_list_alloc(ctx, M->n_row - 1);
+ for (i = 1; i < M->n_row; ++i) {
+ aff = isl_aff_alloc(isl_local_space_copy(ls));
+ if (aff) {
+ isl_int_set_si(aff->v->el[0], 1);
+ isl_seq_cpy(aff->v->el + 1, M->row[i], M->n_col);
+ }
+ list = isl_aff_list_add(list, aff);
+ }
+ isl_local_space_free(ls);
+
dom = isl_basic_set_finalize(dom);
- if (sol->fn(isl_basic_set_copy(dom), isl_mat_copy(M), sol->user) < 0)
+ if (sol->fn(isl_basic_set_copy(dom), list, sol->user) < 0)
goto error;
isl_basic_set_free(dom);
}
static struct isl_sol_for *sol_for_init(struct isl_basic_map *bmap, int max,
- int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_mat *map,
+ int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_aff_list *list,
void *user),
void *user)
{
}
int isl_basic_map_foreach_lexopt(__isl_keep isl_basic_map *bmap, int max,
- int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_mat *map,
+ int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_aff_list *list,
void *user),
void *user)
{
bmap = isl_basic_map_detect_equalities(bmap);
sol_for = sol_for_init(bmap, max, fn, user);
- if (isl_basic_map_fast_is_empty(bmap))
+ if (isl_basic_map_plain_is_empty(bmap))
/* nothing */;
else {
struct isl_tab *tab;
}
int isl_basic_map_foreach_lexmin(__isl_keep isl_basic_map *bmap,
- int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_mat *map,
+ int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_aff_list *list,
void *user),
void *user)
{
}
int isl_basic_map_foreach_lexmax(__isl_keep isl_basic_map *bmap,
- int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_mat *map,
+ int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_aff_list *list,
void *user),
void *user)
{
return isl_basic_map_foreach_lexopt(bmap, 1, fn, user);
}
+
+int isl_basic_set_foreach_lexmax(__isl_keep isl_basic_set *bset,
+ int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_aff_list *list,
+ void *user),
+ void *user)
+{
+ return isl_basic_map_foreach_lexmax(bset, fn, user);
+}
+
+/* Check if the given sequence of len variables starting at pos
+ * represents a trivial (i.e., zero) solution.
+ * The variables are assumed to be non-negative and to come in pairs,
+ * with each pair representing a variable of unrestricted sign.
+ * The solution is trivial if each such pair in the sequence consists
+ * of two identical values, meaning that the variable being represented
+ * has value zero.
+ */
+static int region_is_trivial(struct isl_tab *tab, int pos, int len)
+{
+ int i;
+
+ if (len == 0)
+ return 0;
+
+ for (i = 0; i < len; i += 2) {
+ int neg_row;
+ int pos_row;
+
+ neg_row = tab->var[pos + i].is_row ?
+ tab->var[pos + i].index : -1;
+ pos_row = tab->var[pos + i + 1].is_row ?
+ tab->var[pos + i + 1].index : -1;
+
+ if ((neg_row < 0 ||
+ isl_int_is_zero(tab->mat->row[neg_row][1])) &&
+ (pos_row < 0 ||
+ isl_int_is_zero(tab->mat->row[pos_row][1])))
+ continue;
+
+ if (neg_row < 0 || pos_row < 0)
+ return 0;
+ if (isl_int_ne(tab->mat->row[neg_row][1],
+ tab->mat->row[pos_row][1]))
+ return 0;
+ }
+
+ return 1;
+}
+
+/* Return the index of the first trivial region or -1 if all regions
+ * are non-trivial.
+ */
+static int first_trivial_region(struct isl_tab *tab,
+ int n_region, struct isl_region *region)
+{
+ int i;
+
+ for (i = 0; i < n_region; ++i) {
+ if (region_is_trivial(tab, region[i].pos, region[i].len))
+ return i;
+ }
+
+ return -1;
+}
+
+/* Check if the solution is optimal, i.e., whether the first
+ * n_op entries are zero.
+ */
+static int is_optimal(__isl_keep isl_vec *sol, int n_op)
+{
+ int i;
+
+ for (i = 0; i < n_op; ++i)
+ if (!isl_int_is_zero(sol->el[1 + i]))
+ return 0;
+ return 1;
+}
+
+/* Add constraints to "tab" that ensure that any solution is significantly
+ * better that that represented by "sol". That is, find the first
+ * relevant (within first n_op) non-zero coefficient and force it (along
+ * with all previous coefficients) to be zero.
+ * If the solution is already optimal (all relevant coefficients are zero),
+ * then just mark the table as empty.
+ */
+static int force_better_solution(struct isl_tab *tab,
+ __isl_keep isl_vec *sol, int n_op)
+{
+ int i;
+ isl_ctx *ctx;
+ isl_vec *v = NULL;
+
+ if (!sol)
+ return -1;
+
+ for (i = 0; i < n_op; ++i)
+ if (!isl_int_is_zero(sol->el[1 + i]))
+ break;
+
+ if (i == n_op) {
+ if (isl_tab_mark_empty(tab) < 0)
+ return -1;
+ return 0;
+ }
+
+ ctx = isl_vec_get_ctx(sol);
+ v = isl_vec_alloc(ctx, 1 + tab->n_var);
+ if (!v)
+ return -1;
+
+ for (; i >= 0; --i) {
+ v = isl_vec_clr(v);
+ isl_int_set_si(v->el[1 + i], -1);
+ if (add_lexmin_eq(tab, v->el) < 0)
+ goto error;
+ }
+
+ isl_vec_free(v);
+ return 0;
+error:
+ isl_vec_free(v);
+ return -1;
+}
+
+struct isl_trivial {
+ int update;
+ int region;
+ int side;
+ struct isl_tab_undo *snap;
+};
+
+/* Return the lexicographically smallest non-trivial solution of the
+ * given ILP problem.
+ *
+ * All variables are assumed to be non-negative.
+ *
+ * n_op is the number of initial coordinates to optimize.
+ * That is, once a solution has been found, we will only continue looking
+ * for solution that result in significantly better values for those
+ * initial coordinates. That is, we only continue looking for solutions
+ * that increase the number of initial zeros in this sequence.
+ *
+ * A solution is non-trivial, if it is non-trivial on each of the
+ * specified regions. Each region represents a sequence of pairs
+ * of variables. A solution is non-trivial on such a region if
+ * at least one of these pairs consists of different values, i.e.,
+ * such that the non-negative variable represented by the pair is non-zero.
+ *
+ * Whenever a conflict is encountered, all constraints involved are
+ * reported to the caller through a call to "conflict".
+ *
+ * We perform a simple branch-and-bound backtracking search.
+ * Each level in the search represents initially trivial region that is forced
+ * to be non-trivial.
+ * At each level we consider n cases, where n is the length of the region.
+ * In terms of the n/2 variables of unrestricted signs being encoded by
+ * the region, we consider the cases
+ * x_0 >= 1
+ * x_0 <= -1
+ * x_0 = 0 and x_1 >= 1
+ * x_0 = 0 and x_1 <= -1
+ * x_0 = 0 and x_1 = 0 and x_2 >= 1
+ * x_0 = 0 and x_1 = 0 and x_2 <= -1
+ * ...
+ * The cases are considered in this order, assuming that each pair
+ * x_i_a x_i_b represents the value x_i_b - x_i_a.
+ * That is, x_0 >= 1 is enforced by adding the constraint
+ * x_0_b - x_0_a >= 1
+ */
+__isl_give isl_vec *isl_tab_basic_set_non_trivial_lexmin(
+ __isl_take isl_basic_set *bset, int n_op, int n_region,
+ struct isl_region *region,
+ int (*conflict)(int con, void *user), void *user)
+{
+ int i, j;
+ int r;
+ isl_ctx *ctx = isl_basic_set_get_ctx(bset);
+ isl_vec *v = NULL;
+ isl_vec *sol = isl_vec_alloc(ctx, 0);
+ struct isl_tab *tab;
+ struct isl_trivial *triv = NULL;
+ int level, init;
+
+ tab = tab_for_lexmin(isl_basic_map_from_range(bset), NULL, 0, 0);
+ if (!tab)
+ goto error;
+ tab->conflict = conflict;
+ tab->conflict_user = user;
+
+ v = isl_vec_alloc(ctx, 1 + tab->n_var);
+ triv = isl_calloc_array(ctx, struct isl_trivial, n_region);
+ if (!v || !triv)
+ goto error;
+
+ level = 0;
+ init = 1;
+
+ while (level >= 0) {
+ int side, base;
+
+ if (init) {
+ tab = cut_to_integer_lexmin(tab);
+ if (!tab)
+ goto error;
+ if (tab->empty)
+ goto backtrack;
+ r = first_trivial_region(tab, n_region, region);
+ if (r < 0) {
+ for (i = 0; i < level; ++i)
+ triv[i].update = 1;
+ isl_vec_free(sol);
+ sol = isl_tab_get_sample_value(tab);
+ if (!sol)
+ goto error;
+ if (is_optimal(sol, n_op))
+ break;
+ goto backtrack;
+ }
+ if (level >= n_region)
+ isl_die(ctx, isl_error_internal,
+ "nesting level too deep", goto error);
+ if (isl_tab_extend_cons(tab,
+ 2 * region[r].len + 2 * n_op) < 0)
+ goto error;
+ triv[level].region = r;
+ triv[level].side = 0;
+ }
+
+ r = triv[level].region;
+ side = triv[level].side;
+ base = 2 * (side/2);
+
+ if (side >= region[r].len) {
+backtrack:
+ level--;
+ init = 0;
+ if (level >= 0)
+ if (isl_tab_rollback(tab, triv[level].snap) < 0)
+ goto error;
+ continue;
+ }
+
+ if (triv[level].update) {
+ if (force_better_solution(tab, sol, n_op) < 0)
+ goto error;
+ triv[level].update = 0;
+ }
+
+ if (side == base && base >= 2) {
+ for (j = base - 2; j < base; ++j) {
+ v = isl_vec_clr(v);
+ isl_int_set_si(v->el[1 + region[r].pos + j], 1);
+ if (add_lexmin_eq(tab, v->el) < 0)
+ goto error;
+ }
+ }
+
+ triv[level].snap = isl_tab_snap(tab);
+ if (isl_tab_push_basis(tab) < 0)
+ goto error;
+
+ v = isl_vec_clr(v);
+ isl_int_set_si(v->el[0], -1);
+ isl_int_set_si(v->el[1 + region[r].pos + side], -1);
+ isl_int_set_si(v->el[1 + region[r].pos + (side ^ 1)], 1);
+ tab = add_lexmin_ineq(tab, v->el);
+
+ triv[level].side++;
+ level++;
+ init = 1;
+ }
+
+ free(triv);
+ isl_vec_free(v);
+ isl_tab_free(tab);
+ isl_basic_set_free(bset);
+
+ return sol;
+error:
+ free(triv);
+ isl_vec_free(v);
+ isl_tab_free(tab);
+ isl_basic_set_free(bset);
+ isl_vec_free(sol);
+ return NULL;
+}
+
+/* Return the lexicographically smallest rational point in "bset",
+ * assuming that all variables are non-negative.
+ * If "bset" is empty, then return a zero-length vector.
+ */
+ __isl_give isl_vec *isl_tab_basic_set_non_neg_lexmin(
+ __isl_take isl_basic_set *bset)
+{
+ struct isl_tab *tab;
+ isl_ctx *ctx = isl_basic_set_get_ctx(bset);
+ isl_vec *sol;
+
+ tab = tab_for_lexmin(isl_basic_map_from_range(bset), NULL, 0, 0);
+ if (!tab)
+ goto error;
+ if (tab->empty)
+ sol = isl_vec_alloc(ctx, 0);
+ else
+ sol = isl_tab_get_sample_value(tab);
+ isl_tab_free(tab);
+ isl_basic_set_free(bset);
+ return sol;
+error:
+ isl_tab_free(tab);
+ isl_basic_set_free(bset);
+ return NULL;
+}