test_bound.c

   1 #include <assert.h>
   2 #include <limits.h>
   3 #include <math.h>
   4 #include <isl/stream.h>
   5 #include <barvinok/options.h>
   6 #include <barvinok/util.h>
   7 #include "verify.h"
   8 #include "config.h"
   9
  10 #ifdef HAVE_SYS_TIMES_H
  11
  12 #include <sys/times.h>
  13
  14 typedef clock_t         my_clock_t;
  15
  16 static my_clock_t time_diff(struct tms *before, struct tms *after)
  17 {
  18         return after->tms_utime - before->tms_utime;
  19 }
  20
  21 #else
  22
  23 typedef int             my_clock_t;
  24
  25 struct tms {};
  26 static void times(struct tms* time)
  27 {
  28 }
  29 static my_clock_t time_diff(struct tms *before, struct tms *after)
  30 {
  31         return 0;
  32 }
  33
  34 #endif
  35
  36 static struct {
  37     int     method;
  38 } methods[] = {
  39     { ISL_BOUND_BERNSTEIN },
  40     { ISL_BOUND_RANGE },
  41 };
  42
  43 #define nr_methods (sizeof(methods)/sizeof(*methods))
  44
  45 struct options {
  46     struct verify_options    *verify;
  47     int quiet;
  48 };
  49
  50 struct isl_arg options_arg[] = {
  51 ISL_ARG_CHILD(struct options, verify, NULL, verify_options_arg, NULL)
  52 ISL_ARG_BOOL(struct options, quiet, 'q', "quiet", 0, NULL)
  53 ISL_ARG_END
  54 };
  55
  56 ISL_ARG_DEF(options, struct options, options_arg)
  57
  58 struct result_data {
  59     isl_int                 n;
  60     double                  RE_sum[nr_methods];
  61
  62     my_clock_t              ticks[nr_methods];
  63     size_t                  size[nr_methods];
  64 };
  65
  66 void result_data_init(struct result_data *result)
  67 {
  68     int i;
  69     for (i = 0; i < nr_methods; ++i) {
  70         result->RE_sum[i] = 0;
  71         result->ticks[i] = 0;
  72         result->size[i] = 0;
  73     }
  74     isl_int_init(result->n);
  75 }
  76
  77 void result_data_clear(struct result_data *result)
  78 {
  79     int i;
  80     isl_int_clear(result->n);
  81 }
  82
  83 void result_data_print(struct result_data *result, int n)
  84 {
  85     int i;
  86
  87     fprintf(stderr, "%g", (double)result->ticks[0]/n);
  88     for (i = 1; i < nr_methods; ++i)
  89         fprintf(stderr, ", %g", (double)result->ticks[i]/n);
  90     fprintf(stderr, "\n");
  91
  92     fprintf(stderr, "%zd/%d", result->size[0], n);
  93     for (i = 1; i < nr_methods; ++i)
  94         fprintf(stderr, ", %zd/%d", result->size[i], n);
  95     fprintf(stderr, "\n");
  96
  97     fprintf(stderr, "%g\n", isl_int_get_d(result->n));
  98     fprintf(stderr, "%g", result->RE_sum[0]/isl_int_get_d(result->n));
  99     for (i = 1; i < nr_methods; ++i)
 100         fprintf(stderr, ", %g", result->RE_sum[i]/isl_int_get_d(result->n));
 101     fprintf(stderr, "\n");
 102 }
 103
 104 struct verify_point_bound {
 105         struct verify_point_data vpd;
 106         struct result_data *result;
 107         isl_pw_qpolynomial *pwqp;
 108         isl_pw_qpolynomial_fold **pwf;
 109 };
 110
 111 static int verify_point(__isl_take isl_point *pnt, void *user)
 112 {
 113         struct verify_point_bound *vpb = (struct verify_point_bound *) user;
 114         const struct verify_options *options = vpb->vpd.options;
 115         int i;
 116         unsigned nparam;
 117         isl_int max, min, exact, approx;
 118         isl_int n, d;
 119         isl_qpolynomial *opt;
 120         isl_pw_qpolynomial *pwqp;
 121         int cst;
 122
 123         vpb->vpd.n--;
 124
 125         isl_int_init(max);
 126         isl_int_init(min);
 127         isl_int_init(exact);
 128         isl_int_init(approx);
 129         isl_int_init(n);
 130         isl_int_init(d);
 131
 132         pwqp = isl_pw_qpolynomial_copy(vpb->pwqp);
 133
 134         nparam = isl_pw_qpolynomial_dim(pwqp, isl_dim_param);
 135         for (i = 0; i < nparam; ++i) {
 136                 isl_point_get_coordinate(pnt, isl_dim_param, i, &n);
 137                 pwqp = isl_pw_qpolynomial_fix_dim(pwqp, isl_dim_param, i, n);
 138         }
 139
 140         opt = isl_pw_qpolynomial_max(isl_pw_qpolynomial_copy(pwqp));
 141         cst = isl_qpolynomial_is_cst(opt, &n, &d);
 142         isl_qpolynomial_free(opt);
 143         assert(cst == 1);
 144         isl_int_fdiv_q(max, n, d);
 145
 146         opt = isl_pw_qpolynomial_min(pwqp);
 147         cst = isl_qpolynomial_is_cst(opt, &n, &d);
 148         isl_qpolynomial_free(opt);
 149         assert(cst == 1);
 150         isl_int_cdiv_q(min, n, d);
 151
 152         isl_int_sub(exact, max, min);
 153         isl_int_add_ui(exact, exact, 1);
 154
 155         if (options->print_all) {
 156                 fprintf(stderr, "max: "); isl_int_print(stderr, max, 0);
 157                 fprintf(stderr, ", min: "); isl_int_print(stderr, min, 0);
 158                 fprintf(stderr, ", range: "); isl_int_print(stderr, exact, 0);
 159         }
 160
 161         for (i = 0; i < nr_methods; ++i) {
 162                 double error;
 163
 164                 opt = isl_pw_qpolynomial_fold_eval(
 165                                 isl_pw_qpolynomial_fold_copy(vpb->pwf[2 * i]),
 166                                 isl_point_copy(pnt));
 167                 cst = isl_qpolynomial_is_cst(opt, &n, &d);
 168                 isl_qpolynomial_free(opt);
 169                 assert(cst == 1);
 170                 isl_int_fdiv_q(max, n, d);
 171
 172                 opt = isl_pw_qpolynomial_fold_eval(
 173                                 isl_pw_qpolynomial_fold_copy(vpb->pwf[2 * i + 1]),
 174                                 isl_point_copy(pnt));
 175                 cst = isl_qpolynomial_is_cst(opt, &n, &d);
 176                 isl_qpolynomial_free(opt);
 177                 assert(cst == 1);
 178                 isl_int_cdiv_q(min, n, d);
 179
 180                 isl_int_sub(approx, max, min);
 181                 isl_int_add_ui(approx, approx, 1);
 182                 if (options->print_all) {
 183                         fprintf(stderr, ", "); isl_int_print(stderr, approx, 0);
 184                 }
 185
 186                 assert(isl_int_ge(approx, exact));
 187                 isl_int_sub(approx, approx, exact);
 188
 189                 error = fabs(isl_int_get_d(approx)) / isl_int_get_d(exact);
 190                 if (options->print_all)
 191                         fprintf(stderr, " (%g)", error);
 192                 vpb->result->RE_sum[i] += error;
 193         }
 194
 195         if (options->print_all) {
 196                 fprintf(stderr, "\n");
 197         } else if ((vpb->vpd.n % vpb->vpd.s) == 0) {
 198                 printf("o");
 199                 fflush(stdout);
 200         }
 201
 202         isl_int_clear(max);
 203         isl_int_clear(min);
 204         isl_int_clear(exact);
 205         isl_int_clear(approx);
 206         isl_int_clear(n);
 207         isl_int_clear(d);
 208
 209         isl_point_free(pnt);
 210
 211         return vpb->vpd.n >= 1 ? 0 : -1;
 212 }
 213
 214 static void test(__isl_keep isl_pw_qpolynomial *pwqp,
 215         __isl_keep isl_pw_qpolynomial_fold **pwf, struct result_data *result,
 216         struct verify_options *options)
 217 {
 218         struct verify_point_bound vpb = { { options }, result };
 219         isl_set *context;
 220         int r;
 221
 222         vpb.pwf = pwf;
 223         vpb.pwqp = pwqp;
 224         context = isl_pw_qpolynomial_domain(isl_pw_qpolynomial_copy(vpb.pwqp));
 225         context = isl_set_params(context);
 226         context = verify_context_set_bounds(context, options);
 227
 228         r = verify_point_data_init(&vpb.vpd, context);
 229         assert(r == 0);
 230         isl_int_set_si(result->n, vpb.vpd.n);
 231
 232         isl_set_foreach_point(context, verify_point, &vpb);
 233         assert(!vpb.vpd.error);
 234
 235         isl_set_free(context);
 236
 237         verify_point_data_fini(&vpb.vpd);
 238 }
 239
 240 void handle(FILE *in, struct result_data *result, struct verify_options *options)
 241 {
 242     int i;
 243     isl_ctx *ctx = isl_ctx_alloc();
 244     isl_pw_qpolynomial *pwqp, *upper, *lower;
 245     isl_pw_qpolynomial_fold *pwf[2*nr_methods];
 246     struct isl_stream *s;
 247
 248     s = isl_stream_new_file(ctx, in);
 249     pwqp = isl_stream_read_pw_qpolynomial(s);
 250
 251     upper = isl_pw_qpolynomial_to_polynomial(isl_pw_qpolynomial_copy(pwqp), 1);
 252     lower = isl_pw_qpolynomial_to_polynomial(isl_pw_qpolynomial_copy(pwqp), -1);
 253
 254     for (i = 0; i < nr_methods; ++i) {
 255         int j;
 256         struct tms st_cpu;
 257         struct tms en_cpu;
 258
 259         times(&st_cpu);
 260         for (j = 0; j < 2; ++j) {
 261             isl_pw_qpolynomial *poly = j == 0 ? upper : lower;
 262             enum isl_fold type = j == 0 ? isl_fold_max : isl_fold_min;
 263             options->barvinok->isl->bound = methods[i].method;
 264             poly = isl_pw_qpolynomial_copy(poly);
 265             pwf[2*i+j] = isl_pw_qpolynomial_bound(poly, type, NULL);
 266         }
 267         times(&en_cpu);
 268         result->ticks[i] = time_diff(&en_cpu, &st_cpu);
 269         result->size[i] = isl_pw_qpolynomial_fold_size(pwf[2*i]);
 270         result->size[i] = isl_pw_qpolynomial_fold_size(pwf[2*i+1]);
 271         if (options->barvinok->verbose) {
 272             isl_printer *p = isl_printer_to_file(ctx, stdout);
 273             for (j = 0; j < 2; ++j) {
 274                 p = isl_printer_print_pw_qpolynomial_fold(p, pwf[2*i+j]);
 275                 p = isl_printer_end_line(p);
 276             }
 277             isl_printer_free(p);
 278         }
 279     }
 280     test(pwqp, pwf, result, options);
 281
 282     for (i = 0; i < 2*nr_methods; ++i)
 283         isl_pw_qpolynomial_fold_free(pwf[i]);
 284     isl_pw_qpolynomial_free(pwqp);
 285     isl_pw_qpolynomial_free(lower);
 286     isl_pw_qpolynomial_free(upper);
 287
 288     isl_stream_free(s);
 289     isl_ctx_free(ctx);
 290 }
 291
 292 int main(int argc, char **argv)
 293 {
 294     char path[PATH_MAX+1];
 295     struct result_data all_result;
 296     int n = 0;
 297     struct options *options = options_new_with_defaults();
 298
 299     argc = options_parse(options, argc, argv, ISL_ARG_ALL);
 300
 301     if (options->verify->M == INT_MIN)
 302         options->verify->M = 100;
 303     if (options->verify->m == INT_MAX)
 304         options->verify->m = -100;
 305
 306     result_data_init(&all_result);
 307
 308     while (fgets(path, sizeof(path), stdin)) {
 309         struct result_data result;
 310         FILE *in;
 311         int i;
 312
 313         ++n;
 314         result_data_init(&result);
 315         fprintf(stderr, "%s", path);
 316         *strchr(path, '\n') = '\0';
 317         in = fopen(path, "r");
 318         assert(in);
 319         handle(in, &result, options->verify);
 320         fclose(in);
 321
 322         if (!options->quiet)
 323             result_data_print(&result, 1);
 324
 325         isl_int_add(all_result.n, all_result.n, result.n);
 326         for (i = 0; i < nr_methods; ++i) {
 327             all_result.RE_sum[i] += result.RE_sum[i];
 328             all_result.ticks[i] += result.ticks[i];
 329             all_result.size[i] += result.size[i];
 330         }
 331
 332         result_data_clear(&result);
 333
 334         if (!options->quiet) {
 335             fprintf(stderr, "average\n");
 336             result_data_print(&all_result, n);
 337         }
 338     }
 339
 340     result_data_clear(&all_result);
 341
 342     options_free(options);
 343
 344     return 0;
 345 }