polysign_isl.c

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