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