add isl_aff_mod_val
[platform/upstream/isl.git] / isl_tab_pip.c
index 412e8ee..2592a28 100644 (file)
@@ -2,7 +2,7 @@
  * Copyright 2008-2009 Katholieke Universiteit Leuven
  * Copyright 2010      INRIA Saclay
  *
- * Use of this software is governed by the GNU LGPLv2.1 license
+ * Use of this software is governed by the MIT license
  *
  * Written by Sven Verdoolaege, K.U.Leuven, Departement
  * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
@@ -17,6 +17,7 @@
 #include "isl_sample.h"
 #include <isl_mat_private.h>
 #include <isl_aff_private.h>
+#include <isl_options_private.h>
 #include <isl_config.h>
 
 /*
@@ -101,6 +102,8 @@ struct isl_context_op {
        void *(*save)(struct isl_context *context);
        /* restore saved context */
        void (*restore)(struct isl_context *context, void *);
+       /* discard saved context */
+       void (*discard)(void *);
        /* invalidate context */
        void (*invalidate)(struct isl_context *context);
        /* free context */
@@ -116,6 +119,14 @@ struct isl_context_lex {
        struct isl_tab *tab;
 };
 
+/* A stack (linked list) of solutions of subtrees of the search space.
+ *
+ * "M" describes the solution in terms of the dimensions of "dom".
+ * The number of columns of "M" is one more than the total number
+ * of dimensions of "dom".
+ *
+ * If "M" is NULL, then there is no solution on "dom".
+ */
 struct isl_partial_sol {
        int level;
        struct isl_basic_set *dom;
@@ -200,6 +211,7 @@ static void sol_push_sol(struct isl_sol *sol,
        return;
 error:
        isl_basic_set_free(dom);
+       isl_mat_free(M);
        sol->error = 1;
 }
 
@@ -303,13 +315,26 @@ static void sol_pop(struct isl_sol *sol)
                        sol_pop_one(sol);
                } else {
                        struct isl_basic_set *bset;
+                       isl_mat *M;
+                       unsigned n;
 
+                       n = isl_basic_set_dim(partial->next->dom, isl_dim_div);
+                       n -= n_div;
                        bset = sol_domain(sol);
-
                        isl_basic_set_free(partial->next->dom);
                        partial->next->dom = bset;
+                       M = partial->next->M;
+                       if (M) {
+                               M = isl_mat_drop_cols(M, M->n_col - n, n);
+                               partial->next->M = M;
+                               if (!M)
+                                       goto error;
+                       }
                        partial->next->level = sol->level;
 
+                       if (!bset)
+                               goto error;
+
                        sol->partial = partial->next;
                        isl_basic_set_free(partial->dom);
                        isl_mat_free(partial->M);
@@ -317,6 +342,9 @@ static void sol_pop(struct isl_sol *sol)
                }
        } else
                sol_pop_one(sol);
+
+       if (0)
+error:         sol->error = 1;
 }
 
 static void sol_dec_level(struct isl_sol *sol)
@@ -430,6 +458,8 @@ static void sol_add(struct isl_sol *sol, struct isl_tab *tab)
 
        if (tab->empty && !sol->add_empty)
                return;
+       if (sol->context->op->is_empty(sol->context))
+               return;
 
        bset = sol_domain(sol);
 
@@ -742,6 +772,30 @@ static struct isl_vec *get_row_parameter_ineq(struct isl_tab *tab, int row)
        return ineq;
 }
 
+/* Normalize a div expression of the form
+ *
+ *     [(g*f(x) + c)/(g * m)]
+ *
+ * with c the constant term and f(x) the remaining coefficients, to
+ *
+ *     [(f(x) + [c/g])/m]
+ */
+static void normalize_div(__isl_keep isl_vec *div)
+{
+       isl_ctx *ctx = isl_vec_get_ctx(div);
+       int len = div->size - 2;
+
+       isl_seq_gcd(div->el + 2, len, &ctx->normalize_gcd);
+       isl_int_gcd(ctx->normalize_gcd, ctx->normalize_gcd, div->el[0]);
+
+       if (isl_int_is_one(ctx->normalize_gcd))
+               return;
+
+       isl_int_divexact(div->el[0], div->el[0], ctx->normalize_gcd);
+       isl_int_fdiv_q(div->el[1], div->el[1], ctx->normalize_gcd);
+       isl_seq_scale_down(div->el + 2, div->el + 2, ctx->normalize_gcd, len);
+}
+
 /* Return a integer division for use in a parametric cut based on the given row.
  * In particular, let the parametric constant of the row be
  *
@@ -762,8 +816,8 @@ static struct isl_vec *get_row_parameter_div(struct isl_tab *tab, int row)
 
        isl_int_set(div->el[0], tab->mat->row[row][0]);
        get_row_parameter_line(tab, row, div->el + 1);
-       div = isl_vec_normalize(div);
        isl_seq_neg(div->el + 1, div->el + 1, div->size - 1);
+       normalize_div(div);
        isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
 
        return div;
@@ -790,7 +844,7 @@ static struct isl_vec *get_row_split_div(struct isl_tab *tab, int row)
 
        isl_int_set(div->el[0], tab->mat->row[row][0]);
        get_row_parameter_line(tab, row, div->el + 1);
-       div = isl_vec_normalize(div);
+       normalize_div(div);
        isl_seq_fdiv_r(div->el + 1, div->el + 1, div->el[0], div->size - 1);
 
        return div;
@@ -828,7 +882,7 @@ static struct isl_vec *ineq_for_div(struct isl_basic_set *bset, unsigned div)
 }
 
 /* Given a row in the tableau and a div that was created
- * using get_row_split_div and that been constrained to equality, i.e.,
+ * using get_row_split_div and that has been constrained to equality, i.e.,
  *
  *             d = floor(\sum_i {a_i} y_i) = \sum_i {a_i} y_i
  *
@@ -859,7 +913,8 @@ static struct isl_tab *set_row_cst_to_div(struct isl_tab *tab, int row, int div)
        } else {
                int dcol = tab->var[tab->n_var - tab->n_div + div].index;
 
-               isl_int_set_si(tab->mat->row[row][2 + tab->M + dcol], 1);
+               isl_int_add_ui(tab->mat->row[row][2 + tab->M + dcol],
+                               tab->mat->row[row][2 + tab->M + dcol], 1);
        }
 
        return tab;
@@ -1149,8 +1204,8 @@ static int report_conflicting_constraint(struct isl_tab *tab, int con)
 
 /* 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.
+ * (if it represents a constraint) and all constraint columns with
+ * non-zero (and therefore negative) coefficients.
  */
 static int report_conflict(struct isl_tab *tab, int row)
 {
@@ -1602,6 +1657,9 @@ static int add_cut(struct isl_tab *tab, int row)
        return tab->con[r].index;
 }
 
+#define CUT_ALL 1
+#define CUT_ONE 0
+
 /* Given a non-parametric tableau, add cuts until an integer
  * sample point is obtained or until the tableau is determined
  * to be integer infeasible.
@@ -1613,8 +1671,12 @@ static int add_cut(struct isl_tab *tab, int row)
  * combination of variables/constraints plus a non-integral constant,
  * then there is no way to obtain an integer point and we return
  * a tableau that is marked empty.
+ * The parameter cutting_strategy controls the strategy used when adding cuts
+ * to remove non-integer points. CUT_ALL adds all possible cuts
+ * before continuing the search. CUT_ONE adds only one cut at a time.
  */
-static struct isl_tab *cut_to_integer_lexmin(struct isl_tab *tab)
+static struct isl_tab *cut_to_integer_lexmin(struct isl_tab *tab,
+       int cutting_strategy)
 {
        int var;
        int row;
@@ -1636,6 +1698,8 @@ static struct isl_tab *cut_to_integer_lexmin(struct isl_tab *tab)
                        row = add_cut(tab, row);
                        if (row < 0)
                                goto error;
+                       if (cutting_strategy == CUT_ONE)
+                               break;
                } while ((var = next_non_integer_var(tab, var, &flags)) != -1);
                if (restore_lexmin(tab) < 0)
                        goto error;
@@ -1726,7 +1790,7 @@ static struct isl_tab *check_integer_feasible(struct isl_tab *tab)
        if (isl_tab_push_basis(tab) < 0)
                goto error;
 
-       tab = cut_to_integer_lexmin(tab);
+       tab = cut_to_integer_lexmin(tab, CUT_ALL);
        if (!tab)
                goto error;
 
@@ -1920,6 +1984,7 @@ static int add_parametric_cut(struct isl_tab *tab, int row,
 
        n = tab->n_div;
        d = context->op->get_div(context, tab, div);
+       isl_vec_free(div);
        if (d < 0)
                return -1;
 
@@ -1985,8 +2050,6 @@ static int add_parametric_cut(struct isl_tab *tab, int row,
        if (tab->row_sign)
                tab->row_sign[tab->con[r].index] = isl_tab_row_neg;
 
-       isl_vec_free(div);
-
        row = tab->con[r].index;
 
        if (d >= n && context->op->detect_equalities(context, tab) < 0)
@@ -2381,6 +2444,10 @@ static void context_lex_restore(struct isl_context *context, void *save)
        }
 }
 
+static void context_lex_discard(void *save)
+{
+}
+
 static int context_lex_is_ok(struct isl_context *context)
 {
        struct isl_context_lex *clex = (struct isl_context_lex *)context;
@@ -2500,6 +2567,7 @@ struct isl_context_op isl_context_lex_op = {
        context_lex_is_ok,
        context_lex_save,
        context_lex_restore,
+       context_lex_discard,
        context_lex_invalidate,
        context_lex_free,
 };
@@ -2508,7 +2576,6 @@ 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);
@@ -2549,6 +2616,14 @@ error:
        return NULL;
 }
 
+/* Representation of the context when using generalized basis reduction.
+ *
+ * "shifted" contains the offsets of the unit hypercubes that lie inside the
+ * context.  Any rational point in "shifted" can therefore be rounded
+ * up to an integer point in the context.
+ * If the context is constrained by any equality, then "shifted" is not used
+ * as it would be empty.
+ */
 struct isl_context_gbr {
        struct isl_context context;
        struct isl_tab *tab;
@@ -2607,7 +2682,7 @@ static void gbr_init_shifted(struct isl_context_gbr *cgbr)
                }
        }
 
-       cgbr->shifted = isl_tab_from_basic_set(bset);
+       cgbr->shifted = isl_tab_from_basic_set(bset, 0);
 
        for (i = 0; i < bset->n_ineq; ++i)
                isl_int_set(bset->ineq[i][0], cst->el[i]);
@@ -2680,7 +2755,8 @@ static struct isl_vec *gbr_get_sample(struct isl_context_gbr *cgbr)
                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)
+               if (isl_tab_track_bset(cgbr->cone,
+                                       isl_basic_set_copy(bset)) < 0)
                        return NULL;
        }
        if (isl_tab_detect_implicit_equalities(cgbr->cone) < 0)
@@ -2771,6 +2847,15 @@ error:
        return NULL;
 }
 
+/* Add the equality described by "eq" to the context.
+ * If "check" is set, then we check if the context is empty after
+ * adding the equality.
+ * If "update" is set, then we check if the samples are still valid.
+ *
+ * We do not explicitly add shifted copies of the equality to
+ * cgbr->shifted since they would conflict with each other.
+ * Instead, we directly mark cgbr->shifted empty.
+ */
 static void context_gbr_add_eq(struct isl_context *context, isl_int *eq,
                int check, int update)
 {
@@ -2778,6 +2863,11 @@ static void context_gbr_add_eq(struct isl_context *context, isl_int *eq,
 
        cgbr->tab = add_gbr_eq(cgbr->tab, eq);
 
+       if (cgbr->shifted && !cgbr->shifted->empty && use_shifted(cgbr)) {
+               if (isl_tab_mark_empty(cgbr->shifted) < 0)
+                       goto error;
+       }
+
        if (cgbr->cone && cgbr->cone->n_col != cgbr->cone->n_dead) {
                if (isl_tab_extend_cons(cgbr->cone, 2) < 0)
                        goto error;
@@ -3042,7 +3132,8 @@ static int context_gbr_detect_equalities(struct isl_context *context,
                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)
+               if (isl_tab_track_bset(cgbr->cone,
+                                       isl_basic_set_copy(bset)) < 0)
                        goto error;
        }
        if (isl_tab_detect_implicit_equalities(cgbr->cone) < 0)
@@ -3050,7 +3141,9 @@ static int context_gbr_detect_equalities(struct isl_context *context,
 
        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)
+       if (!cgbr->tab)
+               return -1;
+       if (cgbr->tab->bmap->n_ineq > n_ineq)
                propagate_equalities(cgbr, tab, n_ineq);
 
        return 0;
@@ -3187,6 +3280,12 @@ error:
        cgbr->tab = NULL;
 }
 
+static void context_gbr_discard(void *save)
+{
+       struct isl_gbr_tab_undo *snap = (struct isl_gbr_tab_undo *)save;
+       free(snap);
+}
+
 static int context_gbr_is_ok(struct isl_context *context)
 {
        struct isl_context_gbr *cgbr = (struct isl_context_gbr *)context;
@@ -3225,6 +3324,7 @@ struct isl_context_op isl_context_gbr_op = {
        context_gbr_is_ok,
        context_gbr_save,
        context_gbr_restore,
+       context_gbr_discard,
        context_gbr_invalidate,
        context_gbr_free,
 };
@@ -3244,13 +3344,10 @@ static struct isl_context *isl_context_gbr_alloc(struct isl_basic_set *dom)
 
        cgbr->shifted = NULL;
        cgbr->cone = NULL;
-       cgbr->tab = isl_tab_from_basic_set(dom);
+       cgbr->tab = isl_tab_from_basic_set(dom, 1);
        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;
@@ -3526,6 +3623,8 @@ static void find_in_pos(struct isl_sol *sol, struct isl_tab *tab, isl_int *ineq)
 
        if (!sol->error)
                sol->context->op->restore(sol->context, saved);
+       else
+               sol->context->op->discard(saved);
        return;
 error:
        sol->error = 1;
@@ -3777,6 +3876,24 @@ error:
        sol->error = 1;
 }
 
+/* Does "sol" contain a pair of partial solutions that could potentially
+ * be merged?
+ *
+ * We currently only check that "sol" is not in an error state
+ * and that there are at least two partial solutions of which the final two
+ * are defined at the same level.
+ */
+static int sol_has_mergeable_solutions(struct isl_sol *sol)
+{
+       if (sol->error)
+               return 0;
+       if (!sol->partial)
+               return 0;
+       if (!sol->partial->next)
+               return 0;
+       return sol->partial->level == sol->partial->next->level;
+}
+
 /* Compute the lexicographic minimum of the set represented by the main
  * tableau "tab" within the context "sol->context_tab".
  *
@@ -3787,10 +3904,20 @@ error:
  * corresponding rows may not be marked as being non-negative.
  * In parts of the context where the added equality does not hold,
  * the main tableau is marked as being empty.
+ *
+ * Before we embark on the actual computation, we save a copy
+ * of the context.  When we return, we check if there are any
+ * partial solutions that can potentially be merged.  If so,
+ * we perform a rollback to the initial state of the context.
+ * The merging of partial solutions happens inside calls to
+ * sol_dec_level that are pushed onto the undo stack of the context.
+ * If there are no partial solutions that can potentially be merged
+ * then the rollback is skipped as it would just be wasted effort.
  */
 static void find_solutions_main(struct isl_sol *sol, struct isl_tab *tab)
 {
        int row;
+       void *saved;
 
        if (!tab)
                goto error;
@@ -3840,8 +3967,15 @@ static void find_solutions_main(struct isl_sol *sol, struct isl_tab *tab)
                row = tab->n_redundant - 1;
        }
 
+       saved = sol->context->op->save(sol->context);
+
        find_solutions(sol, tab);
 
+       if (sol_has_mergeable_solutions(sol))
+               sol->context->op->restore(sol->context, saved);
+       else
+               sol->context->op->discard(saved);
+
        sol->level = 0;
        sol_pop(sol);
 
@@ -4486,7 +4620,7 @@ static union isl_lex_res basic_map_partial_lexopt_symm(
        bmap = isl_basic_map_finalize(bmap);
 
        n_div = isl_basic_set_dim(dom, isl_dim_div);
-       dom = isl_basic_set_add(dom, isl_dim_set, 1);
+       dom = isl_basic_set_add_dims(dom, isl_dim_set, 1);
        dom = isl_basic_set_extend_constraints(dom, 0, n);
        for (i = 0; i < n; ++i) {
                k = isl_basic_set_alloc_inequality(dom);
@@ -4650,6 +4784,7 @@ static void sol_for_add(struct isl_sol_for *sol,
                        isl_int_set(aff->v->el[0], M->row[0][0]);
                        isl_seq_cpy(aff->v->el + 1, M->row[i], M->n_col);
                }
+               aff = isl_aff_normalize(aff);
                list = isl_aff_list_add(list, aff);
        }
        isl_local_space_free(ls);
@@ -4727,11 +4862,13 @@ int isl_basic_map_foreach_lexopt(__isl_keep isl_basic_map *bmap, int max,
        struct isl_sol_for *sol_for = NULL;
 
        bmap = isl_basic_map_copy(bmap);
+       bmap = isl_basic_map_detect_equalities(bmap);
        if (!bmap)
                return -1;
 
-       bmap = isl_basic_map_detect_equalities(bmap);
        sol_for = sol_for_init(bmap, max, fn, user);
+       if (!sol_for)
+               goto error;
 
        if (isl_basic_map_plain_is_empty(bmap))
                /* nothing */;
@@ -4763,30 +4900,6 @@ int isl_basic_set_foreach_lexopt(__isl_keep isl_basic_set *bset, int max,
        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_aff_list *list,
-                 void *user),
-       void *user)
-{
-       return isl_basic_map_foreach_lexopt(bmap, 0, fn, user);
-}
-
-int isl_basic_map_foreach_lexmax(__isl_keep isl_basic_map *bmap,
-       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,
@@ -4954,13 +5067,19 @@ __isl_give isl_vec *isl_tab_basic_set_non_trivial_lexmin(
 {
        int i, j;
        int r;
-       isl_ctx *ctx = isl_basic_set_get_ctx(bset);
+       isl_ctx *ctx;
        isl_vec *v = NULL;
-       isl_vec *sol = isl_vec_alloc(ctx, 0);
+       isl_vec *sol = NULL;
        struct isl_tab *tab;
        struct isl_trivial *triv = NULL;
        int level, init;
 
+       if (!bset)
+               return NULL;
+
+       ctx = isl_basic_set_get_ctx(bset);
+       sol = isl_vec_alloc(ctx, 0);
+
        tab = tab_for_lexmin(bset, NULL, 0, 0);
        if (!tab)
                goto error;
@@ -4979,7 +5098,7 @@ __isl_give isl_vec *isl_tab_basic_set_non_trivial_lexmin(
                int side, base;
 
                if (init) {
-                       tab = cut_to_integer_lexmin(tab);
+                       tab = cut_to_integer_lexmin(tab, CUT_ONE);
                        if (!tab)
                                goto error;
                        if (tab->empty)
@@ -5076,6 +5195,9 @@ __isl_give isl_vec *isl_tab_basic_set_non_neg_lexmin(
        isl_ctx *ctx = isl_basic_set_get_ctx(bset);
        isl_vec *sol;
 
+       if (!bset)
+               return NULL;
+
        tab = tab_for_lexmin(bset, NULL, 0, 0);
        if (!tab)
                goto error;
@@ -5160,6 +5282,7 @@ static void sol_pma_add(struct isl_sol_pma *sol,
        isl_local_space_free(ls);
        isl_mat_free(M);
        dom = isl_basic_set_simplify(dom);
+       dom = isl_basic_set_finalize(dom);
        pma = isl_pw_multi_aff_alloc(isl_set_from_basic_set(dom), maff);
        sol->pma = isl_pw_multi_aff_add_disjoint(sol->pma, pma);
        if (!sol->pma)