basis_reduction_tab.c

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