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  */
SF(basic_map_partial_lexopt_symm,SUFFIX)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(dom) : NULL;
64 
65 	bmap_in = isl_basic_map_dim(bmap, isl_dim_in);
66 	bmap_param = isl_basic_map_dim(bmap, isl_dim_param);
67 	bmap_all = isl_basic_map_dim(bmap, 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, 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(var->el, bmap->ineq[first] + 1 + n_in, n_out);
82 	for (i = second + 1, n = 2; i < bmap->n_ineq; ++i) {
83 		if (isl_seq_eq(var->el, bmap->ineq[i] + 1 + n_in, n_out) &&
84 		    all_single_occurrence(bmap, i, n_in))
85 			list[n++] = i;
86 	}
87 
88 	cst = isl_mat_alloc(ctx, n, 1 + n_in);
89 	if (!cst)
90 		goto error;
91 
92 	for (i = 0; i < n; ++i)
93 		isl_seq_cpy(cst->row[i], bmap->ineq[list[i]], 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, list[i]) < 0)
100 			goto error;
101 
102 	bmap = isl_basic_map_add_dims(bmap, isl_dim_in, 1);
103 	bmap = isl_basic_map_extend_constraints(bmap, 0, 1);
104 	k = isl_basic_map_alloc_inequality(bmap);
105 	if (k < 0)
106 		goto error;
107 	isl_seq_clr(bmap->ineq[k], 1 + n_in);
108 	isl_int_set_si(bmap->ineq[k][1 + n_in], 1);
109 	isl_seq_cpy(bmap->ineq[k] + 1 + n_in + 1, var->el, n_out);
110 	bmap = isl_basic_map_finalize(bmap);
111 
112 	n_div = isl_basic_set_dim(dom, isl_dim_div);
113 	dom = isl_basic_set_add_dims(dom, isl_dim_set, 1);
114 	dom = isl_basic_set_extend_constraints(dom, 0, n);
115 	for (i = 0; i < n; ++i) {
116 		k = isl_basic_set_alloc_inequality(dom);
117 		if (k < 0)
118 			goto error;
119 		isl_seq_cpy(dom->ineq[k], cst->row[i], 1 + n_in);
120 		isl_int_set_si(dom->ineq[k][1 + n_in], -1);
121 		isl_seq_clr(dom->ineq[k] + 1 + n_in + 1, n_div);
122 	}
123 
124 	isl_vec_free(var);
125 	free(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(map_space);
131 	isl_space_free(set_space);
132 	isl_mat_free(cst);
133 	isl_vec_free(var);
134 	free(list);
135 	isl_basic_set_free(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  */
SF(basic_map_partial_lexopt,SUFFIX)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, &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(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  */
SF(isl_tab_basic_map_partial_lexopt,SUFFIX)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, 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(dom, 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 						    isl_basic_set_copy(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(dom);
233 	isl_basic_map_free(bmap);
234 	return NULL;
235 }
236