add isl_qpolynomial_fold_substitute
[isl.git] / pip.c
blobd043a3654b12eeba9ff208355d5475aba8b84b89
1 /*
2 * Copyright 2008-2009 Katholieke Universiteit Leuven
4 * Use of this software is governed by the GNU LGPLv2.1 license
6 * Written by Sven Verdoolaege, K.U.Leuven, Departement
7 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
8 */
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>
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.
25 * Context constraints in PolyLib format
26 * -1
27 * Problem constraints in PolyLib format
28 * Optional list of options
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
37 struct options {
38 struct isl_options *isl;
39 unsigned verify;
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
48 ISL_ARG_DEF(options, struct options, options_arg)
50 static struct isl_basic_set *to_parameter_domain(struct isl_basic_set *context)
52 struct isl_dim *param_dim;
53 struct isl_basic_set *model;
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);
60 return context;
63 static isl_basic_set *construct_cube(isl_basic_set *context)
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;
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);
80 if (dim >= 8)
81 range = 4;
82 else if (dim >= 5)
83 range = 6;
84 else
85 range = 30;
87 isl_int_init(m);
88 isl_int_init(M);
90 isl_int_set_si(m, -range);
91 isl_int_set_si(M, range);
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);
100 isl_int_clear(m);
101 isl_int_clear(M);
103 return cube;
106 isl_basic_set *plug_in_parameters(isl_basic_set *bset, struct isl_vec *params)
108 int i;
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]);
114 bset = isl_basic_set_remove(bset, isl_dim_param, 0, params->size - 1);
116 isl_vec_free(params);
118 return bset;
121 isl_set *set_plug_in_parameters(isl_set *set, struct isl_vec *params)
123 int i;
125 for (i = 0; i < params->size - 1; ++i)
126 set = isl_set_fix(set, isl_dim_param, i, params->el[1 + i]);
128 set = isl_set_remove(set, isl_dim_param, 0, params->size - 1);
130 isl_vec_free(params);
132 return set;
135 /* Compute the lexicographically minimal (or maximal if max is set)
136 * element of bset for the given values of the parameters, by
137 * successively solving an ilp problem in each direction.
139 struct isl_vec *opt_at(struct isl_basic_set *bset,
140 struct isl_vec *params, int max)
142 unsigned dim;
143 struct isl_vec *opt;
144 struct isl_vec *obj;
145 int i;
147 dim = isl_basic_set_dim(bset, isl_dim_set);
149 bset = plug_in_parameters(bset, params);
151 if (isl_basic_set_fast_is_empty(bset)) {
152 opt = isl_vec_alloc(bset->ctx, 0);
153 isl_basic_set_free(bset);
154 return opt;
157 opt = isl_vec_alloc(bset->ctx, 1 + dim);
158 assert(opt);
160 obj = isl_vec_alloc(bset->ctx, 1 + dim);
161 assert(obj);
163 isl_int_set_si(opt->el[0], 1);
164 isl_int_set_si(obj->el[0], 0);
166 for (i = 0; i < dim; ++i) {
167 enum isl_lp_result res;
169 isl_seq_clr(obj->el + 1, dim);
170 isl_int_set_si(obj->el[1 + i], 1);
171 res = isl_basic_set_solve_ilp(bset, max, obj->el,
172 &opt->el[1 + i], NULL);
173 if (res == isl_lp_empty)
174 goto empty;
175 assert(res == isl_lp_ok);
176 bset = isl_basic_set_fix(bset, isl_dim_set, i, opt->el[1 + i]);
179 isl_basic_set_free(bset);
180 isl_vec_free(obj);
182 return opt;
183 empty:
184 isl_vec_free(opt);
185 opt = isl_vec_alloc(bset->ctx, 0);
186 isl_basic_set_free(bset);
187 isl_vec_free(obj);
189 return opt;
192 struct isl_scan_pip {
193 struct isl_scan_callback callback;
194 isl_basic_set *bset;
195 isl_set *sol;
196 isl_set *empty;
197 int stride;
198 int n;
199 int max;
202 /* Check if the "manually" computed optimum of bset at the "sample"
203 * values of the parameters agrees with the solution of pilp problem
204 * represented by the pair (sol, empty).
205 * In particular, if there is no solution for this value of the parameters,
206 * then it should be an element of the parameter domain "empty".
207 * Otherwise, the optimal solution, should be equal to the result of
208 * plugging in the value of the parameters in "sol".
210 static int scan_one(struct isl_scan_callback *callback,
211 __isl_take isl_vec *sample)
213 struct isl_scan_pip *sp = (struct isl_scan_pip *)callback;
214 struct isl_vec *opt;
216 opt = opt_at(isl_basic_set_copy(sp->bset), isl_vec_copy(sample), sp->max);
217 assert(opt);
219 if (opt->size == 0) {
220 isl_point *sample_pnt;
221 sample_pnt = isl_point_alloc(isl_set_get_dim(sp->empty), sample);
222 assert(isl_set_contains_point(sp->empty, sample_pnt));
223 isl_point_free(sample_pnt);
224 isl_vec_free(opt);
225 } else {
226 isl_set *sol;
227 isl_set *opt_set;
228 opt_set = isl_set_from_basic_set(isl_basic_set_from_vec(opt));
229 sol = set_plug_in_parameters(isl_set_copy(sp->sol), sample);
230 assert(isl_set_is_equal(opt_set, sol));
231 isl_set_free(sol);
232 isl_set_free(opt_set);
235 if (!(sp->n++ % sp->stride)) {
236 printf("o");
237 fflush(stdout);
240 return 0;
243 struct isl_scan_count {
244 struct isl_scan_callback callback;
245 int n;
248 static int count_one(struct isl_scan_callback *callback,
249 __isl_take isl_vec *sample)
251 struct isl_scan_count *sc = (struct isl_scan_count *)callback;
252 isl_vec_free(sample);
254 sc->n++;
256 return 0;
259 static void check_solution(isl_basic_set *bset, isl_basic_set *context,
260 isl_set *sol, isl_set *empty, int max)
262 isl_basic_set *cube;
263 struct isl_scan_count sc;
264 struct isl_scan_pip sp;
265 int i;
266 int r;
268 context = isl_basic_set_underlying_set(context);
270 cube = construct_cube(context);
271 context = isl_basic_set_intersect(context, cube);
273 sc.callback.add = count_one;
274 sc.n = 0;
276 r = isl_basic_set_scan(isl_basic_set_copy(context), &sc.callback);
277 assert (r == 0);
279 sp.callback.add = scan_one;
280 sp.bset = bset;
281 sp.sol = sol;
282 sp.empty = empty;
283 sp.n = 0;
284 sp.stride = sc.n > 70 ? 1 + (sc.n + 1)/70 : 1;
285 sp.max = max;
287 for (i = 0; i < sc.n; i += sp.stride)
288 printf(".");
289 printf("\r");
290 fflush(stdout);
292 r = isl_basic_set_scan(context, &sp.callback);
293 assert(r == 0);
295 printf("\n");
297 isl_basic_set_free(bset);
300 int main(int argc, char **argv)
302 struct isl_ctx *ctx;
303 struct isl_basic_set *context, *bset, *copy, *context_copy;
304 struct isl_set *set;
305 struct isl_set *empty;
306 int neg_one;
307 char s[1024];
308 int urs_parms = 0;
309 int urs_unknowns = 0;
310 int max = 0;
311 int rational = 0;
312 int n;
313 struct options *options;
315 options = options_new_with_defaults();
316 assert(options);
317 argc = options_parse(options, argc, argv, ISL_ARG_ALL);
319 ctx = isl_ctx_alloc_with_options(options_arg, options);
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));
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);
336 if (strncasecmp(s, "Urs_parms", 9) == 0)
337 urs_parms = 1;
338 if (strncasecmp(s, "Urs_unknowns", 12) == 0)
339 urs_unknowns = 1;
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)));
349 if (options->verify) {
350 copy = isl_basic_set_copy(bset);
351 context_copy = isl_basic_set_copy(context);
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);
359 if (options->verify) {
360 assert(!rational);
361 check_solution(copy, context_copy, set, empty, max);
362 } else {
363 isl_set_print(set, stdout, 0, ISL_FORMAT_ISL);
364 fprintf(stdout, "\n");
365 fprintf(stdout, "no solution: ");
366 isl_set_print(empty, stdout, 0, ISL_FORMAT_ISL);
367 fprintf(stdout, "\n");
370 isl_set_free(set);
371 isl_set_free(empty);
372 isl_ctx_free(ctx);
374 return 0;