294 lines · c
1/*2 * Copyright 2008-2009 Katholieke Universiteit Leuven3 *4 * Use of this software is governed by the MIT license5 *6 * Written by Sven Verdoolaege, K.U.Leuven, Departement7 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium8 */9 10#include <assert.h>11#include <isl_map_private.h>12#include <isl_seq.h>13#include "isl_tab.h"14#include <isl_int.h>15#include <isl_config.h>16 17struct tab_lp {18 struct isl_ctx *ctx;19 struct isl_vec *row;20 struct isl_tab *tab;21 struct isl_tab_undo **stack;22 isl_int *obj;23 isl_int opt;24 isl_int opt_denom;25 isl_int tmp;26 isl_int tmp2;27 int neq;28 unsigned dim;29 /* number of constraints in initial product tableau */30 int con_offset;31 /* objective function has fixed or no integer value */32 int is_fixed;33};34 35#ifdef USE_GMP_FOR_MP36#define GBR_type mpq_t37#define GBR_init(v) mpq_init(v)38#define GBR_clear(v) mpq_clear(v)39#define GBR_set(a,b) mpq_set(a,b)40#define GBR_set_ui(a,b) mpq_set_ui(a,b,1)41#define GBR_mul(a,b,c) mpq_mul(a,b,c)42#define GBR_lt(a,b) (mpq_cmp(a,b) < 0)43#define GBR_is_zero(a) (mpq_sgn(a) == 0)44#define GBR_numref(a) mpq_numref(a)45#define GBR_denref(a) mpq_denref(a)46#define GBR_floor(a,b) mpz_fdiv_q(a,GBR_numref(b),GBR_denref(b))47#define GBR_ceil(a,b) mpz_cdiv_q(a,GBR_numref(b),GBR_denref(b))48#define GBR_set_num_neg(a, b) mpz_neg(GBR_numref(*a), b);49#define GBR_set_den(a, b) mpz_set(GBR_denref(*a), b);50#endif /* USE_GMP_FOR_MP */51 52#ifdef USE_IMATH_FOR_MP53#include <imrat.h>54 55#define GBR_type mp_rat56#define GBR_init(v) v = mp_rat_alloc()57#define GBR_clear(v) mp_rat_free(v)58#define GBR_set(a,b) mp_rat_copy(b,a)59#define GBR_set_ui(a,b) mp_rat_set_uvalue(a,b,1)60#define GBR_mul(a,b,c) mp_rat_mul(b,c,a)61#define GBR_lt(a,b) (mp_rat_compare(a,b) < 0)62#define GBR_is_zero(a) (mp_rat_compare_zero(a) == 0)63#ifdef USE_SMALL_INT_OPT64#define GBR_numref(a) isl_sioimath_encode_big(mp_rat_numer_ref(a))65#define GBR_denref(a) isl_sioimath_encode_big(mp_rat_denom_ref(a))66#define GBR_floor(a, b) isl_sioimath_fdiv_q((a), GBR_numref(b), GBR_denref(b))67#define GBR_ceil(a, b) isl_sioimath_cdiv_q((a), GBR_numref(b), GBR_denref(b))68#define GBR_set_num_neg(a, b) \69 do { \70 isl_sioimath_scratchspace_t scratch; \71 impz_neg(mp_rat_numer_ref(*a), \72 isl_sioimath_bigarg_src(*b, &scratch));\73 } while (0)74#define GBR_set_den(a, b) \75 do { \76 isl_sioimath_scratchspace_t scratch; \77 impz_set(mp_rat_denom_ref(*a), \78 isl_sioimath_bigarg_src(*b, &scratch));\79 } while (0)80#else /* USE_SMALL_INT_OPT */81#define GBR_numref(a) mp_rat_numer_ref(a)82#define GBR_denref(a) mp_rat_denom_ref(a)83#define GBR_floor(a,b) impz_fdiv_q(a,GBR_numref(b),GBR_denref(b))84#define GBR_ceil(a,b) impz_cdiv_q(a,GBR_numref(b),GBR_denref(b))85#define GBR_set_num_neg(a, b) impz_neg(GBR_numref(*a), b)86#define GBR_set_den(a, b) impz_set(GBR_denref(*a), b)87#endif /* USE_SMALL_INT_OPT */88#endif /* USE_IMATH_FOR_MP */89 90static struct tab_lp *init_lp(struct isl_tab *tab);91static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim);92static int solve_lp(struct tab_lp *lp);93static void get_obj_val(struct tab_lp* lp, GBR_type *F);94static void delete_lp(struct tab_lp *lp);95static int add_lp_row(struct tab_lp *lp, isl_int *row, int dim);96static void get_alpha(struct tab_lp* lp, int row, GBR_type *alpha);97static int del_lp_row(struct tab_lp *lp) WARN_UNUSED;98static int cut_lp_to_hyperplane(struct tab_lp *lp, isl_int *row);99 100#define GBR_LP struct tab_lp101#define GBR_lp_init(P) init_lp(P)102#define GBR_lp_set_obj(lp, obj, dim) set_lp_obj(lp, obj, dim)103#define GBR_lp_solve(lp) solve_lp(lp)104#define GBR_lp_get_obj_val(lp, F) get_obj_val(lp, F)105#define GBR_lp_delete(lp) delete_lp(lp)106#define GBR_lp_next_row(lp) lp->neq107#define GBR_lp_add_row(lp, row, dim) add_lp_row(lp, row, dim)108#define GBR_lp_get_alpha(lp, row, alpha) get_alpha(lp, row, alpha)109#define GBR_lp_del_row(lp) del_lp_row(lp)110#define GBR_lp_is_fixed(lp) (lp)->is_fixed111#define GBR_lp_cut(lp, obj) cut_lp_to_hyperplane(lp, obj)112#include "basis_reduction_templ.c"113 114/* Set up a tableau for the Cartesian product of bset with itself.115 * This could be optimized by first setting up a tableau for bset116 * and then performing the Cartesian product on the tableau.117 */118static struct isl_tab *gbr_tab(struct isl_tab *tab, struct isl_vec *row)119{120 unsigned dim;121 struct isl_tab *prod;122 123 if (!tab || !row)124 return NULL;125 126 dim = tab->n_var;127 prod = isl_tab_product(tab, tab);128 if (isl_tab_extend_cons(prod, 3 * dim + 1) < 0) {129 isl_tab_free(prod);130 return NULL;131 }132 return prod;133}134 135static struct tab_lp *init_lp(struct isl_tab *tab)136{137 struct tab_lp *lp = NULL;138 139 if (!tab)140 return NULL;141 142 lp = isl_calloc_type(tab->mat->ctx, struct tab_lp);143 if (!lp)144 return NULL;145 146 isl_int_init(lp->opt);147 isl_int_init(lp->opt_denom);148 isl_int_init(lp->tmp);149 isl_int_init(lp->tmp2);150 151 lp->dim = tab->n_var;152 153 lp->ctx = tab->mat->ctx;154 isl_ctx_ref(lp->ctx);155 156 lp->stack = isl_alloc_array(lp->ctx, struct isl_tab_undo *, lp->dim);157 158 lp->row = isl_vec_alloc(lp->ctx, 1 + 2 * lp->dim);159 if (!lp->row)160 goto error;161 lp->tab = gbr_tab(tab, lp->row);162 if (!lp->tab)163 goto error;164 lp->con_offset = lp->tab->n_con;165 lp->obj = NULL;166 lp->neq = 0;167 168 return lp;169error:170 delete_lp(lp);171 return NULL;172}173 174static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim)175{176 lp->obj = row;177}178 179static int solve_lp(struct tab_lp *lp)180{181 enum isl_lp_result res;182 unsigned flags = 0;183 184 lp->is_fixed = 0;185 186 isl_int_set_si(lp->row->el[0], 0);187 isl_seq_cpy(lp->row->el + 1, lp->obj, lp->dim);188 isl_seq_neg(lp->row->el + 1 + lp->dim, lp->obj, lp->dim);189 if (lp->neq)190 flags = ISL_TAB_SAVE_DUAL;191 res = isl_tab_min(lp->tab, lp->row->el, lp->ctx->one,192 &lp->opt, &lp->opt_denom, flags);193 isl_int_mul_ui(lp->opt_denom, lp->opt_denom, 2);194 if (isl_int_abs_lt(lp->opt, lp->opt_denom)) {195 struct isl_vec *sample = isl_tab_get_sample_value(lp->tab);196 if (!sample)197 return -1;198 isl_seq_inner_product(lp->obj, sample->el + 1, lp->dim, &lp->tmp);199 isl_seq_inner_product(lp->obj, sample->el + 1 + lp->dim, lp->dim, &lp->tmp2);200 isl_int_cdiv_q(lp->tmp, lp->tmp, sample->el[0]);201 isl_int_fdiv_q(lp->tmp2, lp->tmp2, sample->el[0]);202 if (isl_int_ge(lp->tmp, lp->tmp2))203 lp->is_fixed = 1;204 isl_vec_free(sample);205 }206 isl_int_divexact_ui(lp->opt_denom, lp->opt_denom, 2);207 if (res < 0)208 return -1;209 if (res != isl_lp_ok)210 isl_die(lp->ctx, isl_error_internal,211 "unexpected missing (bounded) solution", return -1);212 return 0;213}214 215/* The current objective function has a fixed (or no) integer value.216 * Cut the tableau to the hyperplane that fixes this value in217 * both halves of the tableau.218 * Return 1 if the resulting tableau is empty.219 */220static int cut_lp_to_hyperplane(struct tab_lp *lp, isl_int *row)221{222 enum isl_lp_result res;223 224 isl_int_set_si(lp->row->el[0], 0);225 isl_seq_cpy(lp->row->el + 1, row, lp->dim);226 isl_seq_clr(lp->row->el + 1 + lp->dim, lp->dim);227 res = isl_tab_min(lp->tab, lp->row->el, lp->ctx->one,228 &lp->tmp, NULL, 0);229 if (res != isl_lp_ok)230 return -1;231 232 isl_int_neg(lp->row->el[0], lp->tmp);233 if (isl_tab_add_eq(lp->tab, lp->row->el) < 0)234 return -1;235 236 isl_seq_cpy(lp->row->el + 1 + lp->dim, row, lp->dim);237 isl_seq_clr(lp->row->el + 1, lp->dim);238 if (isl_tab_add_eq(lp->tab, lp->row->el) < 0)239 return -1;240 241 lp->con_offset += 2;242 243 return lp->tab->empty;244}245 246static void get_obj_val(struct tab_lp* lp, GBR_type *F)247{248 GBR_set_num_neg(F, lp->opt);249 GBR_set_den(F, lp->opt_denom);250}251 252static void delete_lp(struct tab_lp *lp)253{254 if (!lp)255 return;256 257 isl_int_clear(lp->opt);258 isl_int_clear(lp->opt_denom);259 isl_int_clear(lp->tmp);260 isl_int_clear(lp->tmp2);261 isl_vec_free(lp->row);262 free(lp->stack);263 isl_tab_free(lp->tab);264 isl_ctx_deref(lp->ctx);265 free(lp);266}267 268static int add_lp_row(struct tab_lp *lp, isl_int *row, int dim)269{270 lp->stack[lp->neq] = isl_tab_snap(lp->tab);271 272 isl_int_set_si(lp->row->el[0], 0);273 isl_seq_cpy(lp->row->el + 1, row, lp->dim);274 isl_seq_neg(lp->row->el + 1 + lp->dim, row, lp->dim);275 276 if (isl_tab_add_valid_eq(lp->tab, lp->row->el) < 0)277 return -1;278 279 return lp->neq++;280}281 282static void get_alpha(struct tab_lp* lp, int row, GBR_type *alpha)283{284 row += lp->con_offset;285 GBR_set_num_neg(alpha, lp->tab->dual->el[1 + row]);286 GBR_set_den(alpha, lp->tab->dual->el[0]);287}288 289static int del_lp_row(struct tab_lp *lp)290{291 lp->neq--;292 return isl_tab_rollback(lp->tab, lp->stack[lp->neq]);293}294