bound.c

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