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