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