?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
isl_set_project_out: always update dimension, even for empty sets
[isl.git] / isl_map_piplib.c
blob08a3104914ce8d2e514d2813a25e06a6047ba912
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 "isl_set.h"
11 #include "isl_map.h"
12 #include "isl_mat.h"
13 #include "isl_seq.h"
14 #include "isl_piplib.h"
15 #include "isl_map_piplib.h"
16 #include "isl_map_private.h"
17
18 static void copy_values_from(isl_int *dst, Entier *src, unsigned n)
19 {
20 int i;
21
22 for (i = 0; i < n; ++i)
23 entier_assign(dst[i], src[i]);
24 }
25
26 static void add_value(isl_int *dst, Entier *src)
27 {
28 mpz_add(*dst, *dst, *src);
29 }
30
31 static void copy_constraint_from(isl_int *dst, PipVector *src,
32 unsigned nparam, unsigned n_in, unsigned n_out,
33 unsigned extra, int *pos)
34 {
35 int i;
36
37 copy_values_from(dst, src->the_vector+src->nb_elements-1, 1);
38 copy_values_from(dst+1, src->the_vector, nparam+n_in);
39 isl_seq_clr(dst+1+nparam+n_in, n_out);
40 isl_seq_clr(dst+1+nparam+n_in+n_out, extra);
41 for (i = 0; i + n_in + nparam < src->nb_elements-1; ++i) {
42 int p = pos[i];
43 add_value(&dst[1+nparam+n_in+n_out+p],
44 &src->the_vector[n_in+nparam+i]);
45 }
46 }
47
48 static int add_inequality(struct isl_ctx *ctx,
49 struct isl_basic_map *bmap, int *pos, PipVector *vec)
50 {
51 unsigned nparam = isl_basic_map_n_param(bmap);
52 unsigned n_in = isl_basic_map_n_in(bmap);
53 unsigned n_out = isl_basic_map_n_out(bmap);
54 unsigned n_div = isl_basic_map_n_div(bmap);
55 int i = isl_basic_map_alloc_inequality(bmap);
56 if (i < 0)
57 return -1;
58 copy_constraint_from(bmap->ineq[i], vec,
59 nparam, n_in, n_out, n_div, pos);
60
61 return i;
62 }
63
64 /* For a div d = floor(f/m), add the constraints
65 *
66 * f - m d >= 0
67 * -(f-(n-1)) + m d >= 0
68 *
69 * Note that the second constraint is the negation of
70 *
71 * f - m d >= n
72 */
73 static int add_div_constraints(struct isl_ctx *ctx,
74 struct isl_basic_map *bmap, int *pos, PipNewparm *p, unsigned div)
75 {
76 int i, j;
77 unsigned total = isl_basic_map_total_dim(bmap);
78 unsigned div_pos = 1 + total - bmap->n_div + div;
79
80 i = add_inequality(ctx, bmap, pos, p->vector);
81 if (i < 0)
82 return -1;
83 copy_values_from(&bmap->ineq[i][div_pos], &p->deno, 1);
84 isl_int_neg(bmap->ineq[i][div_pos], bmap->ineq[i][div_pos]);
85
86 j = isl_basic_map_alloc_inequality(bmap);
87 if (j < 0)
88 return -1;
89 isl_seq_neg(bmap->ineq[j], bmap->ineq[i], 1 + total);
90 isl_int_add(bmap->ineq[j][0], bmap->ineq[j][0], bmap->ineq[j][div_pos]);
91 isl_int_sub_ui(bmap->ineq[j][0], bmap->ineq[j][0], 1);
92 return j;
93 }
94
95 static int add_equality(struct isl_ctx *ctx,
96 struct isl_basic_map *bmap, int *pos,
97 unsigned var, PipVector *vec)
98 {
99 int i;
100 unsigned nparam = isl_basic_map_n_param(bmap);
101 unsigned n_in = isl_basic_map_n_in(bmap);
102 unsigned n_out = isl_basic_map_n_out(bmap);
103
104 isl_assert(ctx, var < n_out, return -1);
105
106 i = isl_basic_map_alloc_equality(bmap);
107 if (i < 0)
108 return -1;
109 copy_constraint_from(bmap->eq[i], vec,
110 nparam, n_in, n_out, bmap->extra, pos);
111 isl_int_set_si(bmap->eq[i][1+nparam+n_in+var], -1);
112
113 return i;
114 }
115
116 static int find_div(struct isl_ctx *ctx,
117 struct isl_basic_map *bmap, int *pos, PipNewparm *p)
118 {
119 int i, j;
120 unsigned nparam = isl_basic_map_n_param(bmap);
121 unsigned n_in = isl_basic_map_n_in(bmap);
122 unsigned n_out = isl_basic_map_n_out(bmap);
123
124 i = isl_basic_map_alloc_div(bmap);
125 if (i < 0)
126 return -1;
127
128 copy_constraint_from(bmap->div[i]+1, p->vector,
129 nparam, n_in, n_out, bmap->extra, pos);
130
131 copy_values_from(bmap->div[i], &p->deno, 1);
132 for (j = 0; j < i; ++j)
133 if (isl_seq_eq(bmap->div[i], bmap->div[j],
134 1+1+isl_basic_map_total_dim(bmap)+j)) {
135 isl_basic_map_free_div(bmap, 1);
136 return j;
137 }
138
139 if (add_div_constraints(ctx, bmap, pos, p, i) < 0)
140 return -1;
141
142 return i;
143 }
144
145 /* Count some properties of a quast
146 * - maximal number of new parameters
147 * - maximal depth
148 * - total number of solutions
149 * - total number of empty branches
150 */
151 static void quast_count(PipQuast *q, int *maxnew, int depth, int *maxdepth,
152 int *sol, int *nosol)
153 {
154 PipNewparm *p;
155
156 for (p = q->newparm; p; p = p->next)
157 if (p->rank > *maxnew)
158 *maxnew = p->rank;
159 if (q->condition) {
160 if (++depth > *maxdepth)
161 *maxdepth = depth;
162 quast_count(q->next_else, maxnew, depth, maxdepth, sol, nosol);
163 quast_count(q->next_then, maxnew, depth, maxdepth, sol, nosol);
164 } else {
165 if (q->list)
166 ++(*sol);
167 else
168 ++(*nosol);
169 }
170 }
171
172 /*
173 * pos: array of length bmap->set.extra, mapping each of the existential
174 * variables PIP proposes to an existential variable in bmap
175 * bmap: collects the currently active constraints
176 * rest: collects the empty leaves of the quast (if not NULL)
177 */
178 struct scan_data {
179 struct isl_ctx *ctx;
180 struct isl_basic_map *bmap;
181 struct isl_set **rest;
182 int *pos;
183 };
184
185 /*
186 * New existentially quantified variables are places after the existing ones.
187 */
188 static struct isl_map *scan_quast_r(struct scan_data *data, PipQuast *q,
189 struct isl_map *map)
190 {
191 PipNewparm *p;
192 struct isl_basic_map *bmap = data->bmap;
193 unsigned old_n_div = bmap->n_div;
194 unsigned nparam = isl_basic_map_n_param(bmap);
195 unsigned n_in = isl_basic_map_n_in(bmap);
196 unsigned n_out = isl_basic_map_n_out(bmap);
197
198 if (!map)
199 goto error;
200
201 for (p = q->newparm; p; p = p->next) {
202 int pos;
203 unsigned pip_param = nparam + n_in;
204
205 pos = find_div(data->ctx, bmap, data->pos, p);
206 if (pos < 0)
207 goto error;
208 data->pos[p->rank - pip_param] = pos;
209 }
210
211 if (q->condition) {
212 int pos = add_inequality(data->ctx, bmap, data->pos,
213 q->condition);
214 if (pos < 0)
215 goto error;
216 map = scan_quast_r(data, q->next_then, map);
217
218 if (isl_inequality_negate(bmap, pos))
219 goto error;
220 map = scan_quast_r(data, q->next_else, map);
221
222 if (isl_basic_map_free_inequality(bmap, 1))
223 goto error;
224 } else if (q->list) {
225 PipList *l;
226 int j;
227 /* if bmap->n_out is zero, we are only interested in the domains
228 * where a solution exists and not in the actual solution
229 */
230 for (j = 0, l = q->list; j < n_out && l; ++j, l = l->next)
231 if (add_equality(data->ctx, bmap, data->pos, j,
232 l->vector) < 0)
233 goto error;
234 map = isl_map_add(map, isl_basic_map_copy(bmap));
235 if (isl_basic_map_free_equality(bmap, n_out))
236 goto error;
237 } else if (data->rest) {
238 struct isl_basic_set *bset;
239 bset = isl_basic_set_from_basic_map(isl_basic_map_copy(bmap));
240 bset = isl_basic_set_drop_dims(bset, n_in, n_out);
241 if (!bset)
242 goto error;
243 *data->rest = isl_set_add(*data->rest, bset);
244 }
245
246 if (isl_basic_map_free_inequality(bmap, 2*(bmap->n_div - old_n_div)))
247 goto error;
248 if (isl_basic_map_free_div(bmap, bmap->n_div - old_n_div))
249 goto error;
250 return map;
251 error:
252 isl_map_free(map);
253 return NULL;
254 }
255
256 /*
257 * Returns a map of dimension "keep_dim" with "context" as domain and
258 * as range the first "isl_dim_size(keep_dim, isl_dim_out)" variables
259 * in the quast lists.
260 */
261 static struct isl_map *isl_map_from_quast(struct isl_ctx *ctx, PipQuast *q,
262 struct isl_dim *keep_dim,
263 struct isl_basic_set *context,
264 struct isl_set **rest)
265 {
266 int pip_param;
267 int nexist;
268 int max_depth;
269 int n_sol, n_nosol;
270 struct scan_data data;
271 struct isl_map *map = NULL;
272 struct isl_dim *dims;
273 unsigned nparam;
274 unsigned dim;
275 unsigned keep;
276
277 data.ctx = ctx;
278 data.rest = rest;
279 data.bmap = NULL;
280 data.pos = NULL;
281
282 if (!context || !keep_dim)
283 goto error;
284
285 dim = isl_basic_set_n_dim(context);
286 nparam = isl_basic_set_n_param(context);
287 keep = isl_dim_size(keep_dim, isl_dim_out);
288 pip_param = nparam + dim;
289
290 max_depth = 0;
291 n_sol = 0;
292 n_nosol = 0;
293 nexist = pip_param-1;
294 quast_count(q, &nexist, 0, &max_depth, &n_sol, &n_nosol);
295 nexist -= pip_param-1;
296
297 if (rest) {
298 *rest = isl_set_alloc_dim(isl_dim_copy(context->dim), n_nosol,
299 ISL_MAP_DISJOINT);
300 if (!*rest)
301 goto error;
302 }
303 map = isl_map_alloc_dim(isl_dim_copy(keep_dim), n_sol,
304 ISL_MAP_DISJOINT);
305 if (!map)
306 goto error;
307
308 dims = isl_dim_reverse(isl_dim_copy(context->dim));
309 data.bmap = isl_basic_map_from_basic_set(context, dims);
310 data.bmap = isl_basic_map_extend_dim(data.bmap,
311 keep_dim, nexist, keep, max_depth+2*nexist);
312 if (!data.bmap)
313 goto error2;
314
315 if (data.bmap->extra) {
316 int i;
317 data.pos = isl_alloc_array(ctx, int, data.bmap->extra);
318 if (!data.pos)
319 goto error;
320 for (i = 0; i < data.bmap->n_div; ++i)
321 data.pos[i] = i;
322 }
323
324 map = scan_quast_r(&data, q, map);
325 map = isl_map_finalize(map);
326 if (!map)
327 goto error2;
328 if (rest) {
329 *rest = isl_set_finalize(*rest);
330 if (!*rest)
331 goto error2;
332 }
333 isl_basic_map_free(data.bmap);
334 if (data.pos)
335 free(data.pos);
336 return map;
337 error:
338 isl_basic_set_free(context);
339 isl_dim_free(keep_dim);
340 error2:
341 if (data.pos)
342 free(data.pos);
343 isl_basic_map_free(data.bmap);
344 isl_map_free(map);
345 if (rest) {
346 isl_set_free(*rest);
347 *rest = NULL;
348 }
349 return NULL;
350 }
351
352 static void copy_values_to(Entier *dst, isl_int *src, unsigned n)
353 {
354 int i;
355
356 for (i = 0; i < n; ++i)
357 entier_assign(dst[i], src[i]);
358 }
359
360 static void copy_constraint_to(Entier *dst, isl_int *src,
361 unsigned pip_param, unsigned pip_var,
362 unsigned extra_front, unsigned extra_back)
363 {
364 copy_values_to(dst+1+extra_front+pip_var+pip_param+extra_back, src, 1);
365 copy_values_to(dst+1+extra_front+pip_var, src+1, pip_param);
366 copy_values_to(dst+1+extra_front, src+1+pip_param, pip_var);
367 }
368
369 PipMatrix *isl_basic_map_to_pip(struct isl_basic_map *bmap, unsigned pip_param,
370 unsigned extra_front, unsigned extra_back)
371 {
372 int i;
373 unsigned nrow;
374 unsigned ncol;
375 PipMatrix *M;
376 unsigned off;
377 unsigned pip_var = isl_basic_map_total_dim(bmap) - pip_param;
378
379 nrow = extra_front + bmap->n_eq + bmap->n_ineq;
380 ncol = 1 + extra_front + pip_var + pip_param + extra_back + 1;
381 M = pip_matrix_alloc(nrow, ncol);
382 if (!M)
383 return NULL;
384
385 off = extra_front;
386 for (i = 0; i < bmap->n_eq; ++i) {
387 entier_set_si(M->p[off+i][0], 0);
388 copy_constraint_to(M->p[off+i], bmap->eq[i],
389 pip_param, pip_var, extra_front, extra_back);
390 }
391 off += bmap->n_eq;
392 for (i = 0; i < bmap->n_ineq; ++i) {
393 entier_set_si(M->p[off+i][0], 1);
394 copy_constraint_to(M->p[off+i], bmap->ineq[i],
395 pip_param, pip_var, extra_front, extra_back);
396 }
397 return M;
398 }
399
400 PipMatrix *isl_basic_set_to_pip(struct isl_basic_set *bset, unsigned pip_param,
401 unsigned extra_front, unsigned extra_back)
402 {
403 return isl_basic_map_to_pip((struct isl_basic_map *)bset,
404 pip_param, extra_front, extra_back);
405 }
406
407 struct isl_map *isl_pip_basic_map_lexopt(
408 struct isl_basic_map *bmap, struct isl_basic_set *dom,
409 struct isl_set **empty, int max)
410 {
411 PipOptions *options;
412 PipQuast *sol;
413 struct isl_map *map;
414 struct isl_ctx *ctx;
415 PipMatrix *domain = NULL, *context = NULL;
416 unsigned nparam, n_in, n_out;
417
418 bmap = isl_basic_map_detect_equalities(bmap);
419 if (!bmap || !dom)
420 goto error;
421
422 ctx = bmap->ctx;
423 isl_assert(ctx, isl_basic_map_compatible_domain(bmap, dom), goto error);
424 nparam = isl_basic_map_n_param(bmap);
425 n_in = isl_basic_map_n_in(bmap);
426 n_out = isl_basic_map_n_out(bmap);
427
428 domain = isl_basic_map_to_pip(bmap, nparam + n_in, 0, dom->n_div);
429 if (!domain)
430 goto error;
431 context = isl_basic_map_to_pip((struct isl_basic_map *)dom, 0, 0, 0);
432 if (!context)
433 goto error;
434
435 options = pip_options_init();
436 options->Simplify = 1;
437 options->Maximize = max;
438 options->Urs_unknowns = -1;
439 options->Urs_parms = -1;
440 sol = pip_solve(domain, context, -1, options);
441
442 if (sol) {
443 struct isl_basic_set *copy;
444 copy = isl_basic_set_copy(dom);
445 map = isl_map_from_quast(ctx, sol,
446 isl_dim_copy(bmap->dim), copy, empty);
447 } else {
448 map = isl_map_empty_like_basic_map(bmap);
449 if (empty)
450 *empty = NULL;
451 }
452 if (!map)
453 goto error;
454 if (map->n == 0 && empty) {
455 isl_set_free(*empty);
456 *empty = isl_set_from_basic_set(dom);
457 } else
458 isl_basic_set_free(dom);
459 isl_basic_map_free(bmap);
460
461 pip_quast_free(sol);
462 pip_options_free(options);
463 pip_matrix_free(domain);
464 pip_matrix_free(context);
465
466 return map;
467 error:
468 if (domain)
469 pip_matrix_free(domain);
470 if (context)
471 pip_matrix_free(context);
472 isl_basic_map_free(bmap);
473 isl_basic_set_free(dom);
474 return NULL;
475 }
Integer set library