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