Merge branch 'maint'
[isl.git] / isl_map_piplib.c
blob17181539bdc417b5d41fa04f777ca07ea1d452d3
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 <isl_ctx_private.h>
11 #include <isl_map_private.h>
12 #include <isl/set.h>
13 #include <isl/map.h>
14 #include <isl/mat.h>
15 #include <isl/seq.h>
16 #include "isl_piplib.h"
17 #include "isl_map_piplib.h"
19 static void copy_values_from(isl_int *dst, Entier *src, unsigned n)
21 int i;
23 for (i = 0; i < n; ++i)
24 entier_assign(dst[i], src[i]);
27 static void add_value(isl_int *dst, Entier *src)
29 mpz_add(*dst, *dst, *src);
32 static void copy_constraint_from(isl_int *dst, PipVector *src,
33 unsigned nparam, unsigned n_in, unsigned n_out,
34 unsigned extra, int *pos)
36 int i;
38 copy_values_from(dst, src->the_vector+src->nb_elements-1, 1);
39 copy_values_from(dst+1, src->the_vector, nparam+n_in);
40 isl_seq_clr(dst+1+nparam+n_in, n_out);
41 isl_seq_clr(dst+1+nparam+n_in+n_out, extra);
42 for (i = 0; i + n_in + nparam < src->nb_elements-1; ++i) {
43 int p = pos[i];
44 add_value(&dst[1+nparam+n_in+n_out+p],
45 &src->the_vector[n_in+nparam+i]);
49 static int add_inequality(struct isl_ctx *ctx,
50 struct isl_basic_map *bmap, int *pos, PipVector *vec)
52 unsigned nparam = isl_basic_map_n_param(bmap);
53 unsigned n_in = isl_basic_map_n_in(bmap);
54 unsigned n_out = isl_basic_map_n_out(bmap);
55 unsigned n_div = isl_basic_map_n_div(bmap);
56 int i = isl_basic_map_alloc_inequality(bmap);
57 if (i < 0)
58 return -1;
59 copy_constraint_from(bmap->ineq[i], vec,
60 nparam, n_in, n_out, n_div, pos);
62 return i;
65 /* For a div d = floor(f/m), add the constraints
67 * f - m d >= 0
68 * -(f-(n-1)) + m d >= 0
70 * Note that the second constraint is the negation of
72 * f - m d >= n
74 static int add_div_constraints(struct isl_ctx *ctx,
75 struct isl_basic_map *bmap, int *pos, PipNewparm *p, unsigned div)
77 int i, j;
78 unsigned total = isl_basic_map_total_dim(bmap);
79 unsigned div_pos = 1 + total - bmap->n_div + div;
81 i = add_inequality(ctx, bmap, pos, p->vector);
82 if (i < 0)
83 return -1;
84 copy_values_from(&bmap->ineq[i][div_pos], &p->deno, 1);
85 isl_int_neg(bmap->ineq[i][div_pos], bmap->ineq[i][div_pos]);
87 j = isl_basic_map_alloc_inequality(bmap);
88 if (j < 0)
89 return -1;
90 isl_seq_neg(bmap->ineq[j], bmap->ineq[i], 1 + total);
91 isl_int_add(bmap->ineq[j][0], bmap->ineq[j][0], bmap->ineq[j][div_pos]);
92 isl_int_sub_ui(bmap->ineq[j][0], bmap->ineq[j][0], 1);
93 return j;
96 static int add_equality(struct isl_ctx *ctx,
97 struct isl_basic_map *bmap, int *pos,
98 unsigned var, PipVector *vec)
100 int i;
101 unsigned nparam = isl_basic_map_n_param(bmap);
102 unsigned n_in = isl_basic_map_n_in(bmap);
103 unsigned n_out = isl_basic_map_n_out(bmap);
105 isl_assert(ctx, var < n_out, return -1);
107 i = isl_basic_map_alloc_equality(bmap);
108 if (i < 0)
109 return -1;
110 copy_constraint_from(bmap->eq[i], vec,
111 nparam, n_in, n_out, bmap->extra, pos);
112 isl_int_set_si(bmap->eq[i][1+nparam+n_in+var], -1);
114 return i;
117 static int find_div(struct isl_ctx *ctx,
118 struct isl_basic_map *bmap, int *pos, PipNewparm *p)
120 int i, j;
121 unsigned nparam = isl_basic_map_n_param(bmap);
122 unsigned n_in = isl_basic_map_n_in(bmap);
123 unsigned n_out = isl_basic_map_n_out(bmap);
125 i = isl_basic_map_alloc_div(bmap);
126 if (i < 0)
127 return -1;
129 copy_constraint_from(bmap->div[i]+1, p->vector,
130 nparam, n_in, n_out, bmap->extra, pos);
132 copy_values_from(bmap->div[i], &p->deno, 1);
133 for (j = 0; j < i; ++j)
134 if (isl_seq_eq(bmap->div[i], bmap->div[j],
135 1+1+isl_basic_map_total_dim(bmap)+j)) {
136 isl_basic_map_free_div(bmap, 1);
137 return j;
140 if (add_div_constraints(ctx, bmap, pos, p, i) < 0)
141 return -1;
143 return i;
146 /* Count some properties of a quast
147 * - maximal number of new parameters
148 * - maximal depth
149 * - total number of solutions
150 * - total number of empty branches
152 static void quast_count(PipQuast *q, int *maxnew, int depth, int *maxdepth,
153 int *sol, int *nosol)
155 PipNewparm *p;
157 for (p = q->newparm; p; p = p->next)
158 if (p->rank > *maxnew)
159 *maxnew = p->rank;
160 if (q->condition) {
161 if (++depth > *maxdepth)
162 *maxdepth = depth;
163 quast_count(q->next_else, maxnew, depth, maxdepth, sol, nosol);
164 quast_count(q->next_then, maxnew, depth, maxdepth, sol, nosol);
165 } else {
166 if (q->list)
167 ++(*sol);
168 else
169 ++(*nosol);
174 * pos: array of length bmap->set.extra, mapping each of the existential
175 * variables PIP proposes to an existential variable in bmap
176 * bmap: collects the currently active constraints
177 * rest: collects the empty leaves of the quast (if not NULL)
179 struct scan_data {
180 struct isl_ctx *ctx;
181 struct isl_basic_map *bmap;
182 struct isl_set **rest;
183 int *pos;
187 * New existentially quantified variables are places after the existing ones.
189 static struct isl_map *scan_quast_r(struct scan_data *data, PipQuast *q,
190 struct isl_map *map)
192 PipNewparm *p;
193 struct isl_basic_map *bmap = data->bmap;
194 unsigned old_n_div = bmap->n_div;
195 unsigned nparam = isl_basic_map_n_param(bmap);
196 unsigned n_in = isl_basic_map_n_in(bmap);
197 unsigned n_out = isl_basic_map_n_out(bmap);
199 if (!map)
200 goto error;
202 for (p = q->newparm; p; p = p->next) {
203 int pos;
204 unsigned pip_param = nparam + n_in;
206 pos = find_div(data->ctx, bmap, data->pos, p);
207 if (pos < 0)
208 goto error;
209 data->pos[p->rank - pip_param] = pos;
212 if (q->condition) {
213 int pos = add_inequality(data->ctx, bmap, data->pos,
214 q->condition);
215 if (pos < 0)
216 goto error;
217 map = scan_quast_r(data, q->next_then, map);
219 if (isl_inequality_negate(bmap, pos))
220 goto error;
221 map = scan_quast_r(data, q->next_else, map);
223 if (isl_basic_map_free_inequality(bmap, 1))
224 goto error;
225 } else if (q->list) {
226 PipList *l;
227 int j;
228 /* if bmap->n_out is zero, we are only interested in the domains
229 * where a solution exists and not in the actual solution
231 for (j = 0, l = q->list; j < n_out && l; ++j, l = l->next)
232 if (add_equality(data->ctx, bmap, data->pos, j,
233 l->vector) < 0)
234 goto error;
235 map = isl_map_add_basic_map(map, isl_basic_map_copy(bmap));
236 if (isl_basic_map_free_equality(bmap, n_out))
237 goto error;
238 } else if (data->rest) {
239 struct isl_basic_set *bset;
240 bset = isl_basic_set_from_basic_map(isl_basic_map_copy(bmap));
241 bset = isl_basic_set_drop_dims(bset, n_in, n_out);
242 if (!bset)
243 goto error;
244 *data->rest = isl_set_add_basic_set(*data->rest, bset);
247 if (isl_basic_map_free_inequality(bmap, 2*(bmap->n_div - old_n_div)))
248 goto error;
249 if (isl_basic_map_free_div(bmap, bmap->n_div - old_n_div))
250 goto error;
251 return map;
252 error:
253 isl_map_free(map);
254 return NULL;
258 * Returns a map of dimension "keep_dim" with "context" as domain and
259 * as range the first "isl_space_dim(keep_dim, isl_dim_out)" variables
260 * in the quast lists.
262 static struct isl_map *isl_map_from_quast(struct isl_ctx *ctx, PipQuast *q,
263 isl_space *keep_dim,
264 struct isl_basic_set *context,
265 struct isl_set **rest)
267 int pip_param;
268 int nexist;
269 int max_depth;
270 int n_sol, n_nosol;
271 struct scan_data data;
272 struct isl_map *map = NULL;
273 isl_space *dims;
274 unsigned nparam;
275 unsigned dim;
276 unsigned keep;
278 data.ctx = ctx;
279 data.rest = rest;
280 data.bmap = NULL;
281 data.pos = NULL;
283 if (!context || !keep_dim)
284 goto error;
286 dim = isl_basic_set_n_dim(context);
287 nparam = isl_basic_set_n_param(context);
288 keep = isl_space_dim(keep_dim, isl_dim_out);
289 pip_param = nparam + dim;
291 max_depth = 0;
292 n_sol = 0;
293 n_nosol = 0;
294 nexist = pip_param-1;
295 quast_count(q, &nexist, 0, &max_depth, &n_sol, &n_nosol);
296 nexist -= pip_param-1;
298 if (rest) {
299 *rest = isl_set_alloc_space(isl_space_copy(context->dim), n_nosol,
300 ISL_MAP_DISJOINT);
301 if (!*rest)
302 goto error;
304 map = isl_map_alloc_space(isl_space_copy(keep_dim), n_sol,
305 ISL_MAP_DISJOINT);
306 if (!map)
307 goto error;
309 dims = isl_space_reverse(isl_space_copy(context->dim));
310 data.bmap = isl_basic_map_from_basic_set(context, dims);
311 data.bmap = isl_basic_map_extend_space(data.bmap,
312 keep_dim, nexist, keep, max_depth+2*nexist);
313 if (!data.bmap)
314 goto error2;
316 if (data.bmap->extra) {
317 int i;
318 data.pos = isl_alloc_array(ctx, int, data.bmap->extra);
319 if (!data.pos)
320 goto error;
321 for (i = 0; i < data.bmap->n_div; ++i)
322 data.pos[i] = i;
325 map = scan_quast_r(&data, q, map);
326 map = isl_map_finalize(map);
327 if (!map)
328 goto error2;
329 if (rest) {
330 *rest = isl_set_finalize(*rest);
331 if (!*rest)
332 goto error2;
334 isl_basic_map_free(data.bmap);
335 if (data.pos)
336 free(data.pos);
337 return map;
338 error:
339 isl_basic_set_free(context);
340 isl_space_free(keep_dim);
341 error2:
342 if (data.pos)
343 free(data.pos);
344 isl_basic_map_free(data.bmap);
345 isl_map_free(map);
346 if (rest) {
347 isl_set_free(*rest);
348 *rest = NULL;
350 return NULL;
353 static void copy_values_to(Entier *dst, isl_int *src, unsigned n)
355 int i;
357 for (i = 0; i < n; ++i)
358 entier_assign(dst[i], src[i]);
361 static void copy_constraint_to(Entier *dst, isl_int *src,
362 unsigned pip_param, unsigned pip_var,
363 unsigned extra_front, unsigned extra_back)
365 copy_values_to(dst+1+extra_front+pip_var+pip_param+extra_back, src, 1);
366 copy_values_to(dst+1+extra_front+pip_var, src+1, pip_param);
367 copy_values_to(dst+1+extra_front, src+1+pip_param, pip_var);
370 PipMatrix *isl_basic_map_to_pip(struct isl_basic_map *bmap, unsigned pip_param,
371 unsigned extra_front, unsigned extra_back)
373 int i;
374 unsigned nrow;
375 unsigned ncol;
376 PipMatrix *M;
377 unsigned off;
378 unsigned pip_var = isl_basic_map_total_dim(bmap) - pip_param;
380 nrow = extra_front + bmap->n_eq + bmap->n_ineq;
381 ncol = 1 + extra_front + pip_var + pip_param + extra_back + 1;
382 M = pip_matrix_alloc(nrow, ncol);
383 if (!M)
384 return NULL;
386 off = extra_front;
387 for (i = 0; i < bmap->n_eq; ++i) {
388 entier_set_si(M->p[off+i][0], 0);
389 copy_constraint_to(M->p[off+i], bmap->eq[i],
390 pip_param, pip_var, extra_front, extra_back);
392 off += bmap->n_eq;
393 for (i = 0; i < bmap->n_ineq; ++i) {
394 entier_set_si(M->p[off+i][0], 1);
395 copy_constraint_to(M->p[off+i], bmap->ineq[i],
396 pip_param, pip_var, extra_front, extra_back);
398 return M;
401 PipMatrix *isl_basic_set_to_pip(struct isl_basic_set *bset, unsigned pip_param,
402 unsigned extra_front, unsigned extra_back)
404 return isl_basic_map_to_pip((struct isl_basic_map *)bset,
405 pip_param, extra_front, extra_back);
408 struct isl_map *isl_pip_basic_map_lexopt(
409 struct isl_basic_map *bmap, struct isl_basic_set *dom,
410 struct isl_set **empty, int max)
412 PipOptions *options;
413 PipQuast *sol;
414 struct isl_map *map;
415 struct isl_ctx *ctx;
416 PipMatrix *domain = NULL, *context = NULL;
417 unsigned nparam, n_in, n_out;
419 bmap = isl_basic_map_detect_equalities(bmap);
420 if (!bmap || !dom)
421 goto error;
423 ctx = bmap->ctx;
424 isl_assert(ctx, isl_basic_map_compatible_domain(bmap, dom), goto error);
425 nparam = isl_basic_map_n_param(bmap);
426 n_in = isl_basic_map_n_in(bmap);
427 n_out = isl_basic_map_n_out(bmap);
429 domain = isl_basic_map_to_pip(bmap, nparam + n_in, 0, dom->n_div);
430 if (!domain)
431 goto error;
432 context = isl_basic_map_to_pip((struct isl_basic_map *)dom, 0, 0, 0);
433 if (!context)
434 goto error;
436 options = pip_options_init();
437 options->Simplify = 1;
438 options->Maximize = max;
439 options->Urs_unknowns = -1;
440 options->Urs_parms = -1;
441 sol = pip_solve(domain, context, -1, options);
443 if (sol) {
444 struct isl_basic_set *copy;
445 copy = isl_basic_set_copy(dom);
446 map = isl_map_from_quast(ctx, sol,
447 isl_space_copy(bmap->dim), copy, empty);
448 } else {
449 map = isl_map_empty_like_basic_map(bmap);
450 if (empty)
451 *empty = NULL;
453 if (!map)
454 goto error;
455 if (map->n == 0 && empty) {
456 isl_set_free(*empty);
457 *empty = isl_set_from_basic_set(dom);
458 } else
459 isl_basic_set_free(dom);
460 isl_basic_map_free(bmap);
462 pip_quast_free(sol);
463 pip_options_free(options);
464 pip_matrix_free(domain);
465 pip_matrix_free(context);
467 return map;
468 error:
469 if (domain)
470 pip_matrix_free(domain);
471 if (context)
472 pip_matrix_free(context);
473 isl_basic_map_free(bmap);
474 isl_basic_set_free(dom);
475 return NULL;