1 | /* |
2 | * Copyright 2008-2009 Katholieke Universiteit Leuven |
3 | * Copyright 2010 INRIA Saclay |
4 | * Copyright 2011 Sven Verdoolaege |
5 | * |
6 | * Use of this software is governed by the MIT license |
7 | * |
8 | * Written by Sven Verdoolaege, K.U.Leuven, Departement |
9 | * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium |
10 | * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite, |
11 | * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France |
12 | */ |
13 | |
14 | #define xSF(TYPE,SUFFIX) TYPE ## SUFFIX |
15 | #define SF(TYPE,SUFFIX) xSF(TYPE,SUFFIX) |
16 | |
17 | /* Given a basic map with at least two parallel constraints (as found |
18 | * by the function parallel_constraints), first look for more constraints |
19 | * parallel to the two constraint and replace the found list of parallel |
20 | * constraints by a single constraint with as "input" part the minimum |
21 | * of the input parts of the list of constraints. Then, recursively call |
22 | * basic_map_partial_lexopt (possibly finding more parallel constraints) |
23 | * and plug in the definition of the minimum in the result. |
24 | * |
25 | * As in parallel_constraints, only inequality constraints that only |
26 | * involve input variables that do not occur in any other inequality |
27 | * constraints are considered. |
28 | * |
29 | * More specifically, given a set of constraints |
30 | * |
31 | * a x + b_i(p) >= 0 |
32 | * |
33 | * Replace this set by a single constraint |
34 | * |
35 | * a x + u >= 0 |
36 | * |
37 | * with u a new parameter with constraints |
38 | * |
39 | * u <= b_i(p) |
40 | * |
41 | * Any solution to the new system is also a solution for the original system |
42 | * since |
43 | * |
44 | * a x >= -u >= -b_i(p) |
45 | * |
46 | * Moreover, m = min_i(b_i(p)) satisfies the constraints on u and can |
47 | * therefore be plugged into the solution. |
48 | */ |
49 | static TYPE *SF(basic_map_partial_lexopt_symm,SUFFIX)( |
50 | __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom, |
51 | __isl_give isl_set **empty, int max, int first, int second) |
52 | { |
53 | int i, n, k; |
54 | int *list = NULL; |
55 | isl_size bmap_in, bmap_param, bmap_all; |
56 | unsigned n_in, n_out, n_div; |
57 | isl_ctx *ctx; |
58 | isl_vec *var = NULL; |
59 | isl_mat *cst = NULL; |
60 | isl_space *map_space, *set_space; |
61 | |
62 | map_space = isl_basic_map_get_space(bmap); |
63 | set_space = empty ? isl_basic_set_get_space(bset: dom) : NULL; |
64 | |
65 | bmap_in = isl_basic_map_dim(bmap, type: isl_dim_in); |
66 | bmap_param = isl_basic_map_dim(bmap, type: isl_dim_param); |
67 | bmap_all = isl_basic_map_dim(bmap, type: isl_dim_all); |
68 | if (bmap_in < 0 || bmap_param < 0 || bmap_all < 0) |
69 | goto error; |
70 | n_in = bmap_param + bmap_in; |
71 | n_out = bmap_all - n_in; |
72 | |
73 | ctx = isl_basic_map_get_ctx(bmap); |
74 | list = isl_alloc_array(ctx, int, bmap->n_ineq); |
75 | var = isl_vec_alloc(ctx, size: n_out); |
76 | if ((bmap->n_ineq && !list) || (n_out && !var)) |
77 | goto error; |
78 | |
79 | list[0] = first; |
80 | list[1] = second; |
81 | isl_seq_cpy(dst: var->el, src: bmap->ineq[first] + 1 + n_in, len: n_out); |
82 | for (i = second + 1, n = 2; i < bmap->n_ineq; ++i) { |
83 | if (isl_seq_eq(p1: var->el, p2: bmap->ineq[i] + 1 + n_in, len: n_out) && |
84 | all_single_occurrence(bmap, ineq: i, n: n_in)) |
85 | list[n++] = i; |
86 | } |
87 | |
88 | cst = isl_mat_alloc(ctx, n_row: n, n_col: 1 + n_in); |
89 | if (!cst) |
90 | goto error; |
91 | |
92 | for (i = 0; i < n; ++i) |
93 | isl_seq_cpy(dst: cst->row[i], src: bmap->ineq[list[i]], len: 1 + n_in); |
94 | |
95 | bmap = isl_basic_map_cow(bmap); |
96 | if (!bmap) |
97 | goto error; |
98 | for (i = n - 1; i >= 0; --i) |
99 | if (isl_basic_map_drop_inequality(bmap, pos: list[i]) < 0) |
100 | goto error; |
101 | |
102 | bmap = isl_basic_map_add_dims(bmap, type: isl_dim_in, n: 1); |
103 | bmap = isl_basic_map_extend_constraints(base: bmap, n_eq: 0, n_ineq: 1); |
104 | k = isl_basic_map_alloc_inequality(bmap); |
105 | if (k < 0) |
106 | goto error; |
107 | isl_seq_clr(p: bmap->ineq[k], len: 1 + n_in); |
108 | isl_int_set_si(bmap->ineq[k][1 + n_in], 1); |
109 | isl_seq_cpy(dst: bmap->ineq[k] + 1 + n_in + 1, src: var->el, len: n_out); |
110 | bmap = isl_basic_map_finalize(bmap); |
111 | |
112 | n_div = isl_basic_set_dim(bset: dom, type: isl_dim_div); |
113 | dom = isl_basic_set_add_dims(bset: dom, type: isl_dim_set, n: 1); |
114 | dom = isl_basic_set_extend_constraints(base: dom, n_eq: 0, n_ineq: n); |
115 | for (i = 0; i < n; ++i) { |
116 | k = isl_basic_set_alloc_inequality(bset: dom); |
117 | if (k < 0) |
118 | goto error; |
119 | isl_seq_cpy(dst: dom->ineq[k], src: cst->row[i], len: 1 + n_in); |
120 | isl_int_set_si(dom->ineq[k][1 + n_in], -1); |
121 | isl_seq_clr(p: dom->ineq[k] + 1 + n_in + 1, len: n_div); |
122 | } |
123 | |
124 | isl_vec_free(vec: var); |
125 | free(ptr: list); |
126 | |
127 | return SF(basic_map_partial_lexopt_symm_core,SUFFIX)(bmap, dom, empty, |
128 | max, cst, map_space, set_space); |
129 | error: |
130 | isl_space_free(space: map_space); |
131 | isl_space_free(space: set_space); |
132 | isl_mat_free(mat: cst); |
133 | isl_vec_free(vec: var); |
134 | free(ptr: list); |
135 | isl_basic_set_free(bset: dom); |
136 | isl_basic_map_free(bmap); |
137 | return NULL; |
138 | } |
139 | |
140 | /* Recursive part of isl_tab_basic_map_partial_lexopt*, after detecting |
141 | * equalities and removing redundant constraints. |
142 | * |
143 | * We first check if there are any parallel constraints (left). |
144 | * If not, we are in the base case. |
145 | * If there are parallel constraints, we replace them by a single |
146 | * constraint in basic_map_partial_lexopt_symm_pma and then call |
147 | * this function recursively to look for more parallel constraints. |
148 | */ |
149 | static __isl_give TYPE *SF(basic_map_partial_lexopt,SUFFIX)( |
150 | __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom, |
151 | __isl_give isl_set **empty, int max) |
152 | { |
153 | isl_bool par = isl_bool_false; |
154 | int first, second; |
155 | |
156 | if (!bmap) |
157 | goto error; |
158 | |
159 | if (bmap->ctx->opt->pip_symmetry) |
160 | par = parallel_constraints(bmap, first: &first, second: &second); |
161 | if (par < 0) |
162 | goto error; |
163 | if (!par) |
164 | return SF(basic_map_partial_lexopt_base,SUFFIX)(bmap, dom, |
165 | empty, max); |
166 | |
167 | return SF(basic_map_partial_lexopt_symm,SUFFIX)(bmap, dom, empty, max, |
168 | first, second); |
169 | error: |
170 | isl_basic_set_free(bset: dom); |
171 | isl_basic_map_free(bmap); |
172 | return NULL; |
173 | } |
174 | |
175 | /* Compute the lexicographic minimum (or maximum if "flags" includes |
176 | * ISL_OPT_MAX) of "bmap" over the domain "dom" and return the result as |
177 | * either a map or a piecewise multi-affine expression depending on TYPE. |
178 | * If "empty" is not NULL, then *empty is assigned a set that |
179 | * contains those parts of the domain where there is no solution. |
180 | * If "flags" includes ISL_OPT_FULL, then "dom" is NULL and the optimum |
181 | * should be computed over the domain of "bmap". "empty" is also NULL |
182 | * in this case. |
183 | * If "bmap" is marked as rational (ISL_BASIC_MAP_RATIONAL), |
184 | * then we compute the rational optimum. Otherwise, we compute |
185 | * the integral optimum. |
186 | * |
187 | * We perform some preprocessing. As the PILP solver does not |
188 | * handle implicit equalities very well, we first make sure all |
189 | * the equalities are explicitly available. |
190 | * |
191 | * We also add context constraints to the basic map and remove |
192 | * redundant constraints. This is only needed because of the |
193 | * way we handle simple symmetries. In particular, we currently look |
194 | * for symmetries on the constraints, before we set up the main tableau. |
195 | * It is then no good to look for symmetries on possibly redundant constraints. |
196 | * If the domain was extracted from the basic map, then there is |
197 | * no need to add back those constraints again. |
198 | */ |
199 | __isl_give TYPE *SF(isl_tab_basic_map_partial_lexopt,SUFFIX)( |
200 | __isl_take isl_basic_map *bmap, __isl_take isl_basic_set *dom, |
201 | __isl_give isl_set **empty, unsigned flags) |
202 | { |
203 | int max, full; |
204 | isl_bool compatible; |
205 | |
206 | if (empty) |
207 | *empty = NULL; |
208 | |
209 | full = ISL_FL_ISSET(flags, ISL_OPT_FULL); |
210 | if (full) |
211 | dom = extract_domain(bmap, flags); |
212 | compatible = isl_basic_map_compatible_domain(bmap, bset: dom); |
213 | if (compatible < 0) |
214 | goto error; |
215 | if (!compatible) |
216 | isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid, |
217 | "domain does not match input" , goto error); |
218 | |
219 | max = ISL_FL_ISSET(flags, ISL_OPT_MAX); |
220 | if (isl_basic_set_dim(bset: dom, type: isl_dim_all) == 0) |
221 | return SF(basic_map_partial_lexopt,SUFFIX)(bmap, dom, empty, |
222 | max); |
223 | |
224 | if (!full) |
225 | bmap = isl_basic_map_intersect_domain(bmap, |
226 | bset: isl_basic_set_copy(bset: dom)); |
227 | bmap = isl_basic_map_detect_equalities(bmap); |
228 | bmap = isl_basic_map_remove_redundancies(bmap); |
229 | |
230 | return SF(basic_map_partial_lexopt,SUFFIX)(bmap, dom, empty, max); |
231 | error: |
232 | isl_basic_set_free(bset: dom); |
233 | isl_basic_map_free(bmap); |
234 | return NULL; |
235 | } |
236 | |