isl_map_intersect_{domain,range}: allow intersection with parameter constraints
[isl.git] / isl_map_piplib.c
blob257b4607565fe0930e48f1ed45838473e47d2391
1 #include "isl_set.h"
2 #include "isl_map.h"
3 #include "isl_seq.h"
4 #include "isl_piplib.h"
5 #include "isl_map_piplib.h"
6 #include "isl_map_private.h"
8 static void copy_values_from(isl_int *dst, Entier *src, unsigned n)
10 int i;
12 for (i = 0; i < n; ++i)
13 entier_assign(dst[i], src[i]);
16 static void add_value(isl_int *dst, Entier *src)
18 mpz_add(*dst, *dst, *src);
21 static void copy_constraint_from(isl_int *dst, PipVector *src,
22 unsigned nparam, unsigned n_in, unsigned n_out,
23 unsigned extra, int *pos)
25 int i;
27 copy_values_from(dst, src->the_vector+src->nb_elements-1, 1);
28 copy_values_from(dst+1, src->the_vector, nparam+n_in);
29 isl_seq_clr(dst+1+nparam+n_in, n_out);
30 isl_seq_clr(dst+1+nparam+n_in+n_out, extra);
31 for (i = 0; i + n_in + nparam < src->nb_elements-1; ++i) {
32 int p = pos[i];
33 add_value(&dst[1+nparam+n_in+n_out+p],
34 &src->the_vector[n_in+nparam+i]);
38 static int add_inequality(struct isl_ctx *ctx,
39 struct isl_basic_map *bmap, int *pos, PipVector *vec)
41 unsigned nparam = isl_basic_map_n_param(bmap);
42 unsigned n_in = isl_basic_map_n_in(bmap);
43 unsigned n_out = isl_basic_map_n_out(bmap);
44 unsigned n_div = isl_basic_map_n_div(bmap);
45 int i = isl_basic_map_alloc_inequality(bmap);
46 if (i < 0)
47 return -1;
48 copy_constraint_from(bmap->ineq[i], vec,
49 nparam, n_in, n_out, n_div, pos);
51 return i;
54 /* For a div d = floor(f/m), add the constraints
56 * f - m d >= 0
57 * -(f-(n-1)) + m d >= 0
59 * Note that the second constraint is the negation of
61 * f - m d >= n
63 static int add_div_constraints(struct isl_ctx *ctx,
64 struct isl_basic_map *bmap, int *pos, PipNewparm *p, unsigned div)
66 int i, j;
67 unsigned total = isl_basic_map_total_dim(bmap);
68 unsigned div_pos = 1 + total - bmap->n_div + div;
70 i = add_inequality(ctx, bmap, pos, p->vector);
71 if (i < 0)
72 return -1;
73 copy_values_from(&bmap->ineq[i][div_pos], &p->deno, 1);
74 isl_int_neg(bmap->ineq[i][div_pos], bmap->ineq[i][div_pos]);
76 j = isl_basic_map_alloc_inequality(bmap);
77 if (j < 0)
78 return -1;
79 isl_seq_neg(bmap->ineq[j], bmap->ineq[i], 1 + total);
80 isl_int_add(bmap->ineq[j][0], bmap->ineq[j][0], bmap->ineq[j][div_pos]);
81 isl_int_sub_ui(bmap->ineq[j][0], bmap->ineq[j][0], 1);
82 return j;
85 static int add_equality(struct isl_ctx *ctx,
86 struct isl_basic_map *bmap, int *pos,
87 unsigned var, PipVector *vec)
89 int i;
90 unsigned nparam = isl_basic_map_n_param(bmap);
91 unsigned n_in = isl_basic_map_n_in(bmap);
92 unsigned n_out = isl_basic_map_n_out(bmap);
94 isl_assert(ctx, var < n_out, return -1);
96 i = isl_basic_map_alloc_equality(bmap);
97 if (i < 0)
98 return -1;
99 copy_constraint_from(bmap->eq[i], vec,
100 nparam, n_in, n_out, bmap->extra, pos);
101 isl_int_set_si(bmap->eq[i][1+nparam+n_in+var], -1);
103 return i;
106 static int find_div(struct isl_ctx *ctx,
107 struct isl_basic_map *bmap, int *pos, PipNewparm *p)
109 int i, j;
110 unsigned nparam = isl_basic_map_n_param(bmap);
111 unsigned n_in = isl_basic_map_n_in(bmap);
112 unsigned n_out = isl_basic_map_n_out(bmap);
114 i = isl_basic_map_alloc_div(bmap);
115 if (i < 0)
116 return -1;
118 copy_constraint_from(bmap->div[i]+1, p->vector,
119 nparam, n_in, n_out, bmap->extra, pos);
121 copy_values_from(bmap->div[i], &p->deno, 1);
122 for (j = 0; j < i; ++j)
123 if (isl_seq_eq(bmap->div[i], bmap->div[j],
124 1+1+isl_basic_map_total_dim(bmap)+j)) {
125 isl_basic_map_free_div(bmap, 1);
126 return j;
129 if (add_div_constraints(ctx, bmap, pos, p, i) < 0)
130 return -1;
132 return i;
135 /* Count some properties of a quast
136 * - maximal number of new parameters
137 * - maximal depth
138 * - total number of solutions
139 * - total number of empty branches
141 static void quast_count(PipQuast *q, int *maxnew, int depth, int *maxdepth,
142 int *sol, int *nosol)
144 PipNewparm *p;
146 for (p = q->newparm; p; p = p->next)
147 if (p->rank > *maxnew)
148 *maxnew = p->rank;
149 if (q->condition) {
150 if (++depth > *maxdepth)
151 *maxdepth = depth;
152 quast_count(q->next_else, maxnew, depth, maxdepth, sol, nosol);
153 quast_count(q->next_then, maxnew, depth, maxdepth, sol, nosol);
154 } else {
155 if (q->list)
156 ++(*sol);
157 else
158 ++(*nosol);
163 * pos: array of length bmap->set.extra, mapping each of the existential
164 * variables PIP proposes to an existential variable in bmap
165 * bmap: collects the currently active constraints
166 * rest: collects the empty leaves of the quast (if not NULL)
168 struct scan_data {
169 struct isl_ctx *ctx;
170 struct isl_basic_map *bmap;
171 struct isl_set **rest;
172 int *pos;
176 * New existentially quantified variables are places after the existing ones.
178 static struct isl_map *scan_quast_r(struct scan_data *data, PipQuast *q,
179 struct isl_map *map)
181 PipNewparm *p;
182 struct isl_basic_map *bmap = data->bmap;
183 unsigned old_n_div = bmap->n_div;
184 unsigned nparam = isl_basic_map_n_param(bmap);
185 unsigned n_in = isl_basic_map_n_in(bmap);
186 unsigned n_out = isl_basic_map_n_out(bmap);
188 if (!map)
189 goto error;
191 for (p = q->newparm; p; p = p->next) {
192 int pos;
193 unsigned pip_param = nparam + n_in;
195 pos = find_div(data->ctx, bmap, data->pos, p);
196 if (pos < 0)
197 goto error;
198 data->pos[p->rank - pip_param] = pos;
201 if (q->condition) {
202 int pos = add_inequality(data->ctx, bmap, data->pos,
203 q->condition);
204 if (pos < 0)
205 goto error;
206 map = scan_quast_r(data, q->next_then, map);
208 if (isl_inequality_negate(bmap, pos))
209 goto error;
210 map = scan_quast_r(data, q->next_else, map);
212 if (isl_basic_map_free_inequality(bmap, 1))
213 goto error;
214 } else if (q->list) {
215 PipList *l;
216 int j;
217 /* if bmap->n_out is zero, we are only interested in the domains
218 * where a solution exists and not in the actual solution
220 for (j = 0, l = q->list; j < n_out && l; ++j, l = l->next)
221 if (add_equality(data->ctx, bmap, data->pos, j,
222 l->vector) < 0)
223 goto error;
224 map = isl_map_add(map, isl_basic_map_copy(bmap));
225 if (isl_basic_map_free_equality(bmap, n_out))
226 goto error;
227 } else if (data->rest) {
228 struct isl_basic_set *bset;
229 bset = isl_basic_set_from_basic_map(isl_basic_map_copy(bmap));
230 bset = isl_basic_set_drop_dims(bset, n_in, n_out);
231 if (!bset)
232 goto error;
233 *data->rest = isl_set_add(*data->rest, bset);
236 if (isl_basic_map_free_inequality(bmap, 2*(bmap->n_div - old_n_div)))
237 goto error;
238 if (isl_basic_map_free_div(bmap, bmap->n_div - old_n_div))
239 goto error;
240 return map;
241 error:
242 isl_map_free(map);
243 return NULL;
247 * Returns a map with "context" as domain and as range the first
248 * "keep" variables in the quast lists.
250 static struct isl_map *isl_map_from_quast(struct isl_ctx *ctx, PipQuast *q,
251 unsigned keep,
252 struct isl_basic_set *context,
253 struct isl_set **rest)
255 int pip_param;
256 int nexist;
257 int max_depth;
258 int n_sol, n_nosol;
259 struct scan_data data;
260 struct isl_map *map;
261 struct isl_dim *dims;
262 unsigned nparam;
263 unsigned dim;
265 data.ctx = ctx;
266 data.rest = rest;
267 data.bmap = NULL;
268 data.pos = NULL;
270 if (!context)
271 goto error;
273 dim = isl_basic_set_n_dim(context);
274 nparam = isl_basic_set_n_param(context);
275 pip_param = nparam + dim;
277 max_depth = 0;
278 n_sol = 0;
279 n_nosol = 0;
280 nexist = pip_param-1;
281 quast_count(q, &nexist, 0, &max_depth, &n_sol, &n_nosol);
282 nexist -= pip_param-1;
284 if (rest) {
285 *rest = isl_set_alloc(data.ctx, nparam, dim, n_nosol,
286 ISL_MAP_DISJOINT);
287 if (!*rest)
288 goto error;
290 map = isl_map_alloc(data.ctx, nparam, dim, keep, n_sol,
291 ISL_MAP_DISJOINT);
292 if (!map)
293 goto error;
295 dims = isl_dim_reverse(isl_dim_copy(context->dim));
296 data.bmap = isl_basic_map_from_basic_set(context, dims);
297 data.bmap = isl_basic_map_extend(data.bmap,
298 nparam, dim, keep, nexist, keep, max_depth+2*nexist);
299 if (!data.bmap)
300 goto error;
302 if (data.bmap->extra) {
303 int i;
304 data.pos = isl_alloc_array(ctx, int, data.bmap->extra);
305 if (!data.pos)
306 goto error;
307 for (i = 0; i < data.bmap->n_div; ++i)
308 data.pos[i] = i;
311 map = scan_quast_r(&data, q, map);
312 map = isl_map_finalize(map);
313 if (!map)
314 goto error;
315 if (rest) {
316 *rest = isl_set_finalize(*rest);
317 if (!*rest)
318 goto error;
320 isl_basic_map_free(data.bmap);
321 if (data.pos)
322 free(data.pos);
323 return map;
324 error:
325 if (data.pos)
326 free(data.pos);
327 isl_basic_map_free(data.bmap);
328 isl_map_free(map);
329 if (rest) {
330 isl_set_free(*rest);
331 *rest = NULL;
333 return NULL;
336 static void copy_values_to(Entier *dst, isl_int *src, unsigned n)
338 int i;
340 for (i = 0; i < n; ++i)
341 entier_assign(dst[i], src[i]);
344 static void copy_constraint_to(Entier *dst, isl_int *src,
345 unsigned pip_param, unsigned pip_var,
346 unsigned extra_front, unsigned extra_back)
348 copy_values_to(dst+1+extra_front+pip_var+pip_param+extra_back, src, 1);
349 copy_values_to(dst+1+extra_front+pip_var, src+1, pip_param);
350 copy_values_to(dst+1+extra_front, src+1+pip_param, pip_var);
353 PipMatrix *isl_basic_map_to_pip(struct isl_basic_map *bmap, unsigned pip_param,
354 unsigned extra_front, unsigned extra_back)
356 int i;
357 unsigned nrow;
358 unsigned ncol;
359 PipMatrix *M;
360 unsigned off;
361 unsigned pip_var = isl_basic_map_total_dim(bmap) - pip_param;
363 nrow = extra_front + bmap->n_eq + bmap->n_ineq;
364 ncol = 1 + extra_front + pip_var + pip_param + extra_back + 1;
365 M = pip_matrix_alloc(nrow, ncol);
366 if (!M)
367 return NULL;
369 off = extra_front;
370 for (i = 0; i < bmap->n_eq; ++i) {
371 entier_set_si(M->p[off+i][0], 0);
372 copy_constraint_to(M->p[off+i], bmap->eq[i],
373 pip_param, pip_var, extra_front, extra_back);
375 off += bmap->n_eq;
376 for (i = 0; i < bmap->n_ineq; ++i) {
377 entier_set_si(M->p[off+i][0], 1);
378 copy_constraint_to(M->p[off+i], bmap->ineq[i],
379 pip_param, pip_var, extra_front, extra_back);
381 return M;
384 static struct isl_map *extremum_on(
385 struct isl_basic_map *bmap, struct isl_basic_set *dom,
386 struct isl_set **empty, int max)
388 PipOptions *options;
389 PipQuast *sol;
390 struct isl_map *map;
391 struct isl_ctx *ctx;
392 PipMatrix *domain = NULL, *context = NULL;
393 unsigned nparam, n_in, n_out;
395 if (!bmap || !dom)
396 goto error;
398 ctx = bmap->ctx;
399 isl_assert(ctx, isl_basic_map_compatible_domain(bmap, dom), goto error);
400 nparam = isl_basic_map_n_param(bmap);
401 n_in = isl_basic_map_n_in(bmap);
402 n_out = isl_basic_map_n_out(bmap);
404 domain = isl_basic_map_to_pip(bmap, nparam + n_in, 0, dom->n_div);
405 if (!domain)
406 goto error;
407 context = isl_basic_map_to_pip((struct isl_basic_map *)dom, 0, 0, 0);
408 if (!context)
409 goto error;
411 options = pip_options_init();
412 options->Simplify = 1;
413 options->Maximize = max;
414 options->Urs_unknowns = -1;
415 options->Urs_parms = -1;
416 sol = pip_solve(domain, context, -1, options);
418 if (sol) {
419 struct isl_basic_set *copy;
420 copy = isl_basic_set_copy(dom);
421 map = isl_map_from_quast(ctx, sol, n_out, copy, empty);
422 } else {
423 map = isl_map_empty_like_basic_map(bmap);
424 if (empty)
425 *empty = NULL;
427 if (!map)
428 goto error;
429 if (map->n == 0 && empty) {
430 isl_set_free(*empty);
431 *empty = isl_set_from_basic_set(dom);
432 } else
433 isl_basic_set_free(dom);
434 isl_basic_map_free(bmap);
436 pip_quast_free(sol);
437 pip_options_free(options);
438 pip_matrix_free(domain);
439 pip_matrix_free(context);
441 return map;
442 error:
443 if (domain)
444 pip_matrix_free(domain);
445 if (context)
446 pip_matrix_free(context);
447 isl_basic_map_free(bmap);
448 isl_basic_set_free(dom);
449 return NULL;
452 struct isl_map *isl_pip_basic_map_lexmax(
453 struct isl_basic_map *bmap, struct isl_basic_set *dom,
454 struct isl_set **empty)
456 return extremum_on(bmap, dom, empty, 1);
459 struct isl_map *isl_pip_basic_map_lexmin(
460 struct isl_basic_map *bmap, struct isl_basic_set *dom,
461 struct isl_set **empty)
463 return extremum_on(bmap, dom, empty, 0);
466 struct isl_map *isl_pip_basic_map_compute_divs(struct isl_basic_map *bmap)
468 PipMatrix *domain = NULL, *context = NULL;
469 PipOptions *options;
470 PipQuast *sol;
471 struct isl_ctx *ctx;
472 struct isl_dim *dim;
473 struct isl_map *map;
474 struct isl_set *set;
475 struct isl_basic_set *dom;
476 unsigned nparam;
477 unsigned n_in;
478 unsigned n_out;
480 if (!bmap)
481 goto error;
483 ctx = bmap->ctx;
484 nparam = isl_basic_map_n_param(bmap);
485 n_in = isl_basic_map_n_in(bmap);
486 n_out = isl_basic_map_n_out(bmap);
488 domain = isl_basic_map_to_pip(bmap, nparam + n_in + n_out, 0, 0);
489 if (!domain)
490 goto error;
491 context = pip_matrix_alloc(0, nparam + n_in + n_out + 2);
492 if (!context)
493 goto error;
495 options = pip_options_init();
496 options->Simplify = 1;
497 options->Urs_unknowns = -1;
498 options->Urs_parms = -1;
499 sol = pip_solve(domain, context, -1, options);
501 dom = isl_basic_set_alloc(ctx, nparam, n_in + n_out, 0, 0, 0);
502 map = isl_map_from_quast(ctx, sol, 0, dom, NULL);
504 pip_quast_free(sol);
505 pip_options_free(options);
506 pip_matrix_free(domain);
507 pip_matrix_free(context);
509 dim = isl_dim_copy(bmap->dim);
510 isl_basic_map_free(bmap);
512 set = isl_map_domain(map);
514 return isl_map_from_set(set, dim);
515 error:
516 if (domain)
517 pip_matrix_free(domain);
518 if (context)
519 pip_matrix_free(context);
520 isl_basic_set_free(dom);
521 isl_basic_map_free(bmap);
522 return NULL;