basis_reduction_tab.c

   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
  15 struct 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
  33 static struct tab_lp *init_lp(struct isl_tab *tab);
  34 static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim);
  35 static int solve_lp(struct tab_lp *lp);
  36 static void get_obj_val(struct tab_lp* lp, mpq_t *F);
  37 static void delete_lp(struct tab_lp *lp);
  38 static int add_lp_row(struct tab_lp *lp, isl_int *row, int dim);
  39 static void get_alpha(struct tab_lp* lp, int row, mpq_t *alpha);
  40 static int del_lp_row(struct tab_lp *lp) WARN_UNUSED;
  41 static 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  */
  71 static 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
  88 static 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;
 122 error:
 123         delete_lp(lp);
 124         return NULL;
 125 }
 126
 127 static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim)
 128 {
 129         lp->obj = row;
 130 }
 131
 132 static 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  */
 170 static 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
 196 static 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
 202 static 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
 218 static 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
 232 static 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
 239 static int del_lp_row(struct tab_lp *lp)
 240 {
 241         lp->neq--;
 242         return isl_tab_rollback(lp->tab, lp->stack[lp->neq]);
 243 }