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         int              neq;
  14         unsigned         dim;
  15         /* number of constraints in initial product tableau */
  16         int              con_offset;
  17 };
  18
  19 static struct tab_lp *init_lp(struct isl_basic_set *bset);
  20 static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim);
  21 static int solve_lp(struct tab_lp *lp);
  22 static void get_obj_val(struct tab_lp* lp, mpq_t *F);
  23 static void delete_lp(struct tab_lp *lp);
  24 static int add_lp_row(struct tab_lp *lp, isl_int *row, int dim);
  25 static void get_alpha(struct tab_lp* lp, int row, mpq_t *alpha);
  26 static void del_lp_row(struct tab_lp *lp);
  27
  28 #define GBR_LP                              struct tab_lp
  29 #define GBR_type                            mpq_t
  30 #define GBR_init(v)                         mpq_init(v)
  31 #define GBR_clear(v)                        mpq_clear(v)
  32 #define GBR_set(a,b)                        mpq_set(a,b)
  33 #define GBR_set_ui(a,b)                     mpq_set_ui(a,b,1)
  34 #define GBR_mul(a,b,c)                      mpq_mul(a,b,c)
  35 #define GBR_lt(a,b)                         (mpq_cmp(a,b) < 0)
  36 #define GBR_floor(a,b)                      mpz_fdiv_q(a,mpq_numref(b),mpq_denref(b))
  37 #define GBR_ceil(a,b)                       mpz_cdiv_q(a,mpq_numref(b),mpq_denref(b))
  38 #define GBR_lp_init(P)                      init_lp(P)
  39 #define GBR_lp_set_obj(lp, obj, dim)        set_lp_obj(lp, obj, dim)
  40 #define GBR_lp_solve(lp)                    solve_lp(lp)
  41 #define GBR_lp_get_obj_val(lp, F)           get_obj_val(lp, F)
  42 #define GBR_lp_delete(lp)                   delete_lp(lp)
  43 #define GBR_lp_next_row(lp)                 lp->neq
  44 #define GBR_lp_add_row(lp, row, dim)        add_lp_row(lp, row, dim)
  45 #define GBR_lp_get_alpha(lp, row, alpha)    get_alpha(lp, row, alpha)
  46 #define GBR_lp_del_row(lp)                  del_lp_row(lp);
  47 #include "basis_reduction_templ.c"
  48
  49 /* Set up a tableau for the Cartesian product of bset with itself.
  50  * This could be optimized by first setting up a tableau for bset
  51  * and then performing the Cartesian product on the tableau.
  52  */
  53 static struct isl_tab *gbr_tab(struct isl_basic_set *bset,
  54         struct isl_vec *row)
  55 {
  56         int i, j;
  57         unsigned dim;
  58         struct isl_tab *tab;
  59
  60         if (!bset || !row)
  61                 return NULL;
  62
  63         dim = isl_basic_set_total_dim(bset);
  64         tab = isl_tab_alloc(bset->ctx, 2 * bset->n_ineq + dim + 1, 2 * dim, 0);
  65
  66         for (i = 0; i < 2; ++i) {
  67                 isl_seq_clr(row->el + 1 + (1 - i) * dim, dim);
  68                 for (j = 0; j < bset->n_ineq; ++j) {
  69                         isl_int_set(row->el[0], bset->ineq[j][0]);
  70                         isl_seq_cpy(row->el + 1 + i * dim,
  71                                     bset->ineq[j] + 1, dim);
  72                         tab = isl_tab_add_ineq(tab, row->el);
  73                         if (!tab || tab->empty)
  74                                 return tab;
  75                 }
  76         }
  77
  78         return tab;
  79 }
  80
  81 static struct tab_lp *init_lp(struct isl_basic_set *bset)
  82 {
  83         struct tab_lp *lp = NULL;
  84
  85         if (!bset)
  86                 return NULL;
  87
  88         isl_assert(bset->ctx, bset->n_eq == 0, return NULL);
  89
  90         lp = isl_calloc_type(bset->ctx, struct tab_lp);
  91         if (!lp)
  92                 return NULL;
  93
  94         isl_int_init(lp->opt);
  95         isl_int_init(lp->opt_denom);
  96
  97         lp->dim = isl_basic_set_total_dim(bset);
  98
  99         lp->ctx = bset->ctx;
 100         isl_ctx_ref(lp->ctx);
 101
 102         lp->stack = isl_alloc_array(lp->ctx, struct isl_tab_undo *, lp->dim);
 103
 104         lp->row = isl_vec_alloc(lp->ctx, 1 + 2 * lp->dim);
 105         if (!lp->row)
 106                 goto error;
 107         lp->tab = gbr_tab(bset, lp->row);
 108         if (!lp->tab)
 109                 goto error;
 110         lp->con_offset = lp->tab->n_con;
 111         lp->obj = NULL;
 112         lp->neq = 0;
 113
 114         return lp;
 115 error:
 116         delete_lp(lp);
 117         return NULL;
 118 }
 119
 120 static void set_lp_obj(struct tab_lp *lp, isl_int *row, int dim)
 121 {
 122         lp->obj = row;
 123 }
 124
 125 static int solve_lp(struct tab_lp *lp)
 126 {
 127         enum isl_lp_result res;
 128         unsigned flags = 0;
 129
 130         isl_int_set_si(lp->row->el[0], 0);
 131         isl_seq_cpy(lp->row->el + 1, lp->obj, lp->dim);
 132         isl_seq_neg(lp->row->el + 1 + lp->dim, lp->obj, lp->dim);
 133         if (lp->neq)
 134                 flags = ISL_TAB_SAVE_DUAL;
 135         res = isl_tab_min(lp->tab, lp->row->el, lp->ctx->one,
 136                           &lp->opt, &lp->opt_denom, flags);
 137         if (res != isl_lp_ok)
 138                 return -1;
 139         return 0;
 140 }
 141
 142 static void get_obj_val(struct tab_lp* lp, mpq_t *F)
 143 {
 144         isl_int_neg(mpq_numref(*F), lp->opt);
 145         isl_int_set(mpq_denref(*F), lp->opt_denom);
 146 }
 147
 148 static void delete_lp(struct tab_lp *lp)
 149 {
 150         if (!lp)
 151                 return;
 152
 153         isl_int_clear(lp->opt);
 154         isl_int_clear(lp->opt_denom);
 155         isl_vec_free(lp->row);
 156         free(lp->stack);
 157         isl_tab_free(lp->tab);
 158         isl_ctx_deref(lp->ctx);
 159         free(lp);
 160 }
 161
 162 static int add_lp_row(struct tab_lp *lp, isl_int *row, int dim)
 163 {
 164         lp->stack[lp->neq] = isl_tab_snap(lp->tab);
 165
 166         isl_int_set_si(lp->row->el[0], 0);
 167         isl_seq_cpy(lp->row->el + 1, row, lp->dim);
 168         isl_seq_neg(lp->row->el + 1 + lp->dim, row, lp->dim);
 169
 170         lp->tab = isl_tab_add_valid_eq(lp->tab, lp->row->el);
 171
 172         return lp->neq++;
 173 }
 174
 175 static void get_alpha(struct tab_lp* lp, int row, mpq_t *alpha)
 176 {
 177         row += lp->con_offset;
 178         isl_int_neg(mpq_numref(*alpha), lp->tab->dual->el[1 + row]);
 179         isl_int_set(mpq_denref(*alpha), lp->tab->dual->el[0]);
 180 }
 181
 182 static void del_lp_row(struct tab_lp *lp)
 183 {
 184         lp->neq--;
 185         isl_tab_rollback(lp->tab, lp->stack[lp->neq]);
 186 }