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 15struct tab_lp { 16 struct isl_ctx *ctx; 17 struct isl_vec *row; 18 struct isl_tab *tab; 19 struct isl_tab_undo **stack; 20 isl_int *obj; 21 isl_int opt; 22 isl_int opt_denom; 23 isl_int tmp; 24 isl_int tmp2; 25 int neq; 26 unsigned dim; 27 /* number of constraints in initial product tableau */ 28 int con_offset; 29 /* objective function has fixed or no integer value */ 30 int is_fixed; 31}; 32 33static struct tab_lp *init_lp(struct isl_tab *tab); 34static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim); 35static int solve_lp(struct tab_lp *lp); 36static void get_obj_val(struct tab_lp* lp, mpq_t *F); 37static void delete_lp(struct tab_lp *lp); 38static int add_lp_row(struct tab_lp *lp, isl_int *row, int dim); 39static void get_alpha(struct tab_lp* lp, int row, mpq_t *alpha); 40static int del_lp_row(struct tab_lp *lp) WARN_UNUSED; 41static int cut_lp_to_hyperplane(struct tab_lp *lp, isl_int *row); 42 43#define GBR_LP struct tab_lp 44#define GBR_type mpq_t 45#define GBR_init(v) mpq_init(v) 46#define GBR_clear(v) mpq_clear(v) 47#define GBR_set(a,b) mpq_set(a,b) 48#define GBR_set_ui(a,b) mpq_set_ui(a,b,1) 49#define GBR_mul(a,b,c) mpq_mul(a,b,c) 50#define GBR_lt(a,b) (mpq_cmp(a,b) < 0) 51#define GBR_is_zero(a) (mpq_sgn(a) == 0) 52#define GBR_floor(a,b) mpz_fdiv_q(a,mpq_numref(b),mpq_denref(b)) 53#define GBR_ceil(a,b) mpz_cdiv_q(a,mpq_numref(b),mpq_denref(b)) 54#define GBR_lp_init(P) init_lp(P) 55#define GBR_lp_set_obj(lp, obj, dim) set_lp_obj(lp, obj, dim) 56#define GBR_lp_solve(lp) solve_lp(lp) 57#define GBR_lp_get_obj_val(lp, F) get_obj_val(lp, F) 58#define GBR_lp_delete(lp) delete_lp(lp) 59#define GBR_lp_next_row(lp) lp->neq 60#define GBR_lp_add_row(lp, row, dim) add_lp_row(lp, row, dim) 61#define GBR_lp_get_alpha(lp, row, alpha) get_alpha(lp, row, alpha) 62#define GBR_lp_del_row(lp) del_lp_row(lp) 63#define GBR_lp_is_fixed(lp) (lp)->is_fixed 64#define GBR_lp_cut(lp, obj) cut_lp_to_hyperplane(lp, obj) 65#include "basis_reduction_templ.c" 66 67/* Set up a tableau for the Cartesian product of bset with itself. 68 * This could be optimized by first setting up a tableau for bset 69 * and then performing the Cartesian product on the tableau. 70 */ 71static struct isl_tab *gbr_tab(struct isl_tab *tab, struct isl_vec *row) 72{ 73 unsigned dim; 74 struct isl_tab *prod; 75 76 if (!tab || !row) 77 return NULL; 78 79 dim = tab->n_var; 80 prod = isl_tab_product(tab, tab); 81 if (isl_tab_extend_cons(prod, 3 * dim + 1) < 0) { 82 isl_tab_free(prod); 83 return NULL; 84 } 85 return prod; 86} 87 88static struct tab_lp *init_lp(struct isl_tab *tab) 89{ 90 struct tab_lp *lp = NULL; 91 92 if (!tab) 93 return NULL; 94 95 lp = isl_calloc_type(tab->mat->ctx, struct tab_lp); 96 if (!lp) 97 return NULL; 98 99 isl_int_init(lp->opt); 100 isl_int_init(lp->opt_denom); 101 isl_int_init(lp->tmp); 102 isl_int_init(lp->tmp2); 103 104 lp->dim = tab->n_var; 105 106 lp->ctx = tab->mat->ctx; 107 isl_ctx_ref(lp->ctx); 108 109 lp->stack = isl_alloc_array(lp->ctx, struct isl_tab_undo *, lp->dim); 110 111 lp->row = isl_vec_alloc(lp->ctx, 1 + 2 * lp->dim); 112 if (!lp->row) 113 goto error; 114 lp->tab = gbr_tab(tab, lp->row); 115 if (!lp->tab) 116 goto error; 117 lp->con_offset = lp->tab->n_con; 118 lp->obj = NULL; 119 lp->neq = 0; 120 121 return lp; 122error: 123 delete_lp(lp); 124 return NULL; 125} 126 127static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim) 128{ 129 lp->obj = row; 130} 131 132static int solve_lp(struct tab_lp *lp) 133{ 134 enum isl_lp_result res; 135 unsigned flags = 0; 136 137 lp->is_fixed = 0; 138 139 isl_int_set_si(lp->row->el[0], 0); 140 isl_seq_cpy(lp->row->el + 1, lp->obj, lp->dim); 141 isl_seq_neg(lp->row->el + 1 + lp->dim, lp->obj, lp->dim); 142 if (lp->neq) 143 flags = ISL_TAB_SAVE_DUAL; 144 res = isl_tab_min(lp->tab, lp->row->el, lp->ctx->one, 145 &lp->opt, &lp->opt_denom, flags); 146 isl_int_mul_ui(lp->opt_denom, lp->opt_denom, 2); 147 if (isl_int_abs_lt(lp->opt, lp->opt_denom)) { 148 struct isl_vec *sample = isl_tab_get_sample_value(lp->tab); 149 if (!sample) 150 return -1; 151 isl_seq_inner_product(lp->obj, sample->el + 1, lp->dim, &lp->tmp); 152 isl_seq_inner_product(lp->obj, sample->el + 1 + lp->dim, lp->dim, &lp->tmp2); 153 isl_int_cdiv_q(lp->tmp, lp->tmp, sample->el[0]); 154 isl_int_fdiv_q(lp->tmp2, lp->tmp2, sample->el[0]); 155 if (isl_int_ge(lp->tmp, lp->tmp2)) 156 lp->is_fixed = 1; 157 isl_vec_free(sample); 158 } 159 isl_int_divexact_ui(lp->opt_denom, lp->opt_denom, 2); 160 if (res != isl_lp_ok) 161 return -1; 162 return 0; 163} 164 165/* The current objective function has a fixed (or no) integer value. 166 * Cut the tableau to the hyperplane that fixes this value in 167 * both halves of the tableau. 168 * Return 1 if the resulting tableau is empty. 169 */ 170static int cut_lp_to_hyperplane(struct tab_lp *lp, isl_int *row) 171{ 172 enum isl_lp_result res; 173 174 isl_int_set_si(lp->row->el[0], 0); 175 isl_seq_cpy(lp->row->el + 1, row, lp->dim); 176 isl_seq_clr(lp->row->el + 1 + lp->dim, lp->dim); 177 res = isl_tab_min(lp->tab, lp->row->el, lp->ctx->one, 178 &lp->tmp, NULL, 0); 179 if (res != isl_lp_ok) 180 return -1; 181 182 isl_int_neg(lp->row->el[0], lp->tmp); 183 if (isl_tab_add_eq(lp->tab, lp->row->el) < 0) 184 return -1; 185 186 isl_seq_cpy(lp->row->el + 1 + lp->dim, row, lp->dim); 187 isl_seq_clr(lp->row->el + 1, lp->dim); 188 if (isl_tab_add_eq(lp->tab, lp->row->el) < 0) 189 return -1; 190 191 lp->con_offset += 2; 192 193 return lp->tab->empty; 194} 195 196static void get_obj_val(struct tab_lp* lp, mpq_t *F) 197{ 198 isl_int_neg(mpq_numref(*F), lp->opt); 199 isl_int_set(mpq_denref(*F), lp->opt_denom); 200} 201 202static void delete_lp(struct tab_lp *lp) 203{ 204 if (!lp) 205 return; 206 207 isl_int_clear(lp->opt); 208 isl_int_clear(lp->opt_denom); 209 isl_int_clear(lp->tmp); 210 isl_int_clear(lp->tmp2); 211 isl_vec_free(lp->row); 212 free(lp->stack); 213 isl_tab_free(lp->tab); 214 isl_ctx_deref(lp->ctx); 215 free(lp); 216} 217 218static int add_lp_row(struct tab_lp *lp, isl_int *row, int dim) 219{ 220 lp->stack[lp->neq] = isl_tab_snap(lp->tab); 221 222 isl_int_set_si(lp->row->el[0], 0); 223 isl_seq_cpy(lp->row->el + 1, row, lp->dim); 224 isl_seq_neg(lp->row->el + 1 + lp->dim, row, lp->dim); 225 226 if (isl_tab_add_valid_eq(lp->tab, lp->row->el) < 0) 227 return -1; 228 229 return lp->neq++; 230} 231 232static void get_alpha(struct tab_lp* lp, int row, mpq_t *alpha) 233{ 234 row += lp->con_offset; 235 isl_int_neg(mpq_numref(*alpha), lp->tab->dual->el[1 + row]); 236 isl_int_set(mpq_denref(*alpha), lp->tab->dual->el[0]); 237} 238 239static int del_lp_row(struct tab_lp *lp) 240{ 241 lp->neq--; 242 return isl_tab_rollback(lp->tab, lp->stack[lp->neq]); 243} 244