1/*
2 * Copyright 2010 INRIA Saclay
3 *
4 * Use of this software is governed by the MIT license
5 *
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 */
10
11#include <isl_map_private.h>
12#include <isl/set.h>
13#include <isl_space_private.h>
14#include <isl_seq.h>
15#include <isl_aff_private.h>
16#include <isl_mat_private.h>
17#include <isl_factorization.h>
18
19/*
20 * Let C be a cone and define
21 *
22 * C' := { y | forall x in C : y x >= 0 }
23 *
24 * C' contains the coefficients of all linear constraints
25 * that are valid for C.
26 * Furthermore, C'' = C.
27 *
28 * If C is defined as { x | A x >= 0 }
29 * then any element in C' must be a non-negative combination
30 * of the rows of A, i.e., y = t A with t >= 0. That is,
31 *
32 * C' = { y | exists t >= 0 : y = t A }
33 *
34 * If any of the rows in A actually represents an equality, then
35 * also negative combinations of this row are allowed and so the
36 * non-negativity constraint on the corresponding element of t
37 * can be dropped.
38 *
39 * A polyhedron P = { x | b + A x >= 0 } can be represented
40 * in homogeneous coordinates by the cone
41 * C = { [z,x] | b z + A x >= and z >= 0 }
42 * The valid linear constraints on C correspond to the valid affine
43 * constraints on P.
44 * This is essentially Farkas' lemma.
45 *
46 * Since
47 * [ 1 0 ]
48 * [ w y ] = [t_0 t] [ b A ]
49 *
50 * we have
51 *
52 * C' = { w, y | exists t_0, t >= 0 : y = t A and w = t_0 + t b }
53 * or
54 *
55 * C' = { w, y | exists t >= 0 : y = t A and w - t b >= 0 }
56 *
57 * In practice, we introduce an extra variable (w), shifting all
58 * other variables to the right, and an extra inequality
59 * (w - t b >= 0) corresponding to the positivity constraint on
60 * the homogeneous coordinate.
61 *
62 * When going back from coefficients to solutions, we immediately
63 * plug in 1 for z, which corresponds to shifting all variables
64 * to the left, with the leftmost ending up in the constant position.
65 */
66
67/* Add the given prefix to all named isl_dim_set dimensions in "space".
68 */
69static __isl_give isl_space *isl_space_prefix(__isl_take isl_space *space,
70 const char *prefix)
71{
72 int i;
73 isl_ctx *ctx;
74 isl_size nvar;
75 size_t prefix_len = strlen(s: prefix);
76
77 if (!space)
78 return NULL;
79
80 ctx = isl_space_get_ctx(space);
81 nvar = isl_space_dim(space, type: isl_dim_set);
82 if (nvar < 0)
83 return isl_space_free(space);
84
85 for (i = 0; i < nvar; ++i) {
86 const char *name;
87 char *prefix_name;
88
89 name = isl_space_get_dim_name(space, type: isl_dim_set, pos: i);
90 if (!name)
91 continue;
92
93 prefix_name = isl_alloc_array(ctx, char,
94 prefix_len + strlen(name) + 1);
95 if (!prefix_name)
96 goto error;
97 memcpy(dest: prefix_name, src: prefix, n: prefix_len);
98 strcpy(dest: prefix_name + prefix_len, src: name);
99
100 space = isl_space_set_dim_name(space,
101 type: isl_dim_set, pos: i, name: prefix_name);
102 free(ptr: prefix_name);
103 }
104
105 return space;
106error:
107 isl_space_free(space);
108 return NULL;
109}
110
111/* Given a dimension specification of the solutions space, construct
112 * a dimension specification for the space of coefficients.
113 *
114 * In particular transform
115 *
116 * [params] -> { S }
117 *
118 * to
119 *
120 * { coefficients[[cst, params] -> S] }
121 *
122 * and prefix each dimension name with "c_".
123 */
124static __isl_give isl_space *isl_space_coefficients(__isl_take isl_space *space)
125{
126 isl_space *space_param;
127 isl_size nvar;
128 isl_size nparam;
129
130 nvar = isl_space_dim(space, type: isl_dim_set);
131 nparam = isl_space_dim(space, type: isl_dim_param);
132 if (nvar < 0 || nparam < 0)
133 return isl_space_free(space);
134 space_param = isl_space_copy(space);
135 space_param = isl_space_drop_dims(space: space_param, type: isl_dim_set, first: 0, num: nvar);
136 space_param = isl_space_move_dims(space: space_param, dst_type: isl_dim_set, dst_pos: 0,
137 src_type: isl_dim_param, src_pos: 0, n: nparam);
138 space_param = isl_space_prefix(space: space_param, prefix: "c_");
139 space_param = isl_space_insert_dims(space: space_param, type: isl_dim_set, pos: 0, n: 1);
140 space_param = isl_space_set_dim_name(space: space_param,
141 type: isl_dim_set, pos: 0, name: "c_cst");
142 space = isl_space_drop_dims(space, type: isl_dim_param, first: 0, num: nparam);
143 space = isl_space_prefix(space, prefix: "c_");
144 space = isl_space_join(left: isl_space_from_domain(space: space_param),
145 right: isl_space_from_range(space));
146 space = isl_space_wrap(space);
147 space = isl_space_set_tuple_name(space, type: isl_dim_set, s: "coefficients");
148
149 return space;
150}
151
152/* Drop the given prefix from all named dimensions of type "type" in "space".
153 */
154static __isl_give isl_space *isl_space_unprefix(__isl_take isl_space *space,
155 enum isl_dim_type type, const char *prefix)
156{
157 int i;
158 isl_size n;
159 size_t prefix_len = strlen(s: prefix);
160
161 n = isl_space_dim(space, type);
162 if (n < 0)
163 return isl_space_free(space);
164
165 for (i = 0; i < n; ++i) {
166 const char *name;
167
168 name = isl_space_get_dim_name(space, type, pos: i);
169 if (!name)
170 continue;
171 if (strncmp(s1: name, s2: prefix, n: prefix_len))
172 continue;
173
174 space = isl_space_set_dim_name(space,
175 type, pos: i, name: name + prefix_len);
176 }
177
178 return space;
179}
180
181/* Given a dimension specification of the space of coefficients, construct
182 * a dimension specification for the space of solutions.
183 *
184 * In particular transform
185 *
186 * { coefficients[[cst, params] -> S] }
187 *
188 * to
189 *
190 * [params] -> { S }
191 *
192 * and drop the "c_" prefix from the dimension names.
193 */
194static __isl_give isl_space *isl_space_solutions(__isl_take isl_space *space)
195{
196 isl_size nparam;
197
198 space = isl_space_unwrap(space);
199 space = isl_space_drop_dims(space, type: isl_dim_in, first: 0, num: 1);
200 space = isl_space_unprefix(space, type: isl_dim_in, prefix: "c_");
201 space = isl_space_unprefix(space, type: isl_dim_out, prefix: "c_");
202 nparam = isl_space_dim(space, type: isl_dim_in);
203 if (nparam < 0)
204 return isl_space_free(space);
205 space = isl_space_move_dims(space,
206 dst_type: isl_dim_param, dst_pos: 0, src_type: isl_dim_in, src_pos: 0, n: nparam);
207 space = isl_space_range(space);
208
209 return space;
210}
211
212/* Return the rational universe basic set in the given space.
213 */
214static __isl_give isl_basic_set *rational_universe(__isl_take isl_space *space)
215{
216 isl_basic_set *bset;
217
218 bset = isl_basic_set_universe(space);
219 bset = isl_basic_set_set_rational(bset);
220
221 return bset;
222}
223
224/* Compute the dual of "bset" by applying Farkas' lemma.
225 * As explained above, we add an extra dimension to represent
226 * the coefficient of the constant term when going from solutions
227 * to coefficients (shift == 1) and we drop the extra dimension when going
228 * in the opposite direction (shift == -1).
229 * The dual can be created in an arbitrary space.
230 * The caller is responsible for putting the result in the appropriate space.
231 *
232 * If "bset" is (obviously) empty, then the way this emptiness
233 * is represented by the constraints does not allow for the application
234 * of the standard farkas algorithm. We therefore handle this case
235 * specifically and return the universe basic set.
236 */
237static __isl_give isl_basic_set *farkas(__isl_take isl_basic_set *bset,
238 int shift)
239{
240 int i, j, k;
241 isl_ctx *ctx;
242 isl_space *space;
243 isl_basic_set *dual = NULL;
244 isl_size total;
245
246 total = isl_basic_set_dim(bset, type: isl_dim_all);
247 if (total < 0)
248 return isl_basic_set_free(bset);
249
250 ctx = isl_basic_set_get_ctx(bset);
251 space = isl_space_set_alloc(ctx, nparam: 0, dim: total + shift);
252 if (isl_basic_set_plain_is_empty(bset)) {
253 isl_basic_set_free(bset);
254 return rational_universe(space);
255 }
256
257 dual = isl_basic_set_alloc_space(space, extra: bset->n_eq + bset->n_ineq,
258 n_eq: total, n_ineq: bset->n_ineq + (shift > 0));
259 dual = isl_basic_set_set_rational(bset: dual);
260
261 for (i = 0; i < bset->n_eq + bset->n_ineq; ++i) {
262 k = isl_basic_set_alloc_div(bset: dual);
263 if (k < 0)
264 goto error;
265 isl_int_set_si(dual->div[k][0], 0);
266 }
267
268 for (i = 0; i < total; ++i) {
269 k = isl_basic_set_alloc_equality(bset: dual);
270 if (k < 0)
271 goto error;
272 isl_seq_clr(p: dual->eq[k], len: 1 + shift + total);
273 isl_int_set_si(dual->eq[k][1 + shift + i], -1);
274 for (j = 0; j < bset->n_eq; ++j)
275 isl_int_set(dual->eq[k][1 + shift + total + j],
276 bset->eq[j][1 + i]);
277 for (j = 0; j < bset->n_ineq; ++j)
278 isl_int_set(dual->eq[k][1 + shift + total + bset->n_eq + j],
279 bset->ineq[j][1 + i]);
280 }
281
282 for (i = 0; i < bset->n_ineq; ++i) {
283 k = isl_basic_set_alloc_inequality(bset: dual);
284 if (k < 0)
285 goto error;
286 isl_seq_clr(p: dual->ineq[k],
287 len: 1 + shift + total + bset->n_eq + bset->n_ineq);
288 isl_int_set_si(dual->ineq[k][1 + shift + total + bset->n_eq + i], 1);
289 }
290
291 if (shift > 0) {
292 k = isl_basic_set_alloc_inequality(bset: dual);
293 if (k < 0)
294 goto error;
295 isl_seq_clr(p: dual->ineq[k], len: 2 + total);
296 isl_int_set_si(dual->ineq[k][1], 1);
297 for (j = 0; j < bset->n_eq; ++j)
298 isl_int_neg(dual->ineq[k][2 + total + j],
299 bset->eq[j][0]);
300 for (j = 0; j < bset->n_ineq; ++j)
301 isl_int_neg(dual->ineq[k][2 + total + bset->n_eq + j],
302 bset->ineq[j][0]);
303 }
304
305 dual = isl_basic_set_remove_divs(bset: dual);
306 dual = isl_basic_set_simplify(bset: dual);
307 dual = isl_basic_set_finalize(bset: dual);
308
309 isl_basic_set_free(bset);
310 return dual;
311error:
312 isl_basic_set_free(bset);
313 isl_basic_set_free(bset: dual);
314 return NULL;
315}
316
317/* Construct a basic set containing the tuples of coefficients of all
318 * valid affine constraints on the given basic set, ignoring
319 * the space of input and output and without any further decomposition.
320 */
321static __isl_give isl_basic_set *isl_basic_set_coefficients_base(
322 __isl_take isl_basic_set *bset)
323{
324 return farkas(bset, shift: 1);
325}
326
327/* Return the inverse mapping of "morph".
328 */
329static __isl_give isl_mat *peek_inv(__isl_keep isl_morph *morph)
330{
331 return morph ? morph->inv : NULL;
332}
333
334/* Return a copy of the inverse mapping of "morph".
335 */
336static __isl_give isl_mat *get_inv(__isl_keep isl_morph *morph)
337{
338 return isl_mat_copy(mat: peek_inv(morph));
339}
340
341/* Information about a single factor within isl_basic_set_coefficients_product.
342 *
343 * "start" is the position of the first coefficient (beyond
344 * the one corresponding to the constant term) in this factor.
345 * "dim" is the number of coefficients (other than
346 * the one corresponding to the constant term) in this factor.
347 * "n_line" is the number of lines in "coeff".
348 * "n_ray" is the number of rays (other than lines) in "coeff".
349 * "n_vertex" is the number of vertices in "coeff".
350 *
351 * While iterating over the vertices,
352 * "pos" represents the inequality constraint corresponding
353 * to the current vertex.
354 */
355struct isl_coefficients_factor_data {
356 isl_basic_set *coeff;
357 int start;
358 int dim;
359 int n_line;
360 int n_ray;
361 int n_vertex;
362 int pos;
363};
364
365/* Internal data structure for isl_basic_set_coefficients_product.
366 * "n" is the number of factors in the factorization.
367 * "pos" is the next factor that will be considered.
368 * "start_next" is the position of the first coefficient (beyond
369 * the one corresponding to the constant term) in the next factor.
370 * "factors" contains information about the individual "n" factors.
371 */
372struct isl_coefficients_product_data {
373 int n;
374 int pos;
375 int start_next;
376 struct isl_coefficients_factor_data *factors;
377};
378
379/* Initialize the internal data structure for
380 * isl_basic_set_coefficients_product.
381 */
382static isl_stat isl_coefficients_product_data_init(isl_ctx *ctx,
383 struct isl_coefficients_product_data *data, int n)
384{
385 data->n = n;
386 data->pos = 0;
387 data->start_next = 0;
388 data->factors = isl_calloc_array(ctx,
389 struct isl_coefficients_factor_data, n);
390 if (!data->factors)
391 return isl_stat_error;
392 return isl_stat_ok;
393}
394
395/* Free all memory allocated in "data".
396 */
397static void isl_coefficients_product_data_clear(
398 struct isl_coefficients_product_data *data)
399{
400 int i;
401
402 if (data->factors) {
403 for (i = 0; i < data->n; ++i) {
404 isl_basic_set_free(bset: data->factors[i].coeff);
405 }
406 }
407 free(ptr: data->factors);
408}
409
410/* Does inequality "ineq" in the (dual) basic set "bset" represent a ray?
411 * In particular, does it have a zero denominator
412 * (i.e., a zero coefficient for the constant term)?
413 */
414static int is_ray(__isl_keep isl_basic_set *bset, int ineq)
415{
416 return isl_int_is_zero(bset->ineq[ineq][1]);
417}
418
419/* isl_factorizer_every_factor_basic_set callback that
420 * constructs a basic set containing the tuples of coefficients of all
421 * valid affine constraints on the factor "bset" and
422 * extracts further information that will be used
423 * when combining the results over the different factors.
424 */
425static isl_bool isl_basic_set_coefficients_factor(
426 __isl_keep isl_basic_set *bset, void *user)
427{
428 struct isl_coefficients_product_data *data = user;
429 isl_basic_set *coeff;
430 isl_size n_eq, n_ineq, dim;
431 int i, n_ray, n_vertex;
432
433 coeff = isl_basic_set_coefficients_base(bset: isl_basic_set_copy(bset));
434 data->factors[data->pos].coeff = coeff;
435 if (!coeff)
436 return isl_bool_error;
437
438 dim = isl_basic_set_dim(bset, type: isl_dim_set);
439 n_eq = isl_basic_set_n_equality(bset: coeff);
440 n_ineq = isl_basic_set_n_inequality(bset: coeff);
441 if (dim < 0 || n_eq < 0 || n_ineq < 0)
442 return isl_bool_error;
443 n_ray = n_vertex = 0;
444 for (i = 0; i < n_ineq; ++i) {
445 if (is_ray(bset: coeff, ineq: i))
446 n_ray++;
447 else
448 n_vertex++;
449 }
450 data->factors[data->pos].start = data->start_next;
451 data->factors[data->pos].dim = dim;
452 data->factors[data->pos].n_line = n_eq;
453 data->factors[data->pos].n_ray = n_ray;
454 data->factors[data->pos].n_vertex = n_vertex;
455 data->pos++;
456 data->start_next += dim;
457
458 return isl_bool_true;
459}
460
461/* Clear an entry in the product, given that there is a "total" number
462 * of coefficients (other than that of the constant term).
463 */
464static void clear_entry(isl_int *entry, int total)
465{
466 isl_seq_clr(p: entry, len: 1 + 1 + total);
467}
468
469/* Set the part of the entry corresponding to factor "data",
470 * from the factor coefficients in "src".
471 */
472static void set_factor(isl_int *entry, isl_int *src,
473 struct isl_coefficients_factor_data *data)
474{
475 isl_seq_cpy(dst: entry + 1 + 1 + data->start, src: src + 1 + 1, len: data->dim);
476}
477
478/* Set the part of the entry corresponding to factor "data",
479 * from the factor coefficients in "src" multiplied by "f".
480 */
481static void scale_factor(isl_int *entry, isl_int *src, isl_int f,
482 struct isl_coefficients_factor_data *data)
483{
484 isl_seq_scale(dst: entry + 1 + 1 + data->start, src: src + 1 + 1, f, len: data->dim);
485}
486
487/* Add all lines from the given factor to "bset",
488 * given that there is a "total" number of coefficients
489 * (other than that of the constant term).
490 */
491static __isl_give isl_basic_set *add_lines(__isl_take isl_basic_set *bset,
492 struct isl_coefficients_factor_data *factor, int total)
493{
494 int i;
495
496 for (i = 0; i < factor->n_line; ++i) {
497 int k;
498
499 k = isl_basic_set_alloc_equality(bset);
500 if (k < 0)
501 return isl_basic_set_free(bset);
502 clear_entry(entry: bset->eq[k], total);
503 set_factor(entry: bset->eq[k], src: factor->coeff->eq[i], data: factor);
504 }
505
506 return bset;
507}
508
509/* Add all rays (other than lines) from the given factor to "bset",
510 * given that there is a "total" number of coefficients
511 * (other than that of the constant term).
512 */
513static __isl_give isl_basic_set *add_rays(__isl_take isl_basic_set *bset,
514 struct isl_coefficients_factor_data *data, int total)
515{
516 int i;
517 int n_ineq = data->n_ray + data->n_vertex;
518
519 for (i = 0; i < n_ineq; ++i) {
520 int k;
521
522 if (!is_ray(bset: data->coeff, ineq: i))
523 continue;
524
525 k = isl_basic_set_alloc_inequality(bset);
526 if (k < 0)
527 return isl_basic_set_free(bset);
528 clear_entry(entry: bset->ineq[k], total);
529 set_factor(entry: bset->ineq[k], src: data->coeff->ineq[i], data);
530 }
531
532 return bset;
533}
534
535/* Move to the first vertex of the given factor starting
536 * at inequality constraint "start", setting factor->pos and
537 * returning 1 if a vertex is found.
538 */
539static int factor_first_vertex(struct isl_coefficients_factor_data *factor,
540 int start)
541{
542 int j;
543 int n = factor->n_ray + factor->n_vertex;
544
545 for (j = start; j < n; ++j) {
546 if (is_ray(bset: factor->coeff, ineq: j))
547 continue;
548 factor->pos = j;
549 return 1;
550 }
551
552 return 0;
553}
554
555/* Move to the first constraint in each factor starting at "first"
556 * that represents a vertex.
557 * In particular, skip the initial constraints that correspond to rays.
558 */
559static void first_vertex(struct isl_coefficients_product_data *data, int first)
560{
561 int i;
562
563 for (i = first; i < data->n; ++i)
564 factor_first_vertex(factor: &data->factors[i], start: 0);
565}
566
567/* Move to the next vertex in the product.
568 * In particular, move to the next vertex of the last factor.
569 * If all vertices of this last factor have already been considered,
570 * then move to the next vertex of the previous factor(s)
571 * until a factor is found that still has a next vertex.
572 * Once such a next vertex has been found, the subsequent
573 * factors are reset to the first vertex.
574 * Return 1 if any next vertex was found.
575 */
576static int next_vertex(struct isl_coefficients_product_data *data)
577{
578 int i;
579
580 for (i = data->n - 1; i >= 0; --i) {
581 struct isl_coefficients_factor_data *factor = &data->factors[i];
582
583 if (!factor_first_vertex(factor, start: factor->pos + 1))
584 continue;
585 first_vertex(data, first: i + 1);
586 return 1;
587 }
588
589 return 0;
590}
591
592/* Add a vertex to the product "bset" combining the currently selected
593 * vertices of the factors.
594 *
595 * In the dual representation, the constant term is always zero.
596 * The vertex itself is the sum of the contributions of the factors
597 * with a shared denominator in position 1.
598 *
599 * First compute the shared denominator (lcm) and
600 * then scale the numerators to this shared denominator.
601 */
602static __isl_give isl_basic_set *add_vertex(__isl_take isl_basic_set *bset,
603 struct isl_coefficients_product_data *data)
604{
605 int i;
606 int k;
607 isl_int lcm, f;
608
609 k = isl_basic_set_alloc_inequality(bset);
610 if (k < 0)
611 return isl_basic_set_free(bset);
612
613 isl_int_init(lcm);
614 isl_int_init(f);
615 isl_int_set_si(lcm, 1);
616 for (i = 0; i < data->n; ++i) {
617 struct isl_coefficients_factor_data *factor = &data->factors[i];
618 isl_basic_set *coeff = factor->coeff;
619 int pos = factor->pos;
620 isl_int_lcm(lcm, lcm, coeff->ineq[pos][1]);
621 }
622 isl_int_set_si(bset->ineq[k][0], 0);
623 isl_int_set(bset->ineq[k][1], lcm);
624
625 for (i = 0; i < data->n; ++i) {
626 struct isl_coefficients_factor_data *factor = &data->factors[i];
627 isl_basic_set *coeff = factor->coeff;
628 int pos = factor->pos;
629 isl_int_divexact(f, lcm, coeff->ineq[pos][1]);
630 scale_factor(entry: bset->ineq[k], src: coeff->ineq[pos], f, data: factor);
631 }
632
633 isl_int_clear(f);
634 isl_int_clear(lcm);
635
636 return bset;
637}
638
639/* Combine the duals of the factors in the factorization of a basic set
640 * to form the dual of the entire basic set.
641 * The dual share the coefficient of the constant term.
642 * All other coefficients are specific to a factor.
643 * Any constraint not involving the coefficient of the constant term
644 * can therefor simply be copied into the appropriate position.
645 * This includes all equality constraints since the coefficient
646 * of the constant term can always be increased and therefore
647 * never appears in an equality constraint.
648 * The inequality constraints involving the coefficient of
649 * the constant term need to be combined across factors.
650 * In particular, if this coefficient needs to be greater than or equal
651 * to some linear combination of the other coefficients in each factor,
652 * then it needs to be greater than or equal to the sum of
653 * these linear combinations across the factors.
654 *
655 * Alternatively, the constraints of the dual can be seen
656 * as the vertices, rays and lines of the original basic set.
657 * Clearly, rays and lines can simply be copied,
658 * while vertices needs to be combined across factors.
659 * This means that the number of rays and lines in the product
660 * is equal to the sum of the numbers in the factors,
661 * while the number of vertices is the product
662 * of the number of vertices in the factors. Note that each
663 * factor has at least one vertex.
664 * The only exception is when the factor is the dual of an obviously empty set,
665 * in which case a universe dual is created.
666 * In this case, return a universe dual for the product as well.
667 *
668 * While constructing the vertices, look for the first combination
669 * of inequality constraints that represent a vertex,
670 * construct the corresponding vertex and then move on
671 * to the next combination of inequality constraints until
672 * all combinations have been considered.
673 */
674static __isl_give isl_basic_set *construct_product(isl_ctx *ctx,
675 struct isl_coefficients_product_data *data)
676{
677 int i;
678 int n_line, n_ray, n_vertex;
679 int total;
680 isl_space *space;
681 isl_basic_set *product;
682
683 if (!data->factors)
684 return NULL;
685
686 total = data->start_next;
687
688 n_line = 0;
689 n_ray = 0;
690 n_vertex = 1;
691 for (i = 0; i < data->n; ++i) {
692 n_line += data->factors[i].n_line;
693 n_ray += data->factors[i].n_ray;
694 n_vertex *= data->factors[i].n_vertex;
695 }
696
697 space = isl_space_set_alloc(ctx, nparam: 0, dim: 1 + total);
698 if (n_vertex == 0)
699 return rational_universe(space);
700 product = isl_basic_set_alloc_space(space, extra: 0, n_eq: n_line, n_ineq: n_ray + n_vertex);
701 product = isl_basic_set_set_rational(bset: product);
702
703 for (i = 0; i < data->n; ++i)
704 product = add_lines(bset: product, factor: &data->factors[i], total);
705 for (i = 0; i < data->n; ++i)
706 product = add_rays(bset: product, data: &data->factors[i], total);
707
708 first_vertex(data, first: 0);
709 do {
710 product = add_vertex(bset: product, data);
711 } while (next_vertex(data));
712
713 return product;
714}
715
716/* Given a factorization "f" of a basic set,
717 * construct a basic set containing the tuples of coefficients of all
718 * valid affine constraints on the product of the factors, ignoring
719 * the space of input and output.
720 * Note that this product may not be equal to the original basic set,
721 * if a non-trivial transformation is involved.
722 * This is handled by the caller.
723 *
724 * Compute the tuples of coefficients for each factor separately and
725 * then combine the results.
726 */
727static __isl_give isl_basic_set *isl_basic_set_coefficients_product(
728 __isl_take isl_factorizer *f)
729{
730 struct isl_coefficients_product_data data;
731 isl_ctx *ctx;
732 isl_basic_set *coeff;
733 isl_bool every;
734
735 ctx = isl_factorizer_get_ctx(f);
736 if (isl_coefficients_product_data_init(ctx, data: &data, n: f->n_group) < 0)
737 f = isl_factorizer_free(f);
738 every = isl_factorizer_every_factor_basic_set(f,
739 test: &isl_basic_set_coefficients_factor, user: &data);
740 isl_factorizer_free(f);
741 if (every >= 0)
742 coeff = construct_product(ctx, data: &data);
743 else
744 coeff = NULL;
745 isl_coefficients_product_data_clear(data: &data);
746
747 return coeff;
748}
749
750/* Given a factorization "f" of a basic set,
751 * construct a basic set containing the tuples of coefficients of all
752 * valid affine constraints on the basic set, ignoring
753 * the space of input and output.
754 *
755 * The factorization may involve a linear transformation of the basic set.
756 * In particular, the transformed basic set is formulated
757 * in terms of x' = U x, i.e., x = V x', with V = U^{-1}.
758 * The dual is then computed in terms of y' with y'^t [z; x'] >= 0.
759 * Plugging in y' = [1 0; 0 V^t] y yields
760 * y^t [1 0; 0 V] [z; x'] >= 0, i.e., y^t [z; x] >= 0, which is
761 * the desired set of coefficients y.
762 * Note that this transformation to y' only needs to be applied
763 * if U is not the identity matrix.
764 */
765static __isl_give isl_basic_set *isl_basic_set_coefficients_morphed_product(
766 __isl_take isl_factorizer *f)
767{
768 isl_bool is_identity;
769 isl_space *space;
770 isl_mat *inv;
771 isl_multi_aff *ma;
772 isl_basic_set *coeff;
773
774 if (!f)
775 goto error;
776 is_identity = isl_mat_is_scaled_identity(mat: peek_inv(morph: f->morph));
777 if (is_identity < 0)
778 goto error;
779 if (is_identity)
780 return isl_basic_set_coefficients_product(f);
781
782 inv = get_inv(morph: f->morph);
783 inv = isl_mat_transpose(mat: inv);
784 inv = isl_mat_lin_to_aff(mat: inv);
785
786 coeff = isl_basic_set_coefficients_product(f);
787 space = isl_space_map_from_set(space: isl_basic_set_get_space(bset: coeff));
788 ma = isl_multi_aff_from_aff_mat(space, mat: inv);
789 coeff = isl_basic_set_preimage_multi_aff(bset: coeff, ma);
790
791 return coeff;
792error:
793 isl_factorizer_free(f);
794 return NULL;
795}
796
797/* Construct a basic set containing the tuples of coefficients of all
798 * valid affine constraints on the given basic set, ignoring
799 * the space of input and output.
800 *
801 * The caller has already checked that "bset" does not involve
802 * any local variables. It may have parameters, though.
803 * Treat them as regular variables internally.
804 * This is especially important for the factorization,
805 * since the (original) parameters should be taken into account
806 * explicitly in this factorization.
807 *
808 * Check if the basic set can be factorized.
809 * If so, compute constraints on the coefficients of the factors
810 * separately and combine the results.
811 * Otherwise, compute the results for the input basic set as a whole.
812 */
813static __isl_give isl_basic_set *basic_set_coefficients(
814 __isl_take isl_basic_set *bset)
815{
816 isl_factorizer *f;
817 isl_size nparam;
818
819 nparam = isl_basic_set_dim(bset, type: isl_dim_param);
820 if (nparam < 0)
821 return isl_basic_set_free(bset);
822 bset = isl_basic_set_move_dims(bset, dst_type: isl_dim_set, dst_pos: 0,
823 src_type: isl_dim_param, src_pos: 0, n: nparam);
824
825 f = isl_basic_set_factorizer(bset);
826 if (!f)
827 return isl_basic_set_free(bset);
828 if (f->n_group > 0) {
829 isl_basic_set_free(bset);
830 return isl_basic_set_coefficients_morphed_product(f);
831 }
832 isl_factorizer_free(f);
833 return isl_basic_set_coefficients_base(bset);
834}
835
836/* Construct a basic set containing the tuples of coefficients of all
837 * valid affine constraints on the given basic set.
838 */
839__isl_give isl_basic_set *isl_basic_set_coefficients(
840 __isl_take isl_basic_set *bset)
841{
842 isl_space *space;
843
844 if (!bset)
845 return NULL;
846 if (bset->n_div)
847 isl_die(bset->ctx, isl_error_invalid,
848 "input set not allowed to have local variables",
849 goto error);
850
851 space = isl_basic_set_get_space(bset);
852 space = isl_space_coefficients(space);
853
854 bset = basic_set_coefficients(bset);
855 bset = isl_basic_set_reset_space(bset, space);
856 return bset;
857error:
858 isl_basic_set_free(bset);
859 return NULL;
860}
861
862/* Construct a basic set containing the elements that satisfy all
863 * affine constraints whose coefficient tuples are
864 * contained in the given basic set.
865 */
866__isl_give isl_basic_set *isl_basic_set_solutions(
867 __isl_take isl_basic_set *bset)
868{
869 isl_space *space;
870
871 if (!bset)
872 return NULL;
873 if (bset->n_div)
874 isl_die(bset->ctx, isl_error_invalid,
875 "input set not allowed to have local variables",
876 goto error);
877
878 space = isl_basic_set_get_space(bset);
879 space = isl_space_solutions(space);
880
881 bset = farkas(bset, shift: -1);
882 bset = isl_basic_set_reset_space(bset, space);
883 return bset;
884error:
885 isl_basic_set_free(bset);
886 return NULL;
887}
888
889/* Construct a basic set containing the tuples of coefficients of all
890 * valid affine constraints on the given set.
891 */
892__isl_give isl_basic_set *isl_set_coefficients(__isl_take isl_set *set)
893{
894 int i;
895 isl_basic_set *coeff;
896
897 if (!set)
898 return NULL;
899 if (set->n == 0) {
900 isl_space *space = isl_set_get_space(set);
901 space = isl_space_coefficients(space);
902 isl_set_free(set);
903 return rational_universe(space);
904 }
905
906 coeff = isl_basic_set_coefficients(bset: isl_basic_set_copy(bset: set->p[0]));
907
908 for (i = 1; i < set->n; ++i) {
909 isl_basic_set *bset, *coeff_i;
910 bset = isl_basic_set_copy(bset: set->p[i]);
911 coeff_i = isl_basic_set_coefficients(bset);
912 coeff = isl_basic_set_intersect(bset1: coeff, bset2: coeff_i);
913 }
914
915 isl_set_free(set);
916 return coeff;
917}
918
919/* Wrapper around isl_basic_set_coefficients for use
920 * as a isl_basic_set_list_map callback.
921 */
922static __isl_give isl_basic_set *coefficients_wrap(
923 __isl_take isl_basic_set *bset, void *user)
924{
925 return isl_basic_set_coefficients(bset);
926}
927
928/* Replace the elements of "list" by the result of applying
929 * isl_basic_set_coefficients to them.
930 */
931__isl_give isl_basic_set_list *isl_basic_set_list_coefficients(
932 __isl_take isl_basic_set_list *list)
933{
934 return isl_basic_set_list_map(list, fn: &coefficients_wrap, NULL);
935}
936
937/* Construct a basic set containing the elements that satisfy all
938 * affine constraints whose coefficient tuples are
939 * contained in the given set.
940 */
941__isl_give isl_basic_set *isl_set_solutions(__isl_take isl_set *set)
942{
943 int i;
944 isl_basic_set *sol;
945
946 if (!set)
947 return NULL;
948 if (set->n == 0) {
949 isl_space *space = isl_set_get_space(set);
950 space = isl_space_solutions(space);
951 isl_set_free(set);
952 return rational_universe(space);
953 }
954
955 sol = isl_basic_set_solutions(bset: isl_basic_set_copy(bset: set->p[0]));
956
957 for (i = 1; i < set->n; ++i) {
958 isl_basic_set *bset, *sol_i;
959 bset = isl_basic_set_copy(bset: set->p[i]);
960 sol_i = isl_basic_set_solutions(bset);
961 sol = isl_basic_set_intersect(bset1: sol, bset2: sol_i);
962 }
963
964 isl_set_free(set);
965 return sol;
966}
967

source code of polly/lib/External/isl/isl_farkas.c