isl_constraint: drop constraint iteration
[isl.git] / isl_farkas.c
blob2dc617eacabab2a33071d204946e65148df89b69
1 /*
2 * Copyright 2010 INRIA Saclay
4 * Use of this software is governed by the GNU LGPLv2.1 license
6 * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
7 * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
8 * 91893 Orsay, France
9 */
11 #include <isl_map_private.h>
12 #include <isl/set.h>
13 #include <isl_dim_private.h>
14 #include <isl/seq.h>
17 * Let C be a cone and define
19 * C' := { y | forall x in C : y x >= 0 }
21 * C' contains the coefficients of all linear constraints
22 * that are valid for C.
23 * Furthermore, C'' = C.
25 * If C is defined as { x | A x >= 0 }
26 * then any element in C' must be a non-negative combination
27 * of the rows of A, i.e., y = t A with t >= 0. That is,
29 * C' = { y | exists t >= 0 : y = t A }
31 * If any of the rows in A actually represents an equality, then
32 * also negative combinations of this row are allowed and so the
33 * non-negativity constraint on the corresponding element of t
34 * can be dropped.
36 * A polyhedron P = { x | b + A x >= 0 } can be represented
37 * in homogeneous coordinates by the cone
38 * C = { [z,x] | b z + A x >= and z >= 0 }
39 * The valid linear constraints on C correspond to the valid affine
40 * constraints on P.
41 * This is essentially Farkas' lemma.
43 * Let A' = [b A], then, since
44 * [ 1 0 ]
45 * [ w y ] = [t_0 t] [ b A ]
47 * we have
49 * C' = { w, y | exists t_0, t >= 0 : y = t A' and w = t_0 + t b }
50 * or
52 * C' = { w, y | exists t >= 0 : y = t A' and w - t b >= 0 }
54 * In practice, we introduce an extra variable (w), shifting all
55 * other variables to the right, and an extra inequality
56 * (w - t b >= 0) corresponding to the positivity constraint on
57 * the homogeneous coordinate.
59 * When going back from coefficients to solutions, we immediately
60 * plug in 1 for z, which corresponds to shifting all variables
61 * to the left, with the leftmost ending up in the constant position.
64 /* Add the given prefix to all named isl_dim_set dimensions in "dim".
66 static __isl_give isl_dim *isl_dim_prefix(__isl_take isl_dim *dim,
67 const char *prefix)
69 int i;
70 isl_ctx *ctx;
71 unsigned nvar;
72 size_t prefix_len = strlen(prefix);
74 if (!dim)
75 return NULL;
77 ctx = isl_dim_get_ctx(dim);
78 nvar = isl_dim_size(dim, isl_dim_set);
80 for (i = 0; i < nvar; ++i) {
81 const char *name;
82 char *prefix_name;
84 name = isl_dim_get_name(dim, isl_dim_set, i);
85 if (!name)
86 continue;
88 prefix_name = isl_alloc_array(ctx, char,
89 prefix_len + strlen(name) + 1);
90 if (!prefix_name)
91 goto error;
92 memcpy(prefix_name, prefix, prefix_len);
93 strcpy(prefix_name + prefix_len, name);
95 dim = isl_dim_set_name(dim, isl_dim_set, i, prefix_name);
96 free(prefix_name);
99 return dim;
100 error:
101 isl_dim_free(dim);
102 return NULL;
105 /* Given a dimension specification of the solutions space, construct
106 * a dimension specification for the space of coefficients.
108 * In particular transform
110 * [params] -> { S }
112 * to
114 * { coefficients[[cst, params] -> S] }
116 * and prefix each dimension name with "c_".
118 static __isl_give isl_dim *isl_dim_coefficients(__isl_take isl_dim *dim)
120 isl_dim *dim_param;
121 unsigned nvar;
122 unsigned nparam;
124 nvar = isl_dim_size(dim, isl_dim_set);
125 nparam = isl_dim_size(dim, isl_dim_param);
126 dim_param = isl_dim_copy(dim);
127 dim_param = isl_dim_drop(dim_param, isl_dim_set, 0, nvar);
128 dim_param = isl_dim_move(dim_param, isl_dim_set, 0,
129 isl_dim_param, 0, nparam);
130 dim_param = isl_dim_prefix(dim_param, "c_");
131 dim_param = isl_dim_insert(dim_param, isl_dim_set, 0, 1);
132 dim_param = isl_dim_set_name(dim_param, isl_dim_set, 0, "c_cst");
133 dim = isl_dim_drop(dim, isl_dim_param, 0, nparam);
134 dim = isl_dim_prefix(dim, "c_");
135 dim = isl_dim_join(isl_dim_from_domain(dim_param),
136 isl_dim_from_range(dim));
137 dim = isl_dim_wrap(dim);
138 dim = isl_dim_set_tuple_name(dim, isl_dim_set, "coefficients");
140 return dim;
143 /* Drop the given prefix from all named dimensions of type "type" in "dim".
145 static __isl_give isl_dim *isl_dim_unprefix(__isl_take isl_dim *dim,
146 enum isl_dim_type type, const char *prefix)
148 int i;
149 unsigned n;
150 size_t prefix_len = strlen(prefix);
152 n = isl_dim_size(dim, type);
154 for (i = 0; i < n; ++i) {
155 const char *name;
157 name = isl_dim_get_name(dim, type, i);
158 if (!name)
159 continue;
160 if (strncmp(name, prefix, prefix_len))
161 continue;
163 dim = isl_dim_set_name(dim, type, i, name + prefix_len);
166 return dim;
169 /* Given a dimension specification of the space of coefficients, construct
170 * a dimension specification for the space of solutions.
172 * In particular transform
174 * { coefficients[[cst, params] -> S] }
176 * to
178 * [params] -> { S }
180 * and drop the "c_" prefix from the dimension names.
182 static __isl_give isl_dim *isl_dim_solutions(__isl_take isl_dim *dim)
184 unsigned nparam;
186 dim = isl_dim_unwrap(dim);
187 dim = isl_dim_drop(dim, isl_dim_in, 0, 1);
188 dim = isl_dim_unprefix(dim, isl_dim_in, "c_");
189 dim = isl_dim_unprefix(dim, isl_dim_out, "c_");
190 nparam = isl_dim_size(dim, isl_dim_in);
191 dim = isl_dim_move(dim, isl_dim_param, 0, isl_dim_in, 0, nparam);
192 dim = isl_dim_range(dim);
194 return dim;
197 /* Compute the dual of "bset" by applying Farkas' lemma.
198 * As explained above, we add an extra dimension to represent
199 * the coefficient of the constant term when going from solutions
200 * to coefficients (shift == 1) and we drop the extra dimension when going
201 * in the opposite direction (shift == -1). "dim" is the space in which
202 * the dual should be created.
204 static __isl_give isl_basic_set *farkas(__isl_take isl_dim *dim,
205 __isl_take isl_basic_set *bset, int shift)
207 int i, j, k;
208 isl_basic_set *dual = NULL;
209 unsigned total;
211 total = isl_basic_set_total_dim(bset);
213 dual = isl_basic_set_alloc_dim(dim, bset->n_eq + bset->n_ineq,
214 total, bset->n_ineq + (shift > 0));
215 dual = isl_basic_set_set_rational(dual);
217 for (i = 0; i < bset->n_eq + bset->n_ineq; ++i) {
218 k = isl_basic_set_alloc_div(dual);
219 if (k < 0)
220 goto error;
221 isl_int_set_si(dual->div[k][0], 0);
224 for (i = 0; i < total; ++i) {
225 k = isl_basic_set_alloc_equality(dual);
226 if (k < 0)
227 goto error;
228 isl_seq_clr(dual->eq[k], 1 + shift + total);
229 isl_int_set_si(dual->eq[k][1 + shift + i], -1);
230 for (j = 0; j < bset->n_eq; ++j)
231 isl_int_set(dual->eq[k][1 + shift + total + j],
232 bset->eq[j][1 + i]);
233 for (j = 0; j < bset->n_ineq; ++j)
234 isl_int_set(dual->eq[k][1 + shift + total + bset->n_eq + j],
235 bset->ineq[j][1 + i]);
238 for (i = 0; i < bset->n_ineq; ++i) {
239 k = isl_basic_set_alloc_inequality(dual);
240 if (k < 0)
241 goto error;
242 isl_seq_clr(dual->ineq[k],
243 1 + shift + total + bset->n_eq + bset->n_ineq);
244 isl_int_set_si(dual->ineq[k][1 + shift + total + bset->n_eq + i], 1);
247 if (shift > 0) {
248 k = isl_basic_set_alloc_inequality(dual);
249 if (k < 0)
250 goto error;
251 isl_seq_clr(dual->ineq[k], 2 + total);
252 isl_int_set_si(dual->ineq[k][1], 1);
253 for (j = 0; j < bset->n_eq; ++j)
254 isl_int_neg(dual->ineq[k][2 + total + j],
255 bset->eq[j][0]);
256 for (j = 0; j < bset->n_ineq; ++j)
257 isl_int_neg(dual->ineq[k][2 + total + bset->n_eq + j],
258 bset->ineq[j][0]);
261 dual = isl_basic_set_remove_divs(dual);
262 isl_basic_set_simplify(dual);
263 isl_basic_set_finalize(dual);
265 isl_basic_set_free(bset);
266 return dual;
267 error:
268 isl_basic_set_free(bset);
269 isl_basic_set_free(dual);
270 return NULL;
273 /* Construct a basic set containing the tuples of coefficients of all
274 * valid affine constraints on the given basic set.
276 __isl_give isl_basic_set *isl_basic_set_coefficients(
277 __isl_take isl_basic_set *bset)
279 isl_dim *dim;
281 if (!bset)
282 return NULL;
283 if (bset->n_div)
284 isl_die(bset->ctx, isl_error_invalid,
285 "input set not allowed to have local variables",
286 goto error);
288 dim = isl_basic_set_get_dim(bset);
289 dim = isl_dim_coefficients(dim);
291 return farkas(dim, bset, 1);
292 error:
293 isl_basic_set_free(bset);
294 return NULL;
297 /* Construct a basic set containing the elements that satisfy all
298 * affine constraints whose coefficient tuples are
299 * contained in the given basic set.
301 __isl_give isl_basic_set *isl_basic_set_solutions(
302 __isl_take isl_basic_set *bset)
304 isl_dim *dim;
306 if (!bset)
307 return NULL;
308 if (bset->n_div)
309 isl_die(bset->ctx, isl_error_invalid,
310 "input set not allowed to have local variables",
311 goto error);
313 dim = isl_basic_set_get_dim(bset);
314 dim = isl_dim_solutions(dim);
316 return farkas(dim, bset, -1);
317 error:
318 isl_basic_set_free(bset);
319 return NULL;
322 /* Construct a basic set containing the tuples of coefficients of all
323 * valid affine constraints on the given set.
325 __isl_give isl_basic_set *isl_set_coefficients(__isl_take isl_set *set)
327 int i;
328 isl_basic_set *coeff;
330 if (!set)
331 return NULL;
332 if (set->n == 0) {
333 isl_dim *dim = isl_set_get_dim(set);
334 dim = isl_dim_coefficients(dim);
335 coeff = isl_basic_set_universe(dim);
336 coeff = isl_basic_set_set_rational(coeff);
337 isl_set_free(set);
338 return coeff;
341 coeff = isl_basic_set_coefficients(isl_basic_set_copy(set->p[0]));
343 for (i = 1; i < set->n; ++i) {
344 isl_basic_set *bset, *coeff_i;
345 bset = isl_basic_set_copy(set->p[i]);
346 coeff_i = isl_basic_set_coefficients(bset);
347 coeff = isl_basic_set_intersect(coeff, coeff_i);
350 isl_set_free(set);
351 return coeff;
354 /* Construct a basic set containing the elements that satisfy all
355 * affine constraints whose coefficient tuples are
356 * contained in the given set.
358 __isl_give isl_basic_set *isl_set_solutions(__isl_take isl_set *set)
360 int i;
361 isl_basic_set *sol;
363 if (!set)
364 return NULL;
365 if (set->n == 0) {
366 isl_dim *dim = isl_set_get_dim(set);
367 dim = isl_dim_solutions(dim);
368 sol = isl_basic_set_universe(dim);
369 sol = isl_basic_set_set_rational(sol);
370 isl_set_free(set);
371 return sol;
374 sol = isl_basic_set_solutions(isl_basic_set_copy(set->p[0]));
376 for (i = 1; i < set->n; ++i) {
377 isl_basic_set *bset, *sol_i;
378 bset = isl_basic_set_copy(set->p[i]);
379 sol_i = isl_basic_set_solutions(bset);
380 sol = isl_basic_set_intersect(sol, sol_i);
383 isl_set_free(set);
384 return sol;