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 #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 }