pip.c

   1 /*
   2  * Copyright 2008-2009 Katholieke Universiteit Leuven
   3  *
   4  * Use of this software is governed by the GNU LGPLv2.1 license
   5  *
   6  * Written by Sven Verdoolaege, K.U.Leuven, Departement
   7  * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
   8  */
   9
  10 #include <assert.h>
  11 #include <string.h>
  12 #include <isl/set.h>
  13 #include "isl_tab.h"
  14 #include "isl_map_private.h"
  15 #include "isl_sample.h"
  16 #include "isl_scan.h"
  17 #include <isl/seq.h>
  18 #include <isl/ilp.h>
  19 #include <isl_point_private.h>
  20
  21 /* The input of this program is the same as that of the "example" program
  22  * from the PipLib distribution, except that the "big parameter column"
  23  * should always be -1.
  24  *
  25  * Context constraints in PolyLib format
  26  * -1
  27  * Problem constraints in PolyLib format
  28  * Optional list of options
  29  *
  30  * The options are
  31  *      Maximize        compute maximum instead of minimum
  32  *      Rational        compute rational optimum instead of integer optimum
  33  *      Urs_parms       don't assume parameters are non-negative
  34  *      Urs_unknowns    don't assume unknowns are non-negative
  35  */
  36
  37 struct options {
  38         struct isl_options      *isl;
  39         unsigned                 verify;
  40 };
  41
  42 struct isl_arg options_arg[] = {
  43 ISL_ARG_CHILD(struct options, isl, "isl", isl_options_arg, "isl options")
  44 ISL_ARG_BOOL(struct options, verify, 'T', "verify", 0, NULL)
  45 ISL_ARG_END
  46 };
  47
  48 ISL_ARG_DEF(options, struct options, options_arg)
  49
  50 static __isl_give isl_basic_set *set_bounds(__isl_take isl_basic_set *bset)
  51 {
  52         unsigned nparam;
  53         int i, r;
  54         isl_point *pt, *pt2;
  55         isl_basic_set *box;
  56
  57         nparam = isl_basic_set_dim(bset, isl_dim_param);
  58         r = nparam >= 8 ? 4 : nparam >= 5 ? 6 : 30;
  59
  60         pt = isl_basic_set_sample_point(isl_basic_set_copy(bset));
  61         pt2 = isl_point_copy(pt);
  62
  63         for (i = 0; i < nparam; ++i) {
  64                 pt = isl_point_add_ui(pt, isl_dim_param, i, r);
  65                 pt2 = isl_point_sub_ui(pt2, isl_dim_param, i, r);
  66         }
  67
  68         box = isl_basic_set_box_from_points(pt, pt2);
  69
  70         return isl_basic_set_intersect(bset, box);
  71 }
  72
  73 static struct isl_basic_set *to_parameter_domain(struct isl_basic_set *context)
  74 {
  75         struct isl_dim *param_dim;
  76         struct isl_basic_set *model;
  77
  78         param_dim = isl_dim_set_alloc(context->ctx,
  79                                         isl_basic_set_n_dim(context), 0);
  80         model = isl_basic_set_empty(param_dim);
  81         context = isl_basic_set_from_underlying_set(context, model);
  82
  83         return context;
  84 }
  85
  86 isl_basic_set *plug_in_parameters(isl_basic_set *bset, struct isl_vec *params)
  87 {
  88         int i;
  89
  90         for (i = 0; i < params->size - 1; ++i)
  91                 bset = isl_basic_set_fix(bset,
  92                                          isl_dim_param, i, params->el[1 + i]);
  93
  94         bset = isl_basic_set_remove_dims(bset,
  95                                          isl_dim_param, 0, params->size - 1);
  96
  97         isl_vec_free(params);
  98
  99         return bset;
 100 }
 101
 102 isl_set *set_plug_in_parameters(isl_set *set, struct isl_vec *params)
 103 {
 104         int i;
 105
 106         for (i = 0; i < params->size - 1; ++i)
 107                 set = isl_set_fix(set, isl_dim_param, i, params->el[1 + i]);
 108
 109         set = isl_set_remove_dims(set, isl_dim_param, 0, params->size - 1);
 110
 111         isl_vec_free(params);
 112
 113         return set;
 114 }
 115
 116 /* Compute the lexicographically minimal (or maximal if max is set)
 117  * element of bset for the given values of the parameters, by
 118  * successively solving an ilp problem in each direction.
 119  */
 120 struct isl_vec *opt_at(struct isl_basic_set *bset,
 121         struct isl_vec *params, int max)
 122 {
 123         unsigned dim;
 124         struct isl_vec *opt;
 125         struct isl_vec *obj;
 126         int i;
 127
 128         dim = isl_basic_set_dim(bset, isl_dim_set);
 129
 130         bset = plug_in_parameters(bset, params);
 131
 132         if (isl_basic_set_fast_is_empty(bset)) {
 133                 opt = isl_vec_alloc(bset->ctx, 0);
 134                 isl_basic_set_free(bset);
 135                 return opt;
 136         }
 137
 138         opt = isl_vec_alloc(bset->ctx, 1 + dim);
 139         assert(opt);
 140
 141         obj = isl_vec_alloc(bset->ctx, 1 + dim);
 142         assert(obj);
 143
 144         isl_int_set_si(opt->el[0], 1);
 145         isl_int_set_si(obj->el[0], 0);
 146
 147         for (i = 0; i < dim; ++i) {
 148                 enum isl_lp_result res;
 149
 150                 isl_seq_clr(obj->el + 1, dim);
 151                 isl_int_set_si(obj->el[1 + i], 1);
 152                 res = isl_basic_set_solve_ilp(bset, max, obj->el,
 153                                                 &opt->el[1 + i], NULL);
 154                 if (res == isl_lp_empty)
 155                         goto empty;
 156                 assert(res == isl_lp_ok);
 157                 bset = isl_basic_set_fix(bset, isl_dim_set, i, opt->el[1 + i]);
 158         }
 159
 160         isl_basic_set_free(bset);
 161         isl_vec_free(obj);
 162
 163         return opt;
 164 empty:
 165         isl_vec_free(opt);
 166         opt = isl_vec_alloc(bset->ctx, 0);
 167         isl_basic_set_free(bset);
 168         isl_vec_free(obj);
 169
 170         return opt;
 171 }
 172
 173 struct isl_scan_pip {
 174         struct isl_scan_callback callback;
 175         isl_basic_set *bset;
 176         isl_set *sol;
 177         isl_set *empty;
 178         int stride;
 179         int n;
 180         int max;
 181 };
 182
 183 /* Check if the "manually" computed optimum of bset at the "sample"
 184  * values of the parameters agrees with the solution of pilp problem
 185  * represented by the pair (sol, empty).
 186  * In particular, if there is no solution for this value of the parameters,
 187  * then it should be an element of the parameter domain "empty".
 188  * Otherwise, the optimal solution, should be equal to the result of
 189  * plugging in the value of the parameters in "sol".
 190  */
 191 static int scan_one(struct isl_scan_callback *callback,
 192         __isl_take isl_vec *sample)
 193 {
 194         struct isl_scan_pip *sp = (struct isl_scan_pip *)callback;
 195         struct isl_vec *opt;
 196
 197         opt = opt_at(isl_basic_set_copy(sp->bset), isl_vec_copy(sample), sp->max);
 198         assert(opt);
 199
 200         if (opt->size == 0) {
 201                 isl_point *sample_pnt;
 202                 sample_pnt = isl_point_alloc(isl_set_get_dim(sp->empty), sample);
 203                 assert(isl_set_contains_point(sp->empty, sample_pnt));
 204                 isl_point_free(sample_pnt);
 205                 isl_vec_free(opt);
 206         } else {
 207                 isl_set *sol;
 208                 isl_set *opt_set;
 209                 opt_set = isl_set_from_basic_set(isl_basic_set_from_vec(opt));
 210                 sol = set_plug_in_parameters(isl_set_copy(sp->sol), sample);
 211                 assert(isl_set_is_equal(opt_set, sol));
 212                 isl_set_free(sol);
 213                 isl_set_free(opt_set);
 214         }
 215
 216         if (!(sp->n++ % sp->stride)) {
 217                 printf("o");
 218                 fflush(stdout);
 219         }
 220
 221         return 0;
 222 }
 223
 224 struct isl_scan_count {
 225         struct isl_scan_callback callback;
 226         int n;
 227 };
 228
 229 static int count_one(struct isl_scan_callback *callback,
 230         __isl_take isl_vec *sample)
 231 {
 232         struct isl_scan_count *sc = (struct isl_scan_count *)callback;
 233         isl_vec_free(sample);
 234
 235         sc->n++;
 236
 237         return 0;
 238 }
 239
 240 static void check_solution(isl_basic_set *bset, isl_basic_set *context,
 241         isl_set *sol, isl_set *empty, int max)
 242 {
 243         struct isl_scan_count sc;
 244         struct isl_scan_pip sp;
 245         int i;
 246         int r;
 247
 248         context = set_bounds(context);
 249         context = isl_basic_set_underlying_set(context);
 250
 251         sc.callback.add = count_one;
 252         sc.n = 0;
 253
 254         r = isl_basic_set_scan(isl_basic_set_copy(context), &sc.callback);
 255         assert (r == 0);
 256
 257         sp.callback.add = scan_one;
 258         sp.bset = bset;
 259         sp.sol = sol;
 260         sp.empty = empty;
 261         sp.n = 0;
 262         sp.stride = sc.n > 70 ? 1 + (sc.n + 1)/70 : 1;
 263         sp.max = max;
 264
 265         for (i = 0; i < sc.n; i += sp.stride)
 266                 printf(".");
 267         printf("\r");
 268         fflush(stdout);
 269
 270         r = isl_basic_set_scan(context, &sp.callback);
 271         assert(r == 0);
 272
 273         printf("\n");
 274
 275         isl_basic_set_free(bset);
 276 }
 277
 278 int main(int argc, char **argv)
 279 {
 280         struct isl_ctx *ctx;
 281         struct isl_basic_set *context, *bset, *copy, *context_copy;
 282         struct isl_set *set;
 283         struct isl_set *empty;
 284         int neg_one;
 285         char s[1024];
 286         int urs_parms = 0;
 287         int urs_unknowns = 0;
 288         int max = 0;
 289         int rational = 0;
 290         int n;
 291         struct options *options;
 292
 293         options = options_new_with_defaults();
 294         assert(options);
 295         argc = options_parse(options, argc, argv, ISL_ARG_ALL);
 296
 297         ctx = isl_ctx_alloc_with_options(options_arg, options);
 298
 299         context = isl_basic_set_read_from_file(ctx, stdin, 0);
 300         assert(context);
 301         n = fscanf(stdin, "%d", &neg_one);
 302         assert(n == 1);
 303         assert(neg_one == -1);
 304         bset = isl_basic_set_read_from_file(ctx, stdin,
 305                 isl_basic_set_dim(context, isl_dim_set));
 306
 307         while (fgets(s, sizeof(s), stdin)) {
 308                 if (strncasecmp(s, "Maximize", 8) == 0)
 309                         max = 1;
 310                 if (strncasecmp(s, "Rational", 8) == 0) {
 311                         rational = 1;
 312                         bset = isl_basic_set_set_rational(bset);
 313                 }
 314                 if (strncasecmp(s, "Urs_parms", 9) == 0)
 315                         urs_parms = 1;
 316                 if (strncasecmp(s, "Urs_unknowns", 12) == 0)
 317                         urs_unknowns = 1;
 318         }
 319         if (!urs_parms)
 320                 context = isl_basic_set_intersect(context,
 321                 isl_basic_set_positive_orthant(isl_basic_set_get_dim(context)));
 322         context = to_parameter_domain(context);
 323         if (!urs_unknowns)
 324                 bset = isl_basic_set_intersect(bset,
 325                 isl_basic_set_positive_orthant(isl_basic_set_get_dim(bset)));
 326
 327         if (options->verify) {
 328                 copy = isl_basic_set_copy(bset);
 329                 context_copy = isl_basic_set_copy(context);
 330         }
 331
 332         if (max)
 333                 set = isl_basic_set_partial_lexmax(bset, context, &empty);
 334         else
 335                 set = isl_basic_set_partial_lexmin(bset, context, &empty);
 336
 337         if (options->verify) {
 338                 assert(!rational);
 339                 check_solution(copy, context_copy, set, empty, max);
 340         } else {
 341                 isl_set_print(set, stdout, 0, ISL_FORMAT_ISL);
 342                 fprintf(stdout, "\n");
 343                 fprintf(stdout, "no solution: ");
 344                 isl_set_print(empty, stdout, 0, ISL_FORMAT_ISL);
 345                 fprintf(stdout, "\n");
 346         }
 347
 348         isl_set_free(set);
 349         isl_set_free(empty);
 350         isl_ctx_free(ctx);
 351
 352         return 0;
 353 }