add isl_basic_set_foreach_lexopt
[platform/upstream/isl.git] / isl_tab_pip.c
index b63d3b7..01657a0 100644 (file)
@@ -16,6 +16,8 @@
 #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
@@ -26,7 +28,7 @@
  * The strategy used for obtaining a feasible solution is different
  * from the one used in isl_tab.c.  In particular, in isl_tab.c,
  * upon finding a constraint that is not yet satisfied, we pivot
- * in a row that increases the constant term of row holding the
+ * in a row that increases the constant term of the row holding the
  * constraint, making sure the sample solution remains feasible
  * for all the constraints it already satisfied.
  * Here, we always pivot in the row holding the constraint,
@@ -420,7 +422,7 @@ 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;
+       int row;
        isl_int m;
 
        if (sol->error || !tab)
@@ -585,7 +587,6 @@ static void sol_map_add(struct isl_sol_map *sol,
 {
        int i;
        struct isl_basic_map *bmap = NULL;
-       isl_basic_set *context_bset;
        unsigned n_eq;
        unsigned n_ineq;
        unsigned nparam;
@@ -1113,37 +1114,117 @@ static int first_neg(struct isl_tab *tab)
        return -1;
 }
 
+/* Check whether the invariant that all columns are lexico-positive
+ * 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;
+       int col;
+       int var;
+       int row;
+
+       for (col = tab->n_dead; col < tab->n_col; ++col) {
+               if (tab->col_var[col] >= 0 &&
+                   (tab->col_var[col] < tab->n_param ||
+                    tab->col_var[col] >= tab->n_var - tab->n_div))
+                       continue;
+               for (var = tab->n_param; var < tab->n_var - tab->n_div; ++var) {
+                       if (!tab->var[var].is_row) {
+                               if (tab->var[var].index == col)
+                                       break;
+                               else
+                                       continue;
+                       }
+                       row = tab->var[var].index;
+                       if (isl_int_is_zero(tab->mat->row[row][off + col]))
+                               continue;
+                       if (isl_int_is_pos(tab->mat->row[row][off + col]))
+                               break;
+                       fprintf(stderr, "lexneg column %d (row %d)\n",
+                               col, row);
+               }
+               if (var >= tab->n_var - tab->n_div)
+                       fprintf(stderr, "zero column %d\n", col);
+       }
+}
+
+/* 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
@@ -1254,90 +1335,77 @@ 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)
+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;
-       } 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) {
-                       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;
-               }
+                       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
@@ -1368,8 +1436,9 @@ static struct isl_tab *add_lexmin_ineq(struct isl_tab *tab, isl_int *ineq)
                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;
@@ -1580,8 +1649,9 @@ static struct isl_tab *cut_to_integer_lexmin(struct isl_tab *tab)
                        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;
@@ -1660,7 +1730,6 @@ static int sample_is_finite(struct isl_tab *tab)
 static struct isl_tab *check_integer_feasible(struct isl_tab *tab)
 {
        struct isl_tab_undo *snap;
-       int feasible;
 
        if (!tab)
                return NULL;
@@ -1993,8 +2062,8 @@ static struct isl_tab *tab_for_lexmin(struct isl_basic_map *bmap,
                if (!tab || tab->empty)
                        return tab;
        }
-       if (bmap->n_eq)
-               tab = restore_lexmin(tab);
+       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,
@@ -2114,24 +2183,14 @@ static struct isl_tab *context_lex_peek_tab(struct isl_context *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 (add_lexmin_eq(clex->tab, eq) < 0)
+               goto error;
        if (check) {
                int v = tab_has_valid_sample(clex->tab, eq, 1);
                if (v < 0)
@@ -2490,7 +2549,8 @@ static struct isl_context *isl_context_lex_alloc(struct isl_basic_set *dom)
        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;
@@ -2708,8 +2768,6 @@ error:
 
 static struct isl_tab *add_gbr_eq(struct isl_tab *tab, isl_int *eq)
 {
-       int r;
-
        if (!tab)
                return NULL;
 
@@ -2896,7 +2954,6 @@ 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;
@@ -2970,7 +3027,8 @@ static void propagate_equalities(struct isl_context_gbr *cgbr,
                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);
@@ -2987,8 +3045,6 @@ static int context_gbr_detect_equalities(struct isl_context *context,
 {
        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;
@@ -3616,6 +3672,7 @@ error:
 static void find_solutions(struct isl_sol *sol, struct isl_tab *tab)
 {
        struct isl_context *context;
+       int r;
 
        if (!tab || sol->error)
                goto error;
@@ -3627,7 +3684,7 @@ static void find_solutions(struct isl_sol *sol, struct isl_tab *tab)
        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;
@@ -3721,6 +3778,8 @@ static void find_solutions(struct isl_sol *sol, struct isl_tab *tab)
                if (row < 0)
                        goto error;
        }
+       if (r < 0)
+               goto error;
 done:
        sol_add(sol, tab);
        isl_tab_free(tab);
@@ -3911,9 +3970,9 @@ static __isl_give isl_map *basic_map_partial_lexopt_base(
                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 {
@@ -4432,6 +4491,9 @@ struct isl_map *isl_tab_basic_map_partial_lexopt(
        isl_assert(bmap->ctx,
            isl_basic_map_compatible_domain(bmap, dom), goto error);
 
+       if (isl_basic_set_dim(dom, isl_dim_all) == 0)
+               return basic_map_partial_lexopt(bmap, dom, empty, max);
+
        bmap = isl_basic_map_intersect_domain(bmap, isl_basic_set_copy(dom));
        bmap = isl_basic_map_detect_equalities(bmap);
        bmap = isl_basic_map_remove_redundancies(bmap);
@@ -4446,7 +4508,7 @@ error:
 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;
 };
 
@@ -4468,24 +4530,39 @@ static void sol_for_free_wrap(struct isl_sol *sol)
  *
  * Instead of constructing a basic map, this function calls a user
  * defined function with the current context as a basic set and
- * an affine matrix 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);
@@ -4504,7 +4581,7 @@ static void sol_for_add_wrap(struct isl_sol *sol,
 }
 
 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)
 {
@@ -4549,7 +4626,7 @@ static void sol_for_find_solutions(struct isl_sol_for *sol_for,
 }
 
 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)
 {
@@ -4562,7 +4639,7 @@ int isl_basic_map_foreach_lexopt(__isl_keep isl_basic_map *bmap, int max,
        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;
@@ -4584,8 +4661,16 @@ error:
        return -1;
 }
 
+int isl_basic_set_foreach_lexopt(__isl_keep isl_basic_set *bset, int max,
+       int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_aff_list *list,
+                 void *user),
+       void *user)
+{
+       return isl_basic_map_foreach_lexopt(bset, max, fn, user);
+}
+
 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)
 {
@@ -4593,9 +4678,322 @@ int isl_basic_map_foreach_lexmin(__isl_keep isl_basic_map *bmap,
 }
 
 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;
+}