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