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