1 /* 2 * Copyright 2008-2009 Katholieke Universiteit Leuven 3 * 4 * Use of this software is governed by the MIT license 5 * 6 * Written by Sven Verdoolaege, K.U.Leuven, Departement 7 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium 8 */ 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 17 struct 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_MP 36 #define GBR_type mpq_t 37 #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_MP 53 #include <imrat.h> 54 55 #define GBR_type mp_rat 56 #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_OPT 64 #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 90 static struct tab_lp *init_lp(struct isl_tab *tab); 91 static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim); 92 static int solve_lp(struct tab_lp *lp); 93 static void get_obj_val(struct tab_lp* lp, GBR_type *F); 94 static void delete_lp(struct tab_lp *lp); 95 static int add_lp_row(struct tab_lp *lp, isl_int *row, int dim); 96 static void get_alpha(struct tab_lp* lp, int row, GBR_type *alpha); 97 static int del_lp_row(struct tab_lp *lp) WARN_UNUSED; 98 static int cut_lp_to_hyperplane(struct tab_lp *lp, isl_int *row); 99 100 #define GBR_LP struct tab_lp 101 #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->neq 107 #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_fixed 111 #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 bset 116 * and then performing the Cartesian product on the tableau. 117 */ 118 static 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 135 static 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; 169 error: 170 delete_lp(lp); 171 return NULL; 172 } 173 174 static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim) 175 { 176 lp->obj = row; 177 } 178 179 static 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 in 217 * both halves of the tableau. 218 * Return 1 if the resulting tableau is empty. 219 */ 220 static 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 246 static 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 252 static 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 268 static 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 282 static 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 289 static 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
Indexes created Sun Sep 06 00:23:55 UTC 2026