From 77d8153e0094cc3addfacd8847537a79b4ffb8c7 Mon Sep 17 00:00:00 2001 From: Sven Verdoolaege Date: Mon, 15 Jun 2009 12:02:54 +0200 Subject: [PATCH] introduce generalized basis reduction based sampling Also make gbr the default for computing sample integer points. It may be a bit slower than pip on easy problems, but pip can sometimes get into serious problems, especially in cases where the given set does not contain any integer points. The gbr method performs much better in these cases. --- Makefile.am | 6 +- basis_reduction_tab.c | 184 +++++++++++++++ basis_reduction_templ.c | 216 ++++++++++++++++++ include/isl_ctx.h.in | 11 + isl_basis_reduction.h | 17 ++ isl_ctx.c | 8 +- isl_sample.c | 595 +++++++++++++++++++++++++++++++++++++++++++++++- 7 files changed, 1034 insertions(+), 3 deletions(-) create mode 100644 basis_reduction_tab.c create mode 100644 basis_reduction_templ.c create mode 100644 isl_basis_reduction.h diff --git a/Makefile.am b/Makefile.am index db9270b0..172d369d 100644 --- a/Makefile.am +++ b/Makefile.am @@ -42,6 +42,8 @@ libisl_la_SOURCES = \ $(ISL_POLYLIB) \ $(GET_MEMORY_FUNCTIONS) \ isl_affine_hull.c \ + isl_basis_reduction.h \ + basis_reduction_tab.c \ isl_blk.c \ isl_coalesce.c \ isl_constraint.c \ @@ -117,4 +119,6 @@ pkginclude_HEADERS = \ include/isl_set.h \ include/isl_set_polylib.h -EXTRA_DIST = test_inputs +EXTRA_DIST = \ + basis_reduction_templ.c \ + test_inputs diff --git a/basis_reduction_tab.c b/basis_reduction_tab.c new file mode 100644 index 00000000..ecc379e1 --- /dev/null +++ b/basis_reduction_tab.c @@ -0,0 +1,184 @@ +#include +#include "isl_tab.h" + +struct tab_lp { + struct isl_ctx *ctx; + struct isl_vec *row; + struct isl_tab *tab; + struct isl_tab_undo **stack; + isl_int *obj; + isl_int opt; + isl_int opt_denom; + int neq; + unsigned dim; + int n_ineq; +}; + +static struct tab_lp *init_lp(struct isl_basic_set *bset); +static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim); +static int solve_lp(struct tab_lp *lp); +static void get_obj_val(struct tab_lp* lp, mpq_t *F); +static void delete_lp(struct tab_lp *lp); +static int add_lp_row(struct tab_lp *lp, isl_int *row, int dim); +static void get_alpha(struct tab_lp* lp, int row, mpq_t *alpha); +static void del_lp_row(struct tab_lp *lp); + +#define GBR_LP struct tab_lp +#define GBR_type mpq_t +#define GBR_init(v) mpq_init(v) +#define GBR_clear(v) mpq_clear(v) +#define GBR_set(a,b) mpq_set(a,b) +#define GBR_set_ui(a,b) mpq_set_ui(a,b,1) +#define GBR_mul(a,b,c) mpq_mul(a,b,c) +#define GBR_lt(a,b) (mpq_cmp(a,b) < 0) +#define GBR_floor(a,b) mpz_fdiv_q(a,mpq_numref(b),mpq_denref(b)) +#define GBR_ceil(a,b) mpz_cdiv_q(a,mpq_numref(b),mpq_denref(b)) +#define GBR_lp_init(P) init_lp(P) +#define GBR_lp_set_obj(lp, obj, dim) set_lp_obj(lp, obj, dim) +#define GBR_lp_solve(lp) solve_lp(lp) +#define GBR_lp_get_obj_val(lp, F) get_obj_val(lp, F) +#define GBR_lp_delete(lp) delete_lp(lp) +#define GBR_lp_next_row(lp) lp->neq +#define GBR_lp_add_row(lp, row, dim) add_lp_row(lp, row, dim) +#define GBR_lp_get_alpha(lp, row, alpha) get_alpha(lp, row, alpha) +#define GBR_lp_del_row(lp) del_lp_row(lp); +#include "basis_reduction_templ.c" + +/* Set up a tableau for the Cartesian product of bset with itself. + * This could be optimized by first setting up a tableau for bset + * and then performing the Cartesian product on the tableau. + */ +static struct isl_tab *gbr_tab(struct isl_basic_set *bset, + struct isl_vec *row) +{ + int i, j; + unsigned dim; + struct isl_tab *tab; + + if (!bset || !row) + return NULL; + + dim = isl_basic_set_total_dim(bset); + tab = isl_tab_alloc(bset->ctx, 2 * bset->n_ineq + dim + 1, 2 * dim); + + for (i = 0; i < 2; ++i) { + isl_seq_clr(row->el + 1 + (1 - i) * dim, dim); + for (j = 0; j < bset->n_ineq; ++j) { + isl_int_set(row->el[0], bset->ineq[j][0]); + isl_seq_cpy(row->el + 1 + i * dim, + bset->ineq[j] + 1, dim); + tab = isl_tab_add_ineq(bset->ctx, tab, row->el); + if (!tab || tab->empty) + return tab; + } + } + + return tab; +} + +static struct tab_lp *init_lp(struct isl_basic_set *bset) +{ + struct tab_lp *lp = NULL; + + if (!bset) + return NULL; + + isl_assert(bset->ctx, bset->n_eq == 0, return NULL); + + lp = isl_calloc_type(bset->ctx, struct tab_lp); + if (!lp) + return NULL; + + isl_int_init(lp->opt); + isl_int_init(lp->opt_denom); + + lp->dim = isl_basic_set_total_dim(bset); + lp->n_ineq = bset->n_ineq; + + lp->ctx = bset->ctx; + isl_ctx_ref(lp->ctx); + + lp->stack = isl_alloc_array(lp->ctx, struct isl_tab_undo *, lp->dim); + + lp->row = isl_vec_alloc(lp->ctx, 1 + 2 * lp->dim); + if (!lp->row) + goto error; + lp->tab = gbr_tab(bset, lp->row); + if (!lp->tab) + goto error; + lp->obj = NULL; + lp->neq = 0; + + return lp; +error: + delete_lp(lp); + return NULL; +} + +static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim) +{ + lp->obj = row; +} + +static int solve_lp(struct tab_lp *lp) +{ + enum isl_lp_result res; + unsigned flags = 0; + + isl_int_set_si(lp->row->el[0], 0); + isl_seq_cpy(lp->row->el + 1, lp->obj, lp->dim); + isl_seq_neg(lp->row->el + 1 + lp->dim, lp->obj, lp->dim); + if (lp->neq) + flags = ISL_TAB_SAVE_DUAL; + res = isl_tab_min(lp->ctx, lp->tab, lp->row->el, lp->ctx->one, + &lp->opt, &lp->opt_denom, flags); + if (res != isl_lp_ok) + return -1; + return 0; +} + +static void get_obj_val(struct tab_lp* lp, mpq_t *F) +{ + isl_int_neg(mpq_numref(*F), lp->opt); + isl_int_set(mpq_denref(*F), lp->opt_denom); +} + +static void delete_lp(struct tab_lp *lp) +{ + if (!lp) + return; + + isl_int_clear(lp->opt); + isl_int_clear(lp->opt_denom); + isl_vec_free(lp->row); + free(lp->stack); + isl_tab_free(lp->ctx, lp->tab); + isl_ctx_deref(lp->ctx); + free(lp); +} + +static int add_lp_row(struct tab_lp *lp, isl_int *row, int dim) +{ + lp->stack[lp->neq] = isl_tab_snap(lp->ctx, lp->tab); + + isl_int_set_si(lp->row->el[0], 0); + isl_seq_cpy(lp->row->el + 1, row, lp->dim); + isl_seq_neg(lp->row->el + 1 + lp->dim, row, lp->dim); + + lp->tab = isl_tab_add_valid_eq(lp->ctx, lp->tab, lp->row->el); + + return lp->neq++; +} + +static void get_alpha(struct tab_lp* lp, int row, mpq_t *alpha) +{ + row += 2 * lp->n_ineq; + isl_int_neg(mpq_numref(*alpha), lp->tab->dual->el[1 + row]); + isl_int_set(mpq_denref(*alpha), lp->tab->dual->el[0]); +} + +static void del_lp_row(struct tab_lp *lp) +{ + lp->neq--; + isl_tab_rollback(lp->ctx, lp->tab, lp->stack[lp->neq]); +} diff --git a/basis_reduction_templ.c b/basis_reduction_templ.c new file mode 100644 index 00000000..8aeb6c46 --- /dev/null +++ b/basis_reduction_templ.c @@ -0,0 +1,216 @@ +#include +#include "isl_basis_reduction.h" + +static void save_alpha(GBR_LP *lp, int first, int n, GBR_type *alpha) +{ + int i; + + for (i = 0; i < n; ++i) + GBR_lp_get_alpha(lp, first + i, &alpha[i]); +} + +/* This function implements the algorithm described in + * "An Implementation of the Generalized Basis Reduction Algorithm + * for Integer Programming" of Cook el al. to compute a reduced basis. + * We use \epsilon = 1/4. + * + * If options->gbr_only_first is set, the user is only interested + * in the first direction. In this case we stop the basis reduction when + * - the width in the first direction becomes smaller than 2 + * or + * - we have moved forward all the way to the last direction + * and then back again all the way to the first direction. + */ +struct isl_mat *isl_basic_set_reduced_basis(struct isl_basic_set *bset) +{ + unsigned dim; + struct isl_mat *basis; + int unbounded; + int i; + GBR_LP *lp = NULL; + GBR_type F_old, alpha, F_new; + int row; + isl_int tmp; + struct isl_vec *b_tmp; + GBR_type *F = NULL; + GBR_type *alpha_buffer[2] = { NULL, NULL }; + GBR_type *alpha_saved; + GBR_type F_saved; + int use_saved = 0; + isl_int mu[2]; + GBR_type mu_F[2]; + GBR_type two; + + if (!bset) + return NULL; + + dim = isl_basic_set_total_dim(bset); + basis = isl_mat_identity(bset->ctx, dim); + if (!basis) + return NULL; + + if (dim == 1) + return basis; + + isl_int_init(tmp); + isl_int_init(mu[0]); + isl_int_init(mu[1]); + + GBR_init(alpha); + GBR_init(F_old); + GBR_init(F_new); + GBR_init(F_saved); + GBR_init(mu_F[0]); + GBR_init(mu_F[1]); + GBR_init(two); + + b_tmp = isl_vec_alloc(bset->ctx, dim); + if (!b_tmp) + goto error; + + F = isl_alloc_array(bset->ctx, GBR_type, dim); + alpha_buffer[0] = isl_alloc_array(bset->ctx, GBR_type, dim); + alpha_buffer[1] = isl_alloc_array(bset->ctx, GBR_type, dim); + alpha_saved = alpha_buffer[0]; + + if (!F || !alpha_buffer[0] || !alpha_buffer[1]) + goto error; + + for (i = 0; i < dim; ++i) { + GBR_init(F[i]); + GBR_init(alpha_buffer[0][i]); + GBR_init(alpha_buffer[1][i]); + } + + GBR_set_ui(two, 2); + + lp = GBR_lp_init(bset); + if (!lp) + goto error; + + i = 0; + + GBR_lp_set_obj(lp, basis->row[0], dim); + bset->ctx->stats->gbr_solved_lps++; + unbounded = GBR_lp_solve(lp); + isl_assert(bset->ctx, !unbounded, goto error); + GBR_lp_get_obj_val(lp, &F[0]); + + do { + if (use_saved) { + row = GBR_lp_next_row(lp); + GBR_set(F_new, F_saved); + GBR_set(alpha, alpha_saved[i]); + } else { + row = GBR_lp_add_row(lp, basis->row[i], dim); + GBR_lp_set_obj(lp, basis->row[i+1], dim); + bset->ctx->stats->gbr_solved_lps++; + unbounded = GBR_lp_solve(lp); + isl_assert(bset->ctx, !unbounded, goto error); + GBR_lp_get_obj_val(lp, &F_new); + + GBR_lp_get_alpha(lp, row, &alpha); + + if (i > 0) + save_alpha(lp, row-i, i, alpha_saved); + + GBR_lp_del_row(lp); + } + GBR_set(F[i+1], F_new); + + GBR_floor(mu[0], alpha); + GBR_ceil(mu[1], alpha); + + if (isl_int_eq(mu[0], mu[1])) + isl_int_set(tmp, mu[0]); + else { + int j; + + for (j = 0; j <= 1; ++j) { + isl_int_set(tmp, mu[j]); + isl_seq_combine(b_tmp->el, + bset->ctx->one, basis->row[i+1], + tmp, basis->row[i], dim); + GBR_lp_set_obj(lp, b_tmp->el, dim); + bset->ctx->stats->gbr_solved_lps++; + unbounded = GBR_lp_solve(lp); + isl_assert(bset->ctx, !unbounded, goto error); + GBR_lp_get_obj_val(lp, &mu_F[j]); + if (i > 0) + save_alpha(lp, row-i, i, alpha_buffer[j]); + } + + if (GBR_lt(mu_F[0], mu_F[1])) + j = 0; + else + j = 1; + + isl_int_set(tmp, mu[j]); + GBR_set(F_new, mu_F[j]); + alpha_saved = alpha_buffer[j]; + } + isl_seq_combine(basis->row[i+1], + bset->ctx->one, basis->row[i+1], + tmp, basis->row[i], dim); + + GBR_set(F_old, F[i]); + + use_saved = 0; + /* mu_F[0] = 4 * F_new; mu_F[1] = 3 * F_old */ + GBR_set_ui(mu_F[0], 4); + GBR_mul(mu_F[0], mu_F[0], F_new); + GBR_set_ui(mu_F[1], 3); + GBR_mul(mu_F[1], mu_F[1], F_old); + if (GBR_lt(mu_F[0], mu_F[1])) { + basis = isl_mat_swap_rows(bset->ctx, basis, i, i + 1); + if (i > 0) { + use_saved = 1; + GBR_set(F_saved, F_new); + GBR_lp_del_row(lp); + --i; + } else { + GBR_set(F[0], F_new); + if (bset->ctx->gbr_only_first && + GBR_lt(F[0], two)) + break; + } + } else { + GBR_lp_add_row(lp, basis->row[i], dim); + ++i; + } + } while (i < dim-1); + + if (0) { +error: + isl_mat_free(bset->ctx, basis); + basis = NULL; + } + + GBR_lp_delete(lp); + + if (alpha_buffer[1]) + for (i = 0; i < dim; ++i) { + GBR_clear(F[i]); + GBR_clear(alpha_buffer[0][i]); + GBR_clear(alpha_buffer[1][i]); + } + free(F); + free(alpha_buffer[0]); + free(alpha_buffer[1]); + + isl_vec_free(b_tmp); + + GBR_clear(alpha); + GBR_clear(F_old); + GBR_clear(F_new); + GBR_clear(F_saved); + GBR_clear(mu_F[0]); + GBR_clear(mu_F[1]); + GBR_clear(two); + + isl_int_clear(tmp); + isl_int_clear(mu[0]); + isl_int_clear(mu[1]); + + return basis; +} diff --git a/include/isl_ctx.h.in b/include/isl_ctx.h.in index e59a6c26..aa86fe15 100644 --- a/include/isl_ctx.h.in +++ b/include/isl_ctx.h.in @@ -39,9 +39,14 @@ extern "C" { * (in case of pointer return type). * The only exception is the isl_ctx argument, which shoud never be NULL. */ +struct isl_stats { + long gbr_solved_lps; +}; struct isl_ctx { int ref; + struct isl_stats *stats; + isl_int one; isl_int negone; @@ -55,6 +60,12 @@ struct isl_ctx { #define ISL_LP_TAB 0 #define ISL_LP_PIP 1 unsigned lp_solver; + + #define ISL_ILP_GBR 0 + #define ISL_ILP_PIP 1 + unsigned ilp_solver; + + unsigned gbr_only_first; }; /* Some helper macros */ diff --git a/isl_basis_reduction.h b/isl_basis_reduction.h new file mode 100644 index 00000000..cc28bae4 --- /dev/null +++ b/isl_basis_reduction.h @@ -0,0 +1,17 @@ +#ifndef ISL_BASIS_REDUCTION_H +#define ISL_BASIS_REDUCTION_H + +#include "isl_set.h" +#include "isl_mat.h" + +#if defined(__cplusplus) +extern "C" { +#endif + +struct isl_mat *isl_basic_set_reduced_basis(struct isl_basic_set *bset); + +#if defined(__cplusplus) +} +#endif + +#endif diff --git a/isl_ctx.c b/isl_ctx.c index 9ee369fe..e9e8cff9 100644 --- a/isl_ctx.c +++ b/isl_ctx.c @@ -8,13 +8,17 @@ struct isl_ctx *isl_ctx_alloc() { struct isl_ctx *ctx = NULL; - ctx = isl_alloc_type(NULL, struct isl_ctx); + ctx = isl_calloc_type(NULL, struct isl_ctx); if (!ctx) goto error; if (isl_hash_table_init(ctx, &ctx->name_hash, 0)) goto error; + ctx->stats = isl_calloc_type(ctx, struct isl_stats); + if (!ctx->stats) + goto error; + ctx->ref = 0; isl_int_init(ctx->one); @@ -30,6 +34,7 @@ struct isl_ctx *isl_ctx_alloc() #endif ctx->lp_solver = ISL_LP_TAB; + ctx->ilp_solver = ISL_ILP_GBR; return ctx; error: @@ -57,5 +62,6 @@ void isl_ctx_free(struct isl_ctx *ctx) isl_blk_clear_cache(ctx); isl_int_clear(ctx->one); isl_int_clear(ctx->negone); + free(ctx->stats); free(ctx); } diff --git a/isl_sample.c b/isl_sample.c index b4bca122..562c4311 100644 --- a/isl_sample.c +++ b/isl_sample.c @@ -5,6 +5,8 @@ #include "isl_seq.h" #include "isl_map_private.h" #include "isl_equalities.h" +#include "isl_tab.h" +#include "isl_basis_reduction.h" static struct isl_vec *empty_sample(struct isl_basic_set *bset) { @@ -171,6 +173,589 @@ static struct isl_vec *sample_eq(struct isl_basic_set *bset, return sample; } +/* Given a basic set "bset" and an affine function "f"/"denom", + * check if bset is bounded and non-empty and if so, return the minimal + * and maximal value attained by the affine function in "min" and "max". + * The minimal value is rounded up to the nearest integer, while the + * maximal value is rounded down. + * The return value indicates whether the set was empty or unbounded. + */ +static enum isl_lp_result basic_set_range(struct isl_basic_set *bset, + isl_int *f, isl_int denom, isl_int *min, isl_int *max) +{ + unsigned dim; + struct isl_tab *tab; + enum isl_lp_result res; + + if (!bset) + return isl_lp_error; + if (isl_basic_set_fast_is_empty(bset)) + return isl_lp_empty; + + tab = isl_tab_from_basic_set(bset); + res = isl_tab_min(bset->ctx, tab, f, denom, min, NULL, 0); + if (res != isl_lp_ok) { + isl_tab_free(bset->ctx, tab); + return res; + } + dim = isl_basic_set_total_dim(bset); + isl_seq_neg(f, f, 1 + dim); + res = isl_tab_min(bset->ctx, tab, f, denom, max, NULL, 0); + isl_seq_neg(f, f, 1 + dim); + isl_int_neg(*max, *max); + + isl_tab_free(bset->ctx, tab); + return res; +} + +/* Perform a basis reduction on "bset" and return the inverse of + * the new basis, i.e., an affine mapping from the new coordinates to the old, + * in *T. + */ +static struct isl_basic_set *basic_set_reduced(struct isl_basic_set *bset, + struct isl_mat **T) +{ + struct isl_ctx *ctx; + unsigned gbr_only_first; + + *T = NULL; + if (!bset) + return NULL; + + ctx = bset->ctx; + + gbr_only_first = ctx->gbr_only_first; + ctx->gbr_only_first = 1; + *T = isl_basic_set_reduced_basis(bset); + ctx->gbr_only_first = gbr_only_first; + + *T = isl_mat_lin_to_aff(bset->ctx, *T); + *T = isl_mat_right_inverse(bset->ctx, *T); + + bset = isl_basic_set_preimage(bset, isl_mat_copy(bset->ctx, *T)); + if (!bset) + goto error; + + return bset; +error: + isl_mat_free(ctx, *T); + *T = NULL; + return NULL; +} + +static struct isl_vec *sample_bounded(struct isl_basic_set *bset); + +/* Given a basic set "bset" whose first coordinate ranges between + * "min" and "max", step through all values from min to max, until + * the slice of bset with the first coordinate fixed to one of these + * values contains an integer point. If such a point is found, return it. + * If none of the slices contains any integer point, then bset itself + * doesn't contain any integer point and an empty sample is returned. + */ +static struct isl_vec *sample_scan(struct isl_basic_set *bset, + isl_int min, isl_int max) +{ + unsigned total; + struct isl_basic_set *slice = NULL; + struct isl_vec *sample = NULL; + isl_int tmp; + + total = isl_basic_set_total_dim(bset); + + isl_int_init(tmp); + for (isl_int_set(tmp, min); isl_int_le(tmp, max); + isl_int_add_ui(tmp, tmp, 1)) { + int k; + + slice = isl_basic_set_copy(bset); + slice = isl_basic_set_cow(slice); + slice = isl_basic_set_extend_constraints(slice, 1, 0); + k = isl_basic_set_alloc_equality(slice); + if (k < 0) + goto error; + isl_int_set(slice->eq[k][0], tmp); + isl_int_set_si(slice->eq[k][1], -1); + isl_seq_clr(slice->eq[k] + 2, total - 1); + slice = isl_basic_set_simplify(slice); + sample = sample_bounded(slice); + slice = NULL; + if (!sample) + goto error; + if (sample->size > 0) + break; + isl_vec_free(sample); + sample = NULL; + } + if (!sample) + sample = empty_sample(bset); + else + isl_basic_set_free(bset); + isl_int_clear(tmp); + return sample; +error: + isl_basic_set_free(bset); + isl_basic_set_free(slice); + isl_int_clear(tmp); + return NULL; +} + +/* Given a basic set that is known to be bounded, find and return + * an integer point in the basic set, if there is any. + * + * After handling some trivial cases, we check the range of the + * first coordinate. If this coordinate can only attain one integer + * value, we are happy. Otherwise, we perform basis reduction and + * determine the new range. + * + * Then we step through all possible values in the range in sample_scan. + * + * If any basis reduction was performed, the sample value found, if any, + * is transformed back to the original space. + */ +static struct isl_vec *sample_bounded(struct isl_basic_set *bset) +{ + unsigned dim; + struct isl_ctx *ctx; + struct isl_vec *sample; + struct isl_vec *obj = NULL; + struct isl_mat *T = NULL; + isl_int min, max; + enum isl_lp_result res; + + if (!bset) + return NULL; + + if (isl_basic_set_fast_is_empty(bset)) + return empty_sample(bset); + + ctx = bset->ctx; + dim = isl_basic_set_total_dim(bset); + if (dim == 0) + return zero_sample(bset); + if (dim == 1) + return interval_sample(bset); + if (bset->n_eq > 0) + return sample_eq(bset, sample_bounded); + + isl_int_init(min); + isl_int_init(max); + obj = isl_vec_alloc(bset->ctx, 1 + dim); + if (!obj) + goto error; + isl_seq_clr(obj->el, 1+ dim); + isl_int_set_si(obj->el[1], 1); + + res = basic_set_range(bset, obj->el, bset->ctx->one, &min, &max); + if (res == isl_lp_error) + goto error; + isl_assert(bset->ctx, res != isl_lp_unbounded, goto error); + if (res == isl_lp_empty || isl_int_lt(max, min)) { + sample = empty_sample(bset); + goto out; + } + + if (isl_int_ne(min, max)) { + bset = basic_set_reduced(bset, &T); + if (!bset) + goto error; + + res = basic_set_range(bset, obj->el, bset->ctx->one, &min, &max); + if (res == isl_lp_error) + goto error; + isl_assert(bset->ctx, res != isl_lp_unbounded, goto error); + if (res == isl_lp_empty || isl_int_lt(max, min)) { + sample = empty_sample(bset); + goto out; + } + } + + sample = sample_scan(bset, min, max); +out: + if (T) { + if (!sample || sample->size == 0) + isl_mat_free(ctx, T); + else + sample = isl_mat_vec_product(ctx, T, sample); + } + isl_vec_free(obj); + isl_int_clear(min); + isl_int_clear(max); + return sample; +error: + isl_mat_free(ctx, T); + isl_basic_set_free(bset); + isl_vec_free(obj); + isl_int_clear(min); + isl_int_clear(max); + return NULL; +} + +/* Given a basic set "bset" and a value "sample" for the first coordinates + * of bset, plug in these values and drop the corresponding coordinates. + * + * We do this by computing the preimage of the transformation + * + * [ 1 0 ] + * x = [ s 0 ] x' + * [ 0 I ] + * + * where [1 s] is the sample value and I is the identity matrix of the + * appropriate dimension. + */ +static struct isl_basic_set *plug_in(struct isl_basic_set *bset, + struct isl_vec *sample) +{ + int i; + unsigned total; + struct isl_mat *T; + + if (!bset || !sample) + goto error; + + total = isl_basic_set_total_dim(bset); + T = isl_mat_alloc(bset->ctx, 1 + total, 1 + total - (sample->size - 1)); + if (!T) + goto error; + + for (i = 0; i < sample->size; ++i) { + isl_int_set(T->row[i][0], sample->el[i]); + isl_seq_clr(T->row[i] + 1, T->n_col - 1); + } + for (i = 0; i < T->n_col - 1; ++i) { + isl_seq_clr(T->row[sample->size + i], T->n_col); + isl_int_set_si(T->row[sample->size + i][1 + i], 1); + } + isl_vec_free(sample); + + bset = isl_basic_set_preimage(bset, T); + return bset; +error: + isl_basic_set_free(bset); + isl_vec_free(sample); + return NULL; +} + +/* Given a basic set "bset", return any (possibly non-integer) point + * in the basic set. + */ +static struct isl_vec *rational_sample(struct isl_basic_set *bset) +{ + struct isl_tab *tab; + struct isl_vec *sample; + + if (!bset) + return NULL; + + tab = isl_tab_from_basic_set(bset); + sample = isl_tab_get_sample_value(bset->ctx, tab); + isl_tab_free(bset->ctx, tab); + + isl_basic_set_free(bset); + + return sample; +} + +/* Given a rational vector, with the denominator in the first element + * of the vector, round up all coordinates. + */ +struct isl_vec *isl_vec_ceil(struct isl_vec *vec) +{ + int i; + + vec = isl_vec_cow(vec); + if (!vec) + return NULL; + + isl_seq_cdiv_q(vec->el + 1, vec->el + 1, vec->el[0], vec->size - 1); + + isl_int_set_si(vec->el[0], 1); + + return vec; +} + +/* Given a linear cone "cone" and a rational point "vec", + * construct a polyhedron with shifted copies of the constraints in "cone", + * i.e., a polyhedron with "cone" as its recession cone, such that each + * point x in this polyhedron is such that the unit box positioned at x + * lies entirely inside the affine cone 'vec + cone'. + * Any rational point in this polyhedron may therefore be rounded up + * to yield an integer point that lies inside said affine cone. + * + * Denote the constraints of cone by " >= 0" and the rational + * point "vec" by v/d. + * Let b_i = . Then the affine cone 'vec + cone' is given + * by - b/d >= 0. + * The polyhedron - ceil{b/d} >= 0 is a subset of this affine cone. + * We prefer this polyhedron over the actual affine cone because it doesn't + * require a scaling of the constraints. + * If each of the vertices of the unit cube positioned at x lies inside + * this polyhedron, then the whole unit cube at x lies inside the affine cone. + * We therefore impose that x' = x + \sum e_i, for any selection of unit + * vectors lies inside the polyhedron, i.e., + * + * - ceil{b/d} = + sum a_i - ceil{b/d} >= 0 + * + * The most stringent of these constraints is the one that selects + * all negative a_i, so the polyhedron we are looking for has constraints + * + * + sum_{a_i < 0} a_i - ceil{b/d} >= 0 + * + * Note that if cone were known to have only non-negative rays + * (which can be accomplished by a unimodular transformation), + * then we would only have to check the points x' = x + e_i + * and we only have to add the smallest negative a_i (if any) + * instead of the sum of all negative a_i. + */ +static struct isl_basic_set *shift_cone(struct isl_basic_set *cone, + struct isl_vec *vec) +{ + int i, j, k; + unsigned total; + + struct isl_basic_set *shift = NULL; + + if (!cone || !vec) + goto error; + + isl_assert(cone->ctx, cone->n_eq == 0, goto error); + + total = isl_basic_set_total_dim(cone); + + shift = isl_basic_set_alloc_dim(isl_basic_set_get_dim(cone), + 0, 0, cone->n_ineq); + + for (i = 0; i < cone->n_ineq; ++i) { + k = isl_basic_set_alloc_inequality(shift); + if (k < 0) + goto error; + isl_seq_cpy(shift->ineq[k] + 1, cone->ineq[i] + 1, total); + isl_seq_inner_product(shift->ineq[k] + 1, vec->el + 1, total, + &shift->ineq[k][0]); + isl_int_cdiv_q(shift->ineq[k][0], + shift->ineq[k][0], vec->el[0]); + isl_int_neg(shift->ineq[k][0], shift->ineq[k][0]); + for (j = 0; j < total; ++j) { + if (isl_int_is_nonneg(shift->ineq[k][1 + j])) + continue; + isl_int_add(shift->ineq[k][0], + shift->ineq[k][0], shift->ineq[k][1 + j]); + } + } + + isl_basic_set_free(cone); + isl_vec_free(vec); + + return isl_basic_set_finalize(shift); +error: + isl_basic_set_free(shift); + isl_basic_set_free(cone); + isl_vec_free(vec); + return NULL; +} + +/* Given a rational point vec in a (transformed) basic set, + * such that cone is the recession cone of the original basic set, + * "round up" the rational point to an integer point. + * + * We first check if the rational point just happens to be integer. + * If not, we transform the cone in the same way as the basic set, + * pick a point x in this cone shifted to the rational point such that + * the whole unit cube at x is also inside this affine cone. + * Then we simply round up the coordinates of x and return the + * resulting integer point. + */ +static struct isl_vec *round_up_in_cone(struct isl_vec *vec, + struct isl_basic_set *cone, struct isl_mat *U) +{ + unsigned total; + + if (!vec || !cone || !U) + goto error; + + isl_assert(vec->ctx, vec->size != 0, goto error); + if (isl_int_is_one(vec->el[0])) { + isl_mat_free(vec->ctx, U); + isl_basic_set_free(cone); + return vec; + } + + total = isl_basic_set_total_dim(cone); + cone = isl_basic_set_preimage(cone, U); + cone = isl_basic_set_remove_dims(cone, 0, total - (vec->size - 1)); + + cone = shift_cone(cone, vec); + + vec = rational_sample(cone); + vec = isl_vec_ceil(vec); + return vec; +error: + isl_mat_free(vec ? vec->ctx : cone ? cone->ctx : NULL, U); + isl_vec_free(vec); + isl_basic_set_free(cone); + return NULL; +} + +/* Concatenate two integer vectors, i.e., two vectors with denominator + * (stored in element 0) equal to 1. + */ +static struct isl_vec *vec_concat(struct isl_vec *vec1, struct isl_vec *vec2) +{ + struct isl_vec *vec; + + if (!vec1 || !vec2) + goto error; + isl_assert(vec1->ctx, vec1->size > 0, goto error); + isl_assert(vec2->ctx, vec2->size > 0, goto error); + isl_assert(vec1->ctx, isl_int_is_one(vec1->el[0]), goto error); + isl_assert(vec2->ctx, isl_int_is_one(vec2->el[0]), goto error); + + vec = isl_vec_alloc(vec1->ctx, vec1->size + vec2->size - 1); + if (!vec) + goto error; + + isl_seq_cpy(vec->el, vec1->el, vec1->size); + isl_seq_cpy(vec->el + vec1->size, vec2->el + 1, vec2->size - 1); + + isl_vec_free(vec1); + isl_vec_free(vec2); + + return vec; +error: + isl_vec_free(vec1); + isl_vec_free(vec2); + return NULL; +} + +/* Drop all constraints in bset that involve any of the dimensions + * first to first+n-1. + */ +static struct isl_basic_set *drop_constraints_involving + (struct isl_basic_set *bset, unsigned first, unsigned n) +{ + int i; + + if (!bset) + return NULL; + + bset = isl_basic_set_cow(bset); + + for (i = bset->n_ineq - 1; i >= 0; --i) { + if (isl_seq_first_non_zero(bset->ineq[i] + 1 + first, n) == -1) + continue; + isl_basic_set_drop_inequality(bset, i); + } + + return bset; +} + +/* Give a basic set "bset" with recession cone "cone", compute and + * return an integer point in bset, if any. + * + * If the recession cone is full-dimensional, then we know that + * bset contains an infinite number of integer points and it is + * fairly easy to pick one of them. + * If the recession cone is not full-dimensional, then we first + * transform bset such that the bounded directions appear as + * the first dimensions of the transformed basic set. + * We do this by using a unimodular transformation that transforms + * the equalities in the recession cone to equalities on the first + * dimensions. + * + * The transformed set is then projected onto its bounded dimensions. + * Note that to compute this projection, we can simply drop all constraints + * involving any of the unbounded dimensions since these constraints + * cannot be combined to produce a constraint on the bounded dimensions. + * To see this, assume that there is such a combination of constraints + * that produces a constraint on the bounded dimensions. This means + * that some combination of the unbounded dimensions has both an upper + * bound and a lower bound in terms of the bounded dimensions, but then + * this combination would be a bounded direction too and would have been + * transformed into a bounded dimensions. + * + * We then compute a sample value in the bounded dimensions. + * If no such value can be found, then the original set did not contain + * any integer points and we are done. + * Otherwise, we plug in the value we found in the bounded dimensions, + * project out these bounded dimensions and end up with a set with + * a full-dimensional recession cone. + * A sample point in this set is computed by "rounding up" any + * rational point in the set. + * + * The sample points in the bounded and unbounded dimensions are + * then combined into a single sample point and transformed back + * to the original space. + */ +static struct isl_vec *sample_with_cone(struct isl_basic_set *bset, + struct isl_basic_set *cone) +{ + unsigned total; + unsigned cone_dim; + struct isl_mat *M, *U; + struct isl_vec *sample; + struct isl_vec *cone_sample; + struct isl_ctx *ctx; + struct isl_basic_set *bounded; + + if (!bset || !cone) + goto error; + + ctx = bset->ctx; + total = isl_basic_set_total_dim(cone); + cone_dim = total - cone->n_eq; + + M = isl_mat_sub_alloc(bset->ctx, cone->eq, 0, cone->n_eq, 1, total); + M = isl_mat_left_hermite(bset->ctx, M, 0, &U, NULL); + if (!M) + goto error; + isl_mat_free(bset->ctx, M); + + U = isl_mat_lin_to_aff(bset->ctx, U); + bset = isl_basic_set_preimage(bset, isl_mat_copy(bset->ctx, U)); + + bounded = isl_basic_set_copy(bset); + bounded = drop_constraints_involving(bounded, total - cone_dim, cone_dim); + bounded = isl_basic_set_drop_dims(bounded, total - cone_dim, cone_dim); + sample = sample_bounded(bounded); + if (!sample || sample->size == 0) { + isl_basic_set_free(bset); + isl_basic_set_free(cone); + isl_mat_free(ctx, U); + return sample; + } + bset = plug_in(bset, isl_vec_copy(sample)); + cone_sample = rational_sample(bset); + cone_sample = round_up_in_cone(cone_sample, cone, isl_mat_copy(ctx, U)); + sample = vec_concat(sample, cone_sample); + sample = isl_mat_vec_product(ctx, U, sample); + return sample; +error: + isl_basic_set_free(cone); + isl_basic_set_free(bset); + return NULL; +} + +/* Compute and return a sample point in bset using generalized basis + * reduction. We first check if the input set has a non-trivial + * recession cone. If so, we perform some extra preprocessing in + * sample_with_cone. Otherwise, we directly perform generalized basis + * reduction. + */ +static struct isl_vec *gbr_sample_no_lineality(struct isl_basic_set *bset) +{ + unsigned dim; + struct isl_basic_set *cone; + + dim = isl_basic_set_total_dim(bset); + + cone = isl_basic_set_recession_cone(isl_basic_set_copy(bset)); + + if (cone->n_eq < dim) + return sample_with_cone(bset, cone); + + isl_basic_set_free(cone); + return sample_bounded(bset); +} + static struct isl_vec *sample_no_lineality(struct isl_basic_set *bset) { unsigned dim; @@ -185,7 +770,15 @@ static struct isl_vec *sample_no_lineality(struct isl_basic_set *bset) if (dim == 1) return interval_sample(bset); - return isl_pip_basic_set_sample(bset); + switch (bset->ctx->ilp_solver) { + case ISL_ILP_PIP: + return isl_pip_basic_set_sample(bset); + case ISL_ILP_GBR: + return gbr_sample_no_lineality(bset); + } + isl_assert(bset->ctx, 0, ); + isl_basic_set_free(bset); + return NULL; } /* Compute an integer point in "bset" with a lineality space that -- 2.11.4.GIT