polysign_isl.c

   1 #include <assert.h>
   2 #include <isl/val_gmp.h>
   3 #include <isl/lp.h>
   4 #include <isl_set_polylib.h>
   5 #include "polysign.h"
   6
   7 static void values2isl(Value *src, isl_int *dst, int len)
   8 {
   9         int i;
  10
  11         for (i = 0; i < len; ++i)
  12                 isl_int_set_gmp(dst[i], src[i]);
  13 }
  14
  15 static __isl_give isl_mat *extract_equalities(isl_ctx *ctx, Matrix *M)
  16 {
  17         int i, j;
  18         int n;
  19         isl_val *v;
  20         isl_mat *eq;
  21
  22         n = 0;
  23         for (i = 0; i < M->NbRows; ++i)
  24                 if (value_zero_p(M->p[i][0]))
  25                         ++n;
  26
  27         eq = isl_mat_alloc(ctx, n, M->NbColumns - 1);
  28         for (i = 0; i < M->NbRows; ++i) {
  29                 if (!value_zero_p(M->p[i][0]))
  30                         continue;
  31                 for (j = 0; j < M->NbColumns - 1; ++j) {
  32                         v = isl_val_int_from_gmp(ctx, M->p[i][1 + j]);
  33                         eq = isl_mat_set_element_val(eq, i, j, v);
  34                 }
  35         }
  36
  37         return eq;
  38 }
  39
  40 static __isl_give isl_mat *extract_inequalities(isl_ctx *ctx, Matrix *M)
  41 {
  42         int i, j;
  43         int n;
  44         isl_val *v;
  45         isl_mat *ineq;
  46
  47         n = 0;
  48         for (i = 0; i < M->NbRows; ++i)
  49                 if (!value_zero_p(M->p[i][0]))
  50                         ++n;
  51
  52         ineq = isl_mat_alloc(ctx, n, M->NbColumns - 1);
  53         for (i = 0; i < M->NbRows; ++i) {
  54                 if (value_zero_p(M->p[i][0]))
  55                         continue;
  56                 for (j = 0; j < M->NbColumns - 1; ++j) {
  57                         v = isl_val_int_from_gmp(ctx, M->p[i][1 + j]);
  58                         ineq = isl_mat_set_element_val(ineq, i, j, v);
  59                 }
  60         }
  61
  62         return ineq;
  63 }
  64
  65 enum order_sign isl_polyhedron_affine_sign(Polyhedron *D, Matrix *T,
  66                                             struct barvinok_options *options)
  67 {
  68         int i;
  69         isl_ctx *ctx = isl_ctx_alloc();
  70         isl_space *dim;
  71         isl_local_space *ls;
  72         isl_aff *aff;
  73         isl_basic_set *bset;
  74         isl_val *min, *max = NULL;
  75         isl_val *v;
  76         enum order_sign sign = order_undefined;
  77
  78         assert(D->Dimension == T->NbColumns - 1);
  79
  80         dim = isl_space_set_alloc(ctx, 0, D->Dimension);
  81         ls = isl_local_space_from_space(isl_space_copy(dim));
  82         bset = isl_basic_set_new_from_polylib(D, dim);
  83         aff = isl_aff_zero_on_domain(ls);
  84         for (i = 0; i < D->Dimension; ++i) {
  85                 v = isl_val_int_from_gmp(ctx, T->p[0][i]);
  86                 aff = isl_aff_set_coefficient_val(aff, isl_dim_in, i, v);
  87         }
  88         v = isl_val_int_from_gmp(ctx, T->p[0][D->Dimension]);
  89         aff = isl_aff_set_constant_val(aff, v);
  90         v = isl_val_int_from_gmp(ctx, T->p[1][D->Dimension]);
  91         aff = isl_aff_scale_down_val(aff, v);
  92
  93         min = isl_basic_set_min_lp_val(bset, aff);
  94         min = isl_val_ceil(min);
  95         assert(min);
  96
  97         if (isl_val_is_nan(min))
  98                 sign = order_undefined;
  99         else if (isl_val_is_pos(min))
 100                 sign = order_gt;
 101         else {
 102                 max = isl_basic_set_max_lp_val(bset, aff);
 103                 max = isl_val_floor(max);
 104                 assert(max);
 105
 106                 if (isl_val_is_neg(max))
 107                         sign = order_lt;
 108                 else if (isl_val_is_zero(min) && isl_val_is_zero(max))
 109                         sign = order_eq;
 110                 else if (isl_val_is_zero(min))
 111                         sign = order_ge;
 112                 else if (isl_val_is_zero(max))
 113                         sign = order_le;
 114                 else
 115                         sign = order_unknown;
 116         }
 117
 118         isl_basic_set_free(bset);
 119         isl_aff_free(aff);
 120         isl_val_free(min);
 121         isl_val_free(max);
 122         isl_ctx_free(ctx);
 123
 124         return sign;
 125 }
 126
 127 static enum lp_result isl_lp_result2lp_result(enum isl_lp_result res)
 128 {
 129         switch (res) {
 130         case isl_lp_ok:
 131                 return lp_ok;
 132         case isl_lp_unbounded:
 133                 return lp_unbounded;
 134         case isl_lp_empty:
 135                 return lp_empty;
 136         default:
 137                 assert(0);
 138         }
 139 }
 140
 141 enum lp_result isl_constraints_opt(Matrix *C, Value *obj, Value denom,
 142                                     enum lp_dir dir, Value *opt)
 143 {
 144         isl_ctx *ctx = isl_ctx_alloc();
 145         isl_space *dim;
 146         isl_mat *eq, *ineq;
 147         isl_basic_set *bset;
 148         isl_int isl_denom, isl_opt;
 149         isl_vec *aff;
 150         enum isl_lp_result res;
 151         int max = dir == lp_max;
 152
 153         isl_int_init(isl_denom);
 154         isl_int_init(isl_opt);
 155
 156         eq = extract_equalities(ctx, C);
 157         ineq = extract_inequalities(ctx, C);
 158         dim = isl_space_set_alloc(ctx, 0, C->NbColumns - 2);
 159         bset = isl_basic_set_from_constraint_matrices(dim, eq, ineq,
 160                         isl_dim_set, isl_dim_div, isl_dim_param, isl_dim_cst);
 161         aff = isl_vec_alloc(ctx, C->NbColumns - 1);
 162         assert(aff);
 163         values2isl(obj, aff->el + 1, aff->size - 1);
 164         values2isl(obj + aff->size - 1, aff->el, 1);
 165         isl_int_set_gmp(isl_denom, denom);
 166
 167         res = isl_basic_set_solve_lp(bset, max, aff->el, isl_denom, &isl_opt,
 168                                         NULL, NULL);
 169         isl_int_get_gmp(isl_opt, *opt);
 170
 171         isl_vec_free(aff);
 172         isl_int_clear(isl_denom);
 173         isl_int_clear(isl_opt);
 174         isl_basic_set_free(bset);
 175         isl_ctx_free(ctx);
 176
 177         return isl_lp_result2lp_result(res);
 178 }