bound.cc

   1 #include <assert.h>
   2 #include <iostream>
   3 #include <barvinok/evalue.h>
   4 #include <barvinok/options.h>
   5 #include <barvinok/util.h>
   6 #include <barvinok/barvinok.h>
   7 #include <barvinok/bernstein.h>
   8 #include <bound_common.h>
   9 #include "bound_options.h"
  10 #include "evalue_convert.h"
  11 #include "evalue_read.h"
  12 #include "verify.h"
  13
  14 using std::cout;
  15 using std::cerr;
  16 using std::endl;
  17 using namespace barvinok;
  18
  19 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  20
  21 struct verify_point_bound {
  22         struct verify_point_data vpd;
  23         isl_pw_qpolynomial *pwqp;
  24         isl_pw_qpolynomial_fold *pwf;
  25 };
  26
  27 static int verify_point(__isl_take isl_point *pnt, void *user)
  28 {
  29         int i;
  30         unsigned nparam;
  31         struct verify_point_bound *vpb = (struct verify_point_bound *) user;
  32         isl_int v, n, d, b, t;
  33         isl_pw_qpolynomial *pwqp;
  34         isl_qpolynomial *bound;
  35         isl_qpolynomial *opt;
  36         const char *minmax;
  37         int sign;
  38         int ok;
  39         int cst;
  40         FILE *out = vpb->vpd.options->print_all ? stdout : stderr;
  41
  42         vpb->vpd.n--;
  43
  44         if (vpb->vpd.options->barvinok->bernstein_optimize == BV_BERNSTEIN_MAX) {
  45                 minmax = "max";
  46                 sign = 1;
  47         } else {
  48                 minmax = "min";
  49                 sign = -1;
  50         }
  51
  52         isl_int_init(b);
  53         isl_int_init(t);
  54         isl_int_init(v);
  55         isl_int_init(n);
  56         isl_int_init(d);
  57
  58         pwqp = isl_pw_qpolynomial_copy(vpb->pwqp);
  59
  60         nparam = isl_pw_qpolynomial_dim(pwqp, isl_dim_param);
  61         for (i = 0; i < nparam; ++i) {
  62                 isl_point_get_coordinate(pnt, isl_dim_param, i, &v);
  63                 pwqp = isl_pw_qpolynomial_fix_dim(pwqp, isl_dim_param, i, v);
  64         }
  65
  66         bound = isl_pw_qpolynomial_fold_eval(isl_pw_qpolynomial_fold_copy(vpb->pwf),
  67                                                 isl_point_copy(pnt));
  68
  69         if (sign > 0)
  70                 opt = isl_pw_qpolynomial_max(pwqp);
  71         else
  72                 opt = isl_pw_qpolynomial_min(pwqp);
  73
  74         cst = isl_qpolynomial_is_cst(opt, &n, &d);
  75         if (cst != 1)
  76                 goto error;
  77         if (sign > 0)
  78                 isl_int_fdiv_q(v, n, d);
  79         else
  80                 isl_int_cdiv_q(v, n, d);
  81
  82         cst = isl_qpolynomial_is_cst(bound, &n, &d);
  83         if (cst != 1)
  84                 goto error;
  85         if (sign > 0)
  86                 isl_int_fdiv_q(b, n, d);
  87         else
  88                 isl_int_cdiv_q(b, n, d);
  89
  90         if (sign > 0)
  91                 ok = value_ge(b, v);
  92         else
  93                 ok = value_le(b, v);
  94
  95         if (vpb->vpd.options->print_all || !ok) {
  96                 fprintf(out, "%s(", minmax);
  97                 for (i = 0; i < nparam; ++i) {
  98                         if (i)
  99                                 fprintf(out, ", ");
 100                         isl_point_get_coordinate(pnt, isl_dim_param, i, &t);
 101                         isl_int_print(out, t, 0);
 102                 }
 103                 fprintf(out, ") = ");
 104                 isl_int_print(out, n, 0);
 105                 if (!isl_int_is_one(d)) {
 106                         fprintf(out, "/");
 107                         isl_int_print(out, d, 0);
 108                         fprintf(out, " (");
 109                         isl_int_print(out, b, 0);
 110                         fprintf(out, ")");
 111                 }
 112                 fprintf(out, ", %s(EP) = ", minmax);
 113                 isl_int_print(out, v, 0);
 114                 if (ok)
 115                         fprintf(out, ". OK\n");
 116                 else
 117                         fprintf(out, ". NOT OK\n");
 118         } else if ((vpb->vpd.n % vpb->vpd.s) == 0) {
 119                 printf("o");
 120                 fflush(stdout);
 121         }
 122
 123         if (0) {
 124 error:
 125                 ok = 0;
 126         }
 127
 128         isl_qpolynomial_free(bound);
 129         isl_qpolynomial_free(opt);
 130         isl_point_free(pnt);
 131
 132         isl_int_clear(t);
 133         isl_int_clear(d);
 134         isl_int_clear(n);
 135         isl_int_clear(v);
 136         isl_int_clear(b);
 137
 138         if (!ok)
 139                 vpb->vpd.error = 1;
 140
 141         if (vpb->vpd.options->continue_on_error)
 142                 ok = 1;
 143
 144         return (vpb->vpd.n >= 1 && ok) ? 0 : -1;
 145 }
 146
 147 static int verify(__isl_keep isl_pw_qpolynomial_fold *pwf, evalue *EP, unsigned nvar,
 148                   struct verify_options *options)
 149 {
 150         struct verify_point_bound vpb = { { options } };
 151         isl_ctx *ctx;
 152         isl_dim *dim;
 153         isl_set *context;
 154         int r;
 155         unsigned nparam;
 156
 157         ctx = isl_pw_qpolynomial_fold_get_ctx(pwf);
 158         nparam = isl_pw_qpolynomial_fold_dim(pwf, isl_dim_param);
 159         vpb.pwf = pwf;
 160         dim = isl_dim_set_alloc(ctx, nvar + nparam, 0);
 161         vpb.pwqp = isl_pw_qpolynomial_from_evalue(dim, EP);
 162         vpb.pwqp = isl_pw_qpolynomial_move(vpb.pwqp, isl_dim_set, 0,
 163                                                 isl_dim_param, 0, nvar);
 164         context = isl_pw_qpolynomial_fold_domain(
 165                                         isl_pw_qpolynomial_fold_copy(vpb.pwf));
 166         context = verify_context_set_bounds(context, options);
 167
 168         r = verify_point_data_init(&vpb.vpd, context);
 169
 170         if (r == 0)
 171                 isl_set_foreach_point(context, verify_point, &vpb);
 172         if (vpb.vpd.error)
 173                 r = -1;
 174
 175         isl_set_free(context);
 176         isl_pw_qpolynomial_free(vpb.pwqp);
 177
 178         verify_point_data_fini(&vpb.vpd);
 179
 180         return r;
 181 }
 182
 183 static __isl_give isl_pw_qpolynomial_fold *iterate(
 184         __isl_take isl_pw_qpolynomial *pwqp, enum isl_fold type)
 185 {
 186         isl_dim *dim = isl_pw_qpolynomial_get_dim(pwqp);
 187         isl_set *set;
 188         isl_qpolynomial *qp;
 189         isl_qpolynomial_fold *fold;
 190         unsigned nvar;
 191
 192         assert(isl_dim_size(dim, isl_dim_param) == 0);
 193         nvar = isl_dim_size(dim, isl_dim_set);
 194
 195         if (type == isl_fold_min)
 196                 qp = isl_pw_qpolynomial_min(pwqp);
 197         else
 198                 qp = isl_pw_qpolynomial_max(pwqp);
 199
 200         qp = isl_qpolynomial_drop_dims(qp, isl_dim_set, 0, nvar);
 201         fold = isl_qpolynomial_fold_alloc(type, qp);
 202         dim = isl_dim_drop(dim, isl_dim_set, 0, nvar);
 203         set = isl_set_universe(dim);
 204
 205         return isl_pw_qpolynomial_fold_alloc(set, fold);
 206 }
 207
 208 /*
 209  * Split (partition) EP into a partition with (sub)domains containing
 210  * size integer points or less and a partition with (sub)domains
 211  * containing more integer points.
 212  */
 213 static void split_on_domain_size(evalue *EP, evalue **EP_less, evalue **EP_more,
 214                                  int size, barvinok_options *options)
 215 {
 216     assert(value_zero_p(EP->d));
 217     assert(EP->x.p->type == partition);
 218     assert(EP->x.p->size >= 2);
 219
 220     struct evalue_section *s_less = new evalue_section[EP->x.p->size/2];
 221     struct evalue_section *s_more = new evalue_section[EP->x.p->size/2];
 222
 223     int n_less = 0;
 224     int n_more = 0;
 225
 226     Value c;
 227     value_init(c);
 228
 229     for (int i = 0; i < EP->x.p->size/2; ++i) {
 230         Polyhedron *D = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
 231         Polyhedron *D_less = NULL;
 232         Polyhedron *D_more = NULL;
 233         Polyhedron **next_less = &D_less;
 234         Polyhedron **next_more = &D_more;
 235
 236         for (Polyhedron *P = D; P; P = P->next) {
 237             Polyhedron *next = P->next;
 238             P->next = NULL;
 239             barvinok_count_with_options(P, &c, options);
 240             P->next = next;
 241
 242             if (value_zero_p(c))
 243                 continue;
 244
 245             if (value_pos_p(c) && value_cmp_si(c, size) <= 0) {
 246                 *next_less = Polyhedron_Copy(P);
 247                 next_less = &(*next_less)->next;
 248             } else {
 249                 *next_more = Polyhedron_Copy(P);
 250                 next_more = &(*next_more)->next;
 251             }
 252         }
 253
 254         if (D_less) {
 255             s_less[n_less].D = D_less;
 256             s_less[n_less].E = evalue_dup(&EP->x.p->arr[2*i+1]);
 257             n_less++;
 258         } else {
 259             s_more[n_more].D = D_more;
 260             s_more[n_more].E = evalue_dup(&EP->x.p->arr[2*i+1]);
 261             n_more++;
 262         }
 263     }
 264
 265     value_clear(c);
 266
 267     *EP_less = evalue_from_section_array(s_less, n_less);
 268     *EP_more = evalue_from_section_array(s_more, n_more);
 269
 270     delete [] s_less;
 271     delete [] s_more;
 272 }
 273
 274 static __isl_give isl_pw_qpolynomial_fold *optimize(evalue *EP, unsigned nvar,
 275         Polyhedron *C, __isl_take isl_dim *dim,
 276         struct options *options)
 277 {
 278     isl_pw_qpolynomial_fold *pwf;
 279     if (options->iterate > 0) {
 280         evalue *EP_less = NULL;
 281         evalue *EP_more = NULL;
 282         isl_pw_qpolynomial_fold *pwf = NULL, *pwf_more = NULL;
 283
 284         split_on_domain_size(EP, &EP_less, &EP_more, options->iterate,
 285                                 options->verify->barvinok);
 286         if (!EVALUE_IS_ZERO(*EP_less)) {
 287             options->iterate = -1;
 288             pwf = optimize(EP_less, nvar, C, isl_dim_copy(dim), options);
 289         }
 290         if (!EVALUE_IS_ZERO(*EP_more)) {
 291             options->iterate = 0;
 292             pwf_more = optimize(EP_more, nvar, C, isl_dim_copy(dim), options);
 293         }
 294         isl_dim_free(dim);
 295         evalue_free(EP_less);
 296         evalue_free(EP_more);
 297         if (!pwf)
 298             return pwf_more;
 299         if (!pwf_more)
 300             return pwf;
 301         return isl_pw_qpolynomial_fold_add(pwf, pwf_more);
 302     }
 303     int method = options->verify->barvinok->bound;
 304     isl_dim *dim_EP;
 305     isl_pw_qpolynomial *pwqp;
 306     enum isl_fold type = options->lower ? isl_fold_min : isl_fold_max;
 307     dim_EP = isl_dim_insert(dim, isl_dim_param, 0, nvar);
 308     pwqp = isl_pw_qpolynomial_from_evalue(dim_EP, EP);
 309     pwqp = isl_pw_qpolynomial_move(pwqp, isl_dim_set, 0, isl_dim_param, 0, nvar);
 310     if (options->iterate)
 311         pwf = iterate(pwqp, type);
 312     else
 313         pwf = isl_pw_qpolynomial_bound(pwqp, type, method);
 314     return pwf;
 315 }
 316
 317 static int optimize(evalue *EP, const char **all_vars,
 318                     unsigned nvar, unsigned nparam, struct options *options)
 319 {
 320     Polyhedron *U;
 321     U = Universe_Polyhedron(nparam);
 322     int print_solution = 1;
 323     int result = 0;
 324     isl_ctx *ctx = isl_ctx_alloc_with_options(options_arg, options);
 325     isl_dim *dim;
 326     isl_pw_qpolynomial_fold *pwf;
 327
 328     dim = isl_dim_set_alloc(ctx, nparam, 0);
 329     for (int i = 0; i < nparam; ++i)
 330         dim = isl_dim_set_name(dim, isl_dim_param, i, all_vars[nvar + i]);
 331
 332     if (options->verify->verify) {
 333         verify_options_set_range(options->verify, nvar+nparam);
 334         if (!options->verify->barvinok->verbose)
 335             print_solution = 0;
 336     }
 337
 338     if (options->lower)
 339         options->verify->barvinok->bernstein_optimize = BV_BERNSTEIN_MIN;
 340     else
 341         options->verify->barvinok->bernstein_optimize = BV_BERNSTEIN_MAX;
 342     pwf = optimize(EP, nvar, U, dim, options);
 343     assert(pwf);
 344     if (print_solution) {
 345         isl_printer *p = isl_printer_to_file(ctx, stdout);
 346         p = isl_printer_print_pw_qpolynomial_fold(p, pwf);
 347         p = isl_printer_end_line(p);
 348         isl_printer_free(p);
 349     }
 350     if (options->verify->verify) {
 351         result = verify(pwf, EP, nvar, options->verify);
 352     }
 353     isl_pw_qpolynomial_fold_free(pwf);
 354
 355     Polyhedron_Free(U);
 356
 357     isl_ctx_free(ctx);
 358
 359     return result;
 360 }
 361
 362 int main(int argc, char **argv)
 363 {
 364     evalue *EP;
 365     const char **all_vars = NULL;
 366     unsigned nvar;
 367     unsigned nparam;
 368     struct options *options = options_new_with_defaults();
 369     int result = 0;
 370
 371     argc = options_parse(options, argc, argv, ISL_ARG_ALL);
 372
 373     EP = evalue_read_from_file(stdin, options->var_list, &all_vars,
 374                                &nvar, &nparam, options->verify->barvinok->MaxRays);
 375     assert(EP);
 376
 377     if (options->split)
 378         evalue_split_periods(EP, options->split, options->verify->barvinok->MaxRays);
 379
 380     evalue_convert(EP, options->convert, options->verify->barvinok->verbose,
 381                         nvar+nparam, all_vars);
 382
 383     if (EVALUE_IS_ZERO(*EP))
 384         print_evalue(stdout, EP, all_vars);
 385     else
 386         result = optimize(EP, all_vars, nvar, nparam, options);
 387
 388     evalue_free(EP);
 389
 390     Free_ParamNames(all_vars, nvar+nparam);
 391     return result;
 392 }