verify.c

   1 #include <assert.h>
   2 #include <limits.h>
   3 #include <stdlib.h>
   4 #include <isl/ctx.h>
   5 #include <isl/val.h>
   6 #include <isl/point.h>
   7 #include <isl/polynomial.h>
   8 #include <isl/set.h>
   9 #include <barvinok/polylib.h>
  10 #include <barvinok/isl.h>
  11 #include <barvinok/options.h>
  12 #include <barvinok/util.h>
  13 #include "verify.h"
  14
  15 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  16
  17 /* RANGE : normal range for evaluations (-RANGE -> RANGE) */
  18 #define RANGE 50
  19
  20 /* SRANGE : small range for evaluations */
  21 #define SRANGE 15
  22
  23 /* if dimension >= BIDDIM, use SRANGE */
  24 #define BIGDIM 5
  25
  26 /* VSRANGE : very small range for evaluations */
  27 #define VSRANGE 5
  28
  29 /* if dimension >= VBIDDIM, use VSRANGE */
  30 #define VBIGDIM 8
  31
  32 static int set_m(void *opt, long val)
  33 {
  34         struct verify_options *options = (struct verify_options *)opt;
  35         options->m = val;
  36         options->verify = 1;
  37         return 0;
  38 }
  39
  40 static int set_M(void *opt, long val)
  41 {
  42         struct verify_options *options = (struct verify_options *)opt;
  43         options->M = val;
  44         options->verify = 1;
  45         return 0;
  46 }
  47
  48 static int set_r(void *opt, long val)
  49 {
  50         struct verify_options *options = (struct verify_options *)opt;
  51         options->r = val;
  52         options->m = -val;
  53         options->M = val;
  54         options->verify = 1;
  55         return 0;
  56 }
  57
  58 ISL_ARGS_START(struct verify_options, verify_options_args)
  59 ISL_ARG_CHILD(struct verify_options, barvinok, NULL, &barvinok_options_args,
  60         NULL)
  61 ISL_ARG_BOOL(struct verify_options, verify, 'T', "verify", 0, NULL)
  62 ISL_ARG_BOOL(struct verify_options, exact, 'E', "exact", 0, NULL)
  63 ISL_ARG_BOOL(struct verify_options, print_all, 'A', "print-all", 0, NULL)
  64 ISL_ARG_BOOL(struct verify_options, continue_on_error, 'C',
  65         "continue-on-error", 0, NULL)
  66 ISL_ARG_USER_LONG(struct verify_options, m, 'm', "min", set_m, INT_MAX, NULL)
  67 ISL_ARG_USER_LONG(struct verify_options, M, 'M', "max", set_M, INT_MIN, NULL)
  68 ISL_ARG_USER_LONG(struct verify_options, r, 'r', "range", set_r, -1, NULL)
  69 ISL_ARGS_END
  70
  71 void verify_options_set_range(struct verify_options *options, int dim)
  72 {
  73     int r;
  74
  75     if (dim >= VBIGDIM)
  76         r = VSRANGE;
  77     else if (dim >= BIGDIM)
  78         r = SRANGE;
  79     else
  80         r = RANGE;
  81     /* If the user didn't set m or M, then we try to adjust the window
  82      * to the context in check_poly_context_scan.
  83      */
  84     if (options->m == INT_MAX && options->M == INT_MIN)
  85         options->r = r;
  86     if (options->M == INT_MIN)
  87         options->M = r;
  88     if (options->m == INT_MAX)
  89         options->m = -r;
  90
  91     if (options->verify && options->m > options->M) {
  92         fprintf(stderr,"Nothing to do: min > max !\n");
  93         exit(0);
  94     }
  95 }
  96
  97 /* Set the range of values for verification based on
  98  * the total dimensionality of "pwqp", if not already set by the user.
  99  */
 100 isl_stat verify_options_set_range_pwqp(struct verify_options *options,
 101         __isl_keep isl_pw_qpolynomial *pwqp)
 102 {
 103         isl_size total = isl_pw_qpolynomial_dim(pwqp, isl_dim_all);
 104         if (total < 0)
 105                 return isl_stat_error;
 106         verify_options_set_range(options, total);
 107         return isl_stat_ok;
 108 }
 109
 110 static Polyhedron *project_on(Polyhedron *P, int i)
 111 {
 112     unsigned dim = P->Dimension;
 113     Matrix *T;
 114     Polyhedron *I;
 115
 116     if (dim == 1)
 117         return Polyhedron_Copy(P);
 118
 119     T = Matrix_Alloc(2, dim+1);
 120     value_set_si(T->p[0][i], 1);
 121     value_set_si(T->p[1][dim], 1);
 122     I = Polyhedron_Image(P, T, P->NbConstraints);
 123     Matrix_Free(T);
 124     return I;
 125 }
 126
 127 static void set_bounds(Polyhedron *C, Value **rows, int i, unsigned nparam,
 128                         const struct verify_options *options)
 129 {
 130     Value min;
 131     Value max;
 132
 133     value_init(min);
 134     value_init(max);
 135     value_set_si(min, options->m);
 136     value_set_si(max, options->M);
 137
 138     if (options->r > 0) {
 139         Polyhedron *I = project_on(C, i);
 140         line_minmax(I, &min, &max);
 141         if (value_cmp_si(min, options->M) >= 0)
 142             value_add_int(max, min, options->r);
 143         else if (value_cmp_si(max, options->m) <= 0)
 144             value_sub_int(min, max, options->r);
 145         else {
 146             value_set_si(min, options->m);
 147             value_set_si(max, options->M);
 148         }
 149     }
 150
 151     value_set_si(rows[0][0], 1);
 152     value_set_si(rows[0][1+i], 1);
 153     value_oppose(rows[0][1+nparam], min);
 154     value_set_si(rows[1][0], 1);
 155     value_set_si(rows[1][1+i], -1);
 156     value_assign(rows[1][1+nparam], max);
 157
 158     value_clear(min);
 159     value_clear(max);
 160 }
 161
 162 Polyhedron *check_poly_context_scan(Polyhedron *P, Polyhedron **C,
 163                                     unsigned nparam,
 164                                     const struct verify_options *options)
 165 {
 166     int i;
 167     Matrix *MM;
 168     Polyhedron *CC, *CC2, *CS, *U;
 169     unsigned MaxRays = options->barvinok->MaxRays;
 170
 171     if (nparam <= 0)
 172         return NULL;
 173
 174     if (!P)
 175         CC = *C;
 176     else {
 177         CC = Polyhedron_Project(P, nparam);
 178         if (*C) {
 179             CC2 = DomainIntersection(*C, CC, MaxRays);
 180             Domain_Free(CC);
 181             CC = CC2;
 182         }
 183     }
 184
 185     /* Intersect context with range */
 186     MM = Matrix_Alloc(2*nparam, nparam+2);
 187     for (i = 0; i < nparam; ++i)
 188         set_bounds(CC, &MM->p[2*i], i, nparam, options);
 189     CC2 = AddConstraints(MM->p[0], 2*nparam, CC, options->barvinok->MaxRays);
 190     if (CC != *C)
 191         Domain_Free(CC);
 192     CC = CC2;
 193     U = Universe_Polyhedron(0);
 194     CS = Polyhedron_Scan(CC, U, MaxRays & POL_NO_DUAL ? 0 : MaxRays);
 195     Polyhedron_Free(U);
 196     *C = CC;
 197     Matrix_Free(MM);
 198     return CS;
 199 }
 200
 201 void check_poly_init(Polyhedron *C, struct verify_options *options)
 202 {
 203     int d, i;
 204
 205     if (options->print_all)
 206         return;
 207     if (C->Dimension <= 0)
 208         return;
 209
 210     d = options->M - options->m;
 211     if (d > 80)
 212         options->st = 1+d/80;
 213     else
 214         options->st = 1;
 215     for (i = options->m; i <= options->M; i += options->st)
 216         printf(".");
 217     printf( "\r" );
 218     fflush(stdout);
 219 }
 220
 221 static void print_rational(FILE *out, Value n, Value d)
 222 {
 223     value_print(out, VALUE_FMT, n);
 224     if (value_notone_p(d)) {
 225         fprintf(out, "/");
 226         value_print(out, VALUE_FMT, d);
 227     }
 228 }
 229
 230 void check_poly_print(int ok, int nparam, Value *z,
 231                       Value want_n, Value want_d,
 232                       Value got_n, Value got_d,
 233                       const char *op, const char *check,
 234                       const char *long_op,
 235                       const struct verify_options *options)
 236 {
 237     int k;
 238
 239     if (options->print_all) {
 240         printf("%s(", op);
 241         value_print(stdout, VALUE_FMT, z[0]);
 242         for (k = 1; k < nparam; ++k) {
 243             printf(", ");
 244             value_print(stdout, VALUE_FMT, z[k]);
 245         }
 246         printf(") = ");
 247         print_rational(stdout, got_n, got_d);
 248         printf(", %s = ", check);
 249         print_rational(stdout, want_n, want_d);
 250         printf(". ");
 251     }
 252
 253     if (!ok) {
 254         printf("\n");
 255         fflush(stdout);
 256         fprintf(stderr, "Error !\n");
 257         fprintf(stderr, "%s(", op);
 258         value_print(stderr, VALUE_FMT, z[0]);
 259         for (k = 1; k < nparam; ++k) {
 260             fprintf(stderr,", ");
 261             value_print(stderr, VALUE_FMT, z[k]);
 262         }
 263         fprintf(stderr, ") should be ");
 264         print_rational(stderr, want_n, want_d);
 265         fprintf(stderr, ", while %s gives ", long_op);
 266         print_rational(stderr, got_n, got_d);
 267         fprintf(stderr, ".\n");
 268     } else if (options->print_all)
 269         printf("OK.\n");
 270 }
 271
 272 /****************************************************/
 273 /* function check_poly :                            */
 274 /* scans the parameter space from Min to Max (all   */
 275 /* directions). Computes the number of points in    */
 276 /* the polytope using both methods, and compare them*/
 277 /* returns 1 on success                             */
 278 /****************************************************/
 279
 280 int check_poly(Polyhedron *CS, const struct check_poly_data *data,
 281                int nparam, int pos, Value *z,
 282                const struct verify_options *options)
 283 {
 284     if (pos == nparam) {
 285         if (!data->check(data, nparam, z, options) && !options->continue_on_error)
 286             return 0;
 287     } else {
 288         Value LB, UB;
 289         int ok;
 290         value_init(LB);
 291         value_init(UB);
 292         ok = !(lower_upper_bounds(1+pos, CS, z-1, &LB, &UB));
 293         assert(ok);
 294         for (; value_le(LB, UB); value_increment(LB, LB)) {
 295             if (!options->print_all) {
 296                 int k = VALUE_TO_INT(LB);
 297                 if (!pos && !(k % options->st)) {
 298                     printf("o");
 299                     fflush(stdout);
 300                 }
 301             }
 302
 303             value_assign(z[pos], LB);
 304             if (!check_poly(CS->next, data, nparam, pos+1, z, options)) {
 305                 value_clear(LB);
 306                 value_clear(UB);
 307                 return 0;
 308             }
 309         }
 310         value_set_si(z[pos], 0);
 311         value_clear(LB);
 312         value_clear(UB);
 313     }
 314     return 1;
 315 } /* check_poly */
 316
 317 void check_EP_set_scan(struct check_EP_data *data, Polyhedron *C,
 318                         unsigned MaxRays)
 319 {
 320     const evalue *EP = data->EP;
 321     int i;
 322     int n_S = 0;
 323
 324     for (i = 0; i < EP->x.p->size/2; ++i) {
 325         Polyhedron *A = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
 326         for (; A; A = A->next)
 327             ++n_S;
 328     }
 329
 330     data->n_S = n_S;
 331     data->S = ALLOCN(Polyhedron *, n_S);
 332     n_S = 0;
 333     for (i = 0; i < EP->x.p->size/2; ++i) {
 334         Polyhedron *A = EVALUE_DOMAIN(EP->x.p->arr[2*i]);
 335         for (; A; A = A->next) {
 336             Polyhedron *next = A->next;
 337             A->next = NULL;
 338             data->S[n_S++] = Polyhedron_Scan(A, C,
 339                                     MaxRays & POL_NO_DUAL ? 0 : MaxRays);
 340             A->next = next;
 341         }
 342     }
 343 }
 344
 345 void check_EP_clear_scan(struct check_EP_data *data)
 346 {
 347     int i;
 348
 349     for (i = 0; i < data->n_S; ++i)
 350         Domain_Free(data->S[i]);
 351     free(data->S);
 352 }
 353
 354 static int check_EP_on_poly(Polyhedron *P,
 355                             struct check_EP_data *data,
 356                             unsigned nvar, unsigned nparam,
 357                             struct verify_options *options)
 358 {
 359     Polyhedron *CS;
 360     unsigned MaxRays = options->barvinok->MaxRays;
 361     int ok = 1;
 362
 363     CS = check_poly_context_scan(NULL, &P, P->Dimension, options);
 364
 365     check_poly_init(P, options);
 366
 367     if (!(CS && emptyQ2(CS))) {
 368         check_EP_set_scan(data, P, MaxRays);
 369         ok = check_poly(CS, &data->cp, nparam, 0, data->cp.z+1+nvar, options);
 370         check_EP_clear_scan(data);
 371     }
 372
 373     if (!options->print_all)
 374         printf("\n");
 375
 376     if (CS) {
 377         Domain_Free(CS);
 378         Domain_Free(P);
 379     }
 380
 381     return ok;
 382 }
 383
 384 /*
 385  * Project on final dim dimensions
 386  */
 387 Polyhedron *DomainProject(Polyhedron *D, unsigned dim, unsigned MaxRays)
 388 {
 389     Polyhedron *P;
 390     Polyhedron *R;
 391
 392     R = Polyhedron_Project(D, dim);
 393     for (P = D->next; P; P = P->next) {
 394         Polyhedron *R2 = Polyhedron_Project(P, dim);
 395         Polyhedron *R3 = DomainUnion(R, R2, MaxRays);
 396         Polyhedron_Free(R2);
 397         Domain_Free(R);
 398         R = R3;
 399     }
 400     return R;
 401 }
 402
 403 static Polyhedron *evalue_parameter_domain(const evalue *e, unsigned nparam,
 404                                              unsigned MaxRays)
 405 {
 406     int i;
 407     Polyhedron *U = NULL;
 408
 409     if (EVALUE_IS_ZERO(*e))
 410         return Universe_Polyhedron(0);
 411
 412     assert(value_zero_p(e->d));
 413     assert(e->x.p->type == partition);
 414     assert(e->x.p->size >= 2);
 415
 416     for (i = 0; i < e->x.p->size/2; ++i) {
 417         Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2*i]);
 418         Polyhedron *P = DomainProject(D, nparam, MaxRays);
 419         if (!U) {
 420             U = P;
 421         } else {
 422             Polyhedron *T = U;
 423             U = DomainUnion(U, P, MaxRays);
 424             Domain_Free(P);
 425             Domain_Free(T);
 426         }
 427     }
 428     return U;
 429 }
 430
 431 int check_EP(struct check_EP_data *data, unsigned nvar, unsigned nparam,
 432              struct verify_options *options)
 433 {
 434     Vector *p;
 435     int ok = 1;
 436     Polyhedron *D, *P;
 437
 438     p = Vector_Alloc(nvar+nparam+2);
 439     value_set_si(p->p[nvar+nparam+1], 1);
 440
 441     data->cp.z = p->p;
 442
 443     D = evalue_parameter_domain(data->EP, nparam, options->barvinok->MaxRays);
 444
 445     for (P = D; P; P = P->next) {
 446         ok = check_EP_on_poly(P, data, nvar, nparam, options);
 447         if (!ok && !options->continue_on_error)
 448             break;
 449     }
 450
 451     Domain_Free(D);
 452     Vector_Free(p);
 453
 454     return ok;
 455 }
 456
 457 __isl_give isl_set *verify_context_set_bounds(__isl_take isl_set *set,
 458         const struct verify_options *options)
 459 {
 460         int i;
 461         unsigned nparam;
 462         isl_point *pt, *pt2;
 463         isl_set *box;
 464
 465         nparam = isl_set_dim(set, isl_dim_param);
 466
 467         if (options->r > 0) {
 468                 pt = isl_set_sample_point(isl_set_copy(set));
 469                 pt2 = isl_point_copy(pt);
 470
 471                 for (i = 0; i < nparam; ++i) {
 472                         pt = isl_point_add_ui(pt, isl_dim_param, i, options->r);
 473                         pt2 = isl_point_sub_ui(pt2, isl_dim_param, i, options->r);
 474                 }
 475         } else {
 476                 isl_ctx *ctx;
 477                 isl_val *v;
 478
 479                 ctx = isl_set_get_ctx(set);
 480                 pt = isl_point_zero(isl_set_get_space(set));
 481                 v = isl_val_int_from_si(ctx, options->m);
 482                 for (i = 0; i < nparam; ++i)
 483                         pt = isl_point_set_coordinate_val(pt, isl_dim_param, i,
 484                                                             isl_val_copy(v));
 485                 isl_val_free(v);
 486
 487                 pt2 = isl_point_zero(isl_set_get_space(set));
 488                 v = isl_val_int_from_si(ctx, options->M);
 489                 for (i = 0; i < nparam; ++i)
 490                         pt2 = isl_point_set_coordinate_val(pt2, isl_dim_param,
 491                                                             i, isl_val_copy(v));
 492                 isl_val_free(v);
 493         }
 494
 495         box = isl_set_box_from_points(pt, pt2);
 496
 497         return isl_set_intersect(set, box);
 498 }
 499
 500 int verify_point_data_init(struct verify_point_data *vpd,
 501         __isl_keep isl_set *context)
 502 {
 503         isl_val *v;
 504         int i;
 505
 506         context = isl_set_copy(context);
 507         context = isl_set_from_params(context);
 508         context = isl_set_move_dims(context, isl_dim_set, 0, isl_dim_param, 0,
 509                                         isl_set_dim(context, isl_dim_param));
 510         v = isl_pw_qpolynomial_max(isl_set_card(context));
 511         vpd->n = isl_val_cmp_si(v, 200) < 0 ? isl_val_get_num_si(v) : 200;
 512         isl_val_free(v);
 513
 514         if (!vpd->options->print_all) {
 515                 vpd->s = vpd->n < 80 ? 1 : 1 + vpd->n/80;
 516                 for (i = 0; i < vpd->n; i += vpd->s)
 517                         printf(".");
 518                 printf("\r");
 519                 fflush(stdout);
 520         }
 521
 522         vpd->error = !v ? -1 : 0;
 523
 524         return vpd->error;
 525 }
 526
 527 void verify_point_data_fini(struct verify_point_data *vpd)
 528 {
 529         if (!vpd->options->print_all)
 530                 printf("\n");
 531
 532         if (vpd->error)
 533                 fprintf(stderr, "Check failed !\n");
 534 }