bound.c

   1 #include <assert.h>
   2 #include <isl/ctx.h>
   3 #include <isl/val.h>
   4 #include <isl/space.h>
   5 #include <isl/point.h>
   6 #include <isl/set.h>
   7 #include <isl/polynomial.h>
   8 #include <isl/stream.h>
   9 #include <isl/printer.h>
  10 #include <barvinok/options.h>
  11 #include <barvinok/util.h>
  12 #include <barvinok/isl.h>
  13 #include <barvinok/barvinok.h>
  14 #include "verify.h"
  15
  16 struct options {
  17         struct verify_options    *verify;
  18         long split;
  19         int lower;
  20         long iterate;
  21 };
  22
  23 ISL_ARGS_START(struct options, options_args)
  24 ISL_ARG_CHILD(struct options, verify, NULL,
  25         &verify_options_args, "verification")
  26 ISL_ARG_LONG(struct options, split, 0, "split", 0, NULL)
  27 ISL_ARG_OPT_LONG(struct options, iterate, 0, "iterate", 0, -1,
  28         "exact result by iterating over domain (of specified maximal size)")
  29 ISL_ARG_BOOL(struct options, lower, 0, "lower", 0,
  30         "compute lower bound instead of upper bound")
  31 ISL_ARGS_END
  32
  33 ISL_ARG_DEF(options, struct options, options_args)
  34
  35 struct verify_point_bound {
  36         struct verify_point_data vpd;
  37         isl_pw_qpolynomial *pwqp;
  38         isl_pw_qpolynomial_fold *pwf;
  39         enum isl_fold type;
  40 };
  41
  42 static isl_stat verify_point(__isl_take isl_point *pnt, void *user)
  43 {
  44         int i;
  45         unsigned nparam;
  46         struct verify_point_bound *vpb = (struct verify_point_bound *) user;
  47         isl_val *v, *b, *q, *t;
  48         isl_pw_qpolynomial *pwqp;
  49         const char *minmax;
  50         int sign;
  51         isl_bool ok;
  52         FILE *out = vpb->vpd.options->print_all ? stdout : stderr;
  53
  54         vpb->vpd.n--;
  55
  56         if (vpb->type == isl_fold_max) {
  57                 minmax = "max";
  58                 sign = 1;
  59         } else {
  60                 minmax = "min";
  61                 sign = -1;
  62         }
  63
  64         pwqp = isl_pw_qpolynomial_copy(vpb->pwqp);
  65
  66         nparam = isl_pw_qpolynomial_dim(pwqp, isl_dim_param);
  67         for (i = 0; i < nparam; ++i) {
  68                 t = isl_point_get_coordinate_val(pnt, isl_dim_param, i);
  69                 pwqp = isl_pw_qpolynomial_fix_val(pwqp, isl_dim_param, i, t);
  70         }
  71
  72         b = isl_pw_qpolynomial_fold_eval(isl_pw_qpolynomial_fold_copy(vpb->pwf),
  73                                                 isl_point_copy(pnt));
  74
  75         if (sign > 0)
  76                 v = isl_pw_qpolynomial_max(pwqp);
  77         else
  78                 v = isl_pw_qpolynomial_min(pwqp);
  79
  80         if (sign > 0)
  81                 v = isl_val_floor(v);
  82         else
  83                 v = isl_val_ceil(v);
  84
  85         q = isl_val_copy(b);
  86         if (sign > 0)
  87                 b = isl_val_floor(b);
  88         else
  89                 b = isl_val_ceil(b);
  90
  91         if (sign > 0)
  92                 ok = isl_val_ge(b, v);
  93         else
  94                 ok = isl_val_le(b, v);
  95
  96         if (vpb->vpd.options->print_all || !ok) {
  97                 isl_ctx *ctx = isl_point_get_ctx(pnt);
  98                 isl_printer *p;
  99                 p = isl_printer_to_file(ctx, out);
 100                 p = isl_printer_print_str(p, minmax);
 101                 p = isl_printer_print_str(p, "(");
 102                 for (i = 0; i < nparam; ++i) {
 103                         if (i)
 104                                 p = isl_printer_print_str(p, ", ");
 105                         t = isl_point_get_coordinate_val(pnt, isl_dim_param, i);
 106                         p = isl_printer_print_val(p, t);
 107                         isl_val_free(t);
 108                 }
 109                 p = isl_printer_print_str(p, ") = ");
 110                 p = isl_printer_print_val(p, q);
 111                 if (!isl_val_is_int(q)) {
 112                         p = isl_printer_print_str(p, " (");
 113                         p = isl_printer_print_val(p, b);
 114                         p = isl_printer_print_str(p, ")");
 115                 }
 116                 p = isl_printer_print_str(p, ", ");
 117                 p = isl_printer_print_str(p, minmax);
 118                 p = isl_printer_print_str(p, "(EP) = ");
 119                 p = isl_printer_print_val(p, v);
 120                 if (ok)
 121                         p = isl_printer_print_str(p, ". OK");
 122                 else
 123                         p = isl_printer_print_str(p, ". NOT OK");
 124                 p = isl_printer_end_line(p);
 125                 isl_printer_free(p);
 126         } else if ((vpb->vpd.n % vpb->vpd.s) == 0) {
 127                 printf("o");
 128                 fflush(stdout);
 129         }
 130
 131         isl_point_free(pnt);
 132
 133         isl_val_free(b);
 134         isl_val_free(q);
 135         isl_val_free(v);
 136
 137         if (ok < 0 || !ok)
 138                 vpb->vpd.error = 1;
 139
 140         if (vpb->vpd.options->continue_on_error)
 141                 ok = isl_bool_true;
 142
 143         if (vpb->vpd.n < 1)
 144                 return isl_stat_error;
 145         if (ok != isl_bool_true)
 146                 return isl_stat_error;
 147         return isl_stat_ok;
 148 }
 149
 150 static int verify(__isl_keep isl_pw_qpolynomial_fold *pwf,
 151         __isl_take isl_pw_qpolynomial *pwqp, enum isl_fold type,
 152         struct verify_options *options)
 153 {
 154         struct verify_point_bound vpb = { { options } };
 155         isl_set *context;
 156         int r;
 157
 158         vpb.pwf = pwf;
 159         vpb.type = type;
 160         vpb.pwqp = pwqp;
 161         context = isl_pw_qpolynomial_fold_domain(
 162                                         isl_pw_qpolynomial_fold_copy(vpb.pwf));
 163         context = verify_context_set_bounds(context, options);
 164
 165         r = verify_point_data_init(&vpb.vpd, context);
 166
 167         if (r == 0)
 168                 isl_set_foreach_point(context, verify_point, &vpb);
 169         if (vpb.vpd.error)
 170                 r = -1;
 171
 172         isl_set_free(context);
 173         isl_pw_qpolynomial_free(pwqp);
 174
 175         verify_point_data_fini(&vpb.vpd);
 176
 177         return r;
 178 }
 179
 180 static __isl_give isl_pw_qpolynomial_fold *iterate(
 181         __isl_take isl_pw_qpolynomial *pwqp, enum isl_fold type)
 182 {
 183         isl_space *dim = isl_pw_qpolynomial_get_space(pwqp);
 184         isl_set *set;
 185         isl_val *v;
 186         isl_qpolynomial *qp;
 187         isl_qpolynomial_fold *fold;
 188         unsigned nvar;
 189
 190         assert(isl_space_dim(dim, isl_dim_param) == 0);
 191         nvar = isl_space_dim(dim, isl_dim_set);
 192
 193         if (type == isl_fold_min)
 194                 v = isl_pw_qpolynomial_min(pwqp);
 195         else
 196                 v = isl_pw_qpolynomial_max(pwqp);
 197
 198         dim = isl_space_drop_dims(dim, isl_dim_set, 0, nvar);
 199         qp = isl_qpolynomial_val_on_domain(isl_space_copy(dim), v);
 200         fold = isl_qpolynomial_fold_alloc(type, qp);
 201         set = isl_set_universe(dim);
 202
 203         return isl_pw_qpolynomial_fold_alloc(type, set, fold);
 204 }
 205
 206 struct bv_split_data {
 207         int size;
 208         __isl_give isl_pw_qpolynomial *pwqp_less;
 209         __isl_give isl_pw_qpolynomial *pwqp_more;
 210 };
 211
 212 static isl_stat split_on_size(__isl_take isl_set *set,
 213         __isl_take isl_qpolynomial *qp, void *user)
 214 {
 215         struct bv_split_data *data = (struct bv_split_data *)user;
 216         int bounded;
 217         isl_set *set_np;
 218         isl_pw_qpolynomial *pwqp;
 219         int nparam;
 220
 221         nparam = isl_set_dim(set, isl_dim_param);
 222         set_np = isl_set_move_dims(isl_set_copy(set), isl_dim_set, 0,
 223                                         isl_dim_param, 0, nparam);
 224         bounded = isl_set_is_bounded(set_np);
 225         assert(bounded >= 0);
 226         if (bounded) {
 227                 isl_pw_qpolynomial *pwqp;
 228                 isl_val *m;
 229
 230                 pwqp = isl_set_card(set_np);
 231                 m = isl_pw_qpolynomial_max(pwqp);
 232                 bounded = isl_val_cmp_si(m, data->size) <= 0;
 233                 isl_val_free(m);
 234         } else
 235                 isl_set_free(set_np);
 236
 237         pwqp = isl_pw_qpolynomial_alloc(set, qp);
 238         if (bounded)
 239                 data->pwqp_less = isl_pw_qpolynomial_add_disjoint(
 240                                                 data->pwqp_less, pwqp);
 241         else
 242                 data->pwqp_more = isl_pw_qpolynomial_add_disjoint(
 243                                                 data->pwqp_more, pwqp);
 244
 245         return isl_stat_ok;
 246 }
 247
 248 /*
 249  * Split (partition) pwpq into a partition with (sub)domains containing
 250  * size integer points or less and a partition with (sub)domains
 251  * containing more integer points.
 252  */
 253 static int split_on_domain_size(__isl_take isl_pw_qpolynomial *pwqp,
 254         __isl_give isl_pw_qpolynomial **pwqp_less,
 255         __isl_give isl_pw_qpolynomial **pwqp_more,
 256         int size)
 257 {
 258         isl_space *dim;
 259         struct bv_split_data data = { size };
 260         int r;
 261
 262         dim = isl_pw_qpolynomial_get_space(pwqp);
 263         data.pwqp_less = isl_pw_qpolynomial_zero(isl_space_copy(dim));
 264         data.pwqp_more = isl_pw_qpolynomial_zero(dim);
 265         r = isl_pw_qpolynomial_foreach_piece(pwqp, &split_on_size, &data);
 266         *pwqp_less = data.pwqp_less;
 267         *pwqp_more = data.pwqp_more;
 268         isl_pw_qpolynomial_free(pwqp);
 269
 270         return r;
 271 }
 272
 273 static __isl_give isl_pw_qpolynomial_fold *optimize(
 274         __isl_take isl_pw_qpolynomial *pwqp, enum isl_fold type,
 275         struct options *options)
 276 {
 277     isl_pw_qpolynomial_fold *pwf;
 278
 279     if (options->iterate > 0) {
 280         isl_pw_qpolynomial *pwqp_less, *pwqp_more;
 281         isl_pw_qpolynomial_fold *pwf_less, *pwf_more;
 282         split_on_domain_size(pwqp, &pwqp_less, &pwqp_more, options->iterate);
 283         pwf_less = iterate(pwqp_less, type);
 284         pwf_more = isl_pw_qpolynomial_bound(pwqp_more, type, NULL);
 285         pwf = isl_pw_qpolynomial_fold_fold(pwf_less, pwf_more);
 286     } else if (options->iterate)
 287         pwf = iterate(pwqp, type);
 288     else
 289         pwf = isl_pw_qpolynomial_bound(pwqp, type, NULL);
 290     return pwf;
 291 }
 292
 293 static int optimize_and_print(__isl_take isl_pw_qpolynomial *pwqp,
 294         struct options *options)
 295 {
 296     int print_solution = 1;
 297     int result = 0;
 298     isl_pw_qpolynomial_fold *pwf;
 299     enum isl_fold type = options->lower ? isl_fold_min : isl_fold_max;
 300
 301     if (options->verify->verify) {
 302         isl_space *dim = isl_pw_qpolynomial_get_space(pwqp);
 303         unsigned total = isl_space_dim(dim, isl_dim_all);
 304         isl_space_free(dim);
 305         verify_options_set_range(options->verify, total);
 306         if (!options->verify->barvinok->verbose)
 307             print_solution = 0;
 308     }
 309
 310     pwf = optimize(isl_pw_qpolynomial_copy(pwqp), type, options);
 311     assert(pwf);
 312     if (print_solution) {
 313         isl_ctx *ctx = isl_pw_qpolynomial_get_ctx(pwqp);
 314         isl_printer *p = isl_printer_to_file(ctx, stdout);
 315         p = isl_printer_print_pw_qpolynomial_fold(p, pwf);
 316         p = isl_printer_end_line(p);
 317         isl_printer_free(p);
 318     }
 319     if (options->verify->verify) {
 320         enum isl_fold type = options->lower ? isl_fold_min : isl_fold_max;
 321         result = verify(pwf, pwqp, type, options->verify);
 322     } else
 323         isl_pw_qpolynomial_free(pwqp);
 324     isl_pw_qpolynomial_fold_free(pwf);
 325
 326     return result;
 327 }
 328
 329 int main(int argc, char **argv)
 330 {
 331     struct options *options = options_new_with_defaults();
 332     isl_ctx *ctx;
 333     isl_pw_qpolynomial *pwqp;
 334     struct isl_stream *s;
 335     int result = 0;
 336
 337     argc = options_parse(options, argc, argv, ISL_ARG_ALL);
 338     ctx = isl_ctx_alloc_with_options(&options_args, options);
 339
 340     s = isl_stream_new_file(ctx, stdin);
 341     pwqp = isl_stream_read_pw_qpolynomial(s);
 342
 343     if (options->split)
 344         pwqp = isl_pw_qpolynomial_split_periods(pwqp, options->split);
 345
 346     result = optimize_and_print(pwqp, options);
 347
 348     isl_stream_free(s);
 349     isl_ctx_free(ctx);
 350     return result;
 351 }