/*
* 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
+ * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite,
+ * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France
*/
+#include <isl_ctx_private.h>
#include "isl_map_private.h"
#include <isl/seq.h>
#include "isl_tab.h"
#include "isl_sample.h"
#include <isl_mat_private.h>
+#include <isl_aff_private.h>
+#include <isl_options_private.h>
+#include <isl_config.h>
/*
* The implementation of parametric integer linear programming in this file
* 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,
* then the initial sample value may be chosen equal to zero.
* However, we will not make this assumption. Instead, we apply
* the "big parameter" trick. Any variable x is then not directly
- * used in the tableau, but instead it its represented by another
+ * used in the tableau, but instead it is represented by another
* variable x' = M + x, where M is an arbitrarily large (positive)
* value. x' is therefore always non-negative, whatever the value of x.
- * Taking as initial smaple value x' = 0 corresponds to x = -M,
+ * Taking as initial sample value x' = 0 corresponds to x = -M,
* which is always smaller than any possible value of x.
*
* The big parameter trick is used in the main tableau and
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 */
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;
return;
error:
isl_basic_set_free(dom);
+ isl_mat_free(M);
sol->error = 1;
}
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);
}
} else
sol_pop_one(sol);
+
+ if (0)
+error: sol->error = 1;
}
static void sol_dec_level(struct isl_sol *sol)
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)
if (tab->empty && !sol->add_empty)
return;
+ if (sol->context->op->is_empty(sol->context))
+ return;
bset = sol_domain(sol);
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
{
int i;
struct isl_basic_map *bmap = NULL;
- isl_basic_set *context_bset;
unsigned n_eq;
unsigned n_ineq;
unsigned nparam;
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),
+ bmap = isl_basic_map_alloc_space(isl_map_get_space(sol->map),
n_div, n_eq, 2 * n_div + n_ineq);
if (!bmap)
goto error;
bmap = isl_basic_map_finalize(bmap);
sol->map = isl_map_grow(sol->map, 1);
sol->map = isl_map_add_basic_map(sol->map, bmap);
- if (!sol->map)
- goto error;
isl_basic_set_free(dom);
isl_mat_free(M);
+ if (!sol->map)
+ sol->sol.error = 1;
return;
error:
isl_basic_set_free(dom);
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
*
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;
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;
}
/* 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
*
} 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;
}
/* Return the first known violated constraint, i.e., a non-negative
- * contraint that currently has an either obviously negative value
+ * constraint that currently has an either obviously negative value
* or a previously determined to be negative value.
*
* If any constraint has a negative coefficient for the big parameter,
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 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)
+{
+ 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 lexicographicallly
+ * look for a pivoting column that would result in the lexicographically
* smallest increment in the sample point. If there is no such column
* then the tableau is infeasible.
*/
-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;
- } 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
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;
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.
* 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;
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);
- 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 (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;
return i < tab->n_sample;
}
-/* Add a div specifed by "div" to the tableau "tab" and return
+/* Add a div specified by "div" to the tableau "tab" and return
* 1 if the div is obviously non-negative.
*/
static int context_tab_add_div(struct isl_tab *tab, struct isl_vec *div,
n = tab->n_div;
d = context->op->get_div(context, tab, div);
+ isl_vec_free(div);
if (d < 0)
return -1;
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)
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)
return get_div(tab, context, div);
}
+/* Add a div specified by "div" to the context tableau and return
+ * 1 if the div is obviously non-negative.
+ * context_tab_add_div will always return 1, because all variables
+ * in a isl_context_lex tableau are non-negative.
+ * However, if we are using a big parameter in the context, then this only
+ * reflects the non-negativity of the variable used to _encode_ the
+ * div, i.e., div' = M + div, so we can't draw any conclusions.
+ */
static int context_lex_add_div(struct isl_context *context, struct isl_vec *div)
{
struct isl_context_lex *clex = (struct isl_context_lex *)context;
- return context_tab_add_div(clex->tab, div,
+ int nonneg;
+ nonneg = context_tab_add_div(clex->tab, div,
context_lex_add_ineq_wrap, context);
+ if (nonneg < 0)
+ return -1;
+ if (clex->tab->M)
+ return 0;
+ return nonneg;
}
static int context_lex_detect_equalities(struct isl_context *context,
}
}
+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;
context_lex_is_ok,
context_lex_save,
context_lex_restore,
+ context_lex_discard,
context_lex_invalidate,
context_lex_free,
};
{
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);
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;
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;
}
}
- 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]);
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)
static struct isl_tab *add_gbr_eq(struct isl_tab *tab, isl_int *eq)
{
- int r;
-
if (!tab)
return NULL;
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)
{
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;
{
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;
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)
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;
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);
+ cgbr->cone->bmap = isl_basic_map_extend_space(cgbr->cone->bmap,
+ isl_basic_map_get_space(cgbr->cone->bmap), 1, 0, 2);
k = isl_basic_map_alloc_div(cgbr->cone->bmap);
if (k < 0)
return -1;
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;
context_gbr_is_ok,
context_gbr_save,
context_gbr_restore,
+ context_gbr_discard,
context_gbr_invalidate,
context_gbr_free,
};
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;
* a minimization problem, which means that the variables in the
* tableau have value "M - x" rather than "M + x".
*/
-static struct isl_sol_map *sol_map_init(struct isl_basic_map *bmap,
+static struct isl_sol *sol_map_init(struct isl_basic_map *bmap,
struct isl_basic_set *dom, int track_empty, int max)
{
struct isl_sol_map *sol_map = NULL;
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,
+ sol_map->map = isl_map_alloc_space(isl_basic_map_get_space(bmap), 1,
ISL_MAP_DISJOINT);
if (!sol_map->map)
goto error;
goto error;
if (track_empty) {
- sol_map->empty = isl_set_alloc_dim(isl_basic_set_get_dim(dom),
+ sol_map->empty = isl_set_alloc_space(isl_basic_set_get_space(dom),
1, ISL_SET_DISJOINT);
if (!sol_map->empty)
goto error;
}
isl_basic_set_free(dom);
- return sol_map;
+ return &sol_map->sol;
error:
isl_basic_set_free(dom);
sol_map_free(sol_map);
if (!sol->error)
sol->context->op->restore(sol->context, saved);
+ else
+ sol->context->op->discard(saved);
return;
error:
sol->error = 1;
* coefficient are integral, then there is nothing that can be done
* and the tableau has no integral solution.
* If, on the other hand, one or more of the other columns have rational
- * coeffcients, but the parameter coefficients are all integral, then
+ * coefficients, but the parameter coefficients are all integral, then
* we can perform a regular (non-parametric) cut.
* Finally, if there is any parameter coefficient that is non-integral,
* then we need to involve the context tableau. There are two cases here.
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);
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".
*
* 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;
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);
sol->error = 1;
}
-static void sol_map_find_solutions(struct isl_sol_map *sol_map,
- struct isl_tab *tab)
-{
- find_solutions_main(&sol_map->sol, tab);
-}
-
/* Check if integer division "div" of "dom" also occurs in "bmap".
* If so, return its position within the divs.
* If not, return -1.
struct isl_basic_set *dom, unsigned div)
{
int i;
- unsigned b_dim = isl_dim_total(bmap->dim);
- unsigned d_dim = isl_dim_total(dom->dim);
+ unsigned b_dim = isl_space_dim(bmap->dim, isl_dim_all);
+ unsigned d_dim = isl_space_dim(dom->dim, isl_dim_all);
if (isl_int_is_zero(dom->div[div][0]))
return -1;
common++;
other = bmap->n_div - common;
if (dom->n_div - common > 0) {
- bmap = isl_basic_map_extend_dim(bmap, isl_dim_copy(bmap->dim),
+ bmap = isl_basic_map_extend_space(bmap, isl_space_copy(bmap->dim),
dom->n_div - common, 0, 0);
if (!bmap)
return NULL;
return NULL;
}
-/* Compute the lexicographic minimum (or maximum if "max" is set)
- * of "bmap" over the domain "dom" and return the result as a map.
- * If "empty" is not NULL, then *empty is assigned a set that
- * contains those parts of the domain where there is no solution.
- * If "bmap" is marked as rational (ISL_BASIC_MAP_RATIONAL),
- * then we compute the rational optimum. Otherwise, we compute
- * the integral optimum.
+/* Base case of isl_tab_basic_map_partial_lexopt, after removing
+ * some obvious symmetries.
*
- * We perform some preprocessing. As the PILP solver does not
- * handle implicit equalities very well, we first make sure all
- * the equalities are explicitly available.
- * We also make sure the divs in the domain are properly order,
+ * We make sure the divs in the domain are properly ordered,
* because they will be added one by one in the given order
* during the construction of the solution map.
*/
-struct isl_map *isl_tab_basic_map_partial_lexopt(
- struct isl_basic_map *bmap, struct isl_basic_set *dom,
- struct isl_set **empty, int max)
+static struct isl_sol *basic_map_partial_lexopt_base(
+ __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
+ __isl_give isl_set **empty, int max,
+ struct isl_sol *(*init)(__isl_keep isl_basic_map *bmap,
+ __isl_take isl_basic_set *dom, int track_empty, int max))
{
struct isl_tab *tab;
- struct isl_map *result = NULL;
- struct isl_sol_map *sol_map = NULL;
+ struct isl_sol *sol = NULL;
struct isl_context *context;
- struct isl_basic_map *eq;
-
- if (empty)
- *empty = NULL;
- if (!bmap || !dom)
- goto error2;
-
- isl_assert(bmap->ctx,
- isl_basic_map_compatible_domain(bmap, dom), goto error2);
-
- 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);
bmap = align_context_divs(bmap, dom);
}
- sol_map = sol_map_init(bmap, dom, !!empty, max);
- if (!sol_map)
+ sol = init(bmap, dom, !!empty, max);
+ if (!sol)
goto error;
- context = sol_map->sol.context;
- if (isl_basic_set_fast_is_empty(context->op->peek_basic_set(context)))
+ context = sol->context;
+ if (isl_basic_set_plain_is_empty(context->op->peek_basic_set(context)))
/* nothing */;
- else if (isl_basic_map_fast_is_empty(bmap))
- sol_map_add_empty_if_needed(sol_map,
+ else if (isl_basic_map_plain_is_empty(bmap)) {
+ if (sol->add_empty)
+ sol->add_empty(sol,
isl_basic_set_copy(context->op->peek_basic_set(context)));
- else {
+ } else {
tab = tab_for_lexmin(bmap,
context->op->peek_basic_set(context), 1, max);
tab = context->op->detect_nonnegative_parameters(context, tab);
- sol_map_find_solutions(sol_map, tab);
+ find_solutions_main(sol, tab);
}
- if (sol_map->sol.error)
+ if (sol->error)
goto error;
- result = isl_map_copy(sol_map->map);
- if (empty)
- *empty = isl_set_copy(sol_map->empty);
- sol_free(&sol_map->sol);
isl_basic_map_free(bmap);
- return result;
-error2:
- isl_basic_set_free(dom);
+ return sol;
error:
- sol_free(&sol_map->sol);
+ sol_free(sol);
isl_basic_map_free(bmap);
return NULL;
}
-struct isl_sol_for {
- struct isl_sol sol;
- int (*fn)(__isl_take isl_basic_set *dom,
- __isl_take isl_mat *map, void *user);
- void *user;
-};
-
+/* Base case of isl_tab_basic_map_partial_lexopt, after removing
+ * some obvious symmetries.
+ *
+ * We call basic_map_partial_lexopt_base and extract the results.
+ */
+static __isl_give isl_map *basic_map_partial_lexopt_base_map(
+ __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
+ __isl_give isl_set **empty, int max)
+{
+ isl_map *result = NULL;
+ struct isl_sol *sol;
+ struct isl_sol_map *sol_map;
+
+ sol = basic_map_partial_lexopt_base(bmap, dom, empty, max,
+ &sol_map_init);
+ if (!sol)
+ return NULL;
+ sol_map = (struct isl_sol_map *) sol;
+
+ result = isl_map_copy(sol_map->map);
+ if (empty)
+ *empty = isl_set_copy(sol_map->empty);
+ sol_free(&sol_map->sol);
+ return result;
+}
+
+/* Structure used during detection of parallel constraints.
+ * n_in: number of "input" variables: isl_dim_param + isl_dim_in
+ * n_out: number of "output" variables: isl_dim_out + isl_dim_div
+ * val: the coefficients of the output variables
+ */
+struct isl_constraint_equal_info {
+ isl_basic_map *bmap;
+ unsigned n_in;
+ unsigned n_out;
+ isl_int *val;
+};
+
+/* Check whether the coefficients of the output variables
+ * of the constraint in "entry" are equal to info->val.
+ */
+static int constraint_equal(const void *entry, const void *val)
+{
+ isl_int **row = (isl_int **)entry;
+ const struct isl_constraint_equal_info *info = val;
+
+ return isl_seq_eq((*row) + 1 + info->n_in, info->val, info->n_out);
+}
+
+/* Check whether "bmap" has a pair of constraints that have
+ * the same coefficients for the output variables.
+ * Note that the coefficients of the existentially quantified
+ * variables need to be zero since the existentially quantified
+ * of the result are usually not the same as those of the input.
+ * the isl_dim_out and isl_dim_div dimensions.
+ * If so, return 1 and return the row indices of the two constraints
+ * in *first and *second.
+ */
+static int parallel_constraints(__isl_keep isl_basic_map *bmap,
+ int *first, int *second)
+{
+ int i;
+ isl_ctx *ctx = isl_basic_map_get_ctx(bmap);
+ struct isl_hash_table *table = NULL;
+ struct isl_hash_table_entry *entry;
+ struct isl_constraint_equal_info info;
+ unsigned n_out;
+ unsigned n_div;
+
+ ctx = isl_basic_map_get_ctx(bmap);
+ table = isl_hash_table_alloc(ctx, bmap->n_ineq);
+ if (!table)
+ goto error;
+
+ info.n_in = isl_basic_map_dim(bmap, isl_dim_param) +
+ isl_basic_map_dim(bmap, isl_dim_in);
+ info.bmap = bmap;
+ n_out = isl_basic_map_dim(bmap, isl_dim_out);
+ n_div = isl_basic_map_dim(bmap, isl_dim_div);
+ info.n_out = n_out + n_div;
+ for (i = 0; i < bmap->n_ineq; ++i) {
+ uint32_t hash;
+
+ info.val = bmap->ineq[i] + 1 + info.n_in;
+ if (isl_seq_first_non_zero(info.val, n_out) < 0)
+ continue;
+ if (isl_seq_first_non_zero(info.val + n_out, n_div) >= 0)
+ continue;
+ hash = isl_seq_get_hash(info.val, info.n_out);
+ entry = isl_hash_table_find(ctx, table, hash,
+ constraint_equal, &info, 1);
+ if (!entry)
+ goto error;
+ if (entry->data)
+ break;
+ entry->data = &bmap->ineq[i];
+ }
+
+ if (i < bmap->n_ineq) {
+ *first = ((isl_int **)entry->data) - bmap->ineq;
+ *second = i;
+ }
+
+ isl_hash_table_free(ctx, table);
+
+ return i < bmap->n_ineq;
+error:
+ isl_hash_table_free(ctx, table);
+ return -1;
+}
+
+/* Given a set of upper bounds in "var", add constraints to "bset"
+ * that make the i-th bound smallest.
+ *
+ * In particular, if there are n bounds b_i, then add the constraints
+ *
+ * b_i <= b_j for j > i
+ * b_i < b_j for j < i
+ */
+static __isl_give isl_basic_set *select_minimum(__isl_take isl_basic_set *bset,
+ __isl_keep isl_mat *var, int i)
+{
+ isl_ctx *ctx;
+ int j, k;
+
+ ctx = isl_mat_get_ctx(var);
+
+ for (j = 0; j < var->n_row; ++j) {
+ if (j == i)
+ continue;
+ k = isl_basic_set_alloc_inequality(bset);
+ if (k < 0)
+ goto error;
+ isl_seq_combine(bset->ineq[k], ctx->one, var->row[j],
+ ctx->negone, var->row[i], var->n_col);
+ isl_int_set_si(bset->ineq[k][var->n_col], 0);
+ if (j < i)
+ isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
+ }
+
+ bset = isl_basic_set_finalize(bset);
+
+ return bset;
+error:
+ isl_basic_set_free(bset);
+ return NULL;
+}
+
+/* Given a set of upper bounds on the last "input" variable m,
+ * construct a set that assigns the minimal upper bound to m, i.e.,
+ * construct a set that divides the space into cells where one
+ * of the upper bounds is smaller than all the others and assign
+ * this upper bound to m.
+ *
+ * In particular, if there are n bounds b_i, then the result
+ * consists of n basic sets, each one of the form
+ *
+ * m = b_i
+ * b_i <= b_j for j > i
+ * b_i < b_j for j < i
+ */
+static __isl_give isl_set *set_minimum(__isl_take isl_space *dim,
+ __isl_take isl_mat *var)
+{
+ int i, k;
+ isl_basic_set *bset = NULL;
+ isl_ctx *ctx;
+ isl_set *set = NULL;
+
+ if (!dim || !var)
+ goto error;
+
+ ctx = isl_space_get_ctx(dim);
+ set = isl_set_alloc_space(isl_space_copy(dim),
+ var->n_row, ISL_SET_DISJOINT);
+
+ for (i = 0; i < var->n_row; ++i) {
+ bset = isl_basic_set_alloc_space(isl_space_copy(dim), 0,
+ 1, var->n_row - 1);
+ k = isl_basic_set_alloc_equality(bset);
+ if (k < 0)
+ goto error;
+ isl_seq_cpy(bset->eq[k], var->row[i], var->n_col);
+ isl_int_set_si(bset->eq[k][var->n_col], -1);
+ bset = select_minimum(bset, var, i);
+ set = isl_set_add_basic_set(set, bset);
+ }
+
+ isl_space_free(dim);
+ isl_mat_free(var);
+ return set;
+error:
+ isl_basic_set_free(bset);
+ isl_set_free(set);
+ isl_space_free(dim);
+ isl_mat_free(var);
+ return NULL;
+}
+
+/* Given that the last input variable of "bmap" represents the minimum
+ * of the bounds in "cst", check whether we need to split the domain
+ * based on which bound attains the minimum.
+ *
+ * A split is needed when the minimum appears in an integer division
+ * or in an equality. Otherwise, it is only needed if it appears in
+ * an upper bound that is different from the upper bounds on which it
+ * is defined.
+ */
+static int need_split_basic_map(__isl_keep isl_basic_map *bmap,
+ __isl_keep isl_mat *cst)
+{
+ int i, j;
+ unsigned total;
+ unsigned pos;
+
+ pos = cst->n_col - 1;
+ total = isl_basic_map_dim(bmap, isl_dim_all);
+
+ for (i = 0; i < bmap->n_div; ++i)
+ if (!isl_int_is_zero(bmap->div[i][2 + pos]))
+ return 1;
+
+ for (i = 0; i < bmap->n_eq; ++i)
+ if (!isl_int_is_zero(bmap->eq[i][1 + pos]))
+ return 1;
+
+ for (i = 0; i < bmap->n_ineq; ++i) {
+ if (isl_int_is_nonneg(bmap->ineq[i][1 + pos]))
+ continue;
+ if (!isl_int_is_negone(bmap->ineq[i][1 + pos]))
+ return 1;
+ if (isl_seq_first_non_zero(bmap->ineq[i] + 1 + pos + 1,
+ total - pos - 1) >= 0)
+ return 1;
+
+ for (j = 0; j < cst->n_row; ++j)
+ if (isl_seq_eq(bmap->ineq[i], cst->row[j], cst->n_col))
+ break;
+ if (j >= cst->n_row)
+ return 1;
+ }
+
+ return 0;
+}
+
+/* Given that the last set variable of "bset" represents the minimum
+ * of the bounds in "cst", check whether we need to split the domain
+ * based on which bound attains the minimum.
+ *
+ * We simply call need_split_basic_map here. This is safe because
+ * the position of the minimum is computed from "cst" and not
+ * from "bmap".
+ */
+static int need_split_basic_set(__isl_keep isl_basic_set *bset,
+ __isl_keep isl_mat *cst)
+{
+ return need_split_basic_map((isl_basic_map *)bset, cst);
+}
+
+/* Given that the last set variable of "set" represents the minimum
+ * of the bounds in "cst", check whether we need to split the domain
+ * based on which bound attains the minimum.
+ */
+static int need_split_set(__isl_keep isl_set *set, __isl_keep isl_mat *cst)
+{
+ int i;
+
+ for (i = 0; i < set->n; ++i)
+ if (need_split_basic_set(set->p[i], cst))
+ return 1;
+
+ return 0;
+}
+
+/* Given a set of which the last set variable is the minimum
+ * of the bounds in "cst", split each basic set in the set
+ * in pieces where one of the bounds is (strictly) smaller than the others.
+ * This subdivision is given in "min_expr".
+ * The variable is subsequently projected out.
+ *
+ * We only do the split when it is needed.
+ * For example if the last input variable m = min(a,b) and the only
+ * constraints in the given basic set are lower bounds on m,
+ * i.e., l <= m = min(a,b), then we can simply project out m
+ * to obtain l <= a and l <= b, without having to split on whether
+ * m is equal to a or b.
+ */
+static __isl_give isl_set *split(__isl_take isl_set *empty,
+ __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
+{
+ int n_in;
+ int i;
+ isl_space *dim;
+ isl_set *res;
+
+ if (!empty || !min_expr || !cst)
+ goto error;
+
+ n_in = isl_set_dim(empty, isl_dim_set);
+ dim = isl_set_get_space(empty);
+ dim = isl_space_drop_dims(dim, isl_dim_set, n_in - 1, 1);
+ res = isl_set_empty(dim);
+
+ for (i = 0; i < empty->n; ++i) {
+ isl_set *set;
+
+ set = isl_set_from_basic_set(isl_basic_set_copy(empty->p[i]));
+ if (need_split_basic_set(empty->p[i], cst))
+ set = isl_set_intersect(set, isl_set_copy(min_expr));
+ set = isl_set_remove_dims(set, isl_dim_set, n_in - 1, 1);
+
+ res = isl_set_union_disjoint(res, set);
+ }
+
+ isl_set_free(empty);
+ isl_set_free(min_expr);
+ isl_mat_free(cst);
+ return res;
+error:
+ isl_set_free(empty);
+ isl_set_free(min_expr);
+ isl_mat_free(cst);
+ return NULL;
+}
+
+/* Given a map of which the last input variable is the minimum
+ * of the bounds in "cst", split each basic set in the set
+ * in pieces where one of the bounds is (strictly) smaller than the others.
+ * This subdivision is given in "min_expr".
+ * The variable is subsequently projected out.
+ *
+ * The implementation is essentially the same as that of "split".
+ */
+static __isl_give isl_map *split_domain(__isl_take isl_map *opt,
+ __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
+{
+ int n_in;
+ int i;
+ isl_space *dim;
+ isl_map *res;
+
+ if (!opt || !min_expr || !cst)
+ goto error;
+
+ n_in = isl_map_dim(opt, isl_dim_in);
+ dim = isl_map_get_space(opt);
+ dim = isl_space_drop_dims(dim, isl_dim_in, n_in - 1, 1);
+ res = isl_map_empty(dim);
+
+ for (i = 0; i < opt->n; ++i) {
+ isl_map *map;
+
+ map = isl_map_from_basic_map(isl_basic_map_copy(opt->p[i]));
+ if (need_split_basic_map(opt->p[i], cst))
+ map = isl_map_intersect_domain(map,
+ isl_set_copy(min_expr));
+ map = isl_map_remove_dims(map, isl_dim_in, n_in - 1, 1);
+
+ res = isl_map_union_disjoint(res, map);
+ }
+
+ isl_map_free(opt);
+ isl_set_free(min_expr);
+ isl_mat_free(cst);
+ return res;
+error:
+ isl_map_free(opt);
+ isl_set_free(min_expr);
+ isl_mat_free(cst);
+ return NULL;
+}
+
+static __isl_give isl_map *basic_map_partial_lexopt(
+ __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
+ __isl_give isl_set **empty, int max);
+
+union isl_lex_res {
+ void *p;
+ isl_map *map;
+ isl_pw_multi_aff *pma;
+};
+
+/* This function is called from basic_map_partial_lexopt_symm.
+ * The last variable of "bmap" and "dom" corresponds to the minimum
+ * of the bounds in "cst". "map_space" is the space of the original
+ * input relation (of basic_map_partial_lexopt_symm) and "set_space"
+ * is the space of the original domain.
+ *
+ * We recursively call basic_map_partial_lexopt and then plug in
+ * the definition of the minimum in the result.
+ */
+static __isl_give union isl_lex_res basic_map_partial_lexopt_symm_map_core(
+ __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
+ __isl_give isl_set **empty, int max, __isl_take isl_mat *cst,
+ __isl_take isl_space *map_space, __isl_take isl_space *set_space)
+{
+ isl_map *opt;
+ isl_set *min_expr;
+ union isl_lex_res res;
+
+ min_expr = set_minimum(isl_basic_set_get_space(dom), isl_mat_copy(cst));
+
+ opt = basic_map_partial_lexopt(bmap, dom, empty, max);
+
+ if (empty) {
+ *empty = split(*empty,
+ isl_set_copy(min_expr), isl_mat_copy(cst));
+ *empty = isl_set_reset_space(*empty, set_space);
+ }
+
+ opt = split_domain(opt, min_expr, cst);
+ opt = isl_map_reset_space(opt, map_space);
+
+ res.map = opt;
+ return res;
+}
+
+/* Given a basic map with at least two parallel constraints (as found
+ * by the function parallel_constraints), first look for more constraints
+ * parallel to the two constraint and replace the found list of parallel
+ * constraints by a single constraint with as "input" part the minimum
+ * of the input parts of the list of constraints. Then, recursively call
+ * basic_map_partial_lexopt (possibly finding more parallel constraints)
+ * and plug in the definition of the minimum in the result.
+ *
+ * More specifically, given a set of constraints
+ *
+ * a x + b_i(p) >= 0
+ *
+ * Replace this set by a single constraint
+ *
+ * a x + u >= 0
+ *
+ * with u a new parameter with constraints
+ *
+ * u <= b_i(p)
+ *
+ * Any solution to the new system is also a solution for the original system
+ * since
+ *
+ * a x >= -u >= -b_i(p)
+ *
+ * Moreover, m = min_i(b_i(p)) satisfies the constraints on u and can
+ * therefore be plugged into the solution.
+ */
+static union isl_lex_res basic_map_partial_lexopt_symm(
+ __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
+ __isl_give isl_set **empty, int max, int first, int second,
+ __isl_give union isl_lex_res (*core)(__isl_take isl_basic_map *bmap,
+ __isl_take isl_basic_set *dom,
+ __isl_give isl_set **empty,
+ int max, __isl_take isl_mat *cst,
+ __isl_take isl_space *map_space,
+ __isl_take isl_space *set_space))
+{
+ int i, n, k;
+ int *list = NULL;
+ unsigned n_in, n_out, n_div;
+ isl_ctx *ctx;
+ isl_vec *var = NULL;
+ isl_mat *cst = NULL;
+ isl_space *map_space, *set_space;
+ union isl_lex_res res;
+
+ map_space = isl_basic_map_get_space(bmap);
+ set_space = empty ? isl_basic_set_get_space(dom) : NULL;
+
+ n_in = isl_basic_map_dim(bmap, isl_dim_param) +
+ isl_basic_map_dim(bmap, isl_dim_in);
+ n_out = isl_basic_map_dim(bmap, isl_dim_all) - n_in;
+
+ ctx = isl_basic_map_get_ctx(bmap);
+ list = isl_alloc_array(ctx, int, bmap->n_ineq);
+ var = isl_vec_alloc(ctx, n_out);
+ if (!list || !var)
+ goto error;
+
+ list[0] = first;
+ list[1] = second;
+ isl_seq_cpy(var->el, bmap->ineq[first] + 1 + n_in, n_out);
+ for (i = second + 1, n = 2; i < bmap->n_ineq; ++i) {
+ if (isl_seq_eq(var->el, bmap->ineq[i] + 1 + n_in, n_out))
+ list[n++] = i;
+ }
+
+ cst = isl_mat_alloc(ctx, n, 1 + n_in);
+ if (!cst)
+ goto error;
+
+ for (i = 0; i < n; ++i)
+ isl_seq_cpy(cst->row[i], bmap->ineq[list[i]], 1 + n_in);
+
+ bmap = isl_basic_map_cow(bmap);
+ if (!bmap)
+ goto error;
+ for (i = n - 1; i >= 0; --i)
+ if (isl_basic_map_drop_inequality(bmap, list[i]) < 0)
+ goto error;
+
+ bmap = isl_basic_map_add(bmap, isl_dim_in, 1);
+ bmap = isl_basic_map_extend_constraints(bmap, 0, 1);
+ k = isl_basic_map_alloc_inequality(bmap);
+ if (k < 0)
+ goto error;
+ isl_seq_clr(bmap->ineq[k], 1 + n_in);
+ isl_int_set_si(bmap->ineq[k][1 + n_in], 1);
+ isl_seq_cpy(bmap->ineq[k] + 1 + n_in + 1, var->el, n_out);
+ bmap = isl_basic_map_finalize(bmap);
+
+ n_div = isl_basic_set_dim(dom, isl_dim_div);
+ dom = isl_basic_set_add_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);
+ if (k < 0)
+ goto error;
+ isl_seq_cpy(dom->ineq[k], cst->row[i], 1 + n_in);
+ isl_int_set_si(dom->ineq[k][1 + n_in], -1);
+ isl_seq_clr(dom->ineq[k] + 1 + n_in + 1, n_div);
+ }
+
+ isl_vec_free(var);
+ free(list);
+
+ return core(bmap, dom, empty, max, cst, map_space, set_space);
+error:
+ isl_space_free(map_space);
+ isl_space_free(set_space);
+ isl_mat_free(cst);
+ isl_vec_free(var);
+ free(list);
+ isl_basic_set_free(dom);
+ isl_basic_map_free(bmap);
+ res.p = NULL;
+ return res;
+}
+
+static __isl_give isl_map *basic_map_partial_lexopt_symm_map(
+ __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
+ __isl_give isl_set **empty, int max, int first, int second)
+{
+ return basic_map_partial_lexopt_symm(bmap, dom, empty, max,
+ first, second, &basic_map_partial_lexopt_symm_map_core).map;
+}
+
+/* Recursive part of isl_tab_basic_map_partial_lexopt, after detecting
+ * equalities and removing redundant constraints.
+ *
+ * We first check if there are any parallel constraints (left).
+ * If not, we are in the base case.
+ * If there are parallel constraints, we replace them by a single
+ * constraint in basic_map_partial_lexopt_symm and then call
+ * this function recursively to look for more parallel constraints.
+ */
+static __isl_give isl_map *basic_map_partial_lexopt(
+ __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
+ __isl_give isl_set **empty, int max)
+{
+ int par = 0;
+ int first, second;
+
+ if (!bmap)
+ goto error;
+
+ if (bmap->ctx->opt->pip_symmetry)
+ par = parallel_constraints(bmap, &first, &second);
+ if (par < 0)
+ goto error;
+ if (!par)
+ return basic_map_partial_lexopt_base_map(bmap, dom, empty, max);
+
+ return basic_map_partial_lexopt_symm_map(bmap, dom, empty, max,
+ first, second);
+error:
+ isl_basic_set_free(dom);
+ isl_basic_map_free(bmap);
+ return NULL;
+}
+
+/* Compute the lexicographic minimum (or maximum if "max" is set)
+ * of "bmap" over the domain "dom" and return the result as a map.
+ * If "empty" is not NULL, then *empty is assigned a set that
+ * contains those parts of the domain where there is no solution.
+ * If "bmap" is marked as rational (ISL_BASIC_MAP_RATIONAL),
+ * then we compute the rational optimum. Otherwise, we compute
+ * the integral optimum.
+ *
+ * We perform some preprocessing. As the PILP solver does not
+ * handle implicit equalities very well, we first make sure all
+ * the equalities are explicitly available.
+ *
+ * We also add context constraints to the basic map and remove
+ * redundant constraints. This is only needed because of the
+ * way we handle simple symmetries. In particular, we currently look
+ * for symmetries on the constraints, before we set up the main tableau.
+ * It is then no good to look for symmetries on possibly redundant constraints.
+ */
+struct isl_map *isl_tab_basic_map_partial_lexopt(
+ struct isl_basic_map *bmap, struct isl_basic_set *dom,
+ struct isl_set **empty, int max)
+{
+ if (empty)
+ *empty = NULL;
+ if (!bmap || !dom)
+ goto error;
+
+ isl_assert(bmap->ctx,
+ isl_basic_map_compatible_domain(bmap, dom), goto error);
+
+ 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);
+
+ return basic_map_partial_lexopt(bmap, dom, empty, max);
+error:
+ isl_basic_set_free(dom);
+ isl_basic_map_free(bmap);
+ return NULL;
+}
+
+struct isl_sol_for {
+ struct isl_sol sol;
+ int (*fn)(__isl_take isl_basic_set *dom,
+ __isl_take isl_aff_list *list, void *user);
+ void *user;
+};
+
static void sol_for_free(struct isl_sol_for *sol_for)
{
- if (sol_for->sol.context)
- sol_for->sol.context->op->free(sol_for->sol.context);
- free(sol_for);
+ if (sol_for->sol.context)
+ sol_for->sol.context->op->free(sol_for->sol.context);
+ free(sol_for);
+}
+
+static void sol_for_free_wrap(struct isl_sol *sol)
+{
+ sol_for_free((struct isl_sol_for *)sol);
+}
+
+/* Add the solution identified by the tableau and the context tableau.
+ *
+ * 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
+ * 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;
+
+ 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(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);
+
+ dom = isl_basic_set_finalize(dom);
+
+ if (sol->fn(isl_basic_set_copy(dom), list, sol->user) < 0)
+ goto error;
+
+ isl_basic_set_free(dom);
+ isl_mat_free(M);
+ return;
+error:
+ isl_basic_set_free(dom);
+ isl_mat_free(M);
+ sol->sol.error = 1;
+}
+
+static void sol_for_add_wrap(struct isl_sol *sol,
+ struct isl_basic_set *dom, struct isl_mat *M)
+{
+ 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,
+ int (*fn)(__isl_take isl_basic_set *dom, __isl_take isl_aff_list *list,
+ void *user),
+ void *user)
+{
+ struct isl_sol_for *sol_for = NULL;
+ isl_space *dom_dim;
+ struct isl_basic_set *dom = NULL;
+
+ sol_for = isl_calloc_type(bmap->ctx, struct isl_sol_for);
+ if (!sol_for)
+ goto error;
+
+ dom_dim = isl_space_domain(isl_space_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->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;
+
+ sol_for->sol.context = isl_context_alloc(dom);
+ if (!sol_for->sol.context)
+ goto error;
+
+ isl_basic_set_free(dom);
+ return sol_for;
+error:
+ isl_basic_set_free(dom);
+ sol_for_free(sol_for);
+ return NULL;
+}
+
+static void sol_for_find_solutions(struct isl_sol_for *sol_for,
+ struct isl_tab *tab)
+{
+ find_solutions_main(&sol_for->sol, tab);
+}
+
+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_aff_list *list,
+ void *user),
+ void *user)
+{
+ struct isl_sol_for *sol_for = NULL;
+
+ bmap = isl_basic_map_copy(bmap);
+ bmap = isl_basic_map_detect_equalities(bmap);
+ if (!bmap)
+ return -1;
+
+ sol_for = sol_for_init(bmap, max, fn, user);
+ if (!sol_for)
+ goto error;
+
+ if (isl_basic_map_plain_is_empty(bmap))
+ /* nothing */;
+ else {
+ struct isl_tab *tab;
+ struct isl_context *context = sol_for->sol.context;
+ tab = tab_for_lexmin(bmap,
+ 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_free(&sol_for->sol);
+ isl_basic_map_free(bmap);
+ return 0;
+error:
+ sol_free(&sol_for->sol);
+ isl_basic_map_free(bmap);
+ 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);
+}
+
+/* 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_vec *v = NULL;
+ 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;
+ 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, CUT_ONE);
+ 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;
+
+ if (!bset)
+ return NULL;
+
+ tab = tab_for_lexmin(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;
}
-static void sol_for_free_wrap(struct isl_sol *sol)
+struct isl_sol_pma {
+ struct isl_sol sol;
+ isl_pw_multi_aff *pma;
+ isl_set *empty;
+};
+
+static void sol_pma_free(struct isl_sol_pma *sol_pma)
{
- sol_for_free((struct isl_sol_for *)sol);
+ if (!sol_pma)
+ return;
+ if (sol_pma->sol.context)
+ sol_pma->sol.context->op->free(sol_pma->sol.context);
+ isl_pw_multi_aff_free(sol_pma->pma);
+ isl_set_free(sol_pma->empty);
+ free(sol_pma);
}
-/* Add the solution identified by the tableau and the context tableau.
- *
- * 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
- * an affine matrix reprenting 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.)
+/* 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_for_add(struct isl_sol_for *sol,
- struct isl_basic_set *dom, struct isl_mat *M)
+static void sol_pma_add_empty(struct isl_sol_pma *sol,
+ __isl_take isl_basic_set *bset)
{
- if (sol->sol.error || !dom || !M)
+ 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, bset);
+ if (!sol->empty)
+ sol->sol.error = 1;
+ return;
+error:
+ isl_basic_set_free(bset);
+ sol->sol.error = 1;
+}
+
+/* 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 an isl_pw_multi_aff with a single
+ * cell corresponding to "dom" and affine expressions copied from "M".
+ */
+static void sol_pma_add(struct isl_sol_pma *sol,
+ __isl_take isl_basic_set *dom, __isl_take isl_mat *M)
+{
+ int i;
+ isl_local_space *ls;
+ isl_aff *aff;
+ isl_multi_aff *maff;
+ isl_pw_multi_aff *pma;
+
+ maff = isl_multi_aff_alloc(isl_pw_multi_aff_get_space(sol->pma));
+ ls = isl_basic_set_get_local_space(dom);
+ for (i = 1; i < M->n_row; ++i) {
+ aff = isl_aff_alloc(isl_local_space_copy(ls));
+ if (aff) {
+ 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);
+ maff = isl_multi_aff_set_aff(maff, i - 1, aff);
+ }
+ 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)
+ sol->sol.error = 1;
+}
- if (sol->fn(isl_basic_set_copy(dom), isl_mat_copy(M), sol->user) < 0)
- goto error;
+static void sol_pma_free_wrap(struct isl_sol *sol)
+{
+ sol_pma_free((struct isl_sol_pma *)sol);
+}
- isl_basic_set_free(dom);
- isl_mat_free(M);
- return;
-error:
- isl_basic_set_free(dom);
- isl_mat_free(M);
- sol->sol.error = 1;
+static void sol_pma_add_empty_wrap(struct isl_sol *sol,
+ __isl_take isl_basic_set *bset)
+{
+ sol_pma_add_empty((struct isl_sol_pma *)sol, bset);
}
-static void sol_for_add_wrap(struct isl_sol *sol,
- struct isl_basic_set *dom, struct isl_mat *M)
+static void sol_pma_add_wrap(struct isl_sol *sol,
+ __isl_take isl_basic_set *dom, __isl_take isl_mat *M)
{
- sol_for_add((struct isl_sol_for *)sol, dom, M);
+ sol_pma_add((struct isl_sol_pma *)sol, dom, M);
}
-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,
- void *user),
- void *user)
+/* Construct an isl_sol_pma 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.
+ * If max is set, then we are solving a maximization, rather than
+ * a minimization problem, which means that the variables in the
+ * tableau have value "M - x" rather than "M + x".
+ */
+static struct isl_sol *sol_pma_init(__isl_keep isl_basic_map *bmap,
+ __isl_take isl_basic_set *dom, int track_empty, int max)
{
- struct isl_sol_for *sol_for = NULL;
- struct isl_dim *dom_dim;
- struct isl_basic_set *dom = NULL;
+ struct isl_sol_pma *sol_pma = NULL;
- sol_for = isl_calloc_type(bmap->ctx, struct isl_sol_for);
- if (!sol_for)
+ if (!bmap)
goto error;
- dom_dim = isl_dim_domain(isl_dim_copy(bmap->dim));
- dom = isl_basic_set_universe(dom_dim);
+ sol_pma = isl_calloc_type(bmap->ctx, struct isl_sol_pma);
+ if (!sol_pma)
+ goto error;
- 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->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;
+ sol_pma->sol.rational = ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL);
+ sol_pma->sol.dec_level.callback.run = &sol_dec_level_wrap;
+ sol_pma->sol.dec_level.sol = &sol_pma->sol;
+ sol_pma->sol.max = max;
+ sol_pma->sol.n_out = isl_basic_map_dim(bmap, isl_dim_out);
+ sol_pma->sol.add = &sol_pma_add_wrap;
+ sol_pma->sol.add_empty = track_empty ? &sol_pma_add_empty_wrap : NULL;
+ sol_pma->sol.free = &sol_pma_free_wrap;
+ sol_pma->pma = isl_pw_multi_aff_empty(isl_basic_map_get_space(bmap));
+ if (!sol_pma->pma)
+ goto error;
- sol_for->sol.context = isl_context_alloc(dom);
- if (!sol_for->sol.context)
+ sol_pma->sol.context = isl_context_alloc(dom);
+ if (!sol_pma->sol.context)
goto error;
+ if (track_empty) {
+ sol_pma->empty = isl_set_alloc_space(isl_basic_set_get_space(dom),
+ 1, ISL_SET_DISJOINT);
+ if (!sol_pma->empty)
+ goto error;
+ }
+
isl_basic_set_free(dom);
- return sol_for;
+ return &sol_pma->sol;
error:
isl_basic_set_free(dom);
- sol_for_free(sol_for);
+ sol_pma_free(sol_pma);
return NULL;
}
-static void sol_for_find_solutions(struct isl_sol_for *sol_for,
- struct isl_tab *tab)
+/* Base case of isl_tab_basic_map_partial_lexopt, after removing
+ * some obvious symmetries.
+ *
+ * We call basic_map_partial_lexopt_base and extract the results.
+ */
+static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_base_pma(
+ __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
+ __isl_give isl_set **empty, int max)
{
- find_solutions_main(&sol_for->sol, tab);
+ isl_pw_multi_aff *result = NULL;
+ struct isl_sol *sol;
+ struct isl_sol_pma *sol_pma;
+
+ sol = basic_map_partial_lexopt_base(bmap, dom, empty, max,
+ &sol_pma_init);
+ if (!sol)
+ return NULL;
+ sol_pma = (struct isl_sol_pma *) sol;
+
+ result = isl_pw_multi_aff_copy(sol_pma->pma);
+ if (empty)
+ *empty = isl_set_copy(sol_pma->empty);
+ sol_free(&sol_pma->sol);
+ return result;
}
-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,
- void *user),
- void *user)
+/* Given that the last input variable of "maff" represents the minimum
+ * of some bounds, check whether we need to plug in the expression
+ * of the minimum.
+ *
+ * In particular, check if the last input variable appears in any
+ * of the expressions in "maff".
+ */
+static int need_substitution(__isl_keep isl_multi_aff *maff)
{
- struct isl_sol_for *sol_for = NULL;
+ int i;
+ unsigned pos;
- bmap = isl_basic_map_copy(bmap);
- if (!bmap)
- return -1;
+ pos = isl_multi_aff_dim(maff, isl_dim_in) - 1;
- bmap = isl_basic_map_detect_equalities(bmap);
- sol_for = sol_for_init(bmap, max, fn, user);
+ for (i = 0; i < maff->n; ++i)
+ if (isl_aff_involves_dims(maff->p[i], isl_dim_in, pos, 1))
+ return 1;
- if (isl_basic_map_fast_is_empty(bmap))
- /* nothing */;
- else {
- struct isl_tab *tab;
- struct isl_context *context = sol_for->sol.context;
- tab = tab_for_lexmin(bmap,
- 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)
+ return 0;
+}
+
+/* Given a set of upper bounds on the last "input" variable m,
+ * construct a piecewise affine expression that selects
+ * the minimal upper bound to m, i.e.,
+ * divide the space into cells where one
+ * of the upper bounds is smaller than all the others and select
+ * this upper bound on that cell.
+ *
+ * In particular, if there are n bounds b_i, then the result
+ * consists of n cell, each one of the form
+ *
+ * b_i <= b_j for j > i
+ * b_i < b_j for j < i
+ *
+ * The affine expression on this cell is
+ *
+ * b_i
+ */
+static __isl_give isl_pw_aff *set_minimum_pa(__isl_take isl_space *space,
+ __isl_take isl_mat *var)
+{
+ int i;
+ isl_aff *aff = NULL;
+ isl_basic_set *bset = NULL;
+ isl_ctx *ctx;
+ isl_pw_aff *paff = NULL;
+ isl_space *pw_space;
+ isl_local_space *ls = NULL;
+
+ if (!space || !var)
+ goto error;
+
+ ctx = isl_space_get_ctx(space);
+ ls = isl_local_space_from_space(isl_space_copy(space));
+ pw_space = isl_space_copy(space);
+ pw_space = isl_space_from_domain(pw_space);
+ pw_space = isl_space_add_dims(pw_space, isl_dim_out, 1);
+ paff = isl_pw_aff_alloc_size(pw_space, var->n_row);
+
+ for (i = 0; i < var->n_row; ++i) {
+ isl_pw_aff *paff_i;
+
+ aff = isl_aff_alloc(isl_local_space_copy(ls));
+ bset = isl_basic_set_alloc_space(isl_space_copy(space), 0,
+ 0, var->n_row - 1);
+ if (!aff || !bset)
goto error;
+ isl_int_set_si(aff->v->el[0], 1);
+ isl_seq_cpy(aff->v->el + 1, var->row[i], var->n_col);
+ isl_int_set_si(aff->v->el[1 + var->n_col], 0);
+ bset = select_minimum(bset, var, i);
+ paff_i = isl_pw_aff_alloc(isl_set_from_basic_set(bset), aff);
+ paff = isl_pw_aff_add_disjoint(paff, paff_i);
}
- sol_free(&sol_for->sol);
- isl_basic_map_free(bmap);
- return 0;
+ isl_local_space_free(ls);
+ isl_space_free(space);
+ isl_mat_free(var);
+ return paff;
error:
- sol_free(&sol_for->sol);
- isl_basic_map_free(bmap);
- return -1;
+ isl_aff_free(aff);
+ isl_basic_set_free(bset);
+ isl_pw_aff_free(paff);
+ isl_local_space_free(ls);
+ isl_space_free(space);
+ isl_mat_free(var);
+ return NULL;
}
-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,
- void *user),
- void *user)
+/* Given a piecewise multi-affine expression of which the last input variable
+ * is the minimum of the bounds in "cst", plug in the value of the minimum.
+ * This minimum expression is given in "min_expr_pa".
+ * The set "min_expr" contains the same information, but in the form of a set.
+ * The variable is subsequently projected out.
+ *
+ * The implementation is similar to those of "split" and "split_domain".
+ * If the variable appears in a given expression, then minimum expression
+ * is plugged in. Otherwise, if the variable appears in the constraints
+ * and a split is required, then the domain is split. Otherwise, no split
+ * is performed.
+ */
+static __isl_give isl_pw_multi_aff *split_domain_pma(
+ __isl_take isl_pw_multi_aff *opt, __isl_take isl_pw_aff *min_expr_pa,
+ __isl_take isl_set *min_expr, __isl_take isl_mat *cst)
{
- return isl_basic_map_foreach_lexopt(bmap, 0, fn, user);
+ int n_in;
+ int i;
+ isl_space *space;
+ isl_pw_multi_aff *res;
+
+ if (!opt || !min_expr || !cst)
+ goto error;
+
+ n_in = isl_pw_multi_aff_dim(opt, isl_dim_in);
+ space = isl_pw_multi_aff_get_space(opt);
+ space = isl_space_drop_dims(space, isl_dim_in, n_in - 1, 1);
+ res = isl_pw_multi_aff_empty(space);
+
+ for (i = 0; i < opt->n; ++i) {
+ isl_pw_multi_aff *pma;
+
+ pma = isl_pw_multi_aff_alloc(isl_set_copy(opt->p[i].set),
+ isl_multi_aff_copy(opt->p[i].maff));
+ if (need_substitution(opt->p[i].maff))
+ pma = isl_pw_multi_aff_substitute(pma,
+ isl_dim_in, n_in - 1, min_expr_pa);
+ else if (need_split_set(opt->p[i].set, cst))
+ pma = isl_pw_multi_aff_intersect_domain(pma,
+ isl_set_copy(min_expr));
+ pma = isl_pw_multi_aff_project_out(pma,
+ isl_dim_in, n_in - 1, 1);
+
+ res = isl_pw_multi_aff_add_disjoint(res, pma);
+ }
+
+ isl_pw_multi_aff_free(opt);
+ isl_pw_aff_free(min_expr_pa);
+ isl_set_free(min_expr);
+ isl_mat_free(cst);
+ return res;
+error:
+ isl_pw_multi_aff_free(opt);
+ isl_pw_aff_free(min_expr_pa);
+ isl_set_free(min_expr);
+ isl_mat_free(cst);
+ return NULL;
}
-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,
- void *user),
- void *user)
+static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_pma(
+ __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
+ __isl_give isl_set **empty, int max);
+
+/* This function is called from basic_map_partial_lexopt_symm.
+ * The last variable of "bmap" and "dom" corresponds to the minimum
+ * of the bounds in "cst". "map_space" is the space of the original
+ * input relation (of basic_map_partial_lexopt_symm) and "set_space"
+ * is the space of the original domain.
+ *
+ * We recursively call basic_map_partial_lexopt and then plug in
+ * the definition of the minimum in the result.
+ */
+static __isl_give union isl_lex_res basic_map_partial_lexopt_symm_pma_core(
+ __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
+ __isl_give isl_set **empty, int max, __isl_take isl_mat *cst,
+ __isl_take isl_space *map_space, __isl_take isl_space *set_space)
+{
+ isl_pw_multi_aff *opt;
+ isl_pw_aff *min_expr_pa;
+ isl_set *min_expr;
+ union isl_lex_res res;
+
+ min_expr = set_minimum(isl_basic_set_get_space(dom), isl_mat_copy(cst));
+ min_expr_pa = set_minimum_pa(isl_basic_set_get_space(dom),
+ isl_mat_copy(cst));
+
+ opt = basic_map_partial_lexopt_pma(bmap, dom, empty, max);
+
+ if (empty) {
+ *empty = split(*empty,
+ isl_set_copy(min_expr), isl_mat_copy(cst));
+ *empty = isl_set_reset_space(*empty, set_space);
+ }
+
+ opt = split_domain_pma(opt, min_expr_pa, min_expr, cst);
+ opt = isl_pw_multi_aff_reset_space(opt, map_space);
+
+ res.pma = opt;
+ return res;
+}
+
+static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_symm_pma(
+ __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
+ __isl_give isl_set **empty, int max, int first, int second)
+{
+ return basic_map_partial_lexopt_symm(bmap, dom, empty, max,
+ first, second, &basic_map_partial_lexopt_symm_pma_core).pma;
+}
+
+/* Recursive part of isl_basic_map_partial_lexopt_pw_multi_aff, after detecting
+ * equalities and removing redundant constraints.
+ *
+ * We first check if there are any parallel constraints (left).
+ * If not, we are in the base case.
+ * If there are parallel constraints, we replace them by a single
+ * constraint in basic_map_partial_lexopt_symm_pma and then call
+ * this function recursively to look for more parallel constraints.
+ */
+static __isl_give isl_pw_multi_aff *basic_map_partial_lexopt_pma(
+ __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
+ __isl_give isl_set **empty, int max)
+{
+ int par = 0;
+ int first, second;
+
+ if (!bmap)
+ goto error;
+
+ if (bmap->ctx->opt->pip_symmetry)
+ par = parallel_constraints(bmap, &first, &second);
+ if (par < 0)
+ goto error;
+ if (!par)
+ return basic_map_partial_lexopt_base_pma(bmap, dom, empty, max);
+
+ return basic_map_partial_lexopt_symm_pma(bmap, dom, empty, max,
+ first, second);
+error:
+ isl_basic_set_free(dom);
+ isl_basic_map_free(bmap);
+ return NULL;
+}
+
+/* Compute the lexicographic minimum (or maximum if "max" is set)
+ * of "bmap" over the domain "dom" and return the result as a piecewise
+ * multi-affine expression.
+ * If "empty" is not NULL, then *empty is assigned a set that
+ * contains those parts of the domain where there is no solution.
+ * If "bmap" is marked as rational (ISL_BASIC_MAP_RATIONAL),
+ * then we compute the rational optimum. Otherwise, we compute
+ * the integral optimum.
+ *
+ * We perform some preprocessing. As the PILP solver does not
+ * handle implicit equalities very well, we first make sure all
+ * the equalities are explicitly available.
+ *
+ * We also add context constraints to the basic map and remove
+ * redundant constraints. This is only needed because of the
+ * way we handle simple symmetries. In particular, we currently look
+ * for symmetries on the constraints, before we set up the main tableau.
+ * It is then no good to look for symmetries on possibly redundant constraints.
+ */
+__isl_give isl_pw_multi_aff *isl_basic_map_partial_lexopt_pw_multi_aff(
+ __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom,
+ __isl_give isl_set **empty, int max)
{
- return isl_basic_map_foreach_lexopt(bmap, 1, fn, user);
+ if (empty)
+ *empty = NULL;
+ if (!bmap || !dom)
+ goto error;
+
+ 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_pma(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);
+
+ return basic_map_partial_lexopt_pma(bmap, dom, empty, max);
+error:
+ isl_basic_set_free(dom);
+ isl_basic_map_free(bmap);
+ return NULL;
}
projects / platform / upstream / isl.git / blobdiff